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

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

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

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

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

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

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

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

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

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

 53:   Not Collective

 55:   Level: advanced

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

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

 76:   Not Collective, No Fortran Support

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

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

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

 96:   Level: advanced

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

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

110: extern PetscBool DMGeomModelRegisterAllCalled;

112: PetscErrorCode DMGeomModelRegisterDestroy(void)
113: {
114:   PetscFunctionBegin;
115:   PetscCall(PetscFunctionListDestroy(&DMGeomModelList));
116:   DMGeomModelRegisterAllCalled = PETSC_FALSE;
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

120: /*@
121:   DMSetSnapToGeomModel - Choose a geometry model for this `DM`.

123:   Not Collective

125:   Input Parameters:
126: + dm   - The `DM` object
127: - name - A geometry model name, or `NULL` for the default

129:   Level: intermediate

131: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexCreate()`, `DMSnapToGeomModel()`
132: @*/
133: PetscErrorCode DMSetSnapToGeomModel(DM dm, const char name[])
134: {
135:   char      geomname[PETSC_MAX_PATH_LEN];
136:   PetscBool flg;

138:   PetscFunctionBegin;
139:   if (!name && dm->ops->snaptogeommodel) PetscFunctionReturn(PETSC_SUCCESS);
140:   PetscCall(PetscOptionsGetString(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_geom_model", geomname, sizeof(geomname), &flg));
141:   if (flg) name = geomname;
142:   if (!name) {
143:     PetscObject modelObj;

145:     PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", &modelObj));
146:     if (modelObj) name = "egads";
147:     else {
148:       PetscCall(PetscObjectQuery((PetscObject)dm, "EGADSlite Model", &modelObj));
149:       if (modelObj) name = "egads";
150:     }
151:   }
152:   if (!name) PetscFunctionReturn(PETSC_SUCCESS);

154:   PetscCall(PetscFunctionListFind(DMGeomModelList, name, &dm->ops->snaptogeommodel));
155:   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);
156:   PetscFunctionReturn(PETSC_SUCCESS);
157: }

159: /*@
160:   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.

162:   Not Collective

164:   Input Parameters:
165: + dm      - The `DMPLEX` object
166: . p       - The mesh point
167: . dE      - The coordinate dimension
168: - mcoords - A coordinate point lying on the mesh point

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

173:   Level: intermediate

175:   Note:
176:   Returns the original coordinates if no geometry model is found.

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

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