Actual source code: ex2.c
1: static char help[] = "Tests L2 projection with DMSwarm using delta function particles and deposition of linear shape.\n";
3: #include <petscdmplex.h>
4: #include <petscfe.h>
5: #include <petscdmswarm.h>
6: #include <petscds.h>
7: #include <petscksp.h>
8: #include <petsc/private/petscfeimpl.h>
9: typedef struct {
10: PetscReal L[3]; /* Dimensions of the mesh bounding box */
11: PetscInt particlesPerCell; /* The number of partices per cell */
12: PetscReal particleRelDx; /* Relative particle position perturbation compared to average cell diameter h */
13: PetscReal meshRelDx; /* Relative vertex position perturbation compared to average cell diameter h */
14: PetscInt k; /* Mode number for test function */
15: PetscReal momentTol; /* Tolerance for checking moment conservation */
16: PetscBool shape;
17: PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
18: } AppCtx;
20: /* const char *const ex2FunctionTypes[] = {"linear","x2_x4","sin","ex2FunctionTypes","EX2_FUNCTION_",0}; */
22: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
23: {
24: AppCtx *ctx = (AppCtx *)a_ctx;
25: PetscInt d;
27: u[0] = 0.0;
28: for (d = 0; d < dim; ++d) u[0] += x[d] / (ctx->L[d]);
29: return PETSC_SUCCESS;
30: }
32: static PetscErrorCode x2_x4(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
33: {
34: AppCtx *ctx = (AppCtx *)a_ctx;
35: PetscInt d;
37: u[0] = 1;
38: for (d = 0; d < dim; ++d) u[0] *= PetscSqr(x[d]) * PetscSqr(ctx->L[d]) - PetscPowRealInt(x[d], 4);
39: return PETSC_SUCCESS;
40: }
42: static PetscErrorCode sinx(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *a_ctx)
43: {
44: AppCtx *ctx = (AppCtx *)a_ctx;
46: u[0] = PetscSinScalar(2 * PETSC_PI * ctx->k * x[0] / (ctx->L[0]));
47: return PETSC_SUCCESS;
48: }
50: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
51: {
52: char fstring[PETSC_MAX_PATH_LEN] = "linear";
53: PetscBool flag;
55: PetscFunctionBeginUser;
56: options->particlesPerCell = 1;
57: options->k = 1;
58: options->particleRelDx = 1.e-20;
59: options->meshRelDx = 1.e-20;
60: options->momentTol = 100. * PETSC_MACHINE_EPSILON;
61: options->shape = PETSC_FALSE;
63: PetscOptionsBegin(comm, "", "L2 Projection Options", "DMPLEX");
64: PetscCall(PetscOptionsBool("-shape", "Test linear function projection", "ex2.c", options->shape, &options->shape, NULL));
65: PetscCall(PetscOptionsInt("-k", "Mode number of test", "ex2.c", options->k, &options->k, NULL));
66: PetscCall(PetscOptionsInt("-particlesPerCell", "Number of particles per cell", "ex2.c", options->particlesPerCell, &options->particlesPerCell, NULL));
67: PetscCall(PetscOptionsReal("-particle_perturbation", "Relative perturbation of particles (0,1)", "ex2.c", options->particleRelDx, &options->particleRelDx, NULL));
68: PetscCall(PetscOptionsReal("-mesh_perturbation", "Relative perturbation of mesh points (0,1)", "ex2.c", options->meshRelDx, &options->meshRelDx, NULL));
69: PetscCall(PetscOptionsString("-function", "Name of test function", "ex2.c", fstring, fstring, sizeof(fstring), NULL));
70: PetscCall(PetscStrcmp(fstring, "linear", &flag));
71: if (flag) {
72: options->func = linear;
73: } else {
74: PetscCall(PetscStrcmp(fstring, "sin", &flag));
75: if (flag) {
76: options->func = sinx;
77: } else {
78: PetscCall(PetscStrcmp(fstring, "x2_x4", &flag));
79: options->func = x2_x4;
80: PetscCheck(flag, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown function %s", fstring);
81: }
82: }
83: PetscCall(PetscOptionsReal("-moment_tol", "Tolerance for moment checks", "ex2.c", options->momentTol, &options->momentTol, NULL));
84: PetscOptionsEnd();
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: static PetscErrorCode PerturbVertices(DM dm, AppCtx *user)
89: {
90: PetscRandom rnd;
91: PetscReal interval = user->meshRelDx;
92: Vec coordinates;
93: PetscScalar *coords;
94: PetscReal *hh, low[3], high[3];
95: PetscInt d, cdim, cEnd, N, p, bs;
97: PetscFunctionBeginUser;
98: PetscCall(DMGetBoundingBox(dm, low, high));
99: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
100: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
101: PetscCall(PetscRandomSetInterval(rnd, -interval, interval));
102: PetscCall(PetscRandomSetFromOptions(rnd));
103: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
104: PetscCall(DMGetCoordinateDim(dm, &cdim));
105: PetscCall(PetscCalloc1(cdim, &hh));
106: for (d = 0; d < cdim; ++d) hh[d] = (user->L[d]) / PetscPowReal(cEnd, 1. / cdim);
107: PetscCall(VecGetLocalSize(coordinates, &N));
108: PetscCall(VecGetBlockSize(coordinates, &bs));
109: PetscCheck(bs == cdim, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Coordinate vector has wrong block size %" PetscInt_FMT " != %" PetscInt_FMT, bs, cdim);
110: PetscCall(VecGetArray(coordinates, &coords));
111: for (p = 0; p < N; p += cdim) {
112: PetscScalar *coord = &coords[p], value;
114: for (d = 0; d < cdim; ++d) {
115: PetscCall(PetscRandomGetValue(rnd, &value));
116: coord[d] = PetscMax(low[d], PetscMin(high[d], PetscRealPart(coord[d] + value * hh[d])));
117: }
118: }
119: PetscCall(VecRestoreArray(coordinates, &coords));
120: PetscCall(PetscRandomDestroy(&rnd));
121: PetscCall(PetscFree(hh));
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
126: {
127: PetscReal low[3], high[3];
128: PetscInt cdim, d;
130: PetscFunctionBeginUser;
131: PetscCall(DMCreate(comm, dm));
132: PetscCall(DMSetType(*dm, DMPLEX));
133: PetscCall(DMSetFromOptions(*dm));
134: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
136: PetscCall(DMGetCoordinateDim(*dm, &cdim));
137: PetscCall(DMGetBoundingBox(*dm, low, high));
138: for (d = 0; d < cdim; ++d) user->L[d] = high[d] - low[d];
139: PetscCall(PerturbVertices(*dm, user));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
144: {
145: g0[0] = 1.0;
146: }
148: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
149: {
150: PetscFE fe;
151: PetscDS ds;
152: DMPolytopeType ct;
153: PetscInt dim, cStart;
155: PetscFunctionBeginUser;
156: PetscCall(DMGetDimension(dm, &dim));
157: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
158: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
159: PetscCall(PetscFECreateByCell(PetscObjectComm((PetscObject)dm), dim, 1, ct, NULL, -1, &fe));
160: PetscCall(PetscObjectSetName((PetscObject)fe, "fe"));
161: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
162: PetscCall(DMCreateDS(dm));
163: PetscCall(PetscFEDestroy(&fe));
164: /* Setup to form mass matrix */
165: PetscCall(DMGetDS(dm, &ds));
166: PetscCall(PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
171: {
172: PetscRandom rnd, rndp;
173: PetscReal interval = user->particleRelDx;
174: DMPolytopeType ct;
175: PetscBool simplex;
176: PetscScalar value, *vals;
177: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
178: PetscInt *cellid;
179: PetscInt Ncell, Np = user->particlesPerCell, p, cStart, c, dim, d;
181: PetscFunctionBeginUser;
182: PetscCall(DMGetDimension(dm, &dim));
183: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
184: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
185: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
186: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
187: PetscCall(DMSetType(*sw, DMSWARM));
188: PetscCall(DMSetDimension(*sw, dim));
190: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
191: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
192: PetscCall(PetscRandomSetFromOptions(rnd));
193: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rndp));
194: PetscCall(PetscRandomSetInterval(rndp, -interval, interval));
195: PetscCall(PetscRandomSetFromOptions(rndp));
197: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
198: PetscCall(DMSwarmSetCellDM(*sw, dm));
199: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
200: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
201: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Ncell));
202: PetscCall(DMSwarmSetLocalSizes(*sw, Ncell * Np, 0));
203: PetscCall(DMSetFromOptions(*sw));
204: PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
205: PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
206: PetscCall(DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **)&vals));
208: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
209: for (c = 0; c < Ncell; ++c) {
210: if (Np == 1) {
211: PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL));
212: cellid[c] = c;
213: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
214: } else {
215: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); /* affine */
216: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
217: for (p = 0; p < Np; ++p) {
218: const PetscInt n = c * Np + p;
219: PetscReal sum = 0.0, refcoords[3];
221: cellid[n] = c;
222: for (d = 0; d < dim; ++d) {
223: PetscCall(PetscRandomGetValue(rnd, &value));
224: refcoords[d] = PetscRealPart(value);
225: sum += refcoords[d];
226: }
227: if (simplex && sum > 0.0)
228: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
229: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
230: }
231: }
232: }
233: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
234: for (c = 0; c < Ncell; ++c) {
235: for (p = 0; p < Np; ++p) {
236: const PetscInt n = c * Np + p;
238: for (d = 0; d < dim; ++d) {
239: PetscCall(PetscRandomGetValue(rndp, &value));
240: coords[n * dim + d] += PetscRealPart(value);
241: }
242: PetscCall(user->func(dim, 0.0, &coords[n * dim], 1, &vals[c], user));
243: }
244: }
245: PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
246: PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
247: PetscCall(DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **)&vals));
248: PetscCall(PetscRandomDestroy(&rnd));
249: PetscCall(PetscRandomDestroy(&rndp));
250: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
251: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: static PetscErrorCode CreateParticles_Shape(DM dm, DM *sw, AppCtx *user)
256: {
257: PetscDS prob;
258: PetscFE fe;
259: PetscQuadrature quad;
260: PetscScalar *vals;
261: PetscReal *v0, *J, *invJ, detJ, *coords, *xi0;
262: PetscInt *cellid;
263: const PetscReal *qpoints, *qweights;
264: PetscInt cStart, cEnd, c, Nq, q, dim;
266: PetscFunctionBeginUser;
267: PetscCall(DMGetDimension(dm, &dim));
268: PetscCall(DMGetDS(dm, &prob));
269: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
270: PetscCall(PetscDSGetDiscretization(prob, 0, (PetscObject *)&fe));
271: PetscCall(PetscFEGetQuadrature(fe, &quad));
272: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, &qweights));
273: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
274: PetscCall(DMSetType(*sw, DMSWARM));
275: PetscCall(DMSetDimension(*sw, dim));
276: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
277: PetscCall(DMSwarmSetCellDM(*sw, dm));
278: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_REAL));
279: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
280: PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Nq, 0));
281: PetscCall(DMSetFromOptions(*sw));
282: PetscCall(PetscMalloc4(dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
283: for (c = 0; c < dim; c++) xi0[c] = -1.;
284: PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
285: PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
286: PetscCall(DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **)&vals));
287: for (c = cStart; c < cEnd; ++c) {
288: for (q = 0; q < Nq; ++q) {
289: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
290: cellid[c * Nq + q] = c;
291: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q * dim], &coords[(c * Nq + q) * dim]);
292: PetscCall(linear(dim, 0.0, &coords[(c * Nq + q) * dim], 1, &vals[c * Nq + q], (void *)user));
293: vals[c * Nq + q] *= qweights[q] * detJ;
294: }
295: }
296: PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
297: PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
298: PetscCall(DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **)&vals));
299: PetscCall(PetscFree4(xi0, v0, J, invJ));
300: PetscCall(DMSwarmMigrate(*sw, PETSC_FALSE));
301: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
302: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: static PetscErrorCode computeParticleMoments(DM sw, PetscReal moments[3], AppCtx *user)
307: {
308: DM dm;
309: const PetscReal *coords;
310: const PetscScalar *w;
311: PetscReal mom[3] = {0.0, 0.0, 0.0};
312: PetscInt cell, cStart, cEnd, dim;
314: PetscFunctionBeginUser;
315: PetscCall(DMGetDimension(sw, &dim));
316: PetscCall(DMSwarmGetCellDM(sw, &dm));
317: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
318: PetscCall(DMSwarmSortGetAccess(sw));
319: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
320: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&w));
321: for (cell = cStart; cell < cEnd; ++cell) {
322: PetscInt *pidx;
323: PetscInt Np, p, d;
325: PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx));
326: for (p = 0; p < Np; ++p) {
327: const PetscInt idx = pidx[p];
328: const PetscReal *c = &coords[idx * dim];
330: mom[0] += PetscRealPart(w[idx]);
331: mom[1] += PetscRealPart(w[idx]) * c[0];
332: for (d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx]) * c[d] * c[d];
333: }
334: PetscCall(PetscFree(pidx));
335: }
336: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
337: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
338: PetscCall(DMSwarmSortRestoreAccess(sw));
339: PetscCallMPI(MPIU_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: static void f0_1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
344: {
345: f0[0] = u[0];
346: }
348: static void f0_x(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
349: {
350: f0[0] = x[0] * u[0];
351: }
353: static void f0_r2(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
354: {
355: PetscInt d;
357: f0[0] = 0.0;
358: for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d]) * u[0];
359: }
361: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
362: {
363: PetscDS prob;
364: PetscScalar mom;
366: PetscFunctionBeginUser;
367: PetscCall(DMGetDS(dm, &prob));
368: PetscCall(PetscDSSetObjective(prob, 0, &f0_1));
369: PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
370: moments[0] = PetscRealPart(mom);
371: PetscCall(PetscDSSetObjective(prob, 0, &f0_x));
372: PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
373: moments[1] = PetscRealPart(mom);
374: PetscCall(PetscDSSetObjective(prob, 0, &f0_r2));
375: PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
376: moments[2] = PetscRealPart(mom);
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: static PetscErrorCode TestL2ProjectionParticlesToField(DM dm, DM sw, Vec fhat, AppCtx *user)
381: {
382: const char *fieldnames[1] = {"w_q"};
383: Vec fields[1] = {fhat};
384: PetscReal pmoments[3]; // \int f, \int x f, \int r^2 f
385: PetscReal fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f
387: PetscFunctionBeginUser;
388: PetscCall(DMSwarmProjectFields(sw, dm, 1, fieldnames, fields, SCATTER_FORWARD));
390: /* Check moments of field */
391: PetscCall(computeParticleMoments(sw, pmoments, user));
392: PetscCall(computeFEMMoments(dm, fhat, fmoments, user));
393: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
394: for (PetscInt m = 0; m < 3; ++m) {
395: PetscCheck(PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) <= user->momentTol, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Moment %" PetscInt_FMT " error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]),
396: user->momentTol);
397: }
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: static PetscErrorCode TestL2ProjectionFieldToParticles(DM dm, DM sw, Vec fhat, AppCtx *user)
402: {
403: const char *fieldnames[1] = {"w_q"};
404: Vec fields[1] = {fhat};
405: PetscReal pmoments[3]; // \int f, \int x f, \int r^2 f
406: PetscReal fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f
408: PetscFunctionBeginUser;
409: PetscCall(DMSwarmProjectFields(sw, dm, 1, fieldnames, fields, SCATTER_REVERSE));
411: /* Check moments */
412: PetscCall(computeParticleMoments(sw, pmoments, user));
413: PetscCall(computeFEMMoments(dm, fhat, fmoments, user));
414: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
415: for (PetscInt m = 0; m < 3; ++m) {
416: PetscCheck(PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) <= user->momentTol, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Moment %" PetscInt_FMT " error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]),
417: user->momentTol);
418: }
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /* Interpolate the gradient of an FEM (FVM) field. Code repurposed from DMPlexComputeGradientClementInterpolant */
423: static PetscErrorCode InterpolateGradient(DM dm, Vec locX, Vec locC)
424: {
425: DM_Plex *mesh = (DM_Plex *)dm->data;
426: PetscInt debug = mesh->printFEM;
427: DM dmC;
428: PetscSection section;
429: PetscQuadrature quad = NULL;
430: PetscScalar *interpolant, *gradsum;
431: PetscFEGeom fegeom;
432: PetscReal *coords;
433: const PetscReal *quadPoints, *quadWeights;
434: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset;
436: PetscFunctionBegin;
437: PetscCall(VecGetDM(locC, &dmC));
438: PetscCall(VecSet(locC, 0.0));
439: PetscCall(DMGetDimension(dm, &dim));
440: PetscCall(DMGetCoordinateDim(dm, &coordDim));
441: fegeom.dimEmbed = coordDim;
442: PetscCall(DMGetLocalSection(dm, §ion));
443: PetscCall(PetscSectionGetNumFields(section, &numFields));
444: for (field = 0; field < numFields; ++field) {
445: PetscObject obj;
446: PetscClassId id;
447: PetscInt Nc;
449: PetscCall(DMGetField(dm, field, NULL, &obj));
450: PetscCall(PetscObjectGetClassId(obj, &id));
451: if (id == PETSCFE_CLASSID) {
452: PetscFE fe = (PetscFE)obj;
454: PetscCall(PetscFEGetQuadrature(fe, &quad));
455: PetscCall(PetscFEGetNumComponents(fe, &Nc));
456: } else if (id == PETSCFV_CLASSID) {
457: PetscFV fv = (PetscFV)obj;
459: PetscCall(PetscFVGetQuadrature(fv, &quad));
460: PetscCall(PetscFVGetNumComponents(fv, &Nc));
461: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
462: numComponents += Nc;
463: }
464: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
465: PetscCheck(!(qNc != 1) || !(qNc != numComponents), PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Quadrature components %" PetscInt_FMT " != %" PetscInt_FMT " field components", qNc, numComponents);
466: PetscCall(PetscMalloc6(coordDim * numComponents * 2, &gradsum, coordDim * numComponents, &interpolant, coordDim * Nq, &coords, Nq, &fegeom.detJ, coordDim * coordDim * Nq, &fegeom.J, coordDim * coordDim * Nq, &fegeom.invJ));
467: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
468: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
469: for (v = vStart; v < vEnd; ++v) {
470: PetscInt *star = NULL;
471: PetscInt starSize, st, d, fc;
473: PetscCall(PetscArrayzero(gradsum, coordDim * numComponents));
474: PetscCall(DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
475: for (st = 0; st < starSize * 2; st += 2) {
476: const PetscInt cell = star[st];
477: PetscScalar *grad = &gradsum[coordDim * numComponents];
478: PetscScalar *x = NULL;
480: if ((cell < cStart) || (cell >= cEnd)) continue;
481: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ));
482: PetscCall(DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x));
483: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
484: PetscObject obj;
485: PetscClassId id;
486: PetscInt Nb, Nc, q, qc = 0;
488: PetscCall(PetscArrayzero(grad, coordDim * numComponents));
489: PetscCall(DMGetField(dm, field, NULL, &obj));
490: PetscCall(PetscObjectGetClassId(obj, &id));
491: if (id == PETSCFE_CLASSID) {
492: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
493: PetscCall(PetscFEGetDimension((PetscFE)obj, &Nb));
494: } else if (id == PETSCFV_CLASSID) {
495: PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
496: Nb = 1;
497: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
498: for (q = 0; q < Nq; ++q) {
499: PetscCheck(fegeom.detJ[q] > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT ", quadrature points %" PetscInt_FMT, (double)fegeom.detJ[q], cell, q);
500: if (id == PETSCFE_CLASSID) PetscCall(PetscFEInterpolateGradient_Static((PetscFE)obj, 1, &x[fieldOffset], &fegeom, q, interpolant));
501: else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
502: for (fc = 0; fc < Nc; ++fc) {
503: const PetscReal wt = quadWeights[q * qNc + qc + fc];
505: for (d = 0; d < coordDim; ++d) grad[fc * coordDim + d] += interpolant[fc * dim + d] * wt * fegeom.detJ[q];
506: }
507: }
508: fieldOffset += Nb;
509: }
510: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x));
511: for (fc = 0; fc < numComponents; ++fc) {
512: for (d = 0; d < coordDim; ++d) gradsum[fc * coordDim + d] += grad[fc * coordDim + d];
513: }
514: if (debug) {
515: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT " gradient: [", cell));
516: for (fc = 0; fc < numComponents; ++fc) {
517: for (d = 0; d < coordDim; ++d) {
518: if (fc || d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
519: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc * coordDim + d])));
520: }
521: }
522: PetscCall(PetscPrintf(PETSC_COMM_SELF, "]\n"));
523: }
524: }
525: PetscCall(DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
526: PetscCall(DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES));
527: }
528: PetscCall(PetscFree6(gradsum, interpolant, coords, fegeom.detJ, fegeom.J, fegeom.invJ));
529: PetscFunctionReturn(PETSC_SUCCESS);
530: }
532: static PetscErrorCode TestFieldGradientProjection(DM dm, DM sw, AppCtx *user)
533: {
534: MPI_Comm comm;
535: KSP ksp;
536: Mat M; /* FEM mass matrix */
537: Mat M_p; /* Particle mass matrix */
538: Vec f, rhs, fhat, grad; /* Particle field f, \int phi_i f, FEM field */
539: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
540: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
541: PetscInt m;
543: PetscFunctionBeginUser;
544: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
545: PetscCall(KSPCreate(comm, &ksp));
546: PetscCall(KSPSetOptionsPrefix(ksp, "ptof_"));
547: PetscCall(DMGetGlobalVector(dm, &fhat));
548: PetscCall(DMGetGlobalVector(dm, &rhs));
549: PetscCall(DMGetGlobalVector(dm, &grad));
551: PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
552: PetscCall(MatViewFromOptions(M_p, NULL, "-M_p_view"));
554: /* make particle weight vector */
555: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
557: /* create matrix RHS vector */
558: PetscCall(MatMultTranspose(M_p, f, rhs));
559: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
560: PetscCall(PetscObjectSetName((PetscObject)rhs, "rhs"));
561: PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));
563: PetscCall(DMCreateMatrix(dm, &M));
564: PetscCall(DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user));
566: PetscCall(InterpolateGradient(dm, fhat, grad));
568: PetscCall(MatViewFromOptions(M, NULL, "-M_view"));
569: PetscCall(KSPSetOperators(ksp, M, M));
570: PetscCall(KSPSetFromOptions(ksp));
571: PetscCall(KSPSolve(ksp, rhs, grad));
572: PetscCall(PetscObjectSetName((PetscObject)fhat, "fhat"));
573: PetscCall(VecViewFromOptions(fhat, NULL, "-fhat_view"));
575: /* Check moments of field */
576: PetscCall(computeParticleMoments(sw, pmoments, user));
577: PetscCall(computeFEMMoments(dm, grad, fmoments, user));
578: PetscCall(PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
579: for (m = 0; m < 3; ++m) {
580: PetscCheck(PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]) <= user->momentTol, comm, PETSC_ERR_ARG_WRONG, "Moment %" PetscInt_FMT " error too large %g > %g", m, PetscAbsReal((fmoments[m] - pmoments[m]) / fmoments[m]), user->momentTol);
581: }
583: PetscCall(KSPDestroy(&ksp));
584: PetscCall(MatDestroy(&M));
585: PetscCall(MatDestroy(&M_p));
586: PetscCall(DMRestoreGlobalVector(dm, &fhat));
587: PetscCall(DMRestoreGlobalVector(dm, &rhs));
588: PetscCall(DMRestoreGlobalVector(dm, &grad));
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: static PetscErrorCode TestL2Projection_Shape(DM dm, DM sw, AppCtx *user)
593: {
594: PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
595: KSP ksp;
596: Mat mass;
597: Vec u, rhs, uproj;
598: void *ctxs[1];
599: PetscReal error;
601: PetscFunctionBeginUser;
602: ctxs[0] = (void *)user;
603: funcs[0] = linear;
605: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &u));
606: PetscCall(VecViewFromOptions(u, NULL, "-f_view"));
607: PetscCall(DMGetGlobalVector(dm, &rhs));
608: PetscCall(DMCreateMassMatrix(sw, dm, &mass));
609: PetscCall(MatMultTranspose(mass, u, rhs));
610: PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));
611: PetscCall(MatDestroy(&mass));
612: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &u));
613: PetscCall(DMGetGlobalVector(dm, &uproj));
614: PetscCall(DMCreateMatrix(dm, &mass));
615: PetscCall(DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user));
616: PetscCall(MatViewFromOptions(mass, NULL, "-mass_mat_view"));
617: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
618: PetscCall(KSPSetOperators(ksp, mass, mass));
619: PetscCall(KSPSetFromOptions(ksp));
620: PetscCall(KSPSolve(ksp, rhs, uproj));
621: PetscCall(KSPDestroy(&ksp));
622: PetscCall(PetscObjectSetName((PetscObject)uproj, "Full Projection"));
623: PetscCall(VecViewFromOptions(uproj, NULL, "-proj_vec_view"));
624: PetscCall(DMComputeL2Diff(dm, 0.0, funcs, ctxs, uproj, &error));
625: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error));
626: PetscCall(MatDestroy(&mass));
627: PetscCall(DMRestoreGlobalVector(dm, &rhs));
628: PetscCall(DMRestoreGlobalVector(dm, &uproj));
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: int main(int argc, char *argv[])
633: {
634: MPI_Comm comm;
635: DM dm, sw;
636: AppCtx user;
638: PetscFunctionBeginUser;
639: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
640: comm = PETSC_COMM_WORLD;
641: PetscCall(ProcessOptions(comm, &user));
642: PetscCall(CreateMesh(comm, &dm, &user));
643: PetscCall(CreateFEM(dm, &user));
644: if (!user.shape) {
645: Vec fhat;
647: PetscCall(CreateParticles(dm, &sw, &user));
648: PetscCall(DMGetGlobalVector(dm, &fhat));
649: PetscCall(TestL2ProjectionParticlesToField(dm, sw, fhat, &user));
650: PetscCall(TestL2ProjectionFieldToParticles(dm, sw, fhat, &user));
651: PetscCall(TestFieldGradientProjection(dm, sw, &user));
652: PetscCall(DMRestoreGlobalVector(dm, &fhat));
653: } else {
654: PetscCall(CreateParticles_Shape(dm, &sw, &user));
655: PetscCall(TestL2Projection_Shape(dm, sw, &user));
656: }
657: PetscCall(DMDestroy(&dm));
658: PetscCall(DMDestroy(&sw));
659: PetscCall(PetscFinalize());
660: return 0;
661: }
663: /*TEST
665: # Swarm does not handle complex or quad
666: build:
667: requires: !complex double
669: test:
670: suffix: proj_tri_0
671: requires: triangle
672: args: -dm_plex_box_faces 1,1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
673: filter: grep -v marker | grep -v atomic | grep -v usage
675: test:
676: suffix: proj_tri_2_faces
677: requires: triangle
678: args: -dm_plex_box_faces 2,2 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
679: filter: grep -v marker | grep -v atomic | grep -v usage
681: test:
682: suffix: proj_quad_0
683: requires: triangle
684: args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
685: filter: grep -v marker | grep -v atomic | grep -v usage
687: test:
688: suffix: proj_quad_2_faces
689: requires: triangle
690: args: -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
691: filter: grep -v marker | grep -v atomic | grep -v usage
693: test:
694: suffix: proj_tri_5P
695: requires: triangle
696: args: -dm_plex_box_faces 1,1 -particlesPerCell 5 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
697: filter: grep -v marker | grep -v atomic | grep -v usage
699: test:
700: suffix: proj_quad_5P
701: requires: triangle
702: args: -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -particlesPerCell 5 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
703: filter: grep -v marker | grep -v atomic | grep -v usage
705: test:
706: suffix: proj_tri_mdx
707: requires: triangle
708: args: -dm_plex_box_faces 1,1 -mesh_perturbation 1.0e-1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
709: filter: grep -v marker | grep -v atomic | grep -v usage
711: test:
712: suffix: proj_tri_mdx_5P
713: requires: triangle
714: args: -dm_plex_box_faces 1,1 -particlesPerCell 5 -mesh_perturbation 1.0e-1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
715: filter: grep -v marker | grep -v atomic | grep -v usage
717: test:
718: suffix: proj_tri_3d
719: requires: ctetgen
720: args: -dm_plex_dim 3 -dm_plex_box_faces 1,1,1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
721: filter: grep -v marker | grep -v atomic | grep -v usage
723: test:
724: suffix: proj_tri_3d_2_faces
725: requires: ctetgen
726: args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
727: filter: grep -v marker | grep -v atomic | grep -v usage
729: test:
730: suffix: proj_tri_3d_5P
731: requires: ctetgen
732: args: -dm_plex_dim 3 -dm_plex_box_faces 1,1,1 -particlesPerCell 5 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
733: filter: grep -v marker | grep -v atomic | grep -v usage
735: test:
736: suffix: proj_tri_3d_mdx
737: requires: ctetgen
738: args: -dm_plex_dim 3 -dm_plex_box_faces 1,1,1 -mesh_perturbation 1.0e-1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
739: filter: grep -v marker | grep -v atomic | grep -v usage
741: test:
742: suffix: proj_tri_3d_mdx_5P
743: requires: ctetgen
744: args: -dm_plex_dim 3 -dm_plex_box_faces 1,1,1 -particlesPerCell 5 -mesh_perturbation 1.0e-1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
745: filter: grep -v marker | grep -v atomic | grep -v usage
747: test:
748: suffix: proj_tri_3d_mdx_2_faces
749: requires: ctetgen
750: args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -mesh_perturbation 1.0e-1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
751: filter: grep -v marker | grep -v atomic | grep -v usage
753: test:
754: suffix: proj_tri_3d_mdx_5P_2_faces
755: requires: ctetgen
756: args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -particlesPerCell 5 -mesh_perturbation 1.0e-1 -dm_view -sw_view -petscspace_degree 2 -petscfe_default_quadrature_order {{2 3}} -ptof_pc_type lu -ftop_ksp_rtol 1e-15 -ftop_ksp_type lsqr -ftop_pc_type none
757: filter: grep -v marker | grep -v atomic | grep -v usage
759: test:
760: suffix: proj_quad_lsqr_scale
761: requires: !complex
762: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
763: -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
764: -particlesPerCell 200 \
765: -ptof_pc_type lu \
766: -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type none
768: test:
769: suffix: proj_quad_lsqr_prec_scale
770: requires: !complex
771: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
772: -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
773: -particlesPerCell 200 \
774: -ptof_pc_type lu \
775: -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type bjacobi -ftop_sub_pc_type lu -ftop_sub_pc_factor_shift_type nonzero
776: test:
777: suffix: proj_shape_linear_tri_2d
778: requires: ctetgen triangle
779: args: -shape -dm_plex_box_faces 2,2\
780: -dm_view -sw_view\
781: -petscspace_degree 1 -petscfe_default_quadrature_order 1\
782: -pc_type lu
783: test:
784: suffix: proj_shape_linear_quad_2d
785: requires: ctetgen triangle
786: args: -shape -dm_plex_simplex 0 -dm_plex_box_faces 2,2\
787: -dm_view -sw_view\
788: -petscspace_degree 1 -petscfe_default_quadrature_order 1\
789: -pc_type lu
790: test:
791: suffix: proj_shape_linear_tri_3d
792: requires: ctetgen triangle
793: args: -shape -dm_plex_dim 3 -dm_plex_box_faces 2,2,2\
794: -dm_view -sw_view\
795: -petscspace_degree 1 -petscfe_default_quadrature_order 1\
796: -pc_type lu
797: test:
798: suffix: proj_shape_linear_quad_3d
799: requires: ctetgen triangle
800: args: -shape -dm_plex_dim 3 -dm_plex_simplex 0\
801: -dm_plex_box_faces 2,2,2\
802: -dm_view -sw_view\
803: -petscspace_degree 1 -petscfe_default_quadrature_order 1\
804: -pc_type lu
805: TEST*/