1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #ifdef PETSC_HAVE_EGADS 4 #include <egads.h> 5 #endif 6 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: 9 10 https://github.com/tpaviot/oce/tree/master/src/STEPControl 11 12 The STEP, and inner EXPRESS, formats are ISO standards, so they are documented 13 14 https://stackoverflow.com/questions/26774037/documentation-or-specification-for-step-and-stp-files 15 http://stepmod.sourceforge.net/express_model_spec/ 16 17 but again it seems that there has been a deliberate effort at obfuscation, probably to raise the bar for entrants. 18 */ 19 20 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. 23 24 Not collective 25 26 Input Parameters: 27 + dm - The DMPlex object 28 . p - The mesh point 29 - mcoords - A coordinate point lying on the mesh point 30 31 Output Parameter: 32 . gcoords - The closest coordinate point on the geometry model associated with 'p' to the given point 33 34 Note: Returns the original coordinates if no geometry model is found. Right now the only supported geometry model is EGADS. 35 36 Level: intermediate 37 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; 52 PetscErrorCode ierr; 53 54 PetscFunctionBegin; 55 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 56 #ifdef PETSC_HAVE_EGADS 57 ierr = DMGetLabel(dm, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr); 58 ierr = DMGetLabel(dm, "EGADS Face ID", &faceLabel);CHKERRQ(ierr); 59 ierr = DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);CHKERRQ(ierr); 60 if (!bodyLabel || !faceLabel || !edgeLabel) { 61 for (d = 0; d < cdim; ++d) gcoords[d] = mcoords[d]; 62 PetscFunctionReturn(0); 63 } 64 ierr = DMLabelGetValue(bodyLabel, p, &bodyID);CHKERRQ(ierr); 65 ierr = DMLabelGetValue(faceLabel, p, &faceID);CHKERRQ(ierr); 66 ierr = DMLabelGetValue(edgeLabel, p, &edgeID);CHKERRQ(ierr); 67 ierr = PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr); 68 ierr = PetscContainerGetPointer(modelObj, (void **) &model);CHKERRQ(ierr); 69 ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr); 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 ierr = EG_objectBodyTopo(body, EDGE, edgeID, &edge);CHKERRQ(ierr); 75 ierr = EG_invEvaluate(edge, point, params, result); 76 } else { 77 ierr = EG_objectBodyTopo(body, FACE, faceID, &face);CHKERRQ(ierr); 78 ierr = 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 PetscFunctionReturn(0); 85 } 86