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;
25: PetscQuadratureGetData(quad,&dim,NULL,&Nq,&p,NULL);
26: PetscNew(&g);
27: g->xi = p;
28: g->numCells = numCells;
29: g->numPoints = Nq;
30: g->dim = dim;
31: g->dimEmbed = dimEmbed;
32: g->isCohesive = PETSC_FALSE;
33: N = numCells * Nq;
34: PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ);
35: if (faceData) {
36: PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n);
37: PetscCall(PetscCalloc6(N * dimEmbed * dimEmbed, &(g->suppJ[0]), N * dimEmbed * dimEmbed, &(g->suppJ[1]),
38: N * dimEmbed * dimEmbed, &(g->suppInvJ[0]), N * dimEmbed * dimEmbed, &(g->suppInvJ[1]),
39: N, &(g->suppDetJ[0]), N, &(g->suppDetJ[1])));
40: }
41: PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ);
42: *geom = g;
43: return 0;
44: }
46: /*@C
47: PetscFEGeomDestroy - Destroy a PetscFEGeom object
49: Input Parameter:
50: . geom - PetscFEGeom object
52: Level: beginner
54: .seealso: PetscFEGeomCreate()
55: @*/
56: PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom)
57: {
58: if (!*geom) return 0;
59: PetscFree3((*geom)->v,(*geom)->J,(*geom)->detJ);
60: PetscFree((*geom)->invJ);
61: PetscFree2((*geom)->face,(*geom)->n);
62: PetscFree6((*geom)->suppJ[0],(*geom)->suppJ[1],(*geom)->suppInvJ[0],(*geom)->suppInvJ[1],(*geom)->suppDetJ[0],(*geom)->suppDetJ[1]);
63: PetscFree(*geom);
64: return 0;
65: }
67: /*@C
68: PetscFEGeomGetChunk - Get a chunk of cells in the group as a PetscFEGeom
70: Input Parameters:
71: + geom - PetscFEGeom object
72: . cStart - The first cell in the chunk
73: - cEnd - The first cell not in the chunk
75: Output Parameter:
76: . chunkGeom - The chunk of cells
78: Level: intermediate
80: .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate()
81: @*/
82: PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
83: {
84: PetscInt Nq;
85: PetscInt dE;
89: if (!(*chunkGeom)) {
90: PetscNew(chunkGeom);
91: }
92: Nq = geom->numPoints;
93: dE= geom->dimEmbed;
94: (*chunkGeom)->dim = geom->dim;
95: (*chunkGeom)->dimEmbed = geom->dimEmbed;
96: (*chunkGeom)->numPoints = geom->numPoints;
97: (*chunkGeom)->numCells = cEnd - cStart;
98: (*chunkGeom)->xi = geom->xi;
99: (*chunkGeom)->v = &geom->v[Nq*dE*cStart];
100: (*chunkGeom)->J = &geom->J[Nq*dE*dE*cStart];
101: (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq*dE*dE*cStart] : NULL;
102: (*chunkGeom)->detJ = &geom->detJ[Nq*cStart];
103: (*chunkGeom)->n = geom->n ? &geom->n[Nq*dE*cStart] : NULL;
104: (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL;
105: (*chunkGeom)->suppJ[0] = geom->suppJ[0] ? &geom->suppJ[0][Nq*dE*dE*cStart] : NULL;
106: (*chunkGeom)->suppJ[1] = geom->suppJ[1] ? &geom->suppJ[1][Nq*dE*dE*cStart] : NULL;
107: (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq*dE*dE*cStart] : NULL;
108: (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq*dE*dE*cStart] : NULL;
109: (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq*cStart] : NULL;
110: (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq*cStart] : NULL;
111: (*chunkGeom)->isAffine = geom->isAffine;
112: return 0;
113: }
115: /*@C
116: PetscFEGeomRestoreChunk - Restore the chunk
118: Input Parameters:
119: + geom - PetscFEGeom object
120: . cStart - The first cell in the chunk
121: . cEnd - The first cell not in the chunk
122: - chunkGeom - The chunk of cells
124: Level: intermediate
126: .seealso: PetscFEGeomGetChunk(), PetscFEGeomCreate()
127: @*/
128: PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
129: {
130: PetscFree(*chunkGeom);
131: return 0;
132: }
134: /*@C
135: PetscFEGeomGetPoint - Get the geometry for cell c at point p as a PetscFEGeom
137: Input Parameters:
138: + geom - PetscFEGeom object
139: . c - The cell
140: . p - The point
141: - pcoords - The reference coordinates of point p, or NULL
143: Output Parameter:
144: . pgeom - The geometry of cell c at point p
146: Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom,
147: nothing needs to be done with it afterwards.
149: In the affine case, pgeom must have storage for the integration point coordinates in pgeom->v if pcoords is passed in.
151: Level: intermediate
153: .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate()
154: @*/
155: PetscErrorCode PetscFEGeomGetPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, const PetscReal pcoords[], PetscFEGeom *pgeom)
156: {
157: const PetscInt dim = geom->dim;
158: const PetscInt dE = geom->dimEmbed;
159: const PetscInt Np = geom->numPoints;
162: pgeom->dim = dim;
163: pgeom->dimEmbed = dE;
164: //pgeom->isAffine = geom->isAffine;
165: if (geom->isAffine) {
166: if (!p) {
167: pgeom->xi = geom->xi;
168: pgeom->J = &geom->J[c*Np*dE*dE];
169: pgeom->invJ = &geom->invJ[c*Np*dE*dE];
170: pgeom->detJ = &geom->detJ[c*Np];
171: pgeom->n = geom->n ? &geom->n[c*Np*dE] : NULL;
172: }
173: if (pcoords) {CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c*Np*dE], pgeom->J, pcoords, pgeom->v);}
174: } else {
175: pgeom->v = &geom->v[(c*Np+p)*dE];
176: pgeom->J = &geom->J[(c*Np+p)*dE*dE];
177: pgeom->invJ = &geom->invJ[(c*Np+p)*dE*dE];
178: pgeom->detJ = &geom->detJ[c*Np+p];
179: pgeom->n = geom->n ? &geom->n[(c*Np+p)*dE] : NULL;
180: }
181: return 0;
182: }
184: /*@C
185: PetscFEGeomGetCellPoint - Get the cell geometry for face f at point p as a PetscFEGeom
187: Input Parameters:
188: + geom - PetscFEGeom object
189: . f - The face
190: - p - The point
192: Output Parameter:
193: . pgeom - The cell geometry of face f at point p
195: Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom,
196: nothing needs to be done with it afterwards.
198: Level: intermediate
200: .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate()
201: @*/
202: PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom)
203: {
204: const PetscBool bd = geom->dimEmbed > geom->dim && !geom->isCohesive ? PETSC_TRUE : PETSC_FALSE;
205: const PetscInt dim = bd ? geom->dimEmbed : geom->dim;
206: const PetscInt dE = geom->dimEmbed;
207: const PetscInt Np = geom->numPoints;
210: pgeom->dim = dim;
211: pgeom->dimEmbed = dE;
212: //pgeom->isAffine = geom->isAffine;
213: if (geom->isAffine) {
214: if (!p) {
215: if (bd) {
216: pgeom->J = &geom->suppJ[0][c*Np*dE*dE];
217: pgeom->invJ = &geom->suppInvJ[0][c*Np*dE*dE];
218: pgeom->detJ = &geom->suppDetJ[0][c*Np];
219: } else {
220: pgeom->J = &geom->J[c*Np*dE*dE];
221: pgeom->invJ = &geom->invJ[c*Np*dE*dE];
222: pgeom->detJ = &geom->detJ[c*Np];
223: }
224: }
225: } else {
226: if (bd) {
227: pgeom->J = &geom->suppJ[0][(c*Np+p)*dE*dE];
228: pgeom->invJ = &geom->suppInvJ[0][(c*Np+p)*dE*dE];
229: pgeom->detJ = &geom->suppDetJ[0][c*Np+p];
230: } else {
231: pgeom->J = &geom->J[(c*Np+p)*dE*dE];
232: pgeom->invJ = &geom->invJ[(c*Np+p)*dE*dE];
233: pgeom->detJ = &geom->detJ[c*Np+p];
234: }
235: }
236: return 0;
237: }
239: /*@
240: PetscFEGeomComplete - Calculate derived quntites from base geometry specification
242: Input Parameter:
243: . geom - PetscFEGeom object
245: Level: intermediate
247: .seealso: PetscFEGeomCreate()
248: @*/
249: PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom)
250: {
251: PetscInt i, j, N, dE;
254: N = geom->numPoints * geom->numCells;
255: dE = geom->dimEmbed;
256: switch (dE) {
257: case 3:
258: for (i = 0; i < N; i++) {
259: DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]);
260: if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]);
261: }
262: break;
263: case 2:
264: for (i = 0; i < N; i++) {
265: DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]);
266: if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]);
267: }
268: break;
269: case 1:
270: for (i = 0; i < N; i++) {
271: geom->detJ[i] = PetscAbsReal(geom->J[i]);
272: if (geom->invJ) geom->invJ[i] = 1. / geom->J[i];
273: }
274: break;
275: }
276: if (geom->n) {
277: for (i = 0; i < N; i++) {
278: for (j = 0; j < dE; j++) {
279: geom->n[dE*i + j] = geom->J[dE*dE*i + dE*j + dE-1] * ((dE == 2) ? -1. : 1.);
280: }
281: }
282: }
283: return 0;
284: }