1 static char help[] = "Tests for cell geometry\n\n";
2
3 #include <petscdmplex.h>
4
5 typedef enum {RUN_REFERENCE, RUN_FILE, RUN_DISPLAY} RunType;
6
7 typedef struct {
8 DM dm;
9 RunType runType; /* Type of mesh to use */
10 char filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */
11 PetscBool interpolate; /* Interpolate the mesh */
12 PetscBool transform; /* Use random coordinate transformations */
13 /* Data for input meshes */
14 PetscReal *v0, *J, *invJ, *detJ; /* FEM data */
15 PetscReal *centroid, *normal, *vol; /* FVM data */
16 } AppCtx;
17
ReadMesh(MPI_Comm comm,const char * filename,AppCtx * user,DM * dm)18 static PetscErrorCode ReadMesh(MPI_Comm comm, const char *filename, AppCtx *user, DM *dm)
19 {
20 PetscErrorCode ierr;
21
22 PetscFunctionBegin;
23 ierr = DMPlexCreateFromFile(comm, filename, user->interpolate, dm);CHKERRQ(ierr);
24 ierr = PetscObjectSetName((PetscObject) *dm, "Input Mesh");CHKERRQ(ierr);
25 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
26 PetscFunctionReturn(0);
27 }
28
ProcessOptions(MPI_Comm comm,AppCtx * options)29 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
30 {
31 const char *runTypes[3] = {"reference", "file", "display"};
32 PetscInt run;
33 PetscErrorCode ierr;
34
35 PetscFunctionBeginUser;
36 options->runType = RUN_REFERENCE;
37 options->filename[0] = '\0';
38 options->interpolate = PETSC_FALSE;
39 options->transform = PETSC_FALSE;
40
41 ierr = PetscOptionsBegin(comm, "", "Geometry Test Options", "DMPLEX");CHKERRQ(ierr);
42 run = options->runType;
43 ierr = PetscOptionsEList("-run_type", "The run type", "ex8.c", runTypes, 3, runTypes[options->runType], &run, NULL);CHKERRQ(ierr);
44 options->runType = (RunType) run;
45 ierr = PetscOptionsString("-filename", "The mesh file", "ex8.c", options->filename, options->filename, sizeof(options->filename), NULL);CHKERRQ(ierr);
46 ierr = PetscOptionsBool("-interpolate", "Interpolate the mesh", "ex8.c", options->interpolate, &options->interpolate, NULL);CHKERRQ(ierr);
47 ierr = PetscOptionsBool("-transform", "Use random transforms", "ex8.c", options->transform, &options->transform, NULL);CHKERRQ(ierr);
48
49 if (options->runType == RUN_FILE) {
50 PetscInt dim, cStart, cEnd, numCells, n;
51 PetscBool flag;
52
53 ierr = ReadMesh(PETSC_COMM_WORLD, options->filename, options, &options->dm);CHKERRQ(ierr);
54 ierr = DMGetDimension(options->dm, &dim);CHKERRQ(ierr);
55 ierr = DMPlexGetHeightStratum(options->dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
56 numCells = cEnd-cStart;
57 ierr = PetscMalloc7(numCells*dim,&options->v0,numCells*dim*dim,&options->J,numCells*dim*dim,&options->invJ,numCells,&options->detJ,
58 numCells*dim,&options->centroid,numCells*dim,&options->normal,numCells,&options->vol);CHKERRQ(ierr);
59 n = numCells*dim;
60 ierr = PetscOptionsRealArray("-v0", "Input v0 for each cell", "ex8.c", options->v0, &n, &flag);CHKERRQ(ierr);
61 if (flag && n != numCells*dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of v0 %D should be %D", n, numCells*dim);
62 n = numCells*dim*dim;
63 ierr = PetscOptionsRealArray("-J", "Input Jacobian for each cell", "ex8.c", options->J, &n, &flag);CHKERRQ(ierr);
64 if (flag && n != numCells*dim*dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of J %D should be %D", n, numCells*dim*dim);
65 n = numCells*dim*dim;
66 ierr = PetscOptionsRealArray("-invJ", "Input inverse Jacobian for each cell", "ex8.c", options->invJ, &n, &flag);CHKERRQ(ierr);
67 if (flag && n != numCells*dim*dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of invJ %D should be %D", n, numCells*dim*dim);
68 n = numCells;
69 ierr = PetscOptionsRealArray("-detJ", "Input Jacobian determinant for each cell", "ex8.c", options->detJ, &n, &flag);CHKERRQ(ierr);
70 if (flag && n != numCells) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of detJ %D should be %D", n, numCells);
71 n = numCells*dim;
72 ierr = PetscOptionsRealArray("-centroid", "Input centroid for each cell", "ex8.c", options->centroid, &n, &flag);CHKERRQ(ierr);
73 if (flag && n != numCells*dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of centroid %D should be %D", n, numCells*dim);
74 n = numCells*dim;
75 ierr = PetscOptionsRealArray("-normal", "Input normal for each cell", "ex8.c", options->normal, &n, &flag);CHKERRQ(ierr);
76 if (flag && n != numCells*dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of normal %D should be %D", n, numCells*dim);
77 n = numCells;
78 ierr = PetscOptionsRealArray("-vol", "Input volume for each cell", "ex8.c", options->vol, &n, &flag);CHKERRQ(ierr);
79 if (flag && n != numCells) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of vol %D should be %D", n, numCells);
80 } else if (options->runType == RUN_DISPLAY) {
81 ierr = ReadMesh(PETSC_COMM_WORLD, options->filename, options, &options->dm);CHKERRQ(ierr);
82 }
83 ierr = PetscOptionsEnd();
84
85 if (options->transform) {ierr = PetscPrintf(comm, "Using random transforms\n");CHKERRQ(ierr);}
86 PetscFunctionReturn(0);
87 }
88
ChangeCoordinates(DM dm,PetscInt spaceDim,PetscScalar vertexCoords[])89 static PetscErrorCode ChangeCoordinates(DM dm, PetscInt spaceDim, PetscScalar vertexCoords[])
90 {
91 PetscSection coordSection;
92 Vec coordinates;
93 PetscScalar *coords;
94 PetscInt vStart, vEnd, v, d, coordSize;
95 PetscErrorCode ierr;
96
97 PetscFunctionBegin;
98 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
99 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
100 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
101 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr);
102 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr);
103 for (v = vStart; v < vEnd; ++v) {
104 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr);
105 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr);
106 }
107 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
108 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
109 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
110 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
111 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
112 ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr);
113 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
114 for (v = vStart; v < vEnd; ++v) {
115 PetscInt off;
116
117 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
118 for (d = 0; d < spaceDim; ++d) {
119 coords[off+d] = vertexCoords[(v-vStart)*spaceDim+d];
120 }
121 }
122 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
123 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr);
124 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
125 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
126 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
127 PetscFunctionReturn(0);
128 }
129
130 #define RelativeError(a,b) PetscAbs(a-b)/(1.0+PetscMax(PetscAbs(a),PetscAbs(b)))
131
CheckFEMGeometry(DM dm,PetscInt cell,PetscInt spaceDim,PetscReal v0Ex[],PetscReal JEx[],PetscReal invJEx[],PetscReal detJEx)132 static PetscErrorCode CheckFEMGeometry(DM dm, PetscInt cell, PetscInt spaceDim, PetscReal v0Ex[], PetscReal JEx[], PetscReal invJEx[], PetscReal detJEx)
133 {
134 PetscReal v0[3], J[9], invJ[9], detJ;
135 PetscInt d, i, j;
136 PetscErrorCode ierr;
137
138 PetscFunctionBegin;
139 ierr = DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
140 for (d = 0; d < spaceDim; ++d) {
141 if (v0[d] != v0Ex[d]) {
142 switch (spaceDim) {
143 case 2: SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0 (%g, %g) != (%g, %g)", (double)v0[0], (double)v0[1], (double)v0Ex[0], (double)v0Ex[1]);break;
144 case 3: SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0 (%g, %g, %g) != (%g, %g, %g)", (double)v0[0], (double)v0[1], (double)v0[2], (double)v0Ex[0], (double)v0Ex[1], (double)v0Ex[2]);break;
145 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid space dimension %D", spaceDim);
146 }
147 }
148 }
149 for (i = 0; i < spaceDim; ++i) {
150 for (j = 0; j < spaceDim; ++j) {
151 if (RelativeError(J[i*spaceDim+j],JEx[i*spaceDim+j]) > 10*PETSC_SMALL) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid J[%D,%D]: %g != %g", i, j, (double)J[i*spaceDim+j], (double)JEx[i*spaceDim+j]);
152 if (RelativeError(invJ[i*spaceDim+j],invJEx[i*spaceDim+j]) > 10*PETSC_SMALL) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid invJ[%D,%D]: %g != %g", i, j, (double)invJ[i*spaceDim+j], (double)invJEx[i*spaceDim+j]);
153 }
154 }
155 if (RelativeError(detJ,detJEx) > 10*PETSC_SMALL) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid |J| = %g != %g diff %g", (double)detJ, (double)detJEx,(double)(detJ - detJEx));
156 PetscFunctionReturn(0);
157 }
158
CheckFVMGeometry(DM dm,PetscInt cell,PetscInt spaceDim,PetscReal centroidEx[],PetscReal normalEx[],PetscReal volEx)159 static PetscErrorCode CheckFVMGeometry(DM dm, PetscInt cell, PetscInt spaceDim, PetscReal centroidEx[], PetscReal normalEx[], PetscReal volEx)
160 {
161 PetscReal centroid[3], normal[3], vol;
162 PetscInt d;
163 PetscErrorCode ierr;
164
165 PetscFunctionBegin;
166 ierr = DMPlexComputeCellGeometryFVM(dm, cell, &vol, centroid, normal);CHKERRQ(ierr);
167 for (d = 0; d < spaceDim; ++d) {
168 if (RelativeError(centroid[d],centroidEx[d]) > 10*PETSC_SMALL) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %D, Invalid centroid[%D]: %g != %g diff %g", cell, d, (double)centroid[d], (double)centroidEx[d],(double)(centroid[d]-centroidEx[d]));
169 if (RelativeError(normal[d],normalEx[d]) > 10*PETSC_SMALL) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %D, Invalid normal[%D]: %g != %g", cell, d, (double)normal[d], (double)normalEx[d]);
170 }
171 if (RelativeError(volEx,vol) > 10*PETSC_SMALL) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %D, Invalid volume = %g != %g diff %g", cell, (double)vol, (double)volEx,(double)(vol - volEx));
172 PetscFunctionReturn(0);
173 }
174
CheckCell(DM dm,PetscInt cell,PetscBool transform,PetscReal v0Ex[],PetscReal JEx[],PetscReal invJEx[],PetscReal detJEx,PetscReal centroidEx[],PetscReal normalEx[],PetscReal volEx,PetscReal faceCentroidEx[],PetscReal faceNormalEx[],PetscReal faceVolEx[])175 static PetscErrorCode CheckCell(DM dm, PetscInt cell, PetscBool transform, PetscReal v0Ex[], PetscReal JEx[], PetscReal invJEx[], PetscReal detJEx, PetscReal centroidEx[], PetscReal normalEx[], PetscReal volEx, PetscReal faceCentroidEx[], PetscReal faceNormalEx[], PetscReal faceVolEx[])
176 {
177 const PetscInt *cone;
178 PetscInt coneSize, c;
179 PetscInt dim, depth, cdim;
180 PetscErrorCode ierr;
181
182 PetscFunctionBegin;
183 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
184 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
185 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr);
186 if (v0Ex) {
187 ierr = CheckFEMGeometry(dm, cell, cdim, v0Ex, JEx, invJEx, detJEx);CHKERRQ(ierr);
188 }
189 if (dim == depth && centroidEx) {
190 ierr = CheckFVMGeometry(dm, cell, cdim, centroidEx, normalEx, volEx);CHKERRQ(ierr);
191 if (faceCentroidEx) {
192 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr);
193 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr);
194 for (c = 0; c < coneSize; ++c) {
195 ierr = CheckFVMGeometry(dm, cone[c], dim, &faceCentroidEx[c*dim], &faceNormalEx[c*dim], faceVolEx[c]);CHKERRQ(ierr);
196 }
197 }
198 }
199 if (transform) {
200 Vec coordinates;
201 PetscSection coordSection;
202 PetscScalar *coords = NULL, *origCoords, *newCoords;
203 PetscReal *v0ExT, *JExT, *invJExT, detJExT=0, *centroidExT, *normalExT, volExT=0;
204 PetscReal *faceCentroidExT, *faceNormalExT, faceVolExT;
205 PetscRandom r, ang, ang2;
206 PetscInt coordSize, numCorners, t;
207
208 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
209 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
210 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr);
211 ierr = PetscMalloc2(coordSize, &origCoords, coordSize, &newCoords);CHKERRQ(ierr);
212 ierr = PetscMalloc5(cdim, &v0ExT, cdim*cdim, &JExT, cdim*cdim, &invJExT, cdim, ¢roidExT, cdim, &normalExT);CHKERRQ(ierr);
213 ierr = PetscMalloc2(cdim, &faceCentroidExT, cdim, &faceNormalExT);CHKERRQ(ierr);
214 for (c = 0; c < coordSize; ++c) origCoords[c] = coords[c];
215 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);CHKERRQ(ierr);
216 numCorners = coordSize/cdim;
217
218 ierr = PetscRandomCreate(PETSC_COMM_SELF, &r);CHKERRQ(ierr);
219 ierr = PetscRandomSetFromOptions(r);CHKERRQ(ierr);
220 ierr = PetscRandomSetInterval(r, 0.0, 10.0);CHKERRQ(ierr);
221 ierr = PetscRandomCreate(PETSC_COMM_SELF, &ang);CHKERRQ(ierr);
222 ierr = PetscRandomSetFromOptions(ang);CHKERRQ(ierr);
223 ierr = PetscRandomSetInterval(ang, 0.0, 2*PETSC_PI);CHKERRQ(ierr);
224 ierr = PetscRandomCreate(PETSC_COMM_SELF, &ang2);CHKERRQ(ierr);
225 ierr = PetscRandomSetFromOptions(ang2);CHKERRQ(ierr);
226 ierr = PetscRandomSetInterval(ang2, 0.0, PETSC_PI);CHKERRQ(ierr);
227 for (t = 0; t < 1; ++t) {
228 PetscScalar trans[3];
229 PetscReal R[9], rot[3], rotM[9];
230 PetscReal scale, phi, theta, psi = 0.0, norm;
231 PetscInt d, e, f, p;
232
233 for (c = 0; c < coordSize; ++c) newCoords[c] = origCoords[c];
234 ierr = PetscRandomGetValueReal(r, &scale);CHKERRQ(ierr);
235 ierr = PetscRandomGetValueReal(ang, &phi);CHKERRQ(ierr);
236 ierr = PetscRandomGetValueReal(ang2, &theta);CHKERRQ(ierr);
237 for (d = 0; d < cdim; ++d) {
238 ierr = PetscRandomGetValue(r, &trans[d]);CHKERRQ(ierr);
239 }
240 switch (cdim) {
241 case 2:
242 R[0] = PetscCosReal(phi); R[1] = -PetscSinReal(phi);
243 R[2] = PetscSinReal(phi); R[3] = PetscCosReal(phi);
244 break;
245 case 3:
246 {
247 const PetscReal ct = PetscCosReal(theta), st = PetscSinReal(theta);
248 const PetscReal cp = PetscCosReal(phi), sp = PetscSinReal(phi);
249 const PetscReal cs = PetscCosReal(psi), ss = PetscSinReal(psi);
250 R[0] = ct*cs; R[1] = sp*st*cs - cp*ss; R[2] = sp*ss + cp*st*cs;
251 R[3] = ct*ss; R[4] = cp*cs + sp*st*ss; R[5] = cp*st*ss - sp*cs;
252 R[6] = -st; R[7] = sp*ct; R[8] = cp*ct;
253 break;
254 }
255 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid coordinate dimension %D", cdim);
256 }
257 if (v0Ex) {
258 detJExT = detJEx;
259 for (d = 0; d < cdim; ++d) {
260 v0ExT[d] = v0Ex[d];
261 for (e = 0; e < cdim; ++e) {
262 JExT[d*cdim+e] = JEx[d*cdim+e];
263 invJExT[d*cdim+e] = invJEx[d*cdim+e];
264 }
265 }
266 for (d = 0; d < cdim; ++d) {
267 v0ExT[d] *= scale;
268 v0ExT[d] += PetscRealPart(trans[d]);
269 /* Only scale dimensions in the manifold */
270 for (e = 0; e < dim; ++e) {
271 JExT[d*cdim+e] *= scale;
272 invJExT[d*cdim+e] /= scale;
273 }
274 if (d < dim) detJExT *= scale;
275 }
276 /* Do scaling and translation before rotation, so that we can leave out the normal dimension for lower dimensional manifolds */
277 for (d = 0; d < cdim; ++d) {
278 for (e = 0, rot[d] = 0.0; e < cdim; ++e) {
279 rot[d] += R[d*cdim+e] * v0ExT[e];
280 }
281 }
282 for (d = 0; d < cdim; ++d) v0ExT[d] = rot[d];
283 for (d = 0; d < cdim; ++d) {
284 for (e = 0; e < cdim; ++e) {
285 for (f = 0, rotM[d*cdim+e] = 0.0; f < cdim; ++f) {
286 rotM[d*cdim+e] += R[d*cdim+f] * JExT[f*cdim+e];
287 }
288 }
289 }
290 for (d = 0; d < cdim; ++d) {
291 for (e = 0; e < cdim; ++e) {
292 JExT[d*cdim+e] = rotM[d*cdim+e];
293 }
294 }
295 for (d = 0; d < cdim; ++d) {
296 for (e = 0; e < cdim; ++e) {
297 for (f = 0, rotM[d*cdim+e] = 0.0; f < cdim; ++f) {
298 rotM[d*cdim+e] += invJExT[d*cdim+f] * R[e*cdim+f];
299 }
300 }
301 }
302 for (d = 0; d < cdim; ++d) {
303 for (e = 0; e < cdim; ++e) {
304 invJExT[d*cdim+e] = rotM[d*cdim+e];
305 }
306 }
307 }
308 if (centroidEx) {
309 volExT = volEx;
310 for (d = 0; d < cdim; ++d) {
311 centroidExT[d] = centroidEx[d];
312 normalExT[d] = normalEx[d];
313 }
314 for (d = 0; d < cdim; ++d) {
315 centroidExT[d] *= scale;
316 centroidExT[d] += PetscRealPart(trans[d]);
317 normalExT[d] /= scale;
318 /* Only scale dimensions in the manifold */
319 if (d < dim) volExT *= scale;
320 }
321 /* Do scaling and translation before rotation, so that we can leave out the normal dimension for lower dimensional manifolds */
322 for (d = 0; d < cdim; ++d) {
323 for (e = 0, rot[d] = 0.0; e < cdim; ++e) {
324 rot[d] += R[d*cdim+e] * centroidExT[e];
325 }
326 }
327 for (d = 0; d < cdim; ++d) centroidExT[d] = rot[d];
328 for (d = 0; d < cdim; ++d) {
329 for (e = 0, rot[d] = 0.0; e < cdim; ++e) {
330 rot[d] += R[d*cdim+e] * normalExT[e];
331 }
332 }
333 for (d = 0; d < cdim; ++d) normalExT[d] = rot[d];
334 for (d = 0, norm = 0.0; d < cdim; ++d) norm += PetscSqr(normalExT[d]);
335 norm = PetscSqrtReal(norm);
336 for (d = 0; d < cdim; ++d) normalExT[d] /= norm;
337 }
338 for (d = 0; d < cdim; ++d) {
339 for (p = 0; p < numCorners; ++p) {
340 newCoords[p*cdim+d] *= scale;
341 newCoords[p*cdim+d] += trans[d];
342 }
343 }
344 for (p = 0; p < numCorners; ++p) {
345 for (d = 0; d < cdim; ++d) {
346 for (e = 0, rot[d] = 0.0; e < cdim; ++e) {
347 rot[d] += R[d*cdim+e] * PetscRealPart(newCoords[p*cdim+e]);
348 }
349 }
350 for (d = 0; d < cdim; ++d) newCoords[p*cdim+d] = rot[d];
351 }
352
353 ierr = ChangeCoordinates(dm, cdim, newCoords);CHKERRQ(ierr);
354 if (v0Ex) {
355 ierr = CheckFEMGeometry(dm, 0, cdim, v0ExT, JExT, invJExT, detJExT);CHKERRQ(ierr);
356 }
357 if (dim == depth && centroidEx) {
358 ierr = CheckFVMGeometry(dm, cell, cdim, centroidExT, normalExT, volExT);CHKERRQ(ierr);
359 if (faceCentroidEx) {
360 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr);
361 ierr = DMPlexGetCone(dm, cell, &cone);CHKERRQ(ierr);
362 for (c = 0; c < coneSize; ++c) {
363 PetscInt off = c*cdim;
364
365 faceVolExT = faceVolEx[c];
366 for (d = 0; d < cdim; ++d) {
367 faceCentroidExT[d] = faceCentroidEx[off+d];
368 faceNormalExT[d] = faceNormalEx[off+d];
369 }
370 for (d = 0; d < cdim; ++d) {
371 faceCentroidExT[d] *= scale;
372 faceCentroidExT[d] += PetscRealPart(trans[d]);
373 faceNormalExT[d] /= scale;
374 /* Only scale dimensions in the manifold */
375 if (d < dim-1) {
376 faceVolExT *= scale;
377 }
378 }
379 for (d = 0; d < cdim; ++d) {
380 for (e = 0, rot[d] = 0.0; e < cdim; ++e) {
381 rot[d] += R[d*cdim+e] * faceCentroidExT[e];
382 }
383 }
384 for (d = 0; d < cdim; ++d) faceCentroidExT[d] = rot[d];
385 for (d = 0; d < cdim; ++d) {
386 for (e = 0, rot[d] = 0.0; e < cdim; ++e) {
387 rot[d] += R[d*cdim+e] * faceNormalExT[e];
388 }
389 }
390 for (d = 0; d < cdim; ++d) faceNormalExT[d] = rot[d];
391 for (d = 0, norm = 0.0; d < cdim; ++d) norm += PetscSqr(faceNormalExT[d]);
392 norm = PetscSqrtReal(norm);
393 for (d = 0; d < cdim; ++d) faceNormalExT[d] /= norm;
394
395 ierr = CheckFVMGeometry(dm, cone[c], cdim, faceCentroidExT, faceNormalExT, faceVolExT);CHKERRQ(ierr);
396 }
397 }
398 }
399 }
400 ierr = PetscRandomDestroy(&r);CHKERRQ(ierr);
401 ierr = PetscRandomDestroy(&ang);CHKERRQ(ierr);
402 ierr = PetscRandomDestroy(&ang2);CHKERRQ(ierr);
403 ierr = PetscFree2(origCoords, newCoords);CHKERRQ(ierr);
404 ierr = PetscFree5(v0ExT, JExT, invJExT, centroidExT, normalExT);CHKERRQ(ierr);
405 ierr = PetscFree2(faceCentroidExT, faceNormalExT);CHKERRQ(ierr);
406 }
407 PetscFunctionReturn(0);
408 }
409
TestTriangle(MPI_Comm comm,PetscBool interpolate,PetscBool transform)410 static PetscErrorCode TestTriangle(MPI_Comm comm, PetscBool interpolate, PetscBool transform)
411 {
412 DM dm;
413 PetscInt dim;
414 PetscErrorCode ierr;
415
416 PetscFunctionBegin;
417 /* Create reference triangle */
418 dim = 2;
419 ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
420 ierr = PetscObjectSetName((PetscObject) dm, "triangle");CHKERRQ(ierr);
421 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
422 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr);
423 {
424 PetscInt numPoints[2] = {3, 1};
425 PetscInt coneSize[4] = {3, 0, 0, 0};
426 PetscInt cones[3] = {1, 2, 3};
427 PetscInt coneOrientations[3] = {0, 0, 0};
428 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0};
429
430 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
431 if (interpolate) {
432 DM idm;
433
434 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
435 ierr = DMDestroy(&dm);CHKERRQ(ierr);
436 dm = idm;
437 }
438 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
439 }
440 /* Check reference geometry: determinant is scaled by reference volume (2.0) */
441 {
442 PetscReal v0Ex[2] = {-1.0, -1.0};
443 PetscReal JEx[4] = {1.0, 0.0, 0.0, 1.0};
444 PetscReal invJEx[4] = {1.0, 0.0, 0.0, 1.0};
445 PetscReal detJEx = 1.0;
446 PetscReal centroidEx[2] = {-((PetscReal)1.)/((PetscReal)3.), -((PetscReal)1.)/((PetscReal)3.)};
447 PetscReal normalEx[2] = {0.0, 0.0};
448 PetscReal volEx = 2.0;
449
450 ierr = CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL);CHKERRQ(ierr);
451 }
452 /* Move to 3D: Check reference geometry: determinant is scaled by reference volume (2.0) */
453 dim = 3;
454 {
455 PetscScalar vertexCoords[9] = {-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, -1.0, 1.0, 0.0};
456 PetscReal v0Ex[3] = {-1.0, -1.0, 0.0};
457 PetscReal JEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
458 PetscReal invJEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
459 PetscReal detJEx = 1.0;
460 PetscReal centroidEx[3] = {-((PetscReal)1.)/((PetscReal)3.), -((PetscReal)1.)/((PetscReal)3.), 0.0};
461 PetscReal normalEx[3] = {0.0, 0.0, 1.0};
462 PetscReal volEx = 2.0;
463
464 ierr = ChangeCoordinates(dm, dim, vertexCoords);CHKERRQ(ierr);
465 ierr = CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL);CHKERRQ(ierr);
466 }
467 /* Cleanup */
468 ierr = DMDestroy(&dm);CHKERRQ(ierr);
469 PetscFunctionReturn(0);
470 }
471
TestQuadrilateral(MPI_Comm comm,PetscBool interpolate,PetscBool transform)472 static PetscErrorCode TestQuadrilateral(MPI_Comm comm, PetscBool interpolate, PetscBool transform)
473 {
474 DM dm;
475 PetscInt dim;
476 PetscErrorCode ierr;
477
478 PetscFunctionBegin;
479 /* Create reference quadrilateral */
480 dim = 2;
481 ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
482 ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);
483 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
484 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr);
485 {
486 PetscInt numPoints[2] = {4, 1};
487 PetscInt coneSize[5] = {4, 0, 0, 0, 0};
488 PetscInt cones[4] = {1, 2, 3, 4};
489 PetscInt coneOrientations[4] = {0, 0, 0, 0};
490 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0};
491
492 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
493 if (interpolate) {
494 DM idm;
495
496 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
497 ierr = DMDestroy(&dm);CHKERRQ(ierr);
498 dm = idm;
499 }
500 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
501 }
502 /* Check reference geometry: determinant is scaled by reference volume (2.0) */
503 {
504 PetscReal v0Ex[2] = {-1.0, -1.0};
505 PetscReal JEx[4] = {1.0, 0.0, 0.0, 1.0};
506 PetscReal invJEx[4] = {1.0, 0.0, 0.0, 1.0};
507 PetscReal detJEx = 1.0;
508 PetscReal centroidEx[2] = {0.0, 0.0};
509 PetscReal normalEx[2] = {0.0, 0.0};
510 PetscReal volEx = 4.0;
511
512 ierr = CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL);CHKERRQ(ierr);
513 }
514 /* Move to 3D: Check reference geometry: determinant is scaled by reference volume (4.0) */
515 dim = 3;
516 {
517 PetscScalar vertexCoords[12] = {-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 1.0, 1.0, 0.0, -1.0, 1.0, 0.0};
518 PetscReal v0Ex[3] = {-1.0, -1.0, 0.0};
519 PetscReal JEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
520 PetscReal invJEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
521 PetscReal detJEx = 1.0;
522 PetscReal centroidEx[3] = {0.0, 0.0, 0.0};
523 PetscReal normalEx[3] = {0.0, 0.0, 1.0};
524 PetscReal volEx = 4.0;
525
526 ierr = ChangeCoordinates(dm, dim, vertexCoords);CHKERRQ(ierr);
527 ierr = CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL);CHKERRQ(ierr);
528 }
529 /* Cleanup */
530 ierr = DMDestroy(&dm);CHKERRQ(ierr);
531 PetscFunctionReturn(0);
532 }
533
TestTetrahedron(MPI_Comm comm,PetscBool interpolate,PetscBool transform)534 static PetscErrorCode TestTetrahedron(MPI_Comm comm, PetscBool interpolate, PetscBool transform)
535 {
536 DM dm;
537 PetscInt dim;
538 PetscErrorCode ierr;
539
540 PetscFunctionBegin;
541 /* Create reference tetrahedron */
542 dim = 3;
543 ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
544 ierr = PetscObjectSetName((PetscObject) dm, "tetrahedron");CHKERRQ(ierr);
545 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
546 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr);
547 {
548 PetscInt numPoints[2] = {4, 1};
549 PetscInt coneSize[5] = {4, 0, 0, 0, 0};
550 PetscInt cones[4] = {1, 2, 3, 4};
551 PetscInt coneOrientations[4] = {0, 0, 0, 0};
552 PetscScalar vertexCoords[12] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0};
553
554 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
555 if (interpolate) {
556 DM idm;
557
558 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
559 ierr = DMDestroy(&dm);CHKERRQ(ierr);
560 dm = idm;
561 }
562 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
563 }
564 /* Check reference geometry: determinant is scaled by reference volume (4/3) */
565 {
566 PetscReal v0Ex[3] = {-1.0, -1.0, -1.0};
567 PetscReal JEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
568 PetscReal invJEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
569 PetscReal detJEx = 1.0;
570 PetscReal centroidEx[3] = {-0.5, -0.5, -0.5};
571 PetscReal normalEx[3] = {0.0, 0.0, 0.0};
572 PetscReal volEx = (PetscReal)4.0/(PetscReal)3.0;
573
574 ierr = CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL);CHKERRQ(ierr);
575 }
576 /* Cleanup */
577 ierr = DMDestroy(&dm);CHKERRQ(ierr);
578 PetscFunctionReturn(0);
579 }
580
TestHexahedron(MPI_Comm comm,PetscBool interpolate,PetscBool transform)581 static PetscErrorCode TestHexahedron(MPI_Comm comm, PetscBool interpolate, PetscBool transform)
582 {
583 DM dm;
584 PetscInt dim;
585 PetscErrorCode ierr;
586
587 PetscFunctionBegin;
588 /* Create reference hexahedron */
589 dim = 3;
590 ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
591 ierr = PetscObjectSetName((PetscObject) dm, "hexahedron");CHKERRQ(ierr);
592 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
593 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr);
594 {
595 PetscInt numPoints[2] = {8, 1};
596 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0};
597 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8};
598 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
599 PetscScalar vertexCoords[24] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0,
600 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0};
601
602 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
603 if (interpolate) {
604 DM idm;
605
606 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
607 ierr = DMDestroy(&dm);CHKERRQ(ierr);
608 dm = idm;
609 }
610 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
611 }
612 /* Check reference geometry: determinant is scaled by reference volume 8.0 */
613 {
614 PetscReal v0Ex[3] = {-1.0, -1.0, -1.0};
615 PetscReal JEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
616 PetscReal invJEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
617 PetscReal detJEx = 1.0;
618 PetscReal centroidEx[3] = {0.0, 0.0, 0.0};
619 PetscReal normalEx[3] = {0.0, 0.0, 0.0};
620 PetscReal volEx = 8.0;
621
622 ierr = CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL);CHKERRQ(ierr);
623 }
624 /* Cleanup */
625 ierr = DMDestroy(&dm);CHKERRQ(ierr);
626 PetscFunctionReturn(0);
627 }
628
629 /* This wedge is a tensor product cell, rather than a normal wedge */
TestWedge(MPI_Comm comm,PetscBool interpolate,PetscBool transform)630 static PetscErrorCode TestWedge(MPI_Comm comm, PetscBool interpolate, PetscBool transform)
631 {
632 DM dm;
633 PetscInt dim;
634 PetscErrorCode ierr;
635
636 PetscFunctionBegin;
637 /* Create reference wedge */
638 dim = 3;
639 ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
640 ierr = PetscObjectSetName((PetscObject) dm, "Triangular Prism");CHKERRQ(ierr);
641 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
642 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr);
643 {
644 PetscInt numPoints[2] = {6, 1};
645 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0};
646 PetscInt cones[6] = {1, 2, 3, 4, 5, 6};
647 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0};
648 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0,
649 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0};
650
651 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
652 if (interpolate) {
653 DM idm;
654
655 ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
656 ierr = DMDestroy(&dm);CHKERRQ(ierr);
657 dm = idm;
658 }
659 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
660 }
661 /* Check reference geometry: determinant is scaled by reference volume 4.0 */
662 {
663 #if 0
664 /* FEM geometry is not functional for wedges */
665 PetscReal v0Ex[3] = {-1.0, -1.0, -1.0};
666 PetscReal JEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
667 PetscReal invJEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
668 PetscReal detJEx = 1.0;
669 #endif
670
671 if (interpolate) {
672 PetscReal centroidEx[3] = {-((PetscReal)1.)/((PetscReal)3.), -((PetscReal)1.)/((PetscReal)3.), 0.0};
673 PetscReal normalEx[3] = {0.0, 0.0, 0.0};
674 PetscReal volEx = 4.0;
675 PetscReal faceVolEx[5] = {2.0, 2.0, 4.0, PETSC_SQRT2*4.0, 4.0};
676 PetscReal faceNormalEx[15] = {0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, PETSC_SQRT2/2.0, PETSC_SQRT2/2.0, 0.0, -1.0, 0.0, 0.0};
677 PetscReal faceCentroidEx[15] = {-((PetscReal)1.)/((PetscReal)3.), -((PetscReal)1.)/((PetscReal)3.), -1.0,
678 -((PetscReal)1.)/((PetscReal)3.), -((PetscReal)1.)/((PetscReal)3.), 1.0,
679 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0};
680
681 ierr = CheckCell(dm, 0, transform, NULL, NULL, NULL, 0.0, centroidEx, normalEx, volEx, faceCentroidEx, faceNormalEx, faceVolEx);CHKERRQ(ierr);
682 }
683 }
684 /* Cleanup */
685 ierr = DMDestroy(&dm);CHKERRQ(ierr);
686 PetscFunctionReturn(0);
687 }
688
main(int argc,char ** argv)689 int main(int argc, char **argv)
690 {
691 AppCtx user;
692 PetscErrorCode ierr;
693
694 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
695 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
696 if (user.runType == RUN_REFERENCE) {
697 ierr = TestTriangle(PETSC_COMM_SELF, user.interpolate, user.transform);CHKERRQ(ierr);
698 ierr = TestQuadrilateral(PETSC_COMM_SELF, user.interpolate, user.transform);CHKERRQ(ierr);
699 ierr = TestTetrahedron(PETSC_COMM_SELF, user.interpolate, user.transform);CHKERRQ(ierr);
700 ierr = TestHexahedron(PETSC_COMM_SELF, user.interpolate, user.transform);CHKERRQ(ierr);
701 ierr = TestWedge(PETSC_COMM_SELF, user.interpolate, user.transform);CHKERRQ(ierr);
702 } else if (user.runType == RUN_FILE) {
703 PetscInt dim, cStart, cEnd, c;
704
705 ierr = DMGetDimension(user.dm, &dim);CHKERRQ(ierr);
706 ierr = DMPlexGetHeightStratum(user.dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
707 for (c = 0; c < cEnd-cStart; ++c) {
708 ierr = CheckCell(user.dm, c+cStart, PETSC_FALSE, &user.v0[c*dim], &user.J[c*dim*dim], &user.invJ[c*dim*dim], user.detJ[c], &user.centroid[c*dim], &user.normal[c*dim], user.vol[c], NULL, NULL, NULL);CHKERRQ(ierr);
709 }
710 ierr = PetscFree7(user.v0,user.J,user.invJ,user.detJ,user.centroid,user.normal,user.vol);CHKERRQ(ierr);
711 ierr = DMDestroy(&user.dm);CHKERRQ(ierr);
712 } else if (user.runType == RUN_DISPLAY) {
713 DM gdm, dmCell;
714 Vec cellgeom, facegeom;
715 const PetscScalar *cgeom;
716 PetscInt dim, d, cStart, cEnd, cEndInterior, c;
717
718 ierr = DMGetCoordinateDim(user.dm, &dim);CHKERRQ(ierr);
719 ierr = DMPlexConstructGhostCells(user.dm, NULL, NULL, &gdm);CHKERRQ(ierr);
720 if (gdm) {
721 ierr = DMDestroy(&user.dm);CHKERRQ(ierr);
722 user.dm = gdm;
723 }
724 ierr = DMPlexComputeGeometryFVM(user.dm, &cellgeom, &facegeom);CHKERRQ(ierr);
725 ierr = DMPlexGetHeightStratum(user.dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
726 ierr = DMPlexGetGhostCellStratum(user.dm, &cEndInterior, NULL);CHKERRQ(ierr);
727 if (cEndInterior >= 0) cEnd = cEndInterior;
728 ierr = VecGetDM(cellgeom, &dmCell);CHKERRQ(ierr);
729 ierr = VecGetArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
730 for (c = 0; c < cEnd-cStart; ++c) {
731 PetscFVCellGeom *cg;
732
733 ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr);
734 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell %4D: Centroid (", c);CHKERRQ(ierr);
735 for (d = 0; d < dim; ++d) {
736 if (d > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);}
737 ierr = PetscPrintf(PETSC_COMM_SELF, "%12.2g", cg->centroid[d]);CHKERRQ(ierr);
738 }
739 ierr = PetscPrintf(PETSC_COMM_SELF, ") Vol %12.2g\n", cg->volume);CHKERRQ(ierr);
740 }
741 ierr = VecRestoreArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
742 ierr = VecDestroy(&cellgeom);CHKERRQ(ierr);
743 ierr = VecDestroy(&facegeom);CHKERRQ(ierr);
744 ierr = DMDestroy(&user.dm);CHKERRQ(ierr);
745 }
746 ierr = PetscFinalize();
747 return ierr;
748 }
749
750 /*TEST
751
752 test:
753 suffix: 0
754 args: -dm_view ascii::ascii_info_detail
755 test:
756 suffix: 1
757 args: -interpolate -dm_view ascii::ascii_info_detail
758 test:
759 suffix: 2
760 args: -transform
761 test:
762 suffix: 3
763 args: -interpolate -transform
764 test:
765 suffix: 4
766 requires: exodusii
767 args: -run_type file -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/simpleblock-100.exo -dm_view ascii::ascii_info_detail -v0 -1.5,-0.5,0.5,-0.5,-0.5,0.5,0.5,-0.5,0.5 -J 0.0,0.0,0.5,0.0,0.5,0.0,-0.5,0.0,0.0,0.0,0.0,0.5,0.0,0.5,0.0,-0.5,0.0,0.0,0.0,0.0,0.5,0.0,0.5,0.0,-0.5,0.0,0.0 -invJ 0.0,0.0,-2.0,0.0,2.0,0.0,2.0,0.0,0.0,0.0,0.0,-2.0,0.0,2.0,0.0,2.0,0.0,0.0,0.0,0.0,-2.0,0.0,2.0,0.0,2.0,0.0,0.0 -detJ 0.125,0.125,0.125
768 test:
769 suffix: 5
770 requires: exodusii
771 args: -interpolate -run_type file -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/simpleblock-100.exo -dm_view ascii::ascii_info_detail -v0 -1.5,-0.5,0.5,-0.5,-0.5,0.5,0.5,-0.5,0.5 -J 0.0,0.0,0.5,0.0,0.5,0.0,-0.5,0.0,0.0,0.0,0.0,0.5,0.0,0.5,0.0,-0.5,0.0,0.0,0.0,0.0,0.5,0.0,0.5,0.0,-0.5,0.0,0.0 -invJ 0.0,0.0,-2.0,0.0,2.0,0.0,2.0,0.0,0.0,0.0,0.0,-2.0,0.0,2.0,0.0,2.0,0.0,0.0,0.0,0.0,-2.0,0.0,2.0,0.0,2.0,0.0,0.0 -detJ 0.125,0.125,0.125 -centroid -1.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0 -normal 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 -vol 1.0,1.0,1.0
772
773 TEST*/
774