Actual source code: swarmpic_plex.c
1: #include <petscdm.h>
2: #include <petscdmplex.h>
3: #include <petscdmswarm.h>
4: #include "../src/dm/impls/swarm/data_bucket.h"
6: PetscBool SwarmProjcite = PETSC_FALSE;
7: const char SwarmProjCitation[] = "@article{PusztayKnepleyAdams2022,\n"
8: "title = {Conservative Projection Between FEM and Particle Bases},\n"
9: "author = {Joseph V. Pusztay and Matthew G. Knepley and Mark F. Adams},\n"
10: "journal = {SIAM Journal on Scientific Computing},\n"
11: "volume = {44},\n"
12: "number = {4},\n"
13: "pages = {C310--C319},\n"
14: "doi = {10.1137/21M145407},\n"
15: "year = {2022}\n}\n";
17: PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM, DM, PetscInt, PetscReal *xi);
19: static PetscErrorCode private_PetscFECreateDefault_scalar_pk1(DM dm, PetscInt dim, PetscBool isSimplex, PetscInt qorder, PetscFE *fem)
20: {
21: const PetscInt Nc = 1;
22: PetscQuadrature q, fq;
23: DM K;
24: PetscSpace P;
25: PetscDualSpace Q;
26: PetscInt order, quadPointsPerEdge;
27: PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
29: PetscFunctionBegin;
30: /* Create space */
31: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)dm), &P));
32: /* PetscCall(PetscObjectSetOptionsPrefix((PetscObject) P, prefix)); */
33: PetscCall(PetscSpacePolynomialSetTensor(P, tensor));
34: /* PetscCall(PetscSpaceSetFromOptions(P)); */
35: PetscCall(PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL));
36: PetscCall(PetscSpaceSetDegree(P, 1, PETSC_DETERMINE));
37: PetscCall(PetscSpaceSetNumComponents(P, Nc));
38: PetscCall(PetscSpaceSetNumVariables(P, dim));
39: PetscCall(PetscSpaceSetUp(P));
40: PetscCall(PetscSpaceGetDegree(P, &order, NULL));
41: PetscCall(PetscSpacePolynomialGetTensor(P, &tensor));
42: /* Create dual space */
43: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)dm), &Q));
44: PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE));
45: /* PetscCall(PetscObjectSetOptionsPrefix((PetscObject) Q, prefix)); */
46: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K));
47: PetscCall(PetscDualSpaceSetDM(Q, K));
48: PetscCall(DMDestroy(&K));
49: PetscCall(PetscDualSpaceSetNumComponents(Q, Nc));
50: PetscCall(PetscDualSpaceSetOrder(Q, order));
51: PetscCall(PetscDualSpaceLagrangeSetTensor(Q, tensor));
52: /* PetscCall(PetscDualSpaceSetFromOptions(Q)); */
53: PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE));
54: PetscCall(PetscDualSpaceSetUp(Q));
55: /* Create element */
56: PetscCall(PetscFECreate(PetscObjectComm((PetscObject)dm), fem));
57: /* PetscCall(PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix)); */
58: /* PetscCall(PetscFESetFromOptions(*fem)); */
59: PetscCall(PetscFESetType(*fem, PETSCFEBASIC));
60: PetscCall(PetscFESetBasisSpace(*fem, P));
61: PetscCall(PetscFESetDualSpace(*fem, Q));
62: PetscCall(PetscFESetNumComponents(*fem, Nc));
63: PetscCall(PetscFESetUp(*fem));
64: PetscCall(PetscSpaceDestroy(&P));
65: PetscCall(PetscDualSpaceDestroy(&Q));
66: /* Create quadrature (with specified order if given) */
67: qorder = qorder >= 0 ? qorder : order;
68: quadPointsPerEdge = PetscMax(qorder + 1, 1);
69: if (isSimplex) {
70: PetscCall(PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q));
71: PetscCall(PetscDTStroudConicalQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, &fq));
72: } else {
73: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q));
74: PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, &fq));
75: }
76: PetscCall(PetscFESetQuadrature(*fem, q));
77: PetscCall(PetscFESetFaceQuadrature(*fem, fq));
78: PetscCall(PetscQuadratureDestroy(&q));
79: PetscCall(PetscQuadratureDestroy(&fq));
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: static PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(DM dm, DM dmc, PetscInt nsub)
84: {
85: PetscInt dim, nfaces, nbasis;
86: PetscInt q, npoints_q, e, nel, pcnt, ps, pe, d, k, r, Nfc;
87: DMSwarmCellDM celldm;
88: PetscTabulation T;
89: Vec coorlocal;
90: PetscSection coordSection;
91: PetscScalar *elcoor = NULL;
92: PetscReal *swarm_coor;
93: PetscInt *swarm_cellid;
94: const PetscReal *xiq;
95: PetscQuadrature quadrature;
96: PetscFE fe, feRef;
97: PetscBool is_simplex;
98: const char **coordFields, *cellid;
100: PetscFunctionBegin;
101: PetscCall(DMGetDimension(dmc, &dim));
102: is_simplex = PETSC_FALSE;
103: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
104: PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));
105: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
107: PetscCall(private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe));
109: for (r = 0; r < nsub; r++) {
110: PetscCall(PetscFERefine(fe, &feRef));
111: PetscCall(PetscFECopyQuadrature(feRef, fe));
112: PetscCall(PetscFEDestroy(&feRef));
113: }
115: PetscCall(PetscFEGetQuadrature(fe, &quadrature));
116: PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL));
117: PetscCall(PetscFEGetDimension(fe, &nbasis));
118: PetscCall(PetscFEGetCellTabulation(fe, 1, &T));
120: /* 0->cell, 1->edge, 2->vert */
121: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
122: nel = pe - ps;
124: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
125: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
126: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
127: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
129: PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1));
130: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
131: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
133: PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
134: PetscCall(DMGetCoordinateSection(dmc, &coordSection));
136: pcnt = 0;
137: for (e = 0; e < nel; e++) {
138: PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
140: for (q = 0; q < npoints_q; q++) {
141: for (d = 0; d < dim; d++) {
142: swarm_coor[dim * pcnt + d] = 0.0;
143: for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][q * nbasis + k] * PetscRealPart(elcoor[dim * k + d]);
144: }
145: swarm_cellid[pcnt] = e;
146: pcnt++;
147: }
148: PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
149: }
150: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
151: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
153: PetscCall(PetscFEDestroy(&fe));
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: static PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm, DM dmc, PetscInt npoints)
158: {
159: PetscInt dim;
160: PetscInt ii, jj, q, npoints_q, e, nel, npe, pcnt, ps, pe, d, k, nfaces, Nfc;
161: PetscReal *xi, ds, ds2;
162: PetscReal **basis;
163: DMSwarmCellDM celldm;
164: Vec coorlocal;
165: PetscSection coordSection;
166: PetscScalar *elcoor = NULL;
167: PetscReal *swarm_coor;
168: PetscInt *swarm_cellid;
169: PetscBool is_simplex;
170: const char **coordFields, *cellid;
172: PetscFunctionBegin;
173: PetscCall(DMGetDimension(dmc, &dim));
174: PetscCheck(dim == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only 2D is supported");
175: is_simplex = PETSC_FALSE;
176: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
177: PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));
178: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
179: PetscCheck(is_simplex, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only the simplex is supported");
181: PetscCall(PetscMalloc1(dim * npoints * npoints, &xi));
182: pcnt = 0;
183: ds = 1.0 / (PetscReal)(npoints - 1);
184: ds2 = 1.0 / (PetscReal)npoints;
185: for (jj = 0; jj < npoints; jj++) {
186: for (ii = 0; ii < npoints - jj; ii++) {
187: xi[dim * pcnt + 0] = ii * ds;
188: xi[dim * pcnt + 1] = jj * ds;
190: xi[dim * pcnt + 0] *= (1.0 - 1.2 * ds2);
191: xi[dim * pcnt + 1] *= (1.0 - 1.2 * ds2);
193: xi[dim * pcnt + 0] += 0.35 * ds2;
194: xi[dim * pcnt + 1] += 0.35 * ds2;
195: pcnt++;
196: }
197: }
198: npoints_q = pcnt;
200: npe = 3; /* nodes per element (triangle) */
201: PetscCall(PetscMalloc1(npoints_q, &basis));
202: for (q = 0; q < npoints_q; q++) {
203: PetscCall(PetscMalloc1(npe, &basis[q]));
205: basis[q][0] = 1.0 - xi[dim * q + 0] - xi[dim * q + 1];
206: basis[q][1] = xi[dim * q + 0];
207: basis[q][2] = xi[dim * q + 1];
208: }
210: /* 0->cell, 1->edge, 2->vert */
211: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
212: nel = pe - ps;
214: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
215: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
216: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
217: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
219: PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1));
220: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
221: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
223: PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
224: PetscCall(DMGetCoordinateSection(dmc, &coordSection));
226: pcnt = 0;
227: for (e = 0; e < nel; e++) {
228: PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor));
230: for (q = 0; q < npoints_q; q++) {
231: for (d = 0; d < dim; d++) {
232: swarm_coor[dim * pcnt + d] = 0.0;
233: for (k = 0; k < npe; k++) swarm_coor[dim * pcnt + d] += basis[q][k] * PetscRealPart(elcoor[dim * k + d]);
234: }
235: swarm_cellid[pcnt] = e;
236: pcnt++;
237: }
238: PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor));
239: }
240: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
241: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
243: PetscCall(PetscFree(xi));
244: for (q = 0; q < npoints_q; q++) PetscCall(PetscFree(basis[q]));
245: PetscCall(PetscFree(basis));
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param)
250: {
251: PetscInt dim;
253: PetscFunctionBegin;
254: PetscCall(DMGetDimension(celldm, &dim));
255: switch (layout) {
256: case DMSWARMPIC_LAYOUT_REGULAR:
257: PetscCheck(dim != 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No 3D support for REGULAR+PLEX");
258: PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm, celldm, layout_param));
259: break;
260: case DMSWARMPIC_LAYOUT_GAUSS: {
261: PetscQuadrature quad, facequad;
262: const PetscReal *xi;
263: DMPolytopeType ct;
264: PetscInt cStart, Nq;
266: PetscCall(DMPlexGetHeightStratum(celldm, 0, &cStart, NULL));
267: PetscCall(DMPlexGetCellType(celldm, cStart, &ct));
268: PetscCall(PetscDTCreateDefaultQuadrature(ct, layout_param, &quad, &facequad));
269: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &xi, NULL));
270: PetscCall(private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm, celldm, Nq, (PetscReal *)xi));
271: PetscCall(PetscQuadratureDestroy(&quad));
272: PetscCall(PetscQuadratureDestroy(&facequad));
273: } break;
274: case DMSWARMPIC_LAYOUT_SUBDIVISION:
275: PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm, celldm, layout_param));
276: break;
277: }
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm, DM dmc, PetscInt npoints, PetscReal xi[])
282: {
283: PetscBool is_simplex, is_tensorcell;
284: PetscInt dim, nfaces, ps, pe, p, d, nbasis, pcnt, e, k, nel, Nfc;
285: DMSwarmCellDM celldm;
286: PetscFE fe;
287: PetscQuadrature quadrature;
288: PetscTabulation T;
289: PetscReal *xiq;
290: Vec coorlocal;
291: PetscSection coordSection;
292: PetscScalar *elcoor = NULL;
293: PetscReal *swarm_coor;
294: PetscInt *swarm_cellid;
295: const char **coordFields, *cellid;
297: PetscFunctionBegin;
298: PetscCall(DMGetDimension(dmc, &dim));
300: is_simplex = PETSC_FALSE;
301: is_tensorcell = PETSC_FALSE;
302: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
303: PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));
305: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
307: switch (dim) {
308: case 2:
309: if (nfaces == 4) is_tensorcell = PETSC_TRUE;
310: break;
311: case 3:
312: if (nfaces == 6) is_tensorcell = PETSC_TRUE;
313: break;
314: default:
315: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for 2D, 3D");
316: }
318: /* check points provided fail inside the reference cell */
319: if (is_simplex) {
320: for (p = 0; p < npoints; p++) {
321: PetscReal sum;
322: for (d = 0; d < dim; d++) PetscCheck(xi[dim * p + d] >= -1.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the simplex domain");
323: sum = 0.0;
324: for (d = 0; d < dim; d++) sum += xi[dim * p + d];
325: PetscCheck(sum <= 0.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the simplex domain");
326: }
327: } else if (is_tensorcell) {
328: for (p = 0; p < npoints; p++) {
329: for (d = 0; d < dim; d++) PetscCheck(PetscAbsReal(xi[dim * p + d]) <= 1.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the tensor domain [-1,1]^d");
330: }
331: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for d-simplex and d-tensorcell");
333: PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)dm), &quadrature));
334: PetscCall(PetscMalloc1(npoints * dim, &xiq));
335: PetscCall(PetscArraycpy(xiq, xi, npoints * dim));
336: PetscCall(PetscQuadratureSetData(quadrature, dim, 1, npoints, (const PetscReal *)xiq, NULL));
337: PetscCall(private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe));
338: PetscCall(PetscFESetQuadrature(fe, quadrature));
339: PetscCall(PetscFEGetDimension(fe, &nbasis));
340: PetscCall(PetscFEGetCellTabulation(fe, 1, &T));
342: /* for each cell, interpolate coordaintes and insert the interpolated points coordinates into swarm */
343: /* 0->cell, 1->edge, 2->vert */
344: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
345: nel = pe - ps;
347: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
348: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
349: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
350: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
352: PetscCall(DMSwarmSetLocalSizes(dm, npoints * nel, -1));
353: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
354: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
356: PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
357: PetscCall(DMGetCoordinateSection(dmc, &coordSection));
359: pcnt = 0;
360: for (e = 0; e < nel; e++) {
361: PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
363: for (p = 0; p < npoints; p++) {
364: for (d = 0; d < dim; d++) {
365: swarm_coor[dim * pcnt + d] = 0.0;
366: for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][p * nbasis + k] * PetscRealPart(elcoor[dim * k + d]);
367: }
368: swarm_cellid[pcnt] = e;
369: pcnt++;
370: }
371: PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
372: }
373: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
374: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
376: PetscCall(PetscQuadratureDestroy(&quadrature));
377: PetscCall(PetscFEDestroy(&fe));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }