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;
87: PetscTabulation T;
88: Vec coorlocal;
89: PetscSection coordSection;
90: PetscScalar *elcoor = NULL;
91: PetscReal *swarm_coor;
92: PetscInt *swarm_cellid;
93: const PetscReal *xiq;
94: PetscQuadrature quadrature;
95: PetscFE fe, feRef;
96: PetscBool is_simplex;
98: PetscFunctionBegin;
99: PetscCall(DMGetDimension(dmc, &dim));
100: is_simplex = PETSC_FALSE;
101: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
102: PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));
103: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
105: PetscCall(private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe));
107: for (r = 0; r < nsub; r++) {
108: PetscCall(PetscFERefine(fe, &feRef));
109: PetscCall(PetscFECopyQuadrature(feRef, fe));
110: PetscCall(PetscFEDestroy(&feRef));
111: }
113: PetscCall(PetscFEGetQuadrature(fe, &quadrature));
114: PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL));
115: PetscCall(PetscFEGetDimension(fe, &nbasis));
116: PetscCall(PetscFEGetCellTabulation(fe, 1, &T));
118: /* 0->cell, 1->edge, 2->vert */
119: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
120: nel = pe - ps;
122: PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1));
123: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
124: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
126: PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
127: PetscCall(DMGetCoordinateSection(dmc, &coordSection));
129: pcnt = 0;
130: for (e = 0; e < nel; e++) {
131: PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
133: for (q = 0; q < npoints_q; q++) {
134: for (d = 0; d < dim; d++) {
135: swarm_coor[dim * pcnt + d] = 0.0;
136: for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][q * nbasis + k] * PetscRealPart(elcoor[dim * k + d]);
137: }
138: swarm_cellid[pcnt] = e;
139: pcnt++;
140: }
141: PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
142: }
143: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
144: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
146: PetscCall(PetscFEDestroy(&fe));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm, DM dmc, PetscInt npoints)
151: {
152: PetscInt dim;
153: PetscInt ii, jj, q, npoints_q, e, nel, npe, pcnt, ps, pe, d, k, nfaces;
154: PetscReal *xi, ds, ds2;
155: PetscReal **basis;
156: Vec coorlocal;
157: PetscSection coordSection;
158: PetscScalar *elcoor = NULL;
159: PetscReal *swarm_coor;
160: PetscInt *swarm_cellid;
161: PetscBool is_simplex;
163: PetscFunctionBegin;
164: PetscCall(DMGetDimension(dmc, &dim));
165: PetscCheck(dim == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only 2D is supported");
166: is_simplex = PETSC_FALSE;
167: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
168: PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));
169: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
170: PetscCheck(is_simplex, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only the simplex is supported");
172: PetscCall(PetscMalloc1(dim * npoints * npoints, &xi));
173: pcnt = 0;
174: ds = 1.0 / (PetscReal)(npoints - 1);
175: ds2 = 1.0 / (PetscReal)npoints;
176: for (jj = 0; jj < npoints; jj++) {
177: for (ii = 0; ii < npoints - jj; ii++) {
178: xi[dim * pcnt + 0] = ii * ds;
179: xi[dim * pcnt + 1] = jj * ds;
181: xi[dim * pcnt + 0] *= (1.0 - 1.2 * ds2);
182: xi[dim * pcnt + 1] *= (1.0 - 1.2 * ds2);
184: xi[dim * pcnt + 0] += 0.35 * ds2;
185: xi[dim * pcnt + 1] += 0.35 * ds2;
186: pcnt++;
187: }
188: }
189: npoints_q = pcnt;
191: npe = 3; /* nodes per element (triangle) */
192: PetscCall(PetscMalloc1(npoints_q, &basis));
193: for (q = 0; q < npoints_q; q++) {
194: PetscCall(PetscMalloc1(npe, &basis[q]));
196: basis[q][0] = 1.0 - xi[dim * q + 0] - xi[dim * q + 1];
197: basis[q][1] = xi[dim * q + 0];
198: basis[q][2] = xi[dim * q + 1];
199: }
201: /* 0->cell, 1->edge, 2->vert */
202: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
203: nel = pe - ps;
205: PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1));
206: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
207: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
209: PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
210: PetscCall(DMGetCoordinateSection(dmc, &coordSection));
212: pcnt = 0;
213: for (e = 0; e < nel; e++) {
214: PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor));
216: for (q = 0; q < npoints_q; q++) {
217: for (d = 0; d < dim; d++) {
218: swarm_coor[dim * pcnt + d] = 0.0;
219: for (k = 0; k < npe; k++) swarm_coor[dim * pcnt + d] += basis[q][k] * PetscRealPart(elcoor[dim * k + d]);
220: }
221: swarm_cellid[pcnt] = e;
222: pcnt++;
223: }
224: PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor));
225: }
226: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
227: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
229: PetscCall(PetscFree(xi));
230: for (q = 0; q < npoints_q; q++) PetscCall(PetscFree(basis[q]));
231: PetscCall(PetscFree(basis));
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param)
236: {
237: PetscInt dim;
239: PetscFunctionBegin;
240: PetscCall(DMGetDimension(celldm, &dim));
241: switch (layout) {
242: case DMSWARMPIC_LAYOUT_REGULAR:
243: PetscCheck(dim != 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No 3D support for REGULAR+PLEX");
244: PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm, celldm, layout_param));
245: break;
246: case DMSWARMPIC_LAYOUT_GAUSS: {
247: PetscQuadrature quad, facequad;
248: const PetscReal *xi;
249: DMPolytopeType ct;
250: PetscInt cStart, Nq;
252: PetscCall(DMPlexGetHeightStratum(celldm, 0, &cStart, NULL));
253: PetscCall(DMPlexGetCellType(celldm, cStart, &ct));
254: PetscCall(PetscDTCreateDefaultQuadrature(ct, layout_param, &quad, &facequad));
255: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &xi, NULL));
256: PetscCall(private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm, celldm, Nq, (PetscReal *)xi));
257: PetscCall(PetscQuadratureDestroy(&quad));
258: PetscCall(PetscQuadratureDestroy(&facequad));
259: } break;
260: case DMSWARMPIC_LAYOUT_SUBDIVISION:
261: PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm, celldm, layout_param));
262: break;
263: }
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm, DM dmc, PetscInt npoints, PetscReal xi[])
268: {
269: PetscBool is_simplex, is_tensorcell;
270: PetscInt dim, nfaces, ps, pe, p, d, nbasis, pcnt, e, k, nel;
271: PetscFE fe;
272: PetscQuadrature quadrature;
273: PetscTabulation T;
274: PetscReal *xiq;
275: Vec coorlocal;
276: PetscSection coordSection;
277: PetscScalar *elcoor = NULL;
278: PetscReal *swarm_coor;
279: PetscInt *swarm_cellid;
281: PetscFunctionBegin;
282: PetscCall(DMGetDimension(dmc, &dim));
284: is_simplex = PETSC_FALSE;
285: is_tensorcell = PETSC_FALSE;
286: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
287: PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));
289: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
291: switch (dim) {
292: case 2:
293: if (nfaces == 4) is_tensorcell = PETSC_TRUE;
294: break;
295: case 3:
296: if (nfaces == 6) is_tensorcell = PETSC_TRUE;
297: break;
298: default:
299: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for 2D, 3D");
300: }
302: /* check points provided fail inside the reference cell */
303: if (is_simplex) {
304: for (p = 0; p < npoints; p++) {
305: PetscReal sum;
306: 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");
307: sum = 0.0;
308: for (d = 0; d < dim; d++) sum += xi[dim * p + d];
309: PetscCheck(sum <= 0.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the simplex domain");
310: }
311: } else if (is_tensorcell) {
312: for (p = 0; p < npoints; p++) {
313: 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");
314: }
315: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for d-simplex and d-tensorcell");
317: PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)dm), &quadrature));
318: PetscCall(PetscMalloc1(npoints * dim, &xiq));
319: PetscCall(PetscArraycpy(xiq, xi, npoints * dim));
320: PetscCall(PetscQuadratureSetData(quadrature, dim, 1, npoints, (const PetscReal *)xiq, NULL));
321: PetscCall(private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe));
322: PetscCall(PetscFESetQuadrature(fe, quadrature));
323: PetscCall(PetscFEGetDimension(fe, &nbasis));
324: PetscCall(PetscFEGetCellTabulation(fe, 1, &T));
326: /* for each cell, interpolate coordaintes and insert the interpolated points coordinates into swarm */
327: /* 0->cell, 1->edge, 2->vert */
328: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
329: nel = pe - ps;
331: PetscCall(DMSwarmSetLocalSizes(dm, npoints * nel, -1));
332: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
333: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
335: PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
336: PetscCall(DMGetCoordinateSection(dmc, &coordSection));
338: pcnt = 0;
339: for (e = 0; e < nel; e++) {
340: PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
342: for (p = 0; p < npoints; p++) {
343: for (d = 0; d < dim; d++) {
344: swarm_coor[dim * pcnt + d] = 0.0;
345: for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][p * nbasis + k] * PetscRealPart(elcoor[dim * k + d]);
346: }
347: swarm_cellid[pcnt] = e;
348: pcnt++;
349: }
350: PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
351: }
352: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
353: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
355: PetscCall(PetscQuadratureDestroy(&quadrature));
356: PetscCall(PetscFEDestroy(&fe));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }