Actual source code: fegeom.c
petsc-3.14.6 2021-03-30
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: N = numCells * Nq;
35: PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ);
36: if (faceData) {
37: PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n);
38: PetscCalloc6(N * dimEmbed * dimEmbed, &(g->suppJ[0]), N * dimEmbed * dimEmbed, &(g->suppJ[1]),
39: N * dimEmbed * dimEmbed, &(g->suppInvJ[0]), N * dimEmbed * dimEmbed, &(g->suppInvJ[1]),
40: N, &(g->suppDetJ[0]), N, &(g->suppDetJ[1]));
41: }
42: PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ);
43: *geom = g;
44: return(0);
45: }
47: /*@C
48: PetscFEGeomDestroy - Destroy a PetscFEGeom object
50: Input Parameter:
51: . geom - PetscFEGeom object
53: Level: beginner
55: .seealso: PetscFEGeomCreate()
56: @*/
57: PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom)
58: {
62: if (!*geom) return(0);
63: PetscFree3((*geom)->v,(*geom)->J,(*geom)->detJ);
64: PetscFree((*geom)->invJ);
65: PetscFree2((*geom)->face,(*geom)->n);
66: PetscFree6((*geom)->suppJ[0],(*geom)->suppJ[1],(*geom)->suppInvJ[0],(*geom)->suppInvJ[1],(*geom)->suppDetJ[0],(*geom)->suppDetJ[1]);
67: PetscFree(*geom);
68: return(0);
69: }
71: /*@C
72: PetscFEGeomGetChunk - Get a chunk of cells in the group as a PetscFEGeom
74: Input Parameters:
75: + geom - PetscFEGeom object
76: . cStart - The first cell in the chunk
77: - cEnd - The first cell not in the chunk
79: Output Parameter:
80: . chunkGeom - The chunk of cells
82: Level: intermediate
84: .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate()
85: @*/
86: PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
87: {
88: PetscInt Nq;
89: PetscInt dE;
95: if (!(*chunkGeom)) {
96: PetscNew(chunkGeom);
97: }
98: Nq = geom->numPoints;
99: dE= geom->dimEmbed;
100: (*chunkGeom)->dim = geom->dim;
101: (*chunkGeom)->dimEmbed = geom->dimEmbed;
102: (*chunkGeom)->numPoints = geom->numPoints;
103: (*chunkGeom)->numCells = cEnd - cStart;
104: (*chunkGeom)->xi = geom->xi;
105: (*chunkGeom)->v = &geom->v[Nq*dE*cStart];
106: (*chunkGeom)->J = &geom->J[Nq*dE*dE*cStart];
107: (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq*dE*dE*cStart] : NULL;
108: (*chunkGeom)->detJ = &geom->detJ[Nq*cStart];
109: (*chunkGeom)->n = geom->n ? &geom->n[Nq*dE*cStart] : NULL;
110: (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL;
111: (*chunkGeom)->suppJ[0] = geom->suppJ[0] ? &geom->suppJ[0][Nq*dE*dE*cStart] : NULL;
112: (*chunkGeom)->suppJ[1] = geom->suppJ[1] ? &geom->suppJ[1][Nq*dE*dE*cStart] : NULL;
113: (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq*dE*dE*cStart] : NULL;
114: (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq*dE*dE*cStart] : NULL;
115: (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq*cStart] : NULL;
116: (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq*cStart] : NULL;
117: (*chunkGeom)->isAffine = geom->isAffine;
118: return(0);
119: }
121: /*@C
122: PetscFEGeomRestoreChunk - Restore the chunk
124: Input Parameters:
125: + geom - PetscFEGeom object
126: . cStart - The first cell in the chunk
127: . cEnd - The first cell not in the chunk
128: - chunkGeom - The chunk of cells
130: Level: intermediate
132: .seealso: PetscFEGeomGetChunk(), PetscFEGeomCreate()
133: @*/
134: PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
135: {
139: PetscFree(*chunkGeom);
140: return(0);
141: }
143: /*@
144: PetscFEGeomComplete - Calculate derived quntites from base geometry specification
146: Input Parameter:
147: . geom - PetscFEGeom object
149: Level: intermediate
151: .seealso: PetscFEGeomCreate()
152: @*/
153: PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom)
154: {
155: PetscInt i, j, N, dE;
158: N = geom->numPoints * geom->numCells;
159: dE = geom->dimEmbed;
160: switch (dE) {
161: case 3:
162: for (i = 0; i < N; i++) {
163: DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]);
164: if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]);
165: }
166: break;
167: case 2:
168: for (i = 0; i < N; i++) {
169: DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]);
170: if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]);
171: }
172: break;
173: case 1:
174: for (i = 0; i < N; i++) {
175: geom->detJ[i] = PetscAbsReal(geom->J[i]);
176: if (geom->invJ) geom->invJ[i] = 1. / geom->J[i];
177: }
178: break;
179: }
180: if (geom->n) {
181: for (i = 0; i < N; i++) {
182: for (j = 0; j < dE; j++) {
183: geom->n[dE*i + j] = geom->J[dE*dE*i + dE*j + dE-1] * ((dE == 2) ? -1. : 1.);
184: }
185: }
186: }
187: return(0);
188: }