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 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; 57 PetscErrorCode ierr; 58 59 PetscFunctionBegin; 60 ierr = DMGetCoordinateDim(dm, &dE);CHKERRQ(ierr); 61 #ifdef PETSC_HAVE_EGADS 62 ierr = DMGetLabel(dm, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr); 63 ierr = DMGetLabel(dm, "EGADS Face ID", &faceLabel);CHKERRQ(ierr); 64 ierr = DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);CHKERRQ(ierr); 65 if (!bodyLabel || !faceLabel || !edgeLabel) { 66 for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 67 PetscFunctionReturn(0); 68 } 69 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 70 ierr = DMGetCoordinatesLocal(dm, &coordinatesLocal);CHKERRQ(ierr); 71 ierr = DMLabelGetValue(bodyLabel, p, &bodyID);CHKERRQ(ierr); 72 ierr = DMLabelGetValue(faceLabel, p, &faceID);CHKERRQ(ierr); 73 ierr = DMLabelGetValue(edgeLabel, p, &edgeID);CHKERRQ(ierr); 74 ierr = PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr); 75 ierr = PetscContainerGetPointer(modelObj, (void **) &model);CHKERRQ(ierr); 76 ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr); 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]; 79 80 if (edgeID >= 0) {ierr = EG_objectBodyTopo(body, EDGE, edgeID, &obj);CHKERRQ(ierr); Np = 1;} 81 else if (faceID >= 0) {ierr = EG_objectBodyTopo(body, FACE, faceID, &obj);CHKERRQ(ierr); 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 ierr = DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);CHKERRQ(ierr); 85 Nv /= dE; 86 if (Nv == 1) { 87 ierr = DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);CHKERRQ(ierr); 88 for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 89 PetscFunctionReturn(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) {ierr = EG_invEvaluate(obj, &coords[v*dE], ¶msV[v*3], &resultV[v*3]);CHKERRQ(ierr);} 93 ierr = DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);CHKERRQ(ierr); 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 ierr = EG_getRange(face, range, &peri);CHKERRQ(ierr); 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 ierr = EG_evaluate(obj, params, result);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 113 } 114