1 #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/ 2 3 PetscFunctionList DMGeomModelList = NULL; 4 PetscBool DMGeomModelRegisterAllCalled = PETSC_FALSE; 5 6 #if defined(PETSC_HAVE_EGADS) 7 PETSC_INTERN PetscErrorCode DMSnapToGeomModel_EGADS(DM, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]); 8 #endif 9 10 static PetscErrorCode DMSnapToGeomModelBall(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[]) 11 { 12 PetscInt val; 13 14 PetscFunctionBeginUser; 15 PetscCall(DMGetLabelValue(dm, "marker", p, &val)); 16 if (val >= 0) { 17 PetscReal norm = 0.; 18 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 } 27 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; 32 33 PetscFunctionBeginUser; 34 PetscCall(DMGetBoundingBox(dm, gmin, gmax)); 35 PetscCall(DMGetLabelValue(dm, "generatrix", p, &val)); 36 if (val >= 0) { 37 PetscReal norm = 0.; 38 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 } 48 49 /*@C 50 DMGeomModelRegisterAll - Registers all of the geometry model methods in the `DM` package. 51 52 Not Collective 53 54 Level: advanced 55 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 } 71 72 /*@C 73 DMGeomModelRegister - Adds a geometry model to `DM` 74 75 Not Collective, No Fortran Support 76 77 Input Parameters: 78 + sname - name of a new user-defined geometry model 79 - fnc - geometry model function 80 81 Example Usage: 82 .vb 83 DMGeomModelRegister("my_geom_model", MySnapToGeomModel); 84 .ve 85 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 94 95 Level: advanced 96 97 Note: 98 `DMGeomModelRegister()` may be called multiple times to add several user-defined generators 99 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 } 108 109 PetscErrorCode DMGeomModelRegisterDestroy(void) 110 { 111 PetscFunctionBegin; 112 PetscCall(PetscFunctionListDestroy(&DMGeomModelList)); 113 DMGeomModelRegisterAllCalled = PETSC_FALSE; 114 PetscFunctionReturn(PETSC_SUCCESS); 115 } 116 117 /*@ 118 DMSetSnapToGeomModel - Choose a geometry model for this `DM`. 119 120 Not Collective 121 122 Input Parameters: 123 + dm - The `DM` object 124 - name - A geometry model name, or `NULL` for the default 125 126 Level: intermediate 127 128 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexCreate()`, `DMSnapToGeomModel()` 129 @*/ 130 PetscErrorCode DMSetSnapToGeomModel(DM dm, const char name[]) 131 { 132 char geomname[PETSC_MAX_PATH_LEN]; 133 PetscBool flg; 134 135 PetscFunctionBegin; 136 if (!name && dm->ops->snaptogeommodel) PetscFunctionReturn(PETSC_SUCCESS); 137 PetscCall(PetscOptionsGetString(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_geom_model", geomname, sizeof(geomname), &flg)); 138 if (flg) name = geomname; 139 if (!name) { 140 PetscObject modelObj; 141 142 PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", &modelObj)); 143 if (modelObj) name = "egads"; 144 else { 145 PetscCall(PetscObjectQuery((PetscObject)dm, "EGADSlite Model", &modelObj)); 146 if (modelObj) name = "egads"; 147 } 148 } 149 if (!name) PetscFunctionReturn(PETSC_SUCCESS); 150 151 PetscCall(PetscFunctionListFind(DMGeomModelList, name, &dm->ops->snaptogeommodel)); 152 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); 153 PetscFunctionReturn(PETSC_SUCCESS); 154 } 155 156 /*@ 157 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. 158 159 Not Collective 160 161 Input Parameters: 162 + dm - The `DMPLEX` object 163 . p - The mesh point 164 . dE - The coordinate dimension 165 - mcoords - A coordinate point lying on the mesh point 166 167 Output Parameter: 168 . gcoords - The closest coordinate point on the geometry model associated with 'p' to the given point 169 170 Level: intermediate 171 172 Note: 173 Returns the original coordinates if no geometry model is found. 174 175 The coordinate dimension may be different from the coordinate dimension of the `dm`, for example if the transformation is extrusion. 176 177 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexCreate()`, `DMPlexSetRefinementUniform()` 178 @*/ 179 PetscErrorCode DMSnapToGeomModel(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[]) 180 { 181 PetscFunctionBegin; 182 if (!dm->ops->snaptogeommodel) 183 for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 184 else PetscUseTypeMethod(dm, snaptogeommodel, p, dE, mcoords, gcoords); 185 PetscFunctionReturn(PETSC_SUCCESS); 186 } 187