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: }