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: /* Create space */
19: PetscSpaceCreate(PetscObjectComm((PetscObject)dm), &P);
20: /* PetscObjectSetOptionsPrefix((PetscObject) P, prefix); */
21: PetscSpacePolynomialSetTensor(P, tensor);
22: /* PetscSpaceSetFromOptions(P); */
23: PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
24: PetscSpaceSetDegree(P, 1, PETSC_DETERMINE);
25: PetscSpaceSetNumComponents(P, Nc);
26: PetscSpaceSetNumVariables(P, dim);
27: PetscSpaceSetUp(P);
28: PetscSpaceGetDegree(P, &order, NULL);
29: PetscSpacePolynomialGetTensor(P, &tensor);
30: /* Create dual space */
31: PetscDualSpaceCreate(PetscObjectComm((PetscObject)dm), &Q);
32: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
33: /* PetscObjectSetOptionsPrefix((PetscObject) Q, prefix); */
34: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K);
35: PetscDualSpaceSetDM(Q, K);
36: DMDestroy(&K);
37: PetscDualSpaceSetNumComponents(Q, Nc);
38: PetscDualSpaceSetOrder(Q, order);
39: PetscDualSpaceLagrangeSetTensor(Q, tensor);
40: /* PetscDualSpaceSetFromOptions(Q); */
41: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
42: PetscDualSpaceSetUp(Q);
43: /* Create element */
44: PetscFECreate(PetscObjectComm((PetscObject)dm), fem);
45: /* PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix); */
46: /* PetscFESetFromOptions(*fem); */
47: PetscFESetType(*fem, PETSCFEBASIC);
48: PetscFESetBasisSpace(*fem, P);
49: PetscFESetDualSpace(*fem, Q);
50: PetscFESetNumComponents(*fem, Nc);
51: PetscFESetUp(*fem);
52: PetscSpaceDestroy(&P);
53: PetscDualSpaceDestroy(&Q);
54: /* Create quadrature (with specified order if given) */
55: qorder = qorder >= 0 ? qorder : order;
56: quadPointsPerEdge = PetscMax(qorder + 1, 1);
57: if (isSimplex) {
58: PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
59: PetscDTStroudConicalQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
60: } else {
61: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
62: PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
63: }
64: PetscFESetQuadrature(*fem, q);
65: PetscFESetFaceQuadrature(*fem, fq);
66: PetscQuadratureDestroy(&q);
67: PetscQuadratureDestroy(&fq);
68: return 0;
69: }
71: PetscErrorCode subdivide_triangle(PetscReal v1[2], PetscReal v2[2], PetscReal v3[2], PetscInt depth, PetscInt max, PetscReal xi[], PetscInt *np)
72: {
73: PetscReal v12[2], v23[2], v31[2];
74: PetscInt i;
76: if (depth == max) {
77: PetscReal cx[2];
79: cx[0] = (v1[0] + v2[0] + v3[0]) / 3.0;
80: cx[1] = (v1[1] + v2[1] + v3[1]) / 3.0;
82: xi[2 * (*np) + 0] = cx[0];
83: xi[2 * (*np) + 1] = cx[1];
84: (*np)++;
85: return 0;
86: }
88: /* calculate midpoints of each side */
89: for (i = 0; i < 2; i++) {
90: v12[i] = (v1[i] + v2[i]) / 2.0;
91: v23[i] = (v2[i] + v3[i]) / 2.0;
92: v31[i] = (v3[i] + v1[i]) / 2.0;
93: }
95: /* recursively subdivide new triangles */
96: subdivide_triangle(v1, v12, v31, depth + 1, max, xi, np);
97: subdivide_triangle(v2, v23, v12, depth + 1, max, xi, np);
98: subdivide_triangle(v3, v31, v23, depth + 1, max, xi, np);
99: subdivide_triangle(v12, v23, v31, depth + 1, max, xi, np);
100: return 0;
101: }
103: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_SubDivide(DM dm, DM dmc, PetscInt nsub)
104: {
105: const PetscInt dim = 2;
106: PetscInt q, npoints_q, e, nel, npe, pcnt, ps, pe, d, k, depth;
107: PetscReal *xi;
108: PetscReal **basis;
109: Vec coorlocal;
110: PetscSection coordSection;
111: PetscScalar *elcoor = NULL;
112: PetscReal *swarm_coor;
113: PetscInt *swarm_cellid;
114: PetscReal v1[2], v2[2], v3[2];
116: npoints_q = 1;
117: for (d = 0; d < nsub; d++) npoints_q *= 4;
118: PetscMalloc1(dim * npoints_q, &xi);
120: v1[0] = 0.0;
121: v1[1] = 0.0;
122: v2[0] = 1.0;
123: v2[1] = 0.0;
124: v3[0] = 0.0;
125: v3[1] = 1.0;
126: depth = 0;
127: pcnt = 0;
128: subdivide_triangle(v1, v2, v3, depth, nsub, xi, &pcnt);
130: npe = 3; /* nodes per element (triangle) */
131: PetscMalloc1(npoints_q, &basis);
132: for (q = 0; q < npoints_q; q++) {
133: PetscMalloc1(npe, &basis[q]);
135: basis[q][0] = 1.0 - xi[dim * q + 0] - xi[dim * q + 1];
136: basis[q][1] = xi[dim * q + 0];
137: basis[q][2] = xi[dim * q + 1];
138: }
140: /* 0->cell, 1->edge, 2->vert */
141: DMPlexGetHeightStratum(dmc, 0, &ps, &pe);
142: nel = pe - ps;
144: DMSwarmSetLocalSizes(dm, npoints_q * nel, -1);
145: DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
146: DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
148: DMGetCoordinatesLocal(dmc, &coorlocal);
149: DMGetCoordinateSection(dmc, &coordSection);
151: pcnt = 0;
152: for (e = 0; e < nel; e++) {
153: DMPlexVecGetClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor);
155: for (q = 0; q < npoints_q; q++) {
156: for (d = 0; d < dim; d++) {
157: swarm_coor[dim * pcnt + d] = 0.0;
158: for (k = 0; k < npe; k++) swarm_coor[dim * pcnt + d] += basis[q][k] * PetscRealPart(elcoor[dim * k + d]);
159: }
160: swarm_cellid[pcnt] = e;
161: pcnt++;
162: }
163: DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor);
164: }
165: DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
166: DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
168: PetscFree(xi);
169: for (q = 0; q < npoints_q; q++) PetscFree(basis[q]);
170: PetscFree(basis);
171: return 0;
172: }
174: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(DM dm, DM dmc, PetscInt nsub)
175: {
176: PetscInt dim, nfaces, nbasis;
177: PetscInt q, npoints_q, e, nel, pcnt, ps, pe, d, k, r;
178: PetscTabulation T;
179: Vec coorlocal;
180: PetscSection coordSection;
181: PetscScalar *elcoor = NULL;
182: PetscReal *swarm_coor;
183: PetscInt *swarm_cellid;
184: const PetscReal *xiq;
185: PetscQuadrature quadrature;
186: PetscFE fe, feRef;
187: PetscBool is_simplex;
189: DMGetDimension(dmc, &dim);
190: is_simplex = PETSC_FALSE;
191: DMPlexGetHeightStratum(dmc, 0, &ps, &pe);
192: DMPlexGetConeSize(dmc, ps, &nfaces);
193: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
195: private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);
197: for (r = 0; r < nsub; r++) {
198: PetscFERefine(fe, &feRef);
199: PetscFECopyQuadrature(feRef, fe);
200: PetscFEDestroy(&feRef);
201: }
203: PetscFEGetQuadrature(fe, &quadrature);
204: PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL);
205: PetscFEGetDimension(fe, &nbasis);
206: PetscFEGetCellTabulation(fe, 1, &T);
208: /* 0->cell, 1->edge, 2->vert */
209: DMPlexGetHeightStratum(dmc, 0, &ps, &pe);
210: nel = pe - ps;
212: DMSwarmSetLocalSizes(dm, npoints_q * nel, -1);
213: DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
214: DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
216: DMGetCoordinatesLocal(dmc, &coorlocal);
217: DMGetCoordinateSection(dmc, &coordSection);
219: pcnt = 0;
220: for (e = 0; e < nel; e++) {
221: DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor);
223: for (q = 0; q < npoints_q; q++) {
224: for (d = 0; d < dim; d++) {
225: swarm_coor[dim * pcnt + d] = 0.0;
226: for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][q * nbasis + k] * PetscRealPart(elcoor[dim * k + d]);
227: }
228: swarm_cellid[pcnt] = e;
229: pcnt++;
230: }
231: DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor);
232: }
233: DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
234: DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
236: PetscFEDestroy(&fe);
237: return 0;
238: }
240: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm, DM dmc, PetscInt npoints)
241: {
242: PetscInt dim;
243: PetscInt ii, jj, q, npoints_q, e, nel, npe, pcnt, ps, pe, d, k, nfaces;
244: PetscReal *xi, ds, ds2;
245: PetscReal **basis;
246: Vec coorlocal;
247: PetscSection coordSection;
248: PetscScalar *elcoor = NULL;
249: PetscReal *swarm_coor;
250: PetscInt *swarm_cellid;
251: PetscBool is_simplex;
253: DMGetDimension(dmc, &dim);
255: is_simplex = PETSC_FALSE;
256: DMPlexGetHeightStratum(dmc, 0, &ps, &pe);
257: DMPlexGetConeSize(dmc, ps, &nfaces);
258: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
261: PetscMalloc1(dim * npoints * npoints, &xi);
262: pcnt = 0;
263: ds = 1.0 / ((PetscReal)(npoints - 1));
264: ds2 = 1.0 / ((PetscReal)(npoints));
265: for (jj = 0; jj < npoints; jj++) {
266: for (ii = 0; ii < npoints - jj; ii++) {
267: xi[dim * pcnt + 0] = ii * ds;
268: xi[dim * pcnt + 1] = jj * ds;
270: xi[dim * pcnt + 0] *= (1.0 - 1.2 * ds2);
271: xi[dim * pcnt + 1] *= (1.0 - 1.2 * ds2);
273: xi[dim * pcnt + 0] += 0.35 * ds2;
274: xi[dim * pcnt + 1] += 0.35 * ds2;
275: pcnt++;
276: }
277: }
278: npoints_q = pcnt;
280: npe = 3; /* nodes per element (triangle) */
281: PetscMalloc1(npoints_q, &basis);
282: for (q = 0; q < npoints_q; q++) {
283: PetscMalloc1(npe, &basis[q]);
285: basis[q][0] = 1.0 - xi[dim * q + 0] - xi[dim * q + 1];
286: basis[q][1] = xi[dim * q + 0];
287: basis[q][2] = xi[dim * q + 1];
288: }
290: /* 0->cell, 1->edge, 2->vert */
291: DMPlexGetHeightStratum(dmc, 0, &ps, &pe);
292: nel = pe - ps;
294: DMSwarmSetLocalSizes(dm, npoints_q * nel, -1);
295: DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
296: DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
298: DMGetCoordinatesLocal(dmc, &coorlocal);
299: DMGetCoordinateSection(dmc, &coordSection);
301: pcnt = 0;
302: for (e = 0; e < nel; e++) {
303: DMPlexVecGetClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor);
305: for (q = 0; q < npoints_q; q++) {
306: for (d = 0; d < dim; d++) {
307: swarm_coor[dim * pcnt + d] = 0.0;
308: for (k = 0; k < npe; k++) swarm_coor[dim * pcnt + d] += basis[q][k] * PetscRealPart(elcoor[dim * k + d]);
309: }
310: swarm_cellid[pcnt] = e;
311: pcnt++;
312: }
313: DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor);
314: }
315: DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
316: DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
318: PetscFree(xi);
319: for (q = 0; q < npoints_q; q++) PetscFree(basis[q]);
320: PetscFree(basis);
321: return 0;
322: }
324: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param)
325: {
326: PetscInt dim;
328: DMGetDimension(celldm, &dim);
329: switch (layout) {
330: case DMSWARMPIC_LAYOUT_REGULAR:
332: private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm, celldm, layout_param);
333: break;
334: case DMSWARMPIC_LAYOUT_GAUSS: {
335: PetscInt npoints, npoints1, ps, pe, nfaces;
336: const PetscReal *xi;
337: PetscBool is_simplex;
338: PetscQuadrature quadrature;
340: is_simplex = PETSC_FALSE;
341: DMPlexGetHeightStratum(celldm, 0, &ps, &pe);
342: DMPlexGetConeSize(celldm, ps, &nfaces);
343: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
345: npoints1 = layout_param;
346: if (is_simplex) {
347: PetscDTStroudConicalQuadrature(dim, 1, npoints1, -1.0, 1.0, &quadrature);
348: } else {
349: PetscDTGaussTensorQuadrature(dim, 1, npoints1, -1.0, 1.0, &quadrature);
350: }
351: PetscQuadratureGetData(quadrature, NULL, NULL, &npoints, &xi, NULL);
352: private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm, celldm, npoints, (PetscReal *)xi);
353: PetscQuadratureDestroy(&quadrature);
354: } break;
355: case DMSWARMPIC_LAYOUT_SUBDIVISION:
356: private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm, celldm, layout_param);
357: break;
358: }
359: return 0;
360: }
362: /*
363: typedef struct {
364: PetscReal x,y;
365: } Point2d;
367: static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s)
368: {
369: *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y);
370: return 0;
371: }
372: */
373: /*
374: static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v)
375: {
376: PetscReal s1,s2,s3;
377: PetscBool b1, b2, b3;
379: signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f;
380: signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f;
381: signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f;
383: *v = ((b1 == b2) && (b2 == b3));
384: return 0;
385: }
386: */
387: /*
388: static PetscErrorCode _ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ)
389: {
390: PetscReal x1,y1,x2,y2,x3,y3;
391: PetscReal c,b[2],A[2][2],inv[2][2],detJ,od;
393: x1 = coords[2*0+0];
394: x2 = coords[2*1+0];
395: x3 = coords[2*2+0];
397: y1 = coords[2*0+1];
398: y2 = coords[2*1+1];
399: y3 = coords[2*2+1];
401: c = x1 - 0.5*x1 - 0.5*x1 + 0.5*x2 + 0.5*x3;
402: b[0] = xp[0] - c;
403: c = y1 - 0.5*y1 - 0.5*y1 + 0.5*y2 + 0.5*y3;
404: b[1] = xp[1] - c;
406: A[0][0] = -0.5*x1 + 0.5*x2; A[0][1] = -0.5*x1 + 0.5*x3;
407: A[1][0] = -0.5*y1 + 0.5*y2; A[1][1] = -0.5*y1 + 0.5*y3;
409: detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0];
410: *dJ = PetscAbsReal(detJ);
411: od = 1.0/detJ;
413: inv[0][0] = A[1][1] * od;
414: inv[0][1] = -A[0][1] * od;
415: inv[1][0] = -A[1][0] * od;
416: inv[1][1] = A[0][0] * od;
418: xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1];
419: xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1];
420: return 0;
421: }
422: */
424: static PetscErrorCode ComputeLocalCoordinateAffine2d(PetscReal xp[], PetscScalar coords[], PetscReal xip[], PetscReal *dJ)
425: {
426: PetscReal x1, y1, x2, y2, x3, y3;
427: PetscReal b[2], A[2][2], inv[2][2], detJ, od;
429: x1 = PetscRealPart(coords[2 * 0 + 0]);
430: x2 = PetscRealPart(coords[2 * 1 + 0]);
431: x3 = PetscRealPart(coords[2 * 2 + 0]);
433: y1 = PetscRealPart(coords[2 * 0 + 1]);
434: y2 = PetscRealPart(coords[2 * 1 + 1]);
435: y3 = PetscRealPart(coords[2 * 2 + 1]);
437: b[0] = xp[0] - x1;
438: b[1] = xp[1] - y1;
440: A[0][0] = x2 - x1;
441: A[0][1] = x3 - x1;
442: A[1][0] = y2 - y1;
443: A[1][1] = y3 - y1;
445: detJ = A[0][0] * A[1][1] - A[0][1] * A[1][0];
446: *dJ = PetscAbsReal(detJ);
447: od = 1.0 / detJ;
449: inv[0][0] = A[1][1] * od;
450: inv[0][1] = -A[0][1] * od;
451: inv[1][0] = -A[1][0] * od;
452: inv[1][1] = A[0][0] * od;
454: xip[0] = inv[0][0] * b[0] + inv[0][1] * b[1];
455: xip[1] = inv[1][0] * b[0] + inv[1][1] * b[1];
456: return 0;
457: }
459: PetscErrorCode DMSwarmProjectField_ApproxP1_PLEX_2D(DM swarm, PetscReal *swarm_field, DM dm, Vec v_field)
460: {
461: const PetscReal PLEX_C_EPS = 1.0e-8;
462: Vec v_field_l, denom_l, coor_l, denom;
463: PetscInt k, p, e, npoints;
464: PetscInt *mpfield_cell;
465: PetscReal *mpfield_coor;
466: PetscReal xi_p[2];
467: PetscScalar Ni[3];
468: PetscSection coordSection;
469: PetscScalar *elcoor = NULL;
471: VecZeroEntries(v_field);
473: DMGetLocalVector(dm, &v_field_l);
474: DMGetGlobalVector(dm, &denom);
475: DMGetLocalVector(dm, &denom_l);
476: VecZeroEntries(v_field_l);
477: VecZeroEntries(denom);
478: VecZeroEntries(denom_l);
480: DMGetCoordinatesLocal(dm, &coor_l);
481: DMGetCoordinateSection(dm, &coordSection);
483: DMSwarmGetLocalSize(swarm, &npoints);
484: DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor);
485: DMSwarmGetField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell);
487: for (p = 0; p < npoints; p++) {
488: PetscReal *coor_p, dJ;
489: PetscScalar elfield[3];
490: PetscBool point_located;
492: e = mpfield_cell[p];
493: coor_p = &mpfield_coor[2 * p];
495: DMPlexVecGetClosure(dm, coordSection, coor_l, e, NULL, &elcoor);
497: /*
498: while (!point_located && (failed_counter < 25)) {
499: PointInTriangle(point, coords[0], coords[1], coords[2], &point_located);
500: point.x = coor_p[0];
501: point.y = coor_p[1];
502: point.x += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
503: point.y += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
504: failed_counter++;
505: }
507: if (!point_located) {
508: 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);
509: }
512: else {
513: _ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);
514: xi_p[0] = 0.5*(xi_p[0] + 1.0);
515: xi_p[1] = 0.5*(xi_p[1] + 1.0);
517: 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]);
519: }
520: */
522: ComputeLocalCoordinateAffine2d(coor_p, elcoor, xi_p, &dJ);
523: /*
524: 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]);
525: */
526: /*
527: point_located = PETSC_TRUE;
528: if (xi_p[0] < 0.0) {
529: if (xi_p[0] > -PLEX_C_EPS) {
530: xi_p[0] = 0.0;
531: } else {
532: point_located = PETSC_FALSE;
533: }
534: }
535: if (xi_p[1] < 0.0) {
536: if (xi_p[1] > -PLEX_C_EPS) {
537: xi_p[1] = 0.0;
538: } else {
539: point_located = PETSC_FALSE;
540: }
541: }
542: if (xi_p[1] > (1.0-xi_p[0])) {
543: if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) {
544: xi_p[1] = 1.0 - xi_p[0];
545: } else {
546: point_located = PETSC_FALSE;
547: }
548: }
549: if (!point_located) {
550: PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]);
551: 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]);
552: }
554: */
556: Ni[0] = 1.0 - xi_p[0] - xi_p[1];
557: Ni[1] = xi_p[0];
558: Ni[2] = xi_p[1];
560: point_located = PETSC_TRUE;
561: for (k = 0; k < 3; k++) {
562: if (PetscRealPart(Ni[k]) < -PLEX_C_EPS) point_located = PETSC_FALSE;
563: if (PetscRealPart(Ni[k]) > (1.0 + PLEX_C_EPS)) point_located = PETSC_FALSE;
564: }
565: if (!point_located) {
566: PetscPrintf(PETSC_COMM_SELF, "[Error] xi,eta = %+1.8e, %+1.8e\n", (double)xi_p[0], (double)xi_p[1]);
567: 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]));
568: }
571: for (k = 0; k < 3; k++) {
572: Ni[k] = Ni[k] * dJ;
573: elfield[k] = Ni[k] * swarm_field[p];
574: }
575: DMPlexVecRestoreClosure(dm, coordSection, coor_l, e, NULL, &elcoor);
577: DMPlexVecSetClosure(dm, NULL, v_field_l, e, elfield, ADD_VALUES);
578: DMPlexVecSetClosure(dm, NULL, denom_l, e, Ni, ADD_VALUES);
579: }
581: DMSwarmRestoreField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell);
582: DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor);
584: DMLocalToGlobalBegin(dm, v_field_l, ADD_VALUES, v_field);
585: DMLocalToGlobalEnd(dm, v_field_l, ADD_VALUES, v_field);
586: DMLocalToGlobalBegin(dm, denom_l, ADD_VALUES, denom);
587: DMLocalToGlobalEnd(dm, denom_l, ADD_VALUES, denom);
589: VecPointwiseDivide(v_field, v_field, denom);
591: DMRestoreLocalVector(dm, &v_field_l);
592: DMRestoreLocalVector(dm, &denom_l);
593: DMRestoreGlobalVector(dm, &denom);
594: return 0;
595: }
597: PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm, DM celldm, PetscInt project_type, PetscInt nfields, DMSwarmDataField dfield[], Vec vecs[])
598: {
599: PetscInt f, dim;
601: DMGetDimension(swarm, &dim);
602: switch (dim) {
603: case 2:
604: for (f = 0; f < nfields; f++) {
605: PetscReal *swarm_field;
607: DMSwarmDataFieldGetEntries(dfield[f], (void **)&swarm_field);
608: DMSwarmProjectField_ApproxP1_PLEX_2D(swarm, swarm_field, celldm, vecs[f]);
609: }
610: break;
611: case 3:
612: SETERRQ(PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "No support for 3D");
613: default:
614: break;
615: }
616: return 0;
617: }
619: PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm, DM dmc, PetscInt npoints, PetscReal xi[])
620: {
621: PetscBool is_simplex, is_tensorcell;
622: PetscInt dim, nfaces, ps, pe, p, d, nbasis, pcnt, e, k, nel;
623: PetscFE fe;
624: PetscQuadrature quadrature;
625: PetscTabulation T;
626: PetscReal *xiq;
627: Vec coorlocal;
628: PetscSection coordSection;
629: PetscScalar *elcoor = NULL;
630: PetscReal *swarm_coor;
631: PetscInt *swarm_cellid;
633: DMGetDimension(dmc, &dim);
635: is_simplex = PETSC_FALSE;
636: is_tensorcell = PETSC_FALSE;
637: DMPlexGetHeightStratum(dmc, 0, &ps, &pe);
638: DMPlexGetConeSize(dmc, ps, &nfaces);
640: if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
642: switch (dim) {
643: case 2:
644: if (nfaces == 4) is_tensorcell = PETSC_TRUE;
645: break;
646: case 3:
647: if (nfaces == 6) is_tensorcell = PETSC_TRUE;
648: break;
649: default:
650: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for 2D, 3D");
651: }
653: /* check points provided fail inside the reference cell */
654: if (is_simplex) {
655: for (p = 0; p < npoints; p++) {
656: PetscReal sum;
658: sum = 0.0;
659: for (d = 0; d < dim; d++) sum += xi[dim * p + d];
661: }
662: } else if (is_tensorcell) {
663: for (p = 0; p < npoints; p++) {
665: }
666: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for d-simplex and d-tensorcell");
668: PetscQuadratureCreate(PetscObjectComm((PetscObject)dm), &quadrature);
669: PetscMalloc1(npoints * dim, &xiq);
670: PetscArraycpy(xiq, xi, npoints * dim);
671: PetscQuadratureSetData(quadrature, dim, 1, npoints, (const PetscReal *)xiq, NULL);
672: private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);
673: PetscFESetQuadrature(fe, quadrature);
674: PetscFEGetDimension(fe, &nbasis);
675: PetscFEGetCellTabulation(fe, 1, &T);
677: /* for each cell, interpolate coordaintes and insert the interpolated points coordinates into swarm */
678: /* 0->cell, 1->edge, 2->vert */
679: DMPlexGetHeightStratum(dmc, 0, &ps, &pe);
680: nel = pe - ps;
682: DMSwarmSetLocalSizes(dm, npoints * nel, -1);
683: DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
684: DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
686: DMGetCoordinatesLocal(dmc, &coorlocal);
687: DMGetCoordinateSection(dmc, &coordSection);
689: pcnt = 0;
690: for (e = 0; e < nel; e++) {
691: DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor);
693: for (p = 0; p < npoints; p++) {
694: for (d = 0; d < dim; d++) {
695: swarm_coor[dim * pcnt + d] = 0.0;
696: for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][p * nbasis + k] * PetscRealPart(elcoor[dim * k + d]);
697: }
698: swarm_cellid[pcnt] = e;
699: pcnt++;
700: }
701: DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor);
702: }
703: DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
704: DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
706: PetscQuadratureDestroy(&quadrature);
707: PetscFEDestroy(&fe);
708: return 0;
709: }