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:   DM             cdm;
 44:   DMLabel        bodyLabel, faceLabel, edgeLabel;
 45:   PetscContainer modelObj;
 46:   PetscInt       bodyID, faceID, edgeID;
 47:   ego           *bodies;
 48:   ego            model, geom, body, obj;
 49:   /* result has to hold derviatives, along with the value */
 50:   double         params[3], result[18], paramsV[16*3], resultV[16*3];
 51:   int            Nb, oclass, mtype, *senses;
 52:   Vec            coordinatesLocal;
 53:   PetscScalar   *coords = NULL;
 54:   PetscInt       Nv, v, Np = 0, pm;
 55: #endif
 56:   PetscInt       dE, d;

 60:   DMGetCoordinateDim(dm, &dE);
 61: #ifdef PETSC_HAVE_EGADS
 62:   DMGetLabel(dm, "EGADS Body ID",   &bodyLabel);
 63:   DMGetLabel(dm, "EGADS Face ID",   &faceLabel);
 64:   DMGetLabel(dm, "EGADS Edge ID",   &edgeLabel);
 65:   if (!bodyLabel || !faceLabel || !edgeLabel) {
 66:     for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
 67:     return(0);
 68:   }
 69:   DMGetCoordinateDM(dm, &cdm);
 70:   DMGetCoordinatesLocal(dm, &coordinatesLocal);
 71:   DMLabelGetValue(bodyLabel,   p, &bodyID);
 72:   DMLabelGetValue(faceLabel,   p, &faceID);
 73:   DMLabelGetValue(edgeLabel,   p, &edgeID);
 74:   PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
 75:   PetscContainerGetPointer(modelObj, (void **) &model);
 76:   EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
 77:   if (bodyID < 0 || bodyID >= Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %D is not in [0, %d)", bodyID, Nb);
 78:   body = bodies[bodyID];

 80:   if (edgeID >= 0)      {EG_objectBodyTopo(body, EDGE, edgeID, &obj); Np = 1;}
 81:   else if (faceID >= 0) {EG_objectBodyTopo(body, FACE, faceID, &obj); Np = 2;}
 82:   else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %D is not in edge or face label for EGADS", p);
 83:   /* Calculate parameters (t or u,v) for vertices */
 84:   DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 85:   Nv  /= dE;
 86:   if (Nv == 1) {
 87:     DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 88:     for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
 89:     return(0);
 90:   }
 91:   if (Nv > 16) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %D coordinates associated to point %D", Nv, p);
 92:   for (v = 0; v < Nv; ++v) {EG_invEvaluate(obj, &coords[v*dE], &paramsV[v*3], &resultV[v*3]);}
 93:   DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 94:   /* Calculate parameters (t or u,v) for new vertex at edge midpoint */
 95:   for (pm = 0; pm < Np; ++pm) {
 96:     params[pm] = 0.;
 97:     for (v = 0; v < Nv; ++v) {params[pm] += paramsV[v*3+pm];}
 98:     params[pm] /= Nv;
 99:   }
100:   /* TODO Check
101:     double range[4]; // [umin, umax, vmin, vmax]
102:     int    peri;
103:     EG_getRange(face, range, &peri);
104:     if ((paramsNew[0] < range[0]) || (paramsNew[0] > range[1]) || (paramsNew[1] < range[2]) || (paramsNew[1] > range[3])) SETERRQ();
105:   */
106:   /* Put coordinates for new vertex in result[] */
107:   EG_evaluate(obj, params, result);
108:   for (d = 0; d < dE; ++d) gcoords[d] = result[d];
109: #else
110:   for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
111: #endif
112:   return(0);
113: }