Actual source code: dmgeommodel.c

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

  3: PetscFunctionList DMGeomModelList              = NULL;
  4: PetscBool         DMGeomModelRegisterAllCalled = PETSC_FALSE;

  6: #if defined(PETSC_HAVE_EGADS)
  7: PETSC_INTERN PetscErrorCode DMSnapToGeomModel_EGADS(DM, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
  8: #endif

 10: static PetscErrorCode DMSnapToGeomModelBall(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[])
 11: {
 12:   PetscInt val;

 14:   PetscFunctionBeginUser;
 15:   PetscCall(DMGetLabelValue(dm, "marker", p, &val));
 16:   if (val >= 0) {
 17:     PetscReal norm = 0.;

 19:     for (PetscInt d = 0; d < dE; ++d) norm += PetscSqr(PetscRealPart(mcoords[d]));
 20:     norm = PetscSqrtReal(norm);
 21:     for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d] / norm;
 22:   } else {
 23:     for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
 24:   }
 25:   PetscFunctionReturn(PETSC_SUCCESS);
 26: }

 28: static PetscErrorCode DMSnapToGeomModelCylinder(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[])
 29: {
 30:   PetscReal gmin[3], gmax[3];
 31:   PetscInt  val;

 33:   PetscFunctionBeginUser;
 34:   PetscCall(DMGetBoundingBox(dm, gmin, gmax));
 35:   PetscCall(DMGetLabelValue(dm, "generatrix", p, &val));
 36:   if (val >= 0) {
 37:     PetscReal norm = 0.;

 39:     for (PetscInt d = 0; d < dE - 1; ++d) norm += PetscSqr(PetscRealPart(mcoords[d]));
 40:     norm = PetscSqrtReal(norm);
 41:     for (PetscInt d = 0; d < dE - 1; ++d) gcoords[d] = mcoords[d] * gmax[0] / norm;
 42:     gcoords[dE - 1] = mcoords[dE - 1];
 43:   } else {
 44:     for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
 45:   }
 46:   PetscFunctionReturn(PETSC_SUCCESS);
 47: }

 49: /*@C
 50:   DMGeomModelRegisterAll - Registers all of the geometry model methods in the `DM` package.

 52:   Not Collective

 54:   Level: advanced

 56: .seealso: `DM`, `DMGeomModelRegisterDestroy()`
 57: @*/
 58: PetscErrorCode DMGeomModelRegisterAll(void)
 59: {
 60:   PetscFunctionBegin;
 61:   if (DMGeomModelRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
 62:   DMGeomModelRegisterAllCalled = PETSC_TRUE;
 63:   PetscCall(DMGeomModelRegister("ball", DMSnapToGeomModelBall));
 64:   PetscCall(DMGeomModelRegister("cylinder", DMSnapToGeomModelCylinder));
 65: #if defined(PETSC_HAVE_EGADS)
 66:   // FIXME: Brandon uses DMPlexSnapToGeomModel() here instead
 67:   PetscCall(DMGeomModelRegister("egads", DMSnapToGeomModel_EGADS));
 68: #endif
 69:   PetscFunctionReturn(PETSC_SUCCESS);
 70: }

 72: /*@C
 73:   DMGeomModelRegister -  Adds a geometry model to `DM`

 75:   Not Collective, No Fortran Support

 77:   Input Parameters:
 78: + sname - name of a new user-defined geometry model
 79: - fnc   - geometry model function

 81:   Example Usage:
 82: .vb
 83:    DMGeomModelRegister("my_geom_model", MySnapToGeomModel);
 84: .ve

 86:   Then, your generator can be chosen with the procedural interface via
 87: .vb
 88:   DMSetGeomModel(dm, "my_geom_model",...)
 89: .ve
 90:   or at runtime via the option
 91: .vb
 92:   -dm_geom_model my_geom_model
 93: .ve

 95:   Level: advanced

 97:   Note:
 98:   `DMGeomModelRegister()` may be called multiple times to add several user-defined generators

100: .seealso: `DM`, `DMGeomModelRegisterAll()`, `DMPlexGeomModel()`, `DMGeomModelRegisterDestroy()`
101: @*/
102: PetscErrorCode DMGeomModelRegister(const char sname[], PetscErrorCode (*fnc)(DM, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]))
103: {
104:   PetscFunctionBegin;
105:   PetscCall(PetscFunctionListAdd(&DMGeomModelList, sname, fnc));
106:   PetscFunctionReturn(PETSC_SUCCESS);
107: }

109: /*@
110:   DMGeomModelRegisterDestroy - Frees the list of `DM` geometry models that were registered by `DMGeomModelRegister()` or `DMGeomModelRegisterAll()`.

112:   Not Collective

114:   Level: advanced

116: .seealso: `DM`, `DMGeomModelRegister()`, `DMGeomModelRegisterAll()`
117: @*/
118: PetscErrorCode DMGeomModelRegisterDestroy(void)
119: {
120:   PetscFunctionBegin;
121:   PetscCall(PetscFunctionListDestroy(&DMGeomModelList));
122:   DMGeomModelRegisterAllCalled = PETSC_FALSE;
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: /*@
127:   DMSetSnapToGeomModel - Choose a geometry model for this `DM`.

129:   Not Collective

131:   Input Parameters:
132: + dm   - The `DM` object
133: - name - A geometry model name, or `NULL` for the default

135:   Level: intermediate

137: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexCreate()`, `DMSnapToGeomModel()`
138: @*/
139: PetscErrorCode DMSetSnapToGeomModel(DM dm, const char name[])
140: {
141:   char      geomname[PETSC_MAX_PATH_LEN];
142:   PetscBool flg;

144:   PetscFunctionBegin;
145:   if (!name && dm->ops->snaptogeommodel) PetscFunctionReturn(PETSC_SUCCESS);
146:   PetscCall(PetscOptionsGetString(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_geom_model", geomname, sizeof(geomname), &flg));
147:   if (flg) name = geomname;
148:   if (!name) {
149:     PetscObject modelObj;

151:     PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", &modelObj));
152:     if (modelObj) name = "egads";
153:     else {
154:       PetscCall(PetscObjectQuery((PetscObject)dm, "EGADSlite Model", &modelObj));
155:       if (modelObj) name = "egads";
156:     }
157:   }
158:   if (!name) PetscFunctionReturn(PETSC_SUCCESS);

160:   PetscCall(PetscFunctionListFind(DMGeomModelList, name, &dm->ops->snaptogeommodel));
161:   PetscCheck(dm->ops->snaptogeommodel, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Geometry model %s not registered; you may need to add --download-%s to your ./configure options", name, name);
162:   PetscFunctionReturn(PETSC_SUCCESS);
163: }

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

168:   Not Collective

170:   Input Parameters:
171: + dm      - The `DMPLEX` object
172: . p       - The mesh point
173: . dE      - The coordinate dimension
174: - mcoords - A coordinate point lying on the mesh point

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

179:   Level: intermediate

181:   Note:
182:   Returns the original coordinates if no geometry model is found.

184:   The coordinate dimension may be different from the coordinate dimension of the `dm`, for example if the transformation is extrusion.

186: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexCreate()`, `DMPlexSetRefinementUniform()`
187: @*/
188: PetscErrorCode DMSnapToGeomModel(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[])
189: {
190:   PetscFunctionBegin;
191:   if (!dm->ops->snaptogeommodel)
192:     for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
193:   else PetscUseTypeMethod(dm, snaptogeommodel, p, dE, mcoords, gcoords);
194:   PetscFunctionReturn(PETSC_SUCCESS);
195: }