Actual source code: fegeom.c
1: #include <petsc/private/petscfeimpl.h>
3: /*@C
4: PetscFEGeomCreate - Create a PetscFEGeom object to manage geometry for a group of cells
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: - faceData - Flag to construct geometry data for the faces
12: Output Parameter:
13: . geom - The PetscFEGeom object
15: Level: beginner
17: .seealso: PetscFEGeomDestroy(), PetscFEGeomComplete()
18: @*/
19: PetscErrorCode PetscFEGeomCreate(PetscQuadrature quad, PetscInt numCells, PetscInt dimEmbed, PetscBool faceData, PetscFEGeom **geom)
20: {
21: PetscFEGeom *g;
22: PetscInt dim, Nq, N;
23: const PetscReal *p;
24: PetscErrorCode ierr;
27: PetscQuadratureGetData(quad,&dim,NULL,&Nq,&p,NULL);
28: PetscNew(&g);
29: g->xi = p;
30: g->numCells = numCells;
31: g->numPoints = Nq;
32: g->dim = dim;
33: g->dimEmbed = dimEmbed;
34: g->isHybrid = PETSC_FALSE;
35: N = numCells * Nq;
36: PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ);
37: if (faceData) {
38: PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n);
39: PetscCalloc6(N * dimEmbed * dimEmbed, &(g->suppJ[0]), N * dimEmbed * dimEmbed, &(g->suppJ[1]),
40: N * dimEmbed * dimEmbed, &(g->suppInvJ[0]), N * dimEmbed * dimEmbed, &(g->suppInvJ[1]),
41: N, &(g->suppDetJ[0]), N, &(g->suppDetJ[1]));
42: }
43: PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ);
44: *geom = g;
45: return(0);
46: }
48: /*@C
49: PetscFEGeomDestroy - Destroy a PetscFEGeom object
51: Input Parameter:
52: . geom - PetscFEGeom object
54: Level: beginner
56: .seealso: PetscFEGeomCreate()
57: @*/
58: PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom)
59: {
63: if (!*geom) return(0);
64: PetscFree3((*geom)->v,(*geom)->J,(*geom)->detJ);
65: PetscFree((*geom)->invJ);
66: PetscFree2((*geom)->face,(*geom)->n);
67: PetscFree6((*geom)->suppJ[0],(*geom)->suppJ[1],(*geom)->suppInvJ[0],(*geom)->suppInvJ[1],(*geom)->suppDetJ[0],(*geom)->suppDetJ[1]);
68: PetscFree(*geom);
69: return(0);
70: }
72: /*@C
73: PetscFEGeomGetChunk - Get a chunk of cells in the group as a PetscFEGeom
75: Input Parameters:
76: + geom - PetscFEGeom object
77: . cStart - The first cell in the chunk
78: - cEnd - The first cell not in the chunk
80: Output Parameter:
81: . chunkGeom - The chunk of cells
83: Level: intermediate
85: .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate()
86: @*/
87: PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
88: {
89: PetscInt Nq;
90: PetscInt dE;
96: if (!(*chunkGeom)) {
97: PetscNew(chunkGeom);
98: }
99: Nq = geom->numPoints;
100: dE= geom->dimEmbed;
101: (*chunkGeom)->dim = geom->dim;
102: (*chunkGeom)->dimEmbed = geom->dimEmbed;
103: (*chunkGeom)->numPoints = geom->numPoints;
104: (*chunkGeom)->numCells = cEnd - cStart;
105: (*chunkGeom)->xi = geom->xi;
106: (*chunkGeom)->v = &geom->v[Nq*dE*cStart];
107: (*chunkGeom)->J = &geom->J[Nq*dE*dE*cStart];
108: (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq*dE*dE*cStart] : NULL;
109: (*chunkGeom)->detJ = &geom->detJ[Nq*cStart];
110: (*chunkGeom)->n = geom->n ? &geom->n[Nq*dE*cStart] : NULL;
111: (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL;
112: (*chunkGeom)->suppJ[0] = geom->suppJ[0] ? &geom->suppJ[0][Nq*dE*dE*cStart] : NULL;
113: (*chunkGeom)->suppJ[1] = geom->suppJ[1] ? &geom->suppJ[1][Nq*dE*dE*cStart] : NULL;
114: (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq*dE*dE*cStart] : NULL;
115: (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq*dE*dE*cStart] : NULL;
116: (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq*cStart] : NULL;
117: (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq*cStart] : NULL;
118: (*chunkGeom)->isAffine = geom->isAffine;
119: return(0);
120: }
122: /*@C
123: PetscFEGeomRestoreChunk - Restore the chunk
125: Input Parameters:
126: + geom - PetscFEGeom object
127: . cStart - The first cell in the chunk
128: . cEnd - The first cell not in the chunk
129: - chunkGeom - The chunk of cells
131: Level: intermediate
133: .seealso: PetscFEGeomGetChunk(), PetscFEGeomCreate()
134: @*/
135: PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
136: {
140: PetscFree(*chunkGeom);
141: return(0);
142: }
144: /*@C
145: PetscFEGeomGetPoint - Get the geometry for cell c at point p as a PetscFEGeom
147: Input Parameters:
148: + geom - PetscFEGeom object
149: . c - The cell
150: . p - The point
151: - pcoords - The reference coordinates of point p, or NULL
153: Output Parameter:
154: . pgeom - The geometry of cell c at point p
156: Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom,
157: nothing needs to be done with it afterwards.
159: In the affine case, pgeom must have storage for the integration point coordinates in pgeom->v if pcoords is passed in.
161: Level: intermediate
163: .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate()
164: @*/
165: PetscErrorCode PetscFEGeomGetPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, const PetscReal pcoords[], PetscFEGeom *pgeom)
166: {
167: const PetscInt dim = geom->dim;
168: const PetscInt dE = geom->dimEmbed;
169: const PetscInt Np = geom->numPoints;
172: pgeom->dim = dim;
173: pgeom->dimEmbed = dE;
174: //pgeom->isAffine = geom->isAffine;
175: if (geom->isAffine) {
176: if (!p) {
177: pgeom->xi = geom->xi;
178: pgeom->J = &geom->J[c*Np*dE*dE];
179: pgeom->invJ = &geom->invJ[c*Np*dE*dE];
180: pgeom->detJ = &geom->detJ[c*Np];
181: pgeom->n = geom->n ? &geom->n[c*Np*dE] : NULL;
182: }
183: if (pcoords) {CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c*Np*dE], pgeom->J, pcoords, pgeom->v);}
184: } else {
185: pgeom->v = &geom->v[(c*Np+p)*dE];
186: pgeom->J = &geom->J[(c*Np+p)*dE*dE];
187: pgeom->invJ = &geom->invJ[(c*Np+p)*dE*dE];
188: pgeom->detJ = &geom->detJ[c*Np+p];
189: pgeom->n = geom->n ? &geom->n[(c*Np+p)*dE] : NULL;
190: }
191: return(0);
192: }
194: /*@C
195: PetscFEGeomGetCellPoint - Get the cell geometry for face f at point p as a PetscFEGeom
197: Input Parameters:
198: + geom - PetscFEGeom object
199: . f - The face
200: - p - The point
202: Output Parameter:
203: . pgeom - The cell geometry of face f at point p
205: Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom,
206: nothing needs to be done with it afterwards.
208: Level: intermediate
210: .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate()
211: @*/
212: PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom)
213: {
214: const PetscBool bd = geom->dimEmbed > geom->dim && !geom->isHybrid ? 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;
220: pgeom->dim = dim;
221: pgeom->dimEmbed = dE;
222: //pgeom->isAffine = geom->isAffine;
223: if (geom->isAffine) {
224: if (!p) {
225: if (bd) {
226: pgeom->J = &geom->suppJ[0][c*Np*dE*dE];
227: pgeom->invJ = &geom->suppInvJ[0][c*Np*dE*dE];
228: pgeom->detJ = &geom->suppDetJ[0][c*Np];
229: } else {
230: pgeom->J = &geom->J[c*Np*dE*dE];
231: pgeom->invJ = &geom->invJ[c*Np*dE*dE];
232: pgeom->detJ = &geom->detJ[c*Np];
233: }
234: }
235: } else {
236: if (bd) {
237: pgeom->J = &geom->suppJ[0][(c*Np+p)*dE*dE];
238: pgeom->invJ = &geom->suppInvJ[0][(c*Np+p)*dE*dE];
239: pgeom->detJ = &geom->suppDetJ[0][c*Np+p];
240: } else {
241: pgeom->J = &geom->J[(c*Np+p)*dE*dE];
242: pgeom->invJ = &geom->invJ[(c*Np+p)*dE*dE];
243: pgeom->detJ = &geom->detJ[c*Np+p];
244: }
245: }
246: return(0);
247: }
249: /*@
250: PetscFEGeomComplete - Calculate derived quntites from base geometry specification
252: Input Parameter:
253: . geom - PetscFEGeom object
255: Level: intermediate
257: .seealso: PetscFEGeomCreate()
258: @*/
259: PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom)
260: {
261: PetscInt i, j, N, dE;
264: N = geom->numPoints * geom->numCells;
265: dE = geom->dimEmbed;
266: switch (dE) {
267: case 3:
268: for (i = 0; i < N; i++) {
269: DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]);
270: if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]);
271: }
272: break;
273: case 2:
274: for (i = 0; i < N; i++) {
275: DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]);
276: if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]);
277: }
278: break;
279: case 1:
280: for (i = 0; i < N; i++) {
281: geom->detJ[i] = PetscAbsReal(geom->J[i]);
282: if (geom->invJ) geom->invJ[i] = 1. / geom->J[i];
283: }
284: break;
285: }
286: if (geom->n) {
287: for (i = 0; i < N; i++) {
288: for (j = 0; j < dE; j++) {
289: geom->n[dE*i + j] = geom->J[dE*dE*i + dE*j + dE-1] * ((dE == 2) ? -1. : 1.);
290: }
291: }
292: }
293: return(0);
294: }