1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2
3 /*@C
4 PetscFEGeomCreate - Create a `PetscFEGeom` object to manage geometry for a group of cells
5
6 Input Parameters:
7 + quad - A `PetscQuadrature` determining the tabulation
8 . numCells - The number of cells in the group
9 . dimEmbed - The coordinate dimension
10 - mode - Type of geometry data to store
11
12 Output Parameter:
13 . geom - The `PetscFEGeom` object, which is a struct not a `PetscObject`
14
15 Level: beginner
16
17 .seealso: `PetscFEGeom`, `PetscQuadrature`, `PetscFEGeomDestroy()`, `PetscFEGeomComplete()`
18 @*/
PetscFEGeomCreate(PetscQuadrature quad,PetscInt numCells,PetscInt dimEmbed,PetscFEGeomMode mode,PetscFEGeom ** geom)19 PetscErrorCode PetscFEGeomCreate(PetscQuadrature quad, PetscInt numCells, PetscInt dimEmbed, PetscFEGeomMode mode, PetscFEGeom **geom)
20 {
21 PetscFEGeom *g;
22 PetscInt dim, Nq, N;
23 const PetscReal *p;
24
25 PetscFunctionBegin;
26 PetscCall(PetscQuadratureGetData(quad, &dim, NULL, &Nq, &p, NULL));
27 PetscCall(PetscNew(&g));
28 g->mode = mode;
29 g->xi = p;
30 g->numCells = numCells;
31 g->numPoints = Nq;
32 g->dim = dim;
33 g->dimEmbed = dimEmbed;
34 N = numCells * Nq;
35 PetscCall(PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ));
36 if (mode == PETSC_FEGEOM_BOUNDARY || mode == PETSC_FEGEOM_COHESIVE) {
37 PetscCall(PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n));
38 PetscCall(PetscCalloc6(N * dimEmbed * dimEmbed, &g->suppJ[0], N * dimEmbed * dimEmbed, &g->suppJ[1], N * dimEmbed * dimEmbed, &g->suppInvJ[0], N * dimEmbed * dimEmbed, &g->suppInvJ[1], N, &g->suppDetJ[0], N, &g->suppDetJ[1]));
39 }
40 PetscCall(PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ));
41 *geom = g;
42 PetscFunctionReturn(PETSC_SUCCESS);
43 }
44
45 /*@C
46 PetscFEGeomDestroy - Destroy a `PetscFEGeom` object
47
48 Input Parameter:
49 . geom - `PetscFEGeom` object
50
51 Level: beginner
52
53 .seealso: `PetscFEGeom`, `PetscFEGeomCreate()`
54 @*/
PetscFEGeomDestroy(PetscFEGeom ** geom)55 PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom)
56 {
57 PetscFunctionBegin;
58 if (!*geom) PetscFunctionReturn(PETSC_SUCCESS);
59 PetscCall(PetscFree3((*geom)->v, (*geom)->J, (*geom)->detJ));
60 PetscCall(PetscFree((*geom)->invJ));
61 PetscCall(PetscFree2((*geom)->face, (*geom)->n));
62 PetscCall(PetscFree6((*geom)->suppJ[0], (*geom)->suppJ[1], (*geom)->suppInvJ[0], (*geom)->suppInvJ[1], (*geom)->suppDetJ[0], (*geom)->suppDetJ[1]));
63 PetscCall(PetscFree(*geom));
64 PetscFunctionReturn(PETSC_SUCCESS);
65 }
66
67 /*@C
68 PetscFEGeomGetChunk - Get a chunk of cells in the group as a `PetscFEGeom`
69
70 Input Parameters:
71 + geom - `PetscFEGeom` object
72 . cStart - The first cell in the chunk
73 - cEnd - The first cell not in the chunk
74
75 Output Parameter:
76 . chunkGeom - an array of cells of length `cEnd` - `cStart`
77
78 Level: intermediate
79
80 Note:
81 Use `PetscFEGeomRestoreChunk()` to return the result
82
83 .seealso: `PetscFEGeom`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
84 @*/
PetscFEGeomGetChunk(PetscFEGeom * geom,PetscInt cStart,PetscInt cEnd,PetscFEGeom * chunkGeom[])85 PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom *chunkGeom[])
86 {
87 PetscInt Nq;
88 PetscInt dE;
89
90 PetscFunctionBegin;
91 PetscAssertPointer(geom, 1);
92 PetscAssertPointer(chunkGeom, 4);
93 if (!*chunkGeom) PetscCall(PetscNew(chunkGeom));
94 Nq = geom->numPoints;
95 dE = geom->dimEmbed;
96 (*chunkGeom)->mode = geom->mode;
97 (*chunkGeom)->dim = geom->dim;
98 (*chunkGeom)->dimEmbed = geom->dimEmbed;
99 (*chunkGeom)->numPoints = geom->numPoints;
100 (*chunkGeom)->numCells = cEnd - cStart;
101 (*chunkGeom)->xi = geom->xi;
102 (*chunkGeom)->v = PetscSafePointerPlusOffset(geom->v, Nq * dE * cStart);
103 (*chunkGeom)->J = PetscSafePointerPlusOffset(geom->J, Nq * dE * dE * cStart);
104 (*chunkGeom)->invJ = PetscSafePointerPlusOffset(geom->invJ, Nq * dE * dE * cStart);
105 (*chunkGeom)->detJ = PetscSafePointerPlusOffset(geom->detJ, Nq * cStart);
106 (*chunkGeom)->n = PetscSafePointerPlusOffset(geom->n, Nq * dE * cStart);
107 (*chunkGeom)->face = PetscSafePointerPlusOffset(geom->face, cStart);
108 (*chunkGeom)->suppJ[0] = PetscSafePointerPlusOffset(geom->suppJ[0], Nq * dE * dE * cStart);
109 (*chunkGeom)->suppJ[1] = PetscSafePointerPlusOffset(geom->suppJ[1], Nq * dE * dE * cStart);
110 (*chunkGeom)->suppInvJ[0] = PetscSafePointerPlusOffset(geom->suppInvJ[0], Nq * dE * dE * cStart);
111 (*chunkGeom)->suppInvJ[1] = PetscSafePointerPlusOffset(geom->suppInvJ[1], Nq * dE * dE * cStart);
112 (*chunkGeom)->suppDetJ[0] = PetscSafePointerPlusOffset(geom->suppDetJ[0], Nq * cStart);
113 (*chunkGeom)->suppDetJ[1] = PetscSafePointerPlusOffset(geom->suppDetJ[1], Nq * cStart);
114 (*chunkGeom)->isAffine = geom->isAffine;
115 PetscFunctionReturn(PETSC_SUCCESS);
116 }
117
118 /*@C
119 PetscFEGeomRestoreChunk - Restore the chunk obtained with `PetscFEGeomCreateChunk()`
120
121 Input Parameters:
122 + geom - `PetscFEGeom` object
123 . cStart - The first cell in the chunk
124 . cEnd - The first cell not in the chunk
125 - chunkGeom - The chunk of cells
126
127 Level: intermediate
128
129 .seealso: `PetscFEGeom`, `PetscFEGeomGetChunk()`, `PetscFEGeomCreate()`
130 @*/
PetscFEGeomRestoreChunk(PetscFEGeom * geom,PetscInt cStart,PetscInt cEnd,PetscFEGeom ** chunkGeom)131 PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
132 {
133 PetscFunctionBegin;
134 PetscCall(PetscFree(*chunkGeom));
135 PetscFunctionReturn(PETSC_SUCCESS);
136 }
137
138 /*@C
139 PetscFEGeomGetPoint - Get the geometry for cell `c` at point `p` as a `PetscFEGeom`
140
141 Input Parameters:
142 + geom - `PetscFEGeom` object
143 . c - The cell
144 . p - The point
145 - pcoords - The reference coordinates of point `p`, or `NULL`
146
147 Output Parameter:
148 . pgeom - The geometry of cell `c` at point `p`
149
150 Level: intermediate
151
152 Notes:
153 For affine geometries, this only copies to `pgeom` at point 0. Since we copy pointers into `pgeom`,
154 nothing needs to be done with it afterwards.
155
156 In the affine case, `pgeom` must have storage for the integration point coordinates in pgeom->v if `pcoords` is passed in.
157
158 .seealso: `PetscFEGeom`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
159 @*/
PetscFEGeomGetPoint(PetscFEGeom * geom,PetscInt c,PetscInt p,const PetscReal pcoords[],PetscFEGeom * pgeom)160 PetscErrorCode PetscFEGeomGetPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, const PetscReal pcoords[], PetscFEGeom *pgeom)
161 {
162 const PetscInt dim = geom->dim;
163 const PetscInt dE = geom->dimEmbed;
164 const PetscInt Np = geom->numPoints;
165
166 PetscFunctionBeginHot;
167 pgeom->mode = geom->mode;
168 pgeom->dim = dim;
169 pgeom->dimEmbed = dE;
170 //pgeom->isAffine = geom->isAffine;
171 if (geom->isAffine) {
172 if (!p) {
173 pgeom->xi = geom->xi;
174 pgeom->J = &geom->J[c * Np * dE * dE];
175 pgeom->invJ = &geom->invJ[c * Np * dE * dE];
176 pgeom->detJ = &geom->detJ[c * Np];
177 pgeom->n = PetscSafePointerPlusOffset(geom->n, c * Np * dE);
178 }
179 if (pcoords) CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c * Np * dE], pgeom->J, pcoords, pgeom->v);
180 } else {
181 pgeom->v = &geom->v[(c * Np + p) * dE];
182 pgeom->J = &geom->J[(c * Np + p) * dE * dE];
183 pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE];
184 pgeom->detJ = &geom->detJ[c * Np + p];
185 pgeom->n = PetscSafePointerPlusOffset(geom->n, (c * Np + p) * dE);
186 }
187 PetscFunctionReturn(PETSC_SUCCESS);
188 }
189
190 /*@C
191 PetscFEGeomGetCellPoint - Get the cell geometry for cell `c` at point `p` as a `PetscFEGeom`
192
193 Input Parameters:
194 + geom - `PetscFEGeom` object
195 . c - The cell
196 - p - The point
197
198 Output Parameter:
199 . pgeom - The cell geometry of cell `c` at point `p`
200
201 Level: intermediate
202
203 Notes:
204 For PETSC_FEGEOM_BOUNDARY mode, this gives the geometry for supporting cell 0. For PETSC_FEGEOM_COHESIVE mode,
205 this gives the bulk geometry for that internal face.
206
207 For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into `pgeom`,
208 nothing needs to be done with it afterwards.
209
210 .seealso: `PetscFEGeom`, `PetscFEGeomMode`, `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()`
211 @*/
PetscFEGeomGetCellPoint(PetscFEGeom * geom,PetscInt c,PetscInt p,PetscFEGeom * pgeom)212 PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom)
213 {
214 const PetscBool bd = geom->mode == PETSC_FEGEOM_BOUNDARY ? PETSC_TRUE : PETSC_FALSE;
215 const PetscInt dim = bd ? geom->dimEmbed : geom->dim;
216 const PetscInt dE = geom->dimEmbed;
217 const PetscInt Np = geom->numPoints;
218
219 PetscFunctionBeginHot;
220 pgeom->mode = geom->mode;
221 pgeom->dim = dim;
222 pgeom->dimEmbed = dE;
223 //pgeom->isAffine = geom->isAffine;
224 if (geom->isAffine) {
225 if (!p) {
226 if (bd) {
227 pgeom->J = &geom->suppJ[0][c * Np * dE * dE];
228 pgeom->invJ = &geom->suppInvJ[0][c * Np * dE * dE];
229 pgeom->detJ = &geom->suppDetJ[0][c * Np];
230 } else {
231 pgeom->J = &geom->J[c * Np * dE * dE];
232 pgeom->invJ = &geom->invJ[c * Np * dE * dE];
233 pgeom->detJ = &geom->detJ[c * Np];
234 }
235 }
236 } else {
237 if (bd) {
238 pgeom->J = &geom->suppJ[0][(c * Np + p) * dE * dE];
239 pgeom->invJ = &geom->suppInvJ[0][(c * Np + p) * dE * dE];
240 pgeom->detJ = &geom->suppDetJ[0][c * Np + p];
241 } else {
242 pgeom->J = &geom->J[(c * Np + p) * dE * dE];
243 pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE];
244 pgeom->detJ = &geom->detJ[c * Np + p];
245 }
246 }
247 PetscFunctionReturn(PETSC_SUCCESS);
248 }
249
250 /*@C
251 PetscFEGeomComplete - Calculate derived quantities from a base geometry specification
252
253 Input Parameter:
254 . geom - `PetscFEGeom` object
255
256 Level: intermediate
257
258 .seealso: `PetscFEGeom`, `PetscFEGeomCreate()`
259 @*/
PetscFEGeomComplete(PetscFEGeom * geom)260 PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom)
261 {
262 PetscInt i, j, N, dE;
263
264 PetscFunctionBeginHot;
265 N = geom->numPoints * geom->numCells;
266 dE = geom->dimEmbed;
267 switch (dE) {
268 case 3:
269 for (i = 0; i < N; i++) {
270 DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]);
271 if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]);
272 }
273 break;
274 case 2:
275 for (i = 0; i < N; i++) {
276 DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]);
277 if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]);
278 }
279 break;
280 case 1:
281 for (i = 0; i < N; i++) {
282 geom->detJ[i] = PetscAbsReal(geom->J[i]);
283 if (geom->invJ) geom->invJ[i] = 1. / geom->J[i];
284 }
285 break;
286 }
287 if (geom->n) {
288 for (i = 0; i < N; i++) {
289 for (j = 0; j < dE; j++) geom->n[dE * i + j] = geom->J[dE * dE * i + dE * j + dE - 1] * ((dE == 2) ? -1. : 1.);
290 }
291 }
292 PetscFunctionReturn(PETSC_SUCCESS);
293 }
294