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;
17: PetscErrorCode ierr;
20: /* Create space */
21: PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);
22: /*PetscObjectSetOptionsPrefix((PetscObject) P, prefix);*/
23: PetscSpacePolynomialSetTensor(P, tensor);
24: /*PetscSpaceSetFromOptions(P);*/
25: PetscSpaceSetType(P,PETSCSPACEPOLYNOMIAL);
26: PetscSpaceSetDegree(P,1,PETSC_DETERMINE);
27: PetscSpaceSetNumComponents(P, Nc);
28: PetscSpaceSetNumVariables(P, dim);
29: PetscSpaceSetUp(P);
30: PetscSpaceGetDegree(P, &order, NULL);
31: PetscSpacePolynomialGetTensor(P, &tensor);
32: /* Create dual space */
33: PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
34: PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
35: /*PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);*/
36: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
37: PetscDualSpaceSetDM(Q, K);
38: DMDestroy(&K);
39: PetscDualSpaceSetNumComponents(Q, Nc);
40: PetscDualSpaceSetOrder(Q, order);
41: PetscDualSpaceLagrangeSetTensor(Q, tensor);
42: /*PetscDualSpaceSetFromOptions(Q);*/
43: PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
44: PetscDualSpaceSetUp(Q);
45: /* Create element */
46: PetscFECreate(PetscObjectComm((PetscObject) dm), fem);
47: /*PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);*/
48: /*PetscFESetFromOptions(*fem);*/
49: PetscFESetType(*fem,PETSCFEBASIC);
50: PetscFESetBasisSpace(*fem, P);
51: PetscFESetDualSpace(*fem, Q);
52: PetscFESetNumComponents(*fem, Nc);
53: PetscFESetUp(*fem);
54: PetscSpaceDestroy(&P);
55: PetscDualSpaceDestroy(&Q);
56: /* Create quadrature (with specified order if given) */
57: qorder = qorder >= 0 ? qorder : order;
58: quadPointsPerEdge = PetscMax(qorder + 1,1);
59: if (isSimplex) {
60: PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
61: PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
62: }
63: else {
64: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
65: PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
66: }
67: PetscFESetQuadrature(*fem, q);
68: PetscFESetFaceQuadrature(*fem, fq);
69: PetscQuadratureDestroy(&q);
70: PetscQuadratureDestroy(&fq);
71: return(0);
72: }
74: PetscErrorCode subdivide_triangle(PetscReal v1[2],PetscReal v2[2],PetscReal v3[2],PetscInt depth,PetscInt max,PetscReal xi[],PetscInt *np)
75: {
76: PetscReal v12[2],v23[2],v31[2];
77: PetscInt i;
81: if (depth == max) {
82: PetscReal cx[2];
84: cx[0] = (v1[0] + v2[0] + v3[0])/3.0;
85: cx[1] = (v1[1] + v2[1] + v3[1])/3.0;
87: xi[2*(*np)+0] = cx[0];
88: xi[2*(*np)+1] = cx[1];
89: (*np)++;
90: return(0);
91: }
93: /* calculate midpoints of each side */
94: for (i = 0; i < 2; i++) {
95: v12[i] = (v1[i]+v2[i])/2.0;
96: v23[i] = (v2[i]+v3[i])/2.0;
97: v31[i] = (v3[i]+v1[i])/2.0;
98: }
100: /* recursively subdivide new triangles */
101: subdivide_triangle(v1,v12,v31,depth+1,max,xi,np);
102: subdivide_triangle(v2,v23,v12,depth+1,max,xi,np);
103: subdivide_triangle(v3,v31,v23,depth+1,max,xi,np);
104: subdivide_triangle(v12,v23,v31,depth+1,max,xi,np);
105: return(0);
106: }
108: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_SubDivide(DM dm,DM dmc,PetscInt nsub)
109: {
111: const PetscInt dim = 2;
112: PetscInt q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,depth;
113: PetscReal *xi;
114: PetscReal **basis;
115: Vec coorlocal;
116: PetscSection coordSection;
117: PetscScalar *elcoor = NULL;
118: PetscReal *swarm_coor;
119: PetscInt *swarm_cellid;
120: PetscReal v1[2],v2[2],v3[2];
124: npoints_q = 1;
125: for (d=0; d<nsub; d++) { npoints_q *= 4; }
126: PetscMalloc1(dim*npoints_q,&xi);
128: v1[0] = 0.0; v1[1] = 0.0;
129: v2[0] = 1.0; v2[1] = 0.0;
130: v3[0] = 0.0; v3[1] = 1.0;
131: depth = 0;
132: pcnt = 0;
133: subdivide_triangle(v1,v2,v3,depth,nsub,xi,&pcnt);
135: npe = 3; /* nodes per element (triangle) */
136: PetscMalloc1(npoints_q,&basis);
137: for (q=0; q<npoints_q; q++) {
138: PetscMalloc1(npe,&basis[q]);
140: basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1];
141: basis[q][1] = xi[dim*q+0];
142: basis[q][2] = xi[dim*q+1];
143: }
145: /* 0->cell, 1->edge, 2->vert */
146: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
147: nel = pe - ps;
149: DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
150: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
151: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
153: DMGetCoordinatesLocal(dmc,&coorlocal);
154: DMGetCoordinateSection(dmc,&coordSection);
156: pcnt = 0;
157: for (e=0; e<nel; e++) {
158: DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);
160: for (q=0; q<npoints_q; q++) {
161: for (d=0; d<dim; d++) {
162: swarm_coor[dim*pcnt+d] = 0.0;
163: for (k=0; k<npe; k++) {
164: swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]);
165: }
166: }
167: swarm_cellid[pcnt] = e;
168: pcnt++;
169: }
170: DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);
171: }
172: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
173: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
175: PetscFree(xi);
176: for (q=0; q<npoints_q; q++) {
177: PetscFree(basis[q]);
178: }
179: PetscFree(basis);
181: return(0);
182: }
184: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(DM dm,DM dmc,PetscInt nsub)
185: {
187: PetscInt dim,nfaces,nbasis;
188: PetscInt q,npoints_q,e,nel,pcnt,ps,pe,d,k,r;
189: PetscTabulation T;
190: Vec coorlocal;
191: PetscSection coordSection;
192: PetscScalar *elcoor = NULL;
193: PetscReal *swarm_coor;
194: PetscInt *swarm_cellid;
195: const PetscReal *xiq;
196: PetscQuadrature quadrature;
197: PetscFE fe,feRef;
198: PetscBool is_simplex;
202: DMGetDimension(dmc,&dim);
204: is_simplex = PETSC_FALSE;
205: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
206: DMPlexGetConeSize(dmc, ps, &nfaces);
207: if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }
209: private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);
211: for (r=0; r<nsub; r++) {
212: PetscFERefine(fe,&feRef);
213: PetscFECopyQuadrature(feRef,fe);
214: PetscFEDestroy(&feRef);
215: }
217: PetscFEGetQuadrature(fe,&quadrature);
218: PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL);
219: PetscFEGetDimension(fe,&nbasis);
220: PetscFEGetCellTabulation(fe, 1, &T);
222: /* 0->cell, 1->edge, 2->vert */
223: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
224: nel = pe - ps;
226: DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
227: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
228: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
230: DMGetCoordinatesLocal(dmc,&coorlocal);
231: DMGetCoordinateSection(dmc,&coordSection);
233: pcnt = 0;
234: for (e=0; e<nel; e++) {
235: DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);
237: for (q=0; q<npoints_q; q++) {
238: for (d=0; d<dim; d++) {
239: swarm_coor[dim*pcnt+d] = 0.0;
240: for (k=0; k<nbasis; k++) {
241: swarm_coor[dim*pcnt+d] += T->T[0][q*nbasis + k] * PetscRealPart(elcoor[dim*k+d]);
242: }
243: }
244: swarm_cellid[pcnt] = e;
245: pcnt++;
246: }
247: DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);
248: }
249: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
250: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
252: PetscFEDestroy(&fe);
253: return(0);
254: }
256: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm,DM dmc,PetscInt npoints)
257: {
259: PetscInt dim;
260: PetscInt ii,jj,q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,nfaces;
261: PetscReal *xi,ds,ds2;
262: PetscReal **basis;
263: Vec coorlocal;
264: PetscSection coordSection;
265: PetscScalar *elcoor = NULL;
266: PetscReal *swarm_coor;
267: PetscInt *swarm_cellid;
268: PetscBool is_simplex;
271: DMGetDimension(dmc,&dim);
272: if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only 2D is supported");
273: is_simplex = PETSC_FALSE;
274: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
275: DMPlexGetConeSize(dmc, ps, &nfaces);
276: if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }
277: if (!is_simplex) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only the simplex is supported");
279: PetscMalloc1(dim*npoints*npoints,&xi);
280: pcnt = 0;
281: ds = 1.0/((PetscReal)(npoints-1));
282: ds2 = 1.0/((PetscReal)(npoints));
283: for (jj = 0; jj<npoints; jj++) {
284: for (ii=0; ii<npoints-jj; ii++) {
285: xi[dim*pcnt+0] = ii * ds;
286: xi[dim*pcnt+1] = jj * ds;
288: xi[dim*pcnt+0] *= (1.0 - 1.2*ds2);
289: xi[dim*pcnt+1] *= (1.0 - 1.2*ds2);
291: xi[dim*pcnt+0] += 0.35*ds2;
292: xi[dim*pcnt+1] += 0.35*ds2;
294: pcnt++;
295: }
296: }
297: npoints_q = pcnt;
299: npe = 3; /* nodes per element (triangle) */
300: PetscMalloc1(npoints_q,&basis);
301: for (q=0; q<npoints_q; q++) {
302: PetscMalloc1(npe,&basis[q]);
304: basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1];
305: basis[q][1] = xi[dim*q+0];
306: basis[q][2] = xi[dim*q+1];
307: }
309: /* 0->cell, 1->edge, 2->vert */
310: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
311: nel = pe - ps;
313: DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
314: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
315: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
317: DMGetCoordinatesLocal(dmc,&coorlocal);
318: DMGetCoordinateSection(dmc,&coordSection);
320: pcnt = 0;
321: for (e=0; e<nel; e++) {
322: DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);
324: for (q=0; q<npoints_q; q++) {
325: for (d=0; d<dim; d++) {
326: swarm_coor[dim*pcnt+d] = 0.0;
327: for (k=0; k<npe; k++) {
328: swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]);
329: }
330: }
331: swarm_cellid[pcnt] = e;
332: pcnt++;
333: }
334: DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);
335: }
336: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
337: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
339: PetscFree(xi);
340: for (q=0; q<npoints_q; q++) {
341: PetscFree(basis[q]);
342: }
343: PetscFree(basis);
345: return(0);
346: }
348: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)
349: {
351: PetscInt dim;
354: DMGetDimension(celldm,&dim);
355: switch (layout) {
356: case DMSWARMPIC_LAYOUT_REGULAR:
357: if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No 3D support for REGULAR+PLEX");
358: private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm,celldm,layout_param);
359: break;
360: case DMSWARMPIC_LAYOUT_GAUSS:
361: {
362: PetscInt npoints,npoints1,ps,pe,nfaces;
363: const PetscReal *xi;
364: PetscBool is_simplex;
365: PetscQuadrature quadrature;
367: is_simplex = PETSC_FALSE;
368: DMPlexGetHeightStratum(celldm,0,&ps,&pe);
369: DMPlexGetConeSize(celldm, ps, &nfaces);
370: if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }
372: npoints1 = layout_param;
373: if (is_simplex) {
374: PetscDTStroudConicalQuadrature(dim,1,npoints1,-1.0,1.0,&quadrature);
375: } else {
376: PetscDTGaussTensorQuadrature(dim,1,npoints1,-1.0,1.0,&quadrature);
377: }
378: PetscQuadratureGetData(quadrature,NULL,NULL,&npoints,&xi,NULL);
379: private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm,celldm,npoints,(PetscReal*)xi);
380: PetscQuadratureDestroy(&quadrature);
381: }
382: break;
383: case DMSWARMPIC_LAYOUT_SUBDIVISION:
384: private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm,celldm,layout_param);
385: break;
386: }
387: return(0);
388: }
390: /*
391: typedef struct {
392: PetscReal x,y;
393: } Point2d;
395: static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s)
396: {
398: *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y);
399: return(0);
400: }
401: */
402: /*
403: static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v)
404: {
405: PetscReal s1,s2,s3;
406: PetscBool b1, b2, b3;
409: signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f;
410: signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f;
411: signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f;
413: *v = ((b1 == b2) && (b2 == b3));
414: return(0);
415: }
416: */
417: /*
418: static PetscErrorCode _ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ)
419: {
420: PetscReal x1,y1,x2,y2,x3,y3;
421: PetscReal c,b[2],A[2][2],inv[2][2],detJ,od;
424: x1 = coords[2*0+0];
425: x2 = coords[2*1+0];
426: x3 = coords[2*2+0];
428: y1 = coords[2*0+1];
429: y2 = coords[2*1+1];
430: y3 = coords[2*2+1];
432: c = x1 - 0.5*x1 - 0.5*x1 + 0.5*x2 + 0.5*x3;
433: b[0] = xp[0] - c;
434: c = y1 - 0.5*y1 - 0.5*y1 + 0.5*y2 + 0.5*y3;
435: b[1] = xp[1] - c;
437: A[0][0] = -0.5*x1 + 0.5*x2; A[0][1] = -0.5*x1 + 0.5*x3;
438: A[1][0] = -0.5*y1 + 0.5*y2; A[1][1] = -0.5*y1 + 0.5*y3;
440: detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0];
441: *dJ = PetscAbsReal(detJ);
442: od = 1.0/detJ;
444: inv[0][0] = A[1][1] * od;
445: inv[0][1] = -A[0][1] * od;
446: inv[1][0] = -A[1][0] * od;
447: inv[1][1] = A[0][0] * od;
449: xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1];
450: xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1];
451: return(0);
452: }
453: */
455: static PetscErrorCode ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscScalar coords[],PetscReal xip[],PetscReal *dJ)
456: {
457: PetscReal x1,y1,x2,y2,x3,y3;
458: PetscReal b[2],A[2][2],inv[2][2],detJ,od;
461: x1 = PetscRealPart(coords[2*0+0]);
462: x2 = PetscRealPart(coords[2*1+0]);
463: x3 = PetscRealPart(coords[2*2+0]);
465: y1 = PetscRealPart(coords[2*0+1]);
466: y2 = PetscRealPart(coords[2*1+1]);
467: y3 = PetscRealPart(coords[2*2+1]);
469: b[0] = xp[0] - x1;
470: b[1] = xp[1] - y1;
472: A[0][0] = x2-x1; A[0][1] = x3-x1;
473: A[1][0] = y2-y1; A[1][1] = y3-y1;
475: detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0];
476: *dJ = PetscAbsReal(detJ);
477: od = 1.0/detJ;
479: inv[0][0] = A[1][1] * od;
480: inv[0][1] = -A[0][1] * od;
481: inv[1][0] = -A[1][0] * od;
482: inv[1][1] = A[0][0] * od;
484: xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1];
485: xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1];
486: return(0);
487: }
489: PetscErrorCode DMSwarmProjectField_ApproxP1_PLEX_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field)
490: {
492: const PetscReal PLEX_C_EPS = 1.0e-8;
493: Vec v_field_l,denom_l,coor_l,denom;
494: PetscInt k,p,e,npoints;
495: PetscInt *mpfield_cell;
496: PetscReal *mpfield_coor;
497: PetscReal xi_p[2];
498: PetscScalar Ni[3];
499: PetscSection coordSection;
500: PetscScalar *elcoor = NULL;
503: VecZeroEntries(v_field);
505: DMGetLocalVector(dm,&v_field_l);
506: DMGetGlobalVector(dm,&denom);
507: DMGetLocalVector(dm,&denom_l);
508: VecZeroEntries(v_field_l);
509: VecZeroEntries(denom);
510: VecZeroEntries(denom_l);
512: DMGetCoordinatesLocal(dm,&coor_l);
513: DMGetCoordinateSection(dm,&coordSection);
515: DMSwarmGetLocalSize(swarm,&npoints);
516: DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
517: DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
519: for (p=0; p<npoints; p++) {
520: PetscReal *coor_p,dJ;
521: PetscScalar elfield[3];
522: PetscBool point_located;
524: e = mpfield_cell[p];
525: coor_p = &mpfield_coor[2*p];
527: DMPlexVecGetClosure(dm,coordSection,coor_l,e,NULL,&elcoor);
529: /*
530: while (!point_located && (failed_counter < 25)) {
531: PointInTriangle(point, coords[0], coords[1], coords[2], &point_located);
532: point.x = coor_p[0];
533: point.y = coor_p[1];
534: point.x += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
535: point.y += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
536: failed_counter++;
537: }
539: if (!point_located) {
540: PetscPrintf(PETSC_COMM_SELF,"Failed to locate point (%1.8e,%1.8e) in local mesh (cell %D) with triangle coords (%1.8e,%1.8e) : (%1.8e,%1.8e) : (%1.8e,%1.8e) in %D 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);
541: }
543: if (!point_located) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Failed to locate point (%1.8e,%1.8e) in local mesh (cell %D)\n",point.x,point.y,e);
544: else {
545: _ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);
546: xi_p[0] = 0.5*(xi_p[0] + 1.0);
547: xi_p[1] = 0.5*(xi_p[1] + 1.0);
549: PetscPrintf(PETSC_COMM_SELF,"[p=%D] x(%+1.4e,%+1.4e) -> mapped to element %D xi(%+1.4e,%+1.4e)\n",p,point.x,point.y,e,xi_p[0],xi_p[1]);
551: }
552: */
554: ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);
555: /*
556: PetscPrintf(PETSC_COMM_SELF,"[p=%D] x(%+1.4e,%+1.4e) -> mapped to element %D xi(%+1.4e,%+1.4e)\n",p,point.x,point.y,e,xi_p[0],xi_p[1]);
557: */
558: /*
559: point_located = PETSC_TRUE;
560: if (xi_p[0] < 0.0) {
561: if (xi_p[0] > -PLEX_C_EPS) {
562: xi_p[0] = 0.0;
563: } else {
564: point_located = PETSC_FALSE;
565: }
566: }
567: if (xi_p[1] < 0.0) {
568: if (xi_p[1] > -PLEX_C_EPS) {
569: xi_p[1] = 0.0;
570: } else {
571: point_located = PETSC_FALSE;
572: }
573: }
574: if (xi_p[1] > (1.0-xi_p[0])) {
575: if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) {
576: xi_p[1] = 1.0 - xi_p[0];
577: } else {
578: point_located = PETSC_FALSE;
579: }
580: }
581: if (!point_located) {
582: PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]);
583: PetscPrintf(PETSC_COMM_SELF,"[Error] Failed to locate point (%1.8e,%1.8e) in local mesh (cell %D) 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]);
584: }
585: if (!point_located) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Failed to locate point (%1.8e,%1.8e) in local mesh (cell %D)\n",coor_p[0],coor_p[1],e);
586: */
588: Ni[0] = 1.0 - xi_p[0] - xi_p[1];
589: Ni[1] = xi_p[0];
590: Ni[2] = xi_p[1];
592: point_located = PETSC_TRUE;
593: for (k=0; k<3; k++) {
594: if (PetscRealPart(Ni[k]) < -PLEX_C_EPS) point_located = PETSC_FALSE;
595: if (PetscRealPart(Ni[k]) > (1.0+PLEX_C_EPS)) point_located = PETSC_FALSE;
596: }
597: if (!point_located) {
598: PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",(double)xi_p[0],(double)xi_p[1]);
599: PetscPrintf(PETSC_COMM_SELF,"[Error] Failed to locate point (%1.8e,%1.8e) in local mesh (cell %D) 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]));
600: }
601: if (!point_located) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Failed to locate point (%1.8e,%1.8e) in local mesh (cell %D)\n",(double)coor_p[0],(double)coor_p[1],e);
603: for (k=0; k<3; k++) {
604: Ni[k] = Ni[k] * dJ;
605: elfield[k] = Ni[k] * swarm_field[p];
606: }
608: DMPlexVecRestoreClosure(dm,coordSection,coor_l,e,NULL,&elcoor);
610: DMPlexVecSetClosure(dm, NULL,v_field_l, e, elfield, ADD_VALUES);
611: DMPlexVecSetClosure(dm, NULL,denom_l, e, Ni, ADD_VALUES);
612: }
614: DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
615: DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
617: DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);
618: DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);
619: DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);
620: DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);
622: VecPointwiseDivide(v_field,v_field,denom);
624: DMRestoreLocalVector(dm,&v_field_l);
625: DMRestoreLocalVector(dm,&denom_l);
626: DMRestoreGlobalVector(dm,&denom);
628: return(0);
629: }
631: PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[])
632: {
634: PetscInt f,dim;
637: DMGetDimension(swarm,&dim);
638: switch (dim) {
639: case 2:
640: for (f=0; f<nfields; f++) {
641: PetscReal *swarm_field;
643: DMSwarmDataFieldGetEntries(dfield[f],(void**)&swarm_field);
644: DMSwarmProjectField_ApproxP1_PLEX_2D(swarm,swarm_field,celldm,vecs[f]);
645: }
646: break;
647: case 3:
648: SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D");
649: default:
650: break;
651: }
653: return(0);
654: }
656: PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm,DM dmc,PetscInt npoints,PetscReal xi[])
657: {
658: PetscBool is_simplex,is_tensorcell;
660: PetscInt dim,nfaces,ps,pe,p,d,nbasis,pcnt,e,k,nel;
661: PetscFE fe;
662: PetscQuadrature quadrature;
663: PetscTabulation T;
664: PetscReal *xiq;
665: Vec coorlocal;
666: PetscSection coordSection;
667: PetscScalar *elcoor = NULL;
668: PetscReal *swarm_coor;
669: PetscInt *swarm_cellid;
672: DMGetDimension(dmc,&dim);
674: is_simplex = PETSC_FALSE;
675: is_tensorcell = PETSC_FALSE;
676: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
677: DMPlexGetConeSize(dmc, ps, &nfaces);
679: if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }
681: switch (dim) {
682: case 2:
683: if (nfaces == 4) { is_tensorcell = PETSC_TRUE; }
684: break;
685: case 3:
686: if (nfaces == 6) { is_tensorcell = PETSC_TRUE; }
687: break;
688: default:
689: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only support for 2D, 3D");
690: }
692: /* check points provided fail inside the reference cell */
693: if (is_simplex) {
694: for (p=0; p<npoints; p++) {
695: PetscReal sum;
696: for (d=0; d<dim; d++) {
697: if (xi[dim*p+d] < -1.0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points do not fail inside the simplex domain");
698: }
699: sum = 0.0;
700: for (d=0; d<dim; d++) {
701: sum += xi[dim*p+d];
702: }
703: if (sum > 0.0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points do not fail inside the simplex domain");
704: }
705: } else if (is_tensorcell) {
706: for (p=0; p<npoints; p++) {
707: for (d=0; d<dim; d++) {
708: if (PetscAbsReal(xi[dim*p+d]) > 1.0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points do not fail inside the tensor domain [-1,1]^d");
709: }
710: }
711: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only support for d-simplex and d-tensorcell");
713: PetscQuadratureCreate(PetscObjectComm((PetscObject)dm),&quadrature);
714: PetscMalloc1(npoints*dim,&xiq);
715: PetscArraycpy(xiq,xi,npoints*dim);
716: PetscQuadratureSetData(quadrature,dim,1,npoints,(const PetscReal*)xiq,NULL);
717: private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);
718: PetscFESetQuadrature(fe,quadrature);
719: PetscFEGetDimension(fe,&nbasis);
720: PetscFEGetCellTabulation(fe, 1, &T);
722: /* for each cell, interpolate coordaintes and insert the interpolated points coordinates into swarm */
723: /* 0->cell, 1->edge, 2->vert */
724: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
725: nel = pe - ps;
727: DMSwarmSetLocalSizes(dm,npoints*nel,-1);
728: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
729: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
731: DMGetCoordinatesLocal(dmc,&coorlocal);
732: DMGetCoordinateSection(dmc,&coordSection);
734: pcnt = 0;
735: for (e=0; e<nel; e++) {
736: DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);
738: for (p=0; p<npoints; p++) {
739: for (d=0; d<dim; d++) {
740: swarm_coor[dim*pcnt+d] = 0.0;
741: for (k=0; k<nbasis; k++) {
742: swarm_coor[dim*pcnt+d] += T->T[0][p*nbasis + k] * PetscRealPart(elcoor[dim*k+d]);
743: }
744: }
745: swarm_cellid[pcnt] = e;
746: pcnt++;
747: }
748: DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);
749: }
750: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
751: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
753: PetscQuadratureDestroy(&quadrature);
754: PetscFEDestroy(&fe);
756: return(0);
757: }