1:  #include <petsc/private/dmpleximpl.h>

  3: #ifdef PETSC_HAVE_EGADS
  4: #include <egads.h>
  5: #endif

  7: /* We need to understand how to natively parse STEP files. There seems to be only one open source implementation of
  8:    the STEP parser contained in the OpenCASCADE package. It is enough to make a strong man weep:

 10:      https://github.com/tpaviot/oce/tree/master/src/STEPControl

 12:    The STEP, and inner EXPRESS, formats are ISO standards, so they are documented

 14:      https://stackoverflow.com/questions/26774037/documentation-or-specification-for-step-and-stp-files
 15:      http://stepmod.sourceforge.net/express_model_spec/

 17:    but again it seems that there has been a deliberate effort at obfuscation, probably to raise the bar for entrants.
 18: */


 21: /*@
 22:   DMPlexSnapToGeomModel - Given a coordinate point 'mcoords' on the mesh point 'p', return the closest coordinate point 'gcoords' on the geometry model associated with that point.

 24:   Not collective

 26:   Input Parameters:
 27: + dm      - The DMPlex object
 28: . p       - The mesh point
 29: - mcoords - A coordinate point lying on the mesh point

 31:   Output Parameter:
 32: . gcoords - The closest coordinate point on the geometry model associated with 'p' to the given point

 34:   Note: Returns the original coordinates if no geometry model is found. Right now the only supported geometry model is EGADS.

 36:   Level: intermediate

 38: .seealso: DMRefine(), DMPlexCreate(), DMPlexSetRefinementUniform()
 39: @*/
 40: PetscErrorCode DMPlexSnapToGeomModel(DM dm, PetscInt p, const PetscScalar mcoords[], PetscScalar gcoords[])
 41: {
 42: #ifdef PETSC_HAVE_EGADS
 43:   DMLabel        bodyLabel, faceLabel, edgeLabel;
 44:   PetscContainer modelObj;
 45:   PetscInt       bodyID, faceID, edgeID;
 46:   ego           *bodies;
 47:   ego            model, geom, body, face, edge;
 48:   double         point[3], params[3], result[3];
 49:   int            Nb, oclass, mtype, *senses;
 50: #endif
 51:   PetscInt       cdim, d;

 55:   DMGetCoordinateDim(dm, &cdim);
 56: #ifdef PETSC_HAVE_EGADS
 57:   DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
 58:   DMGetLabel(dm, "EGADS Face ID", &faceLabel);
 59:   DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
 60:   if (!bodyLabel || !faceLabel || !edgeLabel) {
 61:     for (d = 0; d < cdim; ++d) gcoords[d] = mcoords[d];
 62:     return(0);
 63:   }
 64:   DMLabelGetValue(bodyLabel, p, &bodyID);
 65:   DMLabelGetValue(faceLabel, p, &faceID);
 66:   DMLabelGetValue(edgeLabel, p, &edgeID);
 67:   PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
 68:   PetscContainerGetPointer(modelObj, (void **) &model);
 69:   EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
 70:   if (bodyID >= Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %D is not in [0, %d)", bodyID, Nb);
 71:   body = bodies[bodyID];
 72:   for (d = 0; d < cdim; ++d) point[d] = mcoords[d];
 73:   if (edgeID >= 0) {
 74:     EG_objectBodyTopo(body, EDGE, edgeID, &edge);
 75:     EG_invEvaluate(edge, point, params, result);
 76:   } else {
 77:     EG_objectBodyTopo(body, FACE, faceID, &face);
 78:     EG_invEvaluate(face, point, params, result);
 79:   }
 80:   for (d = 0; d < cdim; ++d) gcoords[d] = result[d];
 81: #else
 82:   for (d = 0; d < cdim; ++d) gcoords[d] = mcoords[d];
 83: #endif
 84:   return(0);
 85: }