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();
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: static PetscErrorCode PerturbVertices(DM dm, AppCtx *user)
90: {
91: PetscRandom rnd;
92: PetscReal interval = user->meshRelDx;
93: Vec coordinates;
94: PetscScalar *coords;
95: PetscReal *hh, low[3], high[3];
96: PetscInt d, cdim, cEnd, N, p, bs;
98: PetscFunctionBeginUser;
99: PetscCall(DMGetBoundingBox(dm, low, high));
100: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
101: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
102: PetscCall(PetscRandomSetInterval(rnd, -interval, interval));
103: PetscCall(PetscRandomSetFromOptions(rnd));
104: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
105: PetscCall(DMGetCoordinateDim(dm, &cdim));
106: PetscCall(PetscCalloc1(cdim, &hh));
107: for (d = 0; d < cdim; ++d) hh[d] = (user->L[d]) / PetscPowReal(cEnd, 1. / cdim);
108: PetscCall(VecGetLocalSize(coordinates, &N));
109: PetscCall(VecGetBlockSize(coordinates, &bs));
110: PetscCheck(bs == cdim, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Coordinate vector has wrong block size %" PetscInt_FMT " != %" PetscInt_FMT, bs, cdim);
111: PetscCall(VecGetArray(coordinates, &coords));
112: for (p = 0; p < N; p += cdim) {
113: PetscScalar *coord = &coords[p], value;
115: for (d = 0; d < cdim; ++d) {
116: PetscCall(PetscRandomGetValue(rnd, &value));
117: coord[d] = PetscMax(low[d], PetscMin(high[d], PetscRealPart(coord[d] + value * hh[d])));
118: }
119: }
120: PetscCall(VecRestoreArray(coordinates, &coords));
121: PetscCall(PetscRandomDestroy(&rnd));
122: PetscCall(PetscFree(hh));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
127: {
128: PetscReal low[3], high[3];
129: PetscInt cdim, d;
131: PetscFunctionBeginUser;
132: PetscCall(DMCreate(comm, dm));
133: PetscCall(DMSetType(*dm, DMPLEX));
134: PetscCall(DMSetFromOptions(*dm));
135: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
137: PetscCall(DMGetCoordinateDim(*dm, &cdim));
138: PetscCall(DMGetBoundingBox(*dm, low, high));
139: for (d = 0; d < cdim; ++d) user->L[d] = high[d] - low[d];
140: PetscCall(PerturbVertices(*dm, user));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: 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[])
145: {
146: g0[0] = 1.0;
147: }
149: static PetscErrorCode CreateFEM(DM dm, AppCtx *user)
150: {
151: PetscFE fe;
152: PetscDS ds;
153: DMPolytopeType ct;
154: PetscInt dim, cStart;
156: PetscFunctionBeginUser;
157: PetscCall(DMGetDimension(dm, &dim));
158: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
159: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
160: PetscCall(PetscFECreateByCell(PetscObjectComm((PetscObject)dm), dim, 1, ct, NULL, -1, &fe));
161: PetscCall(PetscObjectSetName((PetscObject)fe, "fe"));
162: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
163: PetscCall(DMCreateDS(dm));
164: PetscCall(PetscFEDestroy(&fe));
165: /* Setup to form mass matrix */
166: PetscCall(DMGetDS(dm, &ds));
167: PetscCall(PetscDSSetJacobian(ds, 0, 0, identity, NULL, NULL, NULL));
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
172: {
173: PetscRandom rnd, rndp;
174: PetscReal interval = user->particleRelDx;
175: DMPolytopeType ct;
176: PetscBool simplex;
177: PetscScalar value, *vals;
178: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
179: PetscInt *cellid;
180: PetscInt Ncell, Np = user->particlesPerCell, p, cStart, c, dim, d;
182: PetscFunctionBeginUser;
183: PetscCall(DMGetDimension(dm, &dim));
184: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
185: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
186: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
187: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
188: PetscCall(DMSetType(*sw, DMSWARM));
189: PetscCall(DMSetDimension(*sw, dim));
191: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
192: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
193: PetscCall(PetscRandomSetFromOptions(rnd));
194: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rndp));
195: PetscCall(PetscRandomSetInterval(rndp, -interval, interval));
196: PetscCall(PetscRandomSetFromOptions(rndp));
198: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
199: PetscCall(DMSwarmSetCellDM(*sw, dm));
200: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
201: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
202: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Ncell));
203: PetscCall(DMSwarmSetLocalSizes(*sw, Ncell * Np, 0));
204: PetscCall(DMSetFromOptions(*sw));
205: PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
206: PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
207: PetscCall(DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **)&vals));
209: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
210: for (c = 0; c < Ncell; ++c) {
211: if (Np == 1) {
212: PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL));
213: cellid[c] = c;
214: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
215: } else {
216: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ)); /* affine */
217: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
218: for (p = 0; p < Np; ++p) {
219: const PetscInt n = c * Np + p;
220: PetscReal sum = 0.0, refcoords[3];
222: cellid[n] = c;
223: for (d = 0; d < dim; ++d) {
224: PetscCall(PetscRandomGetValue(rnd, &value));
225: refcoords[d] = PetscRealPart(value);
226: sum += refcoords[d];
227: }
228: if (simplex && sum > 0.0)
229: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
230: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
231: }
232: }
233: }
234: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
235: for (c = 0; c < Ncell; ++c) {
236: for (p = 0; p < Np; ++p) {
237: const PetscInt n = c * Np + p;
239: for (d = 0; d < dim; ++d) {
240: PetscCall(PetscRandomGetValue(rndp, &value));
241: coords[n * dim + d] += PetscRealPart(value);
242: }
243: PetscCall(user->func(dim, 0.0, &coords[n * dim], 1, &vals[c], user));
244: }
245: }
246: PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
247: PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
248: PetscCall(DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **)&vals));
249: PetscCall(PetscRandomDestroy(&rnd));
250: PetscCall(PetscRandomDestroy(&rndp));
251: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
252: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: static PetscErrorCode CreateParticles_Shape(DM dm, DM *sw, AppCtx *user)
257: {
258: PetscDS prob;
259: PetscFE fe;
260: PetscQuadrature quad;
261: PetscScalar *vals;
262: PetscReal *v0, *J, *invJ, detJ, *coords, *xi0;
263: PetscInt *cellid;
264: const PetscReal *qpoints, *qweights;
265: PetscInt cStart, cEnd, c, Nq, q, dim;
267: PetscFunctionBeginUser;
268: PetscCall(DMGetDimension(dm, &dim));
269: PetscCall(DMGetDS(dm, &prob));
270: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
271: PetscCall(PetscDSGetDiscretization(prob, 0, (PetscObject *)&fe));
272: PetscCall(PetscFEGetQuadrature(fe, &quad));
273: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &qpoints, &qweights));
274: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
275: PetscCall(DMSetType(*sw, DMSWARM));
276: PetscCall(DMSetDimension(*sw, dim));
277: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
278: PetscCall(DMSwarmSetCellDM(*sw, dm));
279: PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_REAL));
280: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
281: PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Nq, 0));
282: PetscCall(DMSetFromOptions(*sw));
283: PetscCall(PetscMalloc4(dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
284: for (c = 0; c < dim; c++) xi0[c] = -1.;
285: PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
286: PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
287: PetscCall(DMSwarmGetField(*sw, "w_q", NULL, NULL, (void **)&vals));
288: for (c = cStart; c < cEnd; ++c) {
289: for (q = 0; q < Nq; ++q) {
290: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
291: cellid[c * Nq + q] = c;
292: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q * dim], &coords[(c * Nq + q) * dim]);
293: PetscCall(linear(dim, 0.0, &coords[(c * Nq + q) * dim], 1, &vals[c * Nq + q], (void *)user));
294: vals[c * Nq + q] *= qweights[q] * detJ;
295: }
296: }
297: PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
298: PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
299: PetscCall(DMSwarmRestoreField(*sw, "w_q", NULL, NULL, (void **)&vals));
300: PetscCall(PetscFree4(xi0, v0, J, invJ));
301: PetscCall(DMSwarmMigrate(*sw, PETSC_FALSE));
302: PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
303: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: static PetscErrorCode computeParticleMoments(DM sw, PetscReal moments[3], AppCtx *user)
308: {
309: DM dm;
310: const PetscReal *coords;
311: const PetscScalar *w;
312: PetscReal mom[3] = {0.0, 0.0, 0.0};
313: PetscInt cell, cStart, cEnd, dim;
315: PetscFunctionBeginUser;
316: PetscCall(DMGetDimension(sw, &dim));
317: PetscCall(DMSwarmGetCellDM(sw, &dm));
318: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
319: PetscCall(DMSwarmSortGetAccess(sw));
320: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
321: PetscCall(DMSwarmGetField(sw, "w_q", NULL, NULL, (void **)&w));
322: for (cell = cStart; cell < cEnd; ++cell) {
323: PetscInt *pidx;
324: PetscInt Np, p, d;
326: PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Np, &pidx));
327: for (p = 0; p < Np; ++p) {
328: const PetscInt idx = pidx[p];
329: const PetscReal *c = &coords[idx * dim];
331: mom[0] += PetscRealPart(w[idx]);
332: mom[1] += PetscRealPart(w[idx]) * c[0];
333: for (d = 0; d < dim; ++d) mom[2] += PetscRealPart(w[idx]) * c[d] * c[d];
334: }
335: PetscCall(PetscFree(pidx));
336: }
337: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
338: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
339: PetscCall(DMSwarmSortRestoreAccess(sw));
340: PetscCall(MPIU_Allreduce(mom, moments, 3, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
344: 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[])
345: {
346: f0[0] = u[0];
347: }
349: 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[])
350: {
351: f0[0] = x[0] * u[0];
352: }
354: 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[])
355: {
356: PetscInt d;
358: f0[0] = 0.0;
359: for (d = 0; d < dim; ++d) f0[0] += PetscSqr(x[d]) * u[0];
360: }
362: static PetscErrorCode computeFEMMoments(DM dm, Vec u, PetscReal moments[3], AppCtx *user)
363: {
364: PetscDS prob;
365: PetscScalar mom;
367: PetscFunctionBeginUser;
368: PetscCall(DMGetDS(dm, &prob));
369: PetscCall(PetscDSSetObjective(prob, 0, &f0_1));
370: PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
371: moments[0] = PetscRealPart(mom);
372: PetscCall(PetscDSSetObjective(prob, 0, &f0_x));
373: PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
374: moments[1] = PetscRealPart(mom);
375: PetscCall(PetscDSSetObjective(prob, 0, &f0_r2));
376: PetscCall(DMPlexComputeIntegralFEM(dm, u, &mom, user));
377: moments[2] = PetscRealPart(mom);
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: static PetscErrorCode TestL2ProjectionParticlesToField(DM dm, DM sw, Vec fhat, AppCtx *user)
382: {
383: const char *fieldnames[1] = {"w_q"};
384: Vec fields[1] = {fhat};
385: PetscReal pmoments[3]; // \int f, \int x f, \int r^2 f
386: PetscReal fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f
388: PetscFunctionBeginUser;
389: PetscCall(DMSwarmProjectFields(sw, 1, fieldnames, fields, SCATTER_FORWARD));
391: /* Check moments of field */
392: PetscCall(computeParticleMoments(sw, pmoments, user));
393: PetscCall(computeFEMMoments(dm, fhat, fmoments, user));
394: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
395: for (PetscInt m = 0; m < 3; ++m) {
396: 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]),
397: user->momentTol);
398: }
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: static PetscErrorCode TestL2ProjectionFieldToParticles(DM dm, DM sw, Vec fhat, AppCtx *user)
403: {
404: const char *fieldnames[1] = {"w_q"};
405: Vec fields[1] = {fhat};
406: PetscReal pmoments[3]; // \int f, \int x f, \int r^2 f
407: PetscReal fmoments[3]; // \int \hat f, \int x \hat f, \int r^2 \hat f
409: PetscFunctionBeginUser;
410: PetscCall(DMSwarmProjectFields(sw, 1, fieldnames, fields, SCATTER_REVERSE));
412: /* Check moments */
413: PetscCall(computeParticleMoments(sw, pmoments, user));
414: PetscCall(computeFEMMoments(dm, fhat, fmoments, user));
415: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
416: for (PetscInt m = 0; m < 3; ++m) {
417: 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]),
418: user->momentTol);
419: }
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: /* Interpolate the gradient of an FEM (FVM) field. Code repurposed from DMPlexComputeGradientClementInterpolant */
424: static PetscErrorCode InterpolateGradient(DM dm, Vec locX, Vec locC)
425: {
426: DM_Plex *mesh = (DM_Plex *)dm->data;
427: PetscInt debug = mesh->printFEM;
428: DM dmC;
429: PetscSection section;
430: PetscQuadrature quad = NULL;
431: PetscScalar *interpolant, *gradsum;
432: PetscFEGeom fegeom;
433: PetscReal *coords;
434: const PetscReal *quadPoints, *quadWeights;
435: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset;
437: PetscFunctionBegin;
438: PetscCall(VecGetDM(locC, &dmC));
439: PetscCall(VecSet(locC, 0.0));
440: PetscCall(DMGetDimension(dm, &dim));
441: PetscCall(DMGetCoordinateDim(dm, &coordDim));
442: fegeom.dimEmbed = coordDim;
443: PetscCall(DMGetLocalSection(dm, §ion));
444: PetscCall(PetscSectionGetNumFields(section, &numFields));
445: for (field = 0; field < numFields; ++field) {
446: PetscObject obj;
447: PetscClassId id;
448: PetscInt Nc;
450: PetscCall(DMGetField(dm, field, NULL, &obj));
451: PetscCall(PetscObjectGetClassId(obj, &id));
452: if (id == PETSCFE_CLASSID) {
453: PetscFE fe = (PetscFE)obj;
455: PetscCall(PetscFEGetQuadrature(fe, &quad));
456: PetscCall(PetscFEGetNumComponents(fe, &Nc));
457: } else if (id == PETSCFV_CLASSID) {
458: PetscFV fv = (PetscFV)obj;
460: PetscCall(PetscFVGetQuadrature(fv, &quad));
461: PetscCall(PetscFVGetNumComponents(fv, &Nc));
462: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
463: numComponents += Nc;
464: }
465: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
466: PetscCheck(!(qNc != 1) || !(qNc != numComponents), PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Quadrature components %" PetscInt_FMT " != %" PetscInt_FMT " field components", qNc, numComponents);
467: 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));
468: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
469: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
470: for (v = vStart; v < vEnd; ++v) {
471: PetscInt *star = NULL;
472: PetscInt starSize, st, d, fc;
474: PetscCall(PetscArrayzero(gradsum, coordDim * numComponents));
475: PetscCall(DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
476: for (st = 0; st < starSize * 2; st += 2) {
477: const PetscInt cell = star[st];
478: PetscScalar *grad = &gradsum[coordDim * numComponents];
479: PetscScalar *x = NULL;
481: if ((cell < cStart) || (cell >= cEnd)) continue;
482: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ));
483: PetscCall(DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x));
484: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
485: PetscObject obj;
486: PetscClassId id;
487: PetscInt Nb, Nc, q, qc = 0;
489: PetscCall(PetscArrayzero(grad, coordDim * numComponents));
490: PetscCall(DMGetField(dm, field, NULL, &obj));
491: PetscCall(PetscObjectGetClassId(obj, &id));
492: if (id == PETSCFE_CLASSID) {
493: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
494: PetscCall(PetscFEGetDimension((PetscFE)obj, &Nb));
495: } else if (id == PETSCFV_CLASSID) {
496: PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
497: Nb = 1;
498: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
499: for (q = 0; q < Nq; ++q) {
500: 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);
501: if (id == PETSCFE_CLASSID) PetscCall(PetscFEInterpolateGradient_Static((PetscFE)obj, 1, &x[fieldOffset], &fegeom, q, interpolant));
502: else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, field);
503: for (fc = 0; fc < Nc; ++fc) {
504: const PetscReal wt = quadWeights[q * qNc + qc + fc];
506: for (d = 0; d < coordDim; ++d) grad[fc * coordDim + d] += interpolant[fc * dim + d] * wt * fegeom.detJ[q];
507: }
508: }
509: fieldOffset += Nb;
510: }
511: PetscCall(DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x));
512: for (fc = 0; fc < numComponents; ++fc) {
513: for (d = 0; d < coordDim; ++d) gradsum[fc * coordDim + d] += grad[fc * coordDim + d];
514: }
515: if (debug) {
516: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT " gradient: [", cell));
517: for (fc = 0; fc < numComponents; ++fc) {
518: for (d = 0; d < coordDim; ++d) {
519: if (fc || d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
520: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc * coordDim + d])));
521: }
522: }
523: PetscCall(PetscPrintf(PETSC_COMM_SELF, "]\n"));
524: }
525: }
526: PetscCall(DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star));
527: PetscCall(DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES));
528: }
529: PetscCall(PetscFree6(gradsum, interpolant, coords, fegeom.detJ, fegeom.J, fegeom.invJ));
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: static PetscErrorCode TestFieldGradientProjection(DM dm, DM sw, AppCtx *user)
534: {
535: MPI_Comm comm;
536: KSP ksp;
537: Mat M; /* FEM mass matrix */
538: Mat M_p; /* Particle mass matrix */
539: Vec f, rhs, fhat, grad; /* Particle field f, \int phi_i f, FEM field */
540: PetscReal pmoments[3]; /* \int f, \int x f, \int r^2 f */
541: PetscReal fmoments[3]; /* \int \hat f, \int x \hat f, \int r^2 \hat f */
542: PetscInt m;
544: PetscFunctionBeginUser;
545: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
546: PetscCall(KSPCreate(comm, &ksp));
547: PetscCall(KSPSetOptionsPrefix(ksp, "ptof_"));
548: PetscCall(DMGetGlobalVector(dm, &fhat));
549: PetscCall(DMGetGlobalVector(dm, &rhs));
550: PetscCall(DMGetGlobalVector(dm, &grad));
552: PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
553: PetscCall(MatViewFromOptions(M_p, NULL, "-M_p_view"));
555: /* make particle weight vector */
556: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f));
558: /* create matrix RHS vector */
559: PetscCall(MatMultTranspose(M_p, f, rhs));
560: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f));
561: PetscCall(PetscObjectSetName((PetscObject)rhs, "rhs"));
562: PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));
564: PetscCall(DMCreateMatrix(dm, &M));
565: PetscCall(DMPlexSNESComputeJacobianFEM(dm, fhat, M, M, user));
567: PetscCall(InterpolateGradient(dm, fhat, grad));
569: PetscCall(MatViewFromOptions(M, NULL, "-M_view"));
570: PetscCall(KSPSetOperators(ksp, M, M));
571: PetscCall(KSPSetFromOptions(ksp));
572: PetscCall(KSPSolve(ksp, rhs, grad));
573: PetscCall(PetscObjectSetName((PetscObject)fhat, "fhat"));
574: PetscCall(VecViewFromOptions(fhat, NULL, "-fhat_view"));
576: /* Check moments of field */
577: PetscCall(computeParticleMoments(sw, pmoments, user));
578: PetscCall(computeFEMMoments(dm, grad, fmoments, user));
579: PetscCall(PetscPrintf(comm, "L2 projection mass: %20.10e, x-momentum: %20.10e, energy: %20.10e.\n", fmoments[0], fmoments[1], fmoments[2]));
580: for (m = 0; m < 3; ++m) {
581: 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);
582: }
584: PetscCall(KSPDestroy(&ksp));
585: PetscCall(MatDestroy(&M));
586: PetscCall(MatDestroy(&M_p));
587: PetscCall(DMRestoreGlobalVector(dm, &fhat));
588: PetscCall(DMRestoreGlobalVector(dm, &rhs));
589: PetscCall(DMRestoreGlobalVector(dm, &grad));
591: PetscFunctionReturn(PETSC_SUCCESS);
592: }
594: static PetscErrorCode TestL2Projection_Shape(DM dm, DM sw, AppCtx *user)
595: {
596: PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);
597: KSP ksp;
598: Mat mass;
599: Vec u, rhs, uproj;
600: void *ctxs[1];
601: PetscReal error;
603: PetscFunctionBeginUser;
604: ctxs[0] = (void *)user;
605: funcs[0] = linear;
607: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &u));
608: PetscCall(VecViewFromOptions(u, NULL, "-f_view"));
609: PetscCall(DMGetGlobalVector(dm, &rhs));
610: PetscCall(DMCreateMassMatrix(sw, dm, &mass));
611: PetscCall(MatMultTranspose(mass, u, rhs));
612: PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));
613: PetscCall(MatDestroy(&mass));
614: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &u));
615: PetscCall(DMGetGlobalVector(dm, &uproj));
616: PetscCall(DMCreateMatrix(dm, &mass));
617: PetscCall(DMPlexSNESComputeJacobianFEM(dm, uproj, mass, mass, user));
618: PetscCall(MatViewFromOptions(mass, NULL, "-mass_mat_view"));
619: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
620: PetscCall(KSPSetOperators(ksp, mass, mass));
621: PetscCall(KSPSetFromOptions(ksp));
622: PetscCall(KSPSolve(ksp, rhs, uproj));
623: PetscCall(KSPDestroy(&ksp));
624: PetscCall(PetscObjectSetName((PetscObject)uproj, "Full Projection"));
625: PetscCall(VecViewFromOptions(uproj, NULL, "-proj_vec_view"));
626: PetscCall(DMComputeL2Diff(dm, 0.0, funcs, ctxs, uproj, &error));
627: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error));
628: PetscCall(MatDestroy(&mass));
629: PetscCall(DMRestoreGlobalVector(dm, &rhs));
630: PetscCall(DMRestoreGlobalVector(dm, &uproj));
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: int main(int argc, char *argv[])
635: {
636: MPI_Comm comm;
637: DM dm, sw;
638: AppCtx user;
640: PetscFunctionBeginUser;
641: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
642: comm = PETSC_COMM_WORLD;
643: PetscCall(ProcessOptions(comm, &user));
644: PetscCall(CreateMesh(comm, &dm, &user));
645: PetscCall(CreateFEM(dm, &user));
646: if (!user.shape) {
647: Vec fhat;
649: PetscCall(CreateParticles(dm, &sw, &user));
650: PetscCall(DMGetGlobalVector(dm, &fhat));
651: PetscCall(TestL2ProjectionParticlesToField(dm, sw, fhat, &user));
652: PetscCall(TestL2ProjectionFieldToParticles(dm, sw, fhat, &user));
653: PetscCall(TestFieldGradientProjection(dm, sw, &user));
654: PetscCall(DMRestoreGlobalVector(dm, &fhat));
655: } else {
656: PetscCall(CreateParticles_Shape(dm, &sw, &user));
657: PetscCall(TestL2Projection_Shape(dm, sw, &user));
658: }
659: PetscCall(DMDestroy(&dm));
660: PetscCall(DMDestroy(&sw));
661: PetscCall(PetscFinalize());
662: return 0;
663: }
665: /*TEST
667: # Swarm does not handle complex or quad
668: build:
669: requires: !complex double
671: test:
672: suffix: proj_tri_0
673: requires: triangle
674: 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
675: filter: grep -v marker | grep -v atomic | grep -v usage
677: test:
678: suffix: proj_tri_2_faces
679: requires: triangle
680: 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
681: filter: grep -v marker | grep -v atomic | grep -v usage
683: test:
684: suffix: proj_quad_0
685: requires: triangle
686: 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
687: filter: grep -v marker | grep -v atomic | grep -v usage
689: test:
690: suffix: proj_quad_2_faces
691: requires: triangle
692: 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
693: filter: grep -v marker | grep -v atomic | grep -v usage
695: test:
696: suffix: proj_tri_5P
697: requires: triangle
698: 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
699: filter: grep -v marker | grep -v atomic | grep -v usage
701: test:
702: suffix: proj_quad_5P
703: requires: triangle
704: 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
705: filter: grep -v marker | grep -v atomic | grep -v usage
707: test:
708: suffix: proj_tri_mdx
709: requires: triangle
710: 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
711: filter: grep -v marker | grep -v atomic | grep -v usage
713: test:
714: suffix: proj_tri_mdx_5P
715: requires: triangle
716: 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
717: filter: grep -v marker | grep -v atomic | grep -v usage
719: test:
720: suffix: proj_tri_3d
721: requires: ctetgen
722: 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
723: filter: grep -v marker | grep -v atomic | grep -v usage
725: test:
726: suffix: proj_tri_3d_2_faces
727: requires: ctetgen
728: 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
729: filter: grep -v marker | grep -v atomic | grep -v usage
731: test:
732: suffix: proj_tri_3d_5P
733: requires: ctetgen
734: 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
735: filter: grep -v marker | grep -v atomic | grep -v usage
737: test:
738: suffix: proj_tri_3d_mdx
739: requires: ctetgen
740: 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
741: filter: grep -v marker | grep -v atomic | grep -v usage
743: test:
744: suffix: proj_tri_3d_mdx_5P
745: requires: ctetgen
746: 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
747: filter: grep -v marker | grep -v atomic | grep -v usage
749: test:
750: suffix: proj_tri_3d_mdx_2_faces
751: requires: ctetgen
752: 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
753: filter: grep -v marker | grep -v atomic | grep -v usage
755: test:
756: suffix: proj_tri_3d_mdx_5P_2_faces
757: requires: ctetgen
758: 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
759: filter: grep -v marker | grep -v atomic | grep -v usage
761: test:
762: suffix: proj_quad_lsqr_scale
763: requires: !complex
764: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
765: -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
766: -particlesPerCell 200 \
767: -ptof_pc_type lu \
768: -ftop_ksp_rtol 1e-17 -ftop_ksp_type lsqr -ftop_pc_type none
770: test:
771: suffix: proj_quad_lsqr_prec_scale
772: requires: !complex
773: args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 \
774: -petscspace_degree 2 -petscfe_default_quadrature_order 3 \
775: -particlesPerCell 200 \
776: -ptof_pc_type lu \
777: -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
778: test:
779: suffix: proj_shape_linear_tri_2d
780: requires: ctetgen triangle
781: args: -shape -dm_plex_box_faces 2,2\
782: -dm_view -sw_view\
783: -petscspace_degree 1 -petscfe_default_quadrature_order 1\
784: -pc_type lu
785: test:
786: suffix: proj_shape_linear_quad_2d
787: requires: ctetgen triangle
788: args: -shape -dm_plex_simplex 0 -dm_plex_box_faces 2,2\
789: -dm_view -sw_view\
790: -petscspace_degree 1 -petscfe_default_quadrature_order 1\
791: -pc_type lu
792: test:
793: suffix: proj_shape_linear_tri_3d
794: requires: ctetgen triangle
795: args: -shape -dm_plex_dim 3 -dm_plex_box_faces 2,2,2\
796: -dm_view -sw_view\
797: -petscspace_degree 1 -petscfe_default_quadrature_order 1\
798: -pc_type lu
799: test:
800: suffix: proj_shape_linear_quad_3d
801: requires: ctetgen triangle
802: args: -shape -dm_plex_dim 3 -dm_plex_simplex 0\
803: -dm_plex_box_faces 2,2,2\
804: -dm_view -sw_view\
805: -petscspace_degree 1 -petscfe_default_quadrature_order 1\
806: -pc_type lu
807: TEST*/