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: PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM, DM, PetscInt, PetscReal *xi);
8: static PetscErrorCode private_PetscFECreateDefault_scalar_pk1(DM dm, PetscInt dim, PetscBool isSimplex, PetscInt qorder, PetscFE *fem)
9: {
10: const PetscInt Nc = 1;
11: PetscQuadrature q, fq;
12: DM K;
13: PetscSpace P;
14: PetscDualSpace Q;
15: PetscInt order, quadPointsPerEdge;
16: PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
18: PetscFunctionBegin;
19: /* Create space */
20: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)dm), &P));
21: /* PetscCall(PetscObjectSetOptionsPrefix((PetscObject) P, prefix)); */
22: PetscCall(PetscSpacePolynomialSetTensor(P, tensor));
23: /* PetscCall(PetscSpaceSetFromOptions(P)); */
24: PetscCall(PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL));
25: PetscCall(PetscSpaceSetDegree(P, 1, PETSC_DETERMINE));
26: PetscCall(PetscSpaceSetNumComponents(P, Nc));
27: PetscCall(PetscSpaceSetNumVariables(P, dim));
28: PetscCall(PetscSpaceSetUp(P));
29: PetscCall(PetscSpaceGetDegree(P, &order, NULL));
30: PetscCall(PetscSpacePolynomialGetTensor(P, &tensor));
31: /* Create dual space */
32: PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)dm), &Q));
33: PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE));
34: /* PetscCall(PetscObjectSetOptionsPrefix((PetscObject) Q, prefix)); */
35: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K));
36: PetscCall(PetscDualSpaceSetDM(Q, K));
37: PetscCall(DMDestroy(&K));
38: PetscCall(PetscDualSpaceSetNumComponents(Q, Nc));
39: PetscCall(PetscDualSpaceSetOrder(Q, order));
40: PetscCall(PetscDualSpaceLagrangeSetTensor(Q, tensor));
41: /* PetscCall(PetscDualSpaceSetFromOptions(Q)); */
42: PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE));
43: PetscCall(PetscDualSpaceSetUp(Q));
44: /* Create element */
45: PetscCall(PetscFECreate(PetscObjectComm((PetscObject)dm), fem));
46: /* PetscCall(PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix)); */
47: /* PetscCall(PetscFESetFromOptions(*fem)); */
48: PetscCall(PetscFESetType(*fem, PETSCFEBASIC));
49: PetscCall(PetscFESetBasisSpace(*fem, P));
50: PetscCall(PetscFESetDualSpace(*fem, Q));
51: PetscCall(PetscFESetNumComponents(*fem, Nc));
52: PetscCall(PetscFESetUp(*fem));
53: PetscCall(PetscSpaceDestroy(&P));
54: PetscCall(PetscDualSpaceDestroy(&Q));
55: /* Create quadrature (with specified order if given) */
56: qorder = qorder >= 0 ? qorder : order;
57: quadPointsPerEdge = PetscMax(qorder + 1, 1);
58: if (isSimplex) {
59: PetscCall(PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q));
60: PetscCall(PetscDTStroudConicalQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, &fq));
61: } else {
62: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q));
63: PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, &fq));
64: }
65: PetscCall(PetscFESetQuadrature(*fem, q));
66: PetscCall(PetscFESetFaceQuadrature(*fem, fq));
67: PetscCall(PetscQuadratureDestroy(&q));
68: PetscCall(PetscQuadratureDestroy(&fq));
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
72: PetscErrorCode subdivide_triangle(PetscReal v1[2], PetscReal v2[2], PetscReal v3[2], PetscInt depth, PetscInt max, PetscReal xi[], PetscInt *np)
73: {
74: PetscReal v12[2], v23[2], v31[2];
75: PetscInt i;
77: PetscFunctionBegin;
78: if (depth == max) {
79: PetscReal cx[2];
81: cx[0] = (v1[0] + v2[0] + v3[0]) / 3.0;
82: cx[1] = (v1[1] + v2[1] + v3[1]) / 3.0;
84: xi[2 * (*np) + 0] = cx[0];
85: xi[2 * (*np) + 1] = cx[1];
86: (*np)++;
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: /* calculate midpoints of each side */
91: for (i = 0; i < 2; i++) {
92: v12[i] = (v1[i] + v2[i]) / 2.0;
93: v23[i] = (v2[i] + v3[i]) / 2.0;
94: v31[i] = (v3[i] + v1[i]) / 2.0;
95: }
97: /* recursively subdivide new triangles */
98: PetscCall(subdivide_triangle(v1, v12, v31, depth + 1, max, xi, np));
99: PetscCall(subdivide_triangle(v2, v23, v12, depth + 1, max, xi, np));
100: PetscCall(subdivide_triangle(v3, v31, v23, depth + 1, max, xi, np));
101: PetscCall(subdivide_triangle(v12, v23, v31, depth + 1, max, xi, np));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_SubDivide(DM dm, DM dmc, PetscInt nsub)
106: {
107: const PetscInt dim = 2;
108: PetscInt q, npoints_q, e, nel, npe, pcnt, ps, pe, d, k, depth;
109: PetscReal *xi;
110: PetscReal **basis;
111: Vec coorlocal;
112: PetscSection coordSection;
113: PetscScalar *elcoor = NULL;
114: PetscReal *swarm_coor;
115: PetscInt *swarm_cellid;
116: PetscReal v1[2], v2[2], v3[2];
118: PetscFunctionBegin;
119: npoints_q = 1;
120: for (d = 0; d < nsub; d++) npoints_q *= 4;
121: PetscCall(PetscMalloc1(dim * npoints_q, &xi));
123: v1[0] = 0.0;
124: v1[1] = 0.0;
125: v2[0] = 1.0;
126: v2[1] = 0.0;
127: v3[0] = 0.0;
128: v3[1] = 1.0;
129: depth = 0;
130: pcnt = 0;
131: PetscCall(subdivide_triangle(v1, v2, v3, depth, nsub, xi, &pcnt));
133: npe = 3; /* nodes per element (triangle) */
134: PetscCall(PetscMalloc1(npoints_q, &basis));
135: for (q = 0; q < npoints_q; q++) {
136: PetscCall(PetscMalloc1(npe, &basis[q]));
138: basis[q][0] = 1.0 - xi[dim * q + 0] - xi[dim * q + 1];
139: basis[q][1] = xi[dim * q + 0];
140: basis[q][2] = xi[dim * q + 1];
141: }
143: /* 0->cell, 1->edge, 2->vert */
144: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
145: nel = pe - ps;
147: PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1));
148: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
149: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
151: PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
152: PetscCall(DMGetCoordinateSection(dmc, &coordSection));
154: pcnt = 0;
155: for (e = 0; e < nel; e++) {
156: PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor));
158: for (q = 0; q < npoints_q; q++) {
159: for (d = 0; d < dim; d++) {
160: swarm_coor[dim * pcnt + d] = 0.0;
161: for (k = 0; k < npe; k++) swarm_coor[dim * pcnt + d] += basis[q][k] * PetscRealPart(elcoor[dim * k + d]);
162: }
163: swarm_cellid[pcnt] = e;
164: pcnt++;
165: }
166: PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor));
167: }
168: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
169: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
171: PetscCall(PetscFree(xi));
172: for (q = 0; q < npoints_q; q++) PetscCall(PetscFree(basis[q]));
173: PetscCall(PetscFree(basis));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(DM dm, DM dmc, PetscInt nsub)
178: {
179: PetscInt dim, nfaces, nbasis;
180: PetscInt q, npoints_q, e, nel, pcnt, ps, pe, d, k, r;
181: PetscTabulation T;
182: Vec coorlocal;
183: PetscSection coordSection;
184: PetscScalar *elcoor = NULL;
185: PetscReal *swarm_coor;
186: PetscInt *swarm_cellid;
187: const PetscReal *xiq;
188: PetscQuadrature quadrature;
189: PetscFE fe, feRef;
190: PetscBool is_simplex;
192: PetscFunctionBegin;
193: PetscCall(DMGetDimension(dmc, &dim));
194: is_simplex = PETSC_FALSE;
195: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
196: PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));
197: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
199: PetscCall(private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe));
201: for (r = 0; r < nsub; r++) {
202: PetscCall(PetscFERefine(fe, &feRef));
203: PetscCall(PetscFECopyQuadrature(feRef, fe));
204: PetscCall(PetscFEDestroy(&feRef));
205: }
207: PetscCall(PetscFEGetQuadrature(fe, &quadrature));
208: PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL));
209: PetscCall(PetscFEGetDimension(fe, &nbasis));
210: PetscCall(PetscFEGetCellTabulation(fe, 1, &T));
212: /* 0->cell, 1->edge, 2->vert */
213: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
214: nel = pe - ps;
216: PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1));
217: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
218: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
220: PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
221: PetscCall(DMGetCoordinateSection(dmc, &coordSection));
223: pcnt = 0;
224: for (e = 0; e < nel; e++) {
225: PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
227: for (q = 0; q < npoints_q; q++) {
228: for (d = 0; d < dim; d++) {
229: swarm_coor[dim * pcnt + d] = 0.0;
230: for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][q * nbasis + k] * PetscRealPart(elcoor[dim * k + d]);
231: }
232: swarm_cellid[pcnt] = e;
233: pcnt++;
234: }
235: PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
236: }
237: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
238: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
240: PetscCall(PetscFEDestroy(&fe));
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm, DM dmc, PetscInt npoints)
245: {
246: PetscInt dim;
247: PetscInt ii, jj, q, npoints_q, e, nel, npe, pcnt, ps, pe, d, k, nfaces;
248: PetscReal *xi, ds, ds2;
249: PetscReal **basis;
250: Vec coorlocal;
251: PetscSection coordSection;
252: PetscScalar *elcoor = NULL;
253: PetscReal *swarm_coor;
254: PetscInt *swarm_cellid;
255: PetscBool is_simplex;
257: PetscFunctionBegin;
258: PetscCall(DMGetDimension(dmc, &dim));
259: PetscCheck(dim == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only 2D is supported");
260: is_simplex = PETSC_FALSE;
261: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
262: PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));
263: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
264: PetscCheck(is_simplex, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only the simplex is supported");
266: PetscCall(PetscMalloc1(dim * npoints * npoints, &xi));
267: pcnt = 0;
268: ds = 1.0 / ((PetscReal)(npoints - 1));
269: ds2 = 1.0 / ((PetscReal)(npoints));
270: for (jj = 0; jj < npoints; jj++) {
271: for (ii = 0; ii < npoints - jj; ii++) {
272: xi[dim * pcnt + 0] = ii * ds;
273: xi[dim * pcnt + 1] = jj * ds;
275: xi[dim * pcnt + 0] *= (1.0 - 1.2 * ds2);
276: xi[dim * pcnt + 1] *= (1.0 - 1.2 * ds2);
278: xi[dim * pcnt + 0] += 0.35 * ds2;
279: xi[dim * pcnt + 1] += 0.35 * ds2;
280: pcnt++;
281: }
282: }
283: npoints_q = pcnt;
285: npe = 3; /* nodes per element (triangle) */
286: PetscCall(PetscMalloc1(npoints_q, &basis));
287: for (q = 0; q < npoints_q; q++) {
288: PetscCall(PetscMalloc1(npe, &basis[q]));
290: basis[q][0] = 1.0 - xi[dim * q + 0] - xi[dim * q + 1];
291: basis[q][1] = xi[dim * q + 0];
292: basis[q][2] = xi[dim * q + 1];
293: }
295: /* 0->cell, 1->edge, 2->vert */
296: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
297: nel = pe - ps;
299: PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1));
300: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
301: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
303: PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
304: PetscCall(DMGetCoordinateSection(dmc, &coordSection));
306: pcnt = 0;
307: for (e = 0; e < nel; e++) {
308: PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor));
310: for (q = 0; q < npoints_q; q++) {
311: for (d = 0; d < dim; d++) {
312: swarm_coor[dim * pcnt + d] = 0.0;
313: for (k = 0; k < npe; k++) swarm_coor[dim * pcnt + d] += basis[q][k] * PetscRealPart(elcoor[dim * k + d]);
314: }
315: swarm_cellid[pcnt] = e;
316: pcnt++;
317: }
318: PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor));
319: }
320: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
321: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
323: PetscCall(PetscFree(xi));
324: for (q = 0; q < npoints_q; q++) PetscCall(PetscFree(basis[q]));
325: PetscCall(PetscFree(basis));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param)
330: {
331: PetscInt dim;
333: PetscFunctionBegin;
334: PetscCall(DMGetDimension(celldm, &dim));
335: switch (layout) {
336: case DMSWARMPIC_LAYOUT_REGULAR:
337: PetscCheck(dim != 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No 3D support for REGULAR+PLEX");
338: PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm, celldm, layout_param));
339: break;
340: case DMSWARMPIC_LAYOUT_GAUSS: {
341: PetscInt npoints, npoints1, ps, pe, nfaces;
342: const PetscReal *xi;
343: PetscBool is_simplex;
344: PetscQuadrature quadrature;
346: is_simplex = PETSC_FALSE;
347: PetscCall(DMPlexGetHeightStratum(celldm, 0, &ps, &pe));
348: PetscCall(DMPlexGetConeSize(celldm, ps, &nfaces));
349: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
351: npoints1 = layout_param;
352: if (is_simplex) {
353: PetscCall(PetscDTStroudConicalQuadrature(dim, 1, npoints1, -1.0, 1.0, &quadrature));
354: } else {
355: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, npoints1, -1.0, 1.0, &quadrature));
356: }
357: PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints, &xi, NULL));
358: PetscCall(private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm, celldm, npoints, (PetscReal *)xi));
359: PetscCall(PetscQuadratureDestroy(&quadrature));
360: } break;
361: case DMSWARMPIC_LAYOUT_SUBDIVISION:
362: PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm, celldm, layout_param));
363: break;
364: }
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: /*
369: typedef struct {
370: PetscReal x,y;
371: } Point2d;
373: static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s)
374: {
375: PetscFunctionBegin;
376: *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y);
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
379: */
380: /*
381: static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v)
382: {
383: PetscReal s1,s2,s3;
384: PetscBool b1, b2, b3;
386: PetscFunctionBegin;
387: signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f;
388: signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f;
389: signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f;
391: *v = ((b1 == b2) && (b2 == b3));
392: PetscFunctionReturn(PETSC_SUCCESS);
393: }
394: */
395: /*
396: static PetscErrorCode _ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ)
397: {
398: PetscReal x1,y1,x2,y2,x3,y3;
399: PetscReal c,b[2],A[2][2],inv[2][2],detJ,od;
401: PetscFunctionBegin;
402: x1 = coords[2*0+0];
403: x2 = coords[2*1+0];
404: x3 = coords[2*2+0];
406: y1 = coords[2*0+1];
407: y2 = coords[2*1+1];
408: y3 = coords[2*2+1];
410: c = x1 - 0.5*x1 - 0.5*x1 + 0.5*x2 + 0.5*x3;
411: b[0] = xp[0] - c;
412: c = y1 - 0.5*y1 - 0.5*y1 + 0.5*y2 + 0.5*y3;
413: b[1] = xp[1] - c;
415: A[0][0] = -0.5*x1 + 0.5*x2; A[0][1] = -0.5*x1 + 0.5*x3;
416: A[1][0] = -0.5*y1 + 0.5*y2; A[1][1] = -0.5*y1 + 0.5*y3;
418: detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0];
419: *dJ = PetscAbsReal(detJ);
420: od = 1.0/detJ;
422: inv[0][0] = A[1][1] * od;
423: inv[0][1] = -A[0][1] * od;
424: inv[1][0] = -A[1][0] * od;
425: inv[1][1] = A[0][0] * od;
427: xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1];
428: xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1];
429: PetscFunctionReturn(PETSC_SUCCESS);
430: }
431: */
433: static PetscErrorCode ComputeLocalCoordinateAffine2d(PetscReal xp[], PetscScalar coords[], PetscReal xip[], PetscReal *dJ)
434: {
435: PetscReal x1, y1, x2, y2, x3, y3;
436: PetscReal b[2], A[2][2], inv[2][2], detJ, od;
438: PetscFunctionBegin;
439: x1 = PetscRealPart(coords[2 * 0 + 0]);
440: x2 = PetscRealPart(coords[2 * 1 + 0]);
441: x3 = PetscRealPart(coords[2 * 2 + 0]);
443: y1 = PetscRealPart(coords[2 * 0 + 1]);
444: y2 = PetscRealPart(coords[2 * 1 + 1]);
445: y3 = PetscRealPart(coords[2 * 2 + 1]);
447: b[0] = xp[0] - x1;
448: b[1] = xp[1] - y1;
450: A[0][0] = x2 - x1;
451: A[0][1] = x3 - x1;
452: A[1][0] = y2 - y1;
453: A[1][1] = y3 - y1;
455: detJ = A[0][0] * A[1][1] - A[0][1] * A[1][0];
456: *dJ = PetscAbsReal(detJ);
457: od = 1.0 / detJ;
459: inv[0][0] = A[1][1] * od;
460: inv[0][1] = -A[0][1] * od;
461: inv[1][0] = -A[1][0] * od;
462: inv[1][1] = A[0][0] * od;
464: xip[0] = inv[0][0] * b[0] + inv[0][1] * b[1];
465: xip[1] = inv[1][0] * b[0] + inv[1][1] * b[1];
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: PetscErrorCode DMSwarmProjectField_ApproxP1_PLEX_2D(DM swarm, PetscReal *swarm_field, DM dm, Vec v_field)
470: {
471: const PetscReal PLEX_C_EPS = 1.0e-8;
472: Vec v_field_l, denom_l, coor_l, denom;
473: PetscInt k, p, e, npoints;
474: PetscInt *mpfield_cell;
475: PetscReal *mpfield_coor;
476: PetscReal xi_p[2];
477: PetscScalar Ni[3];
478: PetscSection coordSection;
479: PetscScalar *elcoor = NULL;
481: PetscFunctionBegin;
482: PetscCall(VecZeroEntries(v_field));
484: PetscCall(DMGetLocalVector(dm, &v_field_l));
485: PetscCall(DMGetGlobalVector(dm, &denom));
486: PetscCall(DMGetLocalVector(dm, &denom_l));
487: PetscCall(VecZeroEntries(v_field_l));
488: PetscCall(VecZeroEntries(denom));
489: PetscCall(VecZeroEntries(denom_l));
491: PetscCall(DMGetCoordinatesLocal(dm, &coor_l));
492: PetscCall(DMGetCoordinateSection(dm, &coordSection));
494: PetscCall(DMSwarmGetLocalSize(swarm, &npoints));
495: PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
496: PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));
498: for (p = 0; p < npoints; p++) {
499: PetscReal *coor_p, dJ;
500: PetscScalar elfield[3];
501: PetscBool point_located;
503: e = mpfield_cell[p];
504: coor_p = &mpfield_coor[2 * p];
506: PetscCall(DMPlexVecGetClosure(dm, coordSection, coor_l, e, NULL, &elcoor));
508: /*
509: while (!point_located && (failed_counter < 25)) {
510: PetscCall(PointInTriangle(point, coords[0], coords[1], coords[2], &point_located));
511: point.x = coor_p[0];
512: point.y = coor_p[1];
513: point.x += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
514: point.y += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
515: failed_counter++;
516: }
518: if (!point_located) {
519: PetscPrintf(PETSC_COMM_SELF,"Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ") with triangle coords (%1.8e,%1.8e) : (%1.8e,%1.8e) : (%1.8e,%1.8e) in %" PetscInt_FMT " iterations\n",point.x,point.y,e,coords[0].x,coords[0].y,coords[1].x,coords[1].y,coords[2].x,coords[2].y,failed_counter);
520: }
522: PetscCheck(point_located,PETSC_COMM_SELF,PETSC_ERR_SUP,"Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ")",point.x,point.y,e);
523: else {
524: PetscCall(_ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ));
525: xi_p[0] = 0.5*(xi_p[0] + 1.0);
526: xi_p[1] = 0.5*(xi_p[1] + 1.0);
528: PetscPrintf(PETSC_COMM_SELF,"[p=%" PetscInt_FMT "] x(%+1.4e,%+1.4e) -> mapped to element %" PetscInt_FMT " xi(%+1.4e,%+1.4e)\n",p,point.x,point.y,e,xi_p[0],xi_p[1]);
530: }
531: */
533: PetscCall(ComputeLocalCoordinateAffine2d(coor_p, elcoor, xi_p, &dJ));
534: /*
535: PetscPrintf(PETSC_COMM_SELF,"[p=%" PetscInt_FMT "] x(%+1.4e,%+1.4e) -> mapped to element %" PetscInt_FMT " xi(%+1.4e,%+1.4e)\n",p,point.x,point.y,e,xi_p[0],xi_p[1]);
536: */
537: /*
538: point_located = PETSC_TRUE;
539: if (xi_p[0] < 0.0) {
540: if (xi_p[0] > -PLEX_C_EPS) {
541: xi_p[0] = 0.0;
542: } else {
543: point_located = PETSC_FALSE;
544: }
545: }
546: if (xi_p[1] < 0.0) {
547: if (xi_p[1] > -PLEX_C_EPS) {
548: xi_p[1] = 0.0;
549: } else {
550: point_located = PETSC_FALSE;
551: }
552: }
553: if (xi_p[1] > (1.0-xi_p[0])) {
554: if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) {
555: xi_p[1] = 1.0 - xi_p[0];
556: } else {
557: point_located = PETSC_FALSE;
558: }
559: }
560: if (!point_located) {
561: PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]);
562: PetscPrintf(PETSC_COMM_SELF,"[Error] Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ") with triangle coords (%1.8e,%1.8e) : (%1.8e,%1.8e) : (%1.8e,%1.8e)\n",coor_p[0],coor_p[1],e,elcoor[0],elcoor[1],elcoor[2],elcoor[3],elcoor[4],elcoor[5]);
563: }
564: PetscCheck(point_located,PETSC_COMM_SELF,PETSC_ERR_SUP,"Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ")",coor_p[0],coor_p[1],e);
565: */
567: Ni[0] = 1.0 - xi_p[0] - xi_p[1];
568: Ni[1] = xi_p[0];
569: Ni[2] = xi_p[1];
571: point_located = PETSC_TRUE;
572: for (k = 0; k < 3; k++) {
573: if (PetscRealPart(Ni[k]) < -PLEX_C_EPS) point_located = PETSC_FALSE;
574: if (PetscRealPart(Ni[k]) > (1.0 + PLEX_C_EPS)) point_located = PETSC_FALSE;
575: }
576: if (!point_located) {
577: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[Error] xi,eta = %+1.8e, %+1.8e\n", (double)xi_p[0], (double)xi_p[1]));
578: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[Error] Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ") with triangle coords (%1.8e,%1.8e) : (%1.8e,%1.8e) : (%1.8e,%1.8e)\n", (double)coor_p[0], (double)coor_p[1], e, (double)PetscRealPart(elcoor[0]), (double)PetscRealPart(elcoor[1]), (double)PetscRealPart(elcoor[2]), (double)PetscRealPart(elcoor[3]), (double)PetscRealPart(elcoor[4]), (double)PetscRealPart(elcoor[5])));
579: }
580: PetscCheck(point_located, PETSC_COMM_SELF, PETSC_ERR_SUP, "Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ")", (double)coor_p[0], (double)coor_p[1], e);
582: for (k = 0; k < 3; k++) {
583: Ni[k] = Ni[k] * dJ;
584: elfield[k] = Ni[k] * swarm_field[p];
585: }
586: PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coor_l, e, NULL, &elcoor));
588: PetscCall(DMPlexVecSetClosure(dm, NULL, v_field_l, e, elfield, ADD_VALUES));
589: PetscCall(DMPlexVecSetClosure(dm, NULL, denom_l, e, Ni, ADD_VALUES));
590: }
592: PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));
593: PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
595: PetscCall(DMLocalToGlobalBegin(dm, v_field_l, ADD_VALUES, v_field));
596: PetscCall(DMLocalToGlobalEnd(dm, v_field_l, ADD_VALUES, v_field));
597: PetscCall(DMLocalToGlobalBegin(dm, denom_l, ADD_VALUES, denom));
598: PetscCall(DMLocalToGlobalEnd(dm, denom_l, ADD_VALUES, denom));
600: PetscCall(VecPointwiseDivide(v_field, v_field, denom));
602: PetscCall(DMRestoreLocalVector(dm, &v_field_l));
603: PetscCall(DMRestoreLocalVector(dm, &denom_l));
604: PetscCall(DMRestoreGlobalVector(dm, &denom));
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm, DM celldm, PetscInt project_type, PetscInt nfields, DMSwarmDataField dfield[], Vec vecs[])
609: {
610: PetscInt f, dim;
612: PetscFunctionBegin;
613: PetscCall(DMGetDimension(swarm, &dim));
614: switch (dim) {
615: case 2:
616: for (f = 0; f < nfields; f++) {
617: PetscReal *swarm_field;
619: PetscCall(DMSwarmDataFieldGetEntries(dfield[f], (void **)&swarm_field));
620: PetscCall(DMSwarmProjectField_ApproxP1_PLEX_2D(swarm, swarm_field, celldm, vecs[f]));
621: }
622: break;
623: case 3:
624: SETERRQ(PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "No support for 3D");
625: default:
626: break;
627: }
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm, DM dmc, PetscInt npoints, PetscReal xi[])
632: {
633: PetscBool is_simplex, is_tensorcell;
634: PetscInt dim, nfaces, ps, pe, p, d, nbasis, pcnt, e, k, nel;
635: PetscFE fe;
636: PetscQuadrature quadrature;
637: PetscTabulation T;
638: PetscReal *xiq;
639: Vec coorlocal;
640: PetscSection coordSection;
641: PetscScalar *elcoor = NULL;
642: PetscReal *swarm_coor;
643: PetscInt *swarm_cellid;
645: PetscFunctionBegin;
646: PetscCall(DMGetDimension(dmc, &dim));
648: is_simplex = PETSC_FALSE;
649: is_tensorcell = PETSC_FALSE;
650: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
651: PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));
653: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
655: switch (dim) {
656: case 2:
657: if (nfaces == 4) is_tensorcell = PETSC_TRUE;
658: break;
659: case 3:
660: if (nfaces == 6) is_tensorcell = PETSC_TRUE;
661: break;
662: default:
663: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for 2D, 3D");
664: }
666: /* check points provided fail inside the reference cell */
667: if (is_simplex) {
668: for (p = 0; p < npoints; p++) {
669: PetscReal sum;
670: 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");
671: sum = 0.0;
672: for (d = 0; d < dim; d++) sum += xi[dim * p + d];
673: PetscCheck(sum <= 0.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the simplex domain");
674: }
675: } else if (is_tensorcell) {
676: for (p = 0; p < npoints; p++) {
677: 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");
678: }
679: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for d-simplex and d-tensorcell");
681: PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)dm), &quadrature));
682: PetscCall(PetscMalloc1(npoints * dim, &xiq));
683: PetscCall(PetscArraycpy(xiq, xi, npoints * dim));
684: PetscCall(PetscQuadratureSetData(quadrature, dim, 1, npoints, (const PetscReal *)xiq, NULL));
685: PetscCall(private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe));
686: PetscCall(PetscFESetQuadrature(fe, quadrature));
687: PetscCall(PetscFEGetDimension(fe, &nbasis));
688: PetscCall(PetscFEGetCellTabulation(fe, 1, &T));
690: /* for each cell, interpolate coordaintes and insert the interpolated points coordinates into swarm */
691: /* 0->cell, 1->edge, 2->vert */
692: PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
693: nel = pe - ps;
695: PetscCall(DMSwarmSetLocalSizes(dm, npoints * nel, -1));
696: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
697: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
699: PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
700: PetscCall(DMGetCoordinateSection(dmc, &coordSection));
702: pcnt = 0;
703: for (e = 0; e < nel; e++) {
704: PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
706: for (p = 0; p < npoints; p++) {
707: for (d = 0; d < dim; d++) {
708: swarm_coor[dim * pcnt + d] = 0.0;
709: for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][p * nbasis + k] * PetscRealPart(elcoor[dim * k + d]);
710: }
711: swarm_cellid[pcnt] = e;
712: pcnt++;
713: }
714: PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor));
715: }
716: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
717: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
719: PetscCall(PetscQuadratureDestroy(&quadrature));
720: PetscCall(PetscFEDestroy(&fe));
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }