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 - faceData - Flag to construct geometry data for the faces
11
12 Output Parameter:
13 . geom - The PetscFEGeom object
14
15 Level: beginner
16
17 .seealso: PetscFEGeomDestroy(), PetscFEGeomComplete()
18 @*/
PetscFEGeomCreate(PetscQuadrature quad,PetscInt numCells,PetscInt dimEmbed,PetscBool faceData,PetscFEGeom ** geom)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;
25
26 PetscFunctionBegin;
27 ierr = PetscQuadratureGetData(quad,&dim,NULL,&Nq,&p,NULL);CHKERRQ(ierr);
28 ierr = PetscNew(&g);CHKERRQ(ierr);
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 ierr = PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ);CHKERRQ(ierr);
36 if (faceData) {
37 ierr = PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n);CHKERRQ(ierr);
38 ierr = 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]));CHKERRQ(ierr);
41 }
42 ierr = PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ);CHKERRQ(ierr);
43 *geom = g;
44 PetscFunctionReturn(0);
45 }
46
47 /*@C
48 PetscFEGeomDestroy - Destroy a PetscFEGeom object
49
50 Input Parameter:
51 . geom - PetscFEGeom object
52
53 Level: beginner
54
55 .seealso: PetscFEGeomCreate()
56 @*/
PetscFEGeomDestroy(PetscFEGeom ** geom)57 PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom)
58 {
59 PetscErrorCode ierr;
60
61 PetscFunctionBegin;
62 if (!*geom) PetscFunctionReturn(0);
63 ierr = PetscFree3((*geom)->v,(*geom)->J,(*geom)->detJ);CHKERRQ(ierr);
64 ierr = PetscFree((*geom)->invJ);CHKERRQ(ierr);
65 ierr = PetscFree2((*geom)->face,(*geom)->n);CHKERRQ(ierr);
66 ierr = PetscFree6((*geom)->suppJ[0],(*geom)->suppJ[1],(*geom)->suppInvJ[0],(*geom)->suppInvJ[1],(*geom)->suppDetJ[0],(*geom)->suppDetJ[1]);CHKERRQ(ierr);
67 ierr = PetscFree(*geom);CHKERRQ(ierr);
68 PetscFunctionReturn(0);
69 }
70
71 /*@C
72 PetscFEGeomGetChunk - Get a chunk of cells in the group as a PetscFEGeom
73
74 Input Parameters:
75 + geom - PetscFEGeom object
76 . cStart - The first cell in the chunk
77 - cEnd - The first cell not in the chunk
78
79 Output Parameter:
80 . chunkGeom - The chunk of cells
81
82 Level: intermediate
83
84 .seealso: PetscFEGeomRestoreChunk(), PetscFEGeomCreate()
85 @*/
PetscFEGeomGetChunk(PetscFEGeom * geom,PetscInt cStart,PetscInt cEnd,PetscFEGeom ** chunkGeom)86 PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
87 {
88 PetscInt Nq;
89 PetscInt dE;
90 PetscErrorCode ierr;
91
92 PetscFunctionBegin;
93 PetscValidPointer(geom,1);
94 PetscValidPointer(chunkGeom,2);
95 if (!(*chunkGeom)) {
96 ierr = PetscNew(chunkGeom);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
119 }
120
121 /*@C
122 PetscFEGeomRestoreChunk - Restore the chunk
123
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
129
130 Level: intermediate
131
132 .seealso: PetscFEGeomGetChunk(), PetscFEGeomCreate()
133 @*/
PetscFEGeomRestoreChunk(PetscFEGeom * geom,PetscInt cStart,PetscInt cEnd,PetscFEGeom ** chunkGeom)134 PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom)
135 {
136 PetscErrorCode ierr;
137
138 PetscFunctionBegin;
139 ierr = PetscFree(*chunkGeom);CHKERRQ(ierr);
140 PetscFunctionReturn(0);
141 }
142
143 /*@
144 PetscFEGeomComplete - Calculate derived quntites from base geometry specification
145
146 Input Parameter:
147 . geom - PetscFEGeom object
148
149 Level: intermediate
150
151 .seealso: PetscFEGeomCreate()
152 @*/
PetscFEGeomComplete(PetscFEGeom * geom)153 PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom)
154 {
155 PetscInt i, j, N, dE;
156
157 PetscFunctionBeginHot;
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 PetscFunctionReturn(0);
188 }
189