1 #include <petscdm.h>
2 #include <petscdmplex.h>
3 #include <petscdmswarm.h>
4 #include "../src/dm/impls/swarm/data_bucket.h"
5
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";
16
17 PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM, DM, PetscInt, PetscReal *xi);
18
private_PetscFECreateDefault_scalar_pk1(DM dm,PetscInt dim,PetscBool isSimplex,PetscInt qorder,PetscFE * fem)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;
28
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 }
82
private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(DM dm,DM dmc,PetscInt nsub)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;
99
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;
106
107 PetscCall(private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe));
108
109 for (r = 0; r < nsub; r++) {
110 PetscCall(PetscFERefine(fe, &feRef));
111 PetscCall(PetscFECopyQuadrature(feRef, fe));
112 PetscCall(PetscFEDestroy(&feRef));
113 }
114
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));
119
120 /* 0->cell, 1->edge, 2->vert */
121 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
122 nel = pe - ps;
123
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));
128
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));
132
133 PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
134 PetscCall(DMGetCoordinateSection(dmc, &coordSection));
135
136 pcnt = 0;
137 for (e = 0; e < nel; e++) {
138 PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
139
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));
152
153 PetscCall(PetscFEDestroy(&fe));
154 PetscFunctionReturn(PETSC_SUCCESS);
155 }
156
private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm,DM dmc,PetscInt npoints)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;
171
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");
180
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;
189
190 xi[dim * pcnt + 0] *= (1.0 - 1.2 * ds2);
191 xi[dim * pcnt + 1] *= (1.0 - 1.2 * ds2);
192
193 xi[dim * pcnt + 0] += 0.35 * ds2;
194 xi[dim * pcnt + 1] += 0.35 * ds2;
195 pcnt++;
196 }
197 }
198 npoints_q = pcnt;
199
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]));
204
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 }
209
210 /* 0->cell, 1->edge, 2->vert */
211 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
212 nel = pe - ps;
213
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));
218
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));
222
223 PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
224 PetscCall(DMGetCoordinateSection(dmc, &coordSection));
225
226 pcnt = 0;
227 for (e = 0; e < nel; e++) {
228 PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor));
229
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));
242
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 }
248
private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)249 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param)
250 {
251 PetscInt dim;
252
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;
265
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 }
280
private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm,DM dmc,PetscInt npoints,PetscReal xi[])281 PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm, DM dmc, PetscInt npoints, PetscReal xi[])
282 {
283 PetscBool is_simplex, is_tensorcell;
284 PetscInt dim, ps, pe, nel, nfaces, Nfc;
285 DMSwarmCellDM celldm;
286 PetscReal *swarm_coor;
287 PetscInt *swarm_cellid;
288 const char **coordFields, *cellid;
289
290 PetscFunctionBegin;
291 PetscCall(DMGetDimension(dmc, &dim));
292
293 is_simplex = PETSC_FALSE;
294 is_tensorcell = PETSC_FALSE;
295 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
296 PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));
297
298 if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
299
300 switch (dim) {
301 case 2:
302 if (nfaces == 4) is_tensorcell = PETSC_TRUE;
303 break;
304 case 3:
305 if (nfaces == 6) is_tensorcell = PETSC_TRUE;
306 break;
307 default:
308 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for 2D, 3D");
309 }
310
311 /* check points provided fail inside the reference cell */
312 if (is_simplex) {
313 for (PetscInt p = 0; p < npoints; p++) {
314 PetscReal sum;
315 for (PetscInt 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");
316 sum = 0.0;
317 for (PetscInt d = 0; d < dim; d++) sum += xi[dim * p + d];
318 PetscCheck(sum <= 0.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the simplex domain");
319 }
320 } else if (is_tensorcell) {
321 for (PetscInt p = 0; p < npoints; p++) {
322 for (PetscInt 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");
323 }
324 } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for d-simplex and d-tensorcell");
325
326 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
327 nel = pe - ps;
328
329 PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
330 PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
331 PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
332 PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
333
334 PetscCall(DMSwarmSetLocalSizes(dm, npoints * nel, PETSC_DECIDE));
335 PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
336 PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
337
338 // Use DMPlexReferenceToCoordinates so that arbitrary discretizations work
339 for (PetscInt e = 0; e < nel; e++) {
340 PetscCall(DMPlexReferenceToCoordinates(dmc, e + ps, npoints, xi, &swarm_coor[npoints * dim * e]));
341 for (PetscInt p = 0; p < npoints; p++) swarm_cellid[e * npoints + p] = e + ps;
342 }
343 PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
344 PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
345 PetscFunctionReturn(PETSC_SUCCESS);
346 }
347