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: }
61: else {
62: PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
63: PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
64: }
65: PetscFESetQuadrature(*fem, q);
66: PetscFESetFaceQuadrature(*fem, fq);
67: PetscQuadratureDestroy(&q);
68: PetscQuadratureDestroy(&fq);
69: return 0;
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: if (depth == max) {
78: PetscReal cx[2];
80: cx[0] = (v1[0] + v2[0] + v3[0])/3.0;
81: cx[1] = (v1[1] + v2[1] + v3[1])/3.0;
83: xi[2*(*np)+0] = cx[0];
84: xi[2*(*np)+1] = cx[1];
85: (*np)++;
86: return 0;
87: }
89: /* calculate midpoints of each side */
90: for (i = 0; i < 2; i++) {
91: v12[i] = (v1[i]+v2[i])/2.0;
92: v23[i] = (v2[i]+v3[i])/2.0;
93: v31[i] = (v3[i]+v1[i])/2.0;
94: }
96: /* recursively subdivide new triangles */
97: subdivide_triangle(v1,v12,v31,depth+1,max,xi,np);
98: subdivide_triangle(v2,v23,v12,depth+1,max,xi,np);
99: subdivide_triangle(v3,v31,v23,depth+1,max,xi,np);
100: subdivide_triangle(v12,v23,v31,depth+1,max,xi,np);
101: return 0;
102: }
104: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_SubDivide(DM dm,DM dmc,PetscInt nsub)
105: {
106: const PetscInt dim = 2;
107: PetscInt q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,depth;
108: PetscReal *xi;
109: PetscReal **basis;
110: Vec coorlocal;
111: PetscSection coordSection;
112: PetscScalar *elcoor = NULL;
113: PetscReal *swarm_coor;
114: PetscInt *swarm_cellid;
115: PetscReal v1[2],v2[2],v3[2];
117: npoints_q = 1;
118: for (d=0; d<nsub; d++) { npoints_q *= 4; }
119: PetscMalloc1(dim*npoints_q,&xi);
121: v1[0] = 0.0; v1[1] = 0.0;
122: v2[0] = 1.0; v2[1] = 0.0;
123: v3[0] = 0.0; v3[1] = 1.0;
124: depth = 0;
125: pcnt = 0;
126: subdivide_triangle(v1,v2,v3,depth,nsub,xi,&pcnt);
128: npe = 3; /* nodes per element (triangle) */
129: PetscMalloc1(npoints_q,&basis);
130: for (q=0; q<npoints_q; q++) {
131: PetscMalloc1(npe,&basis[q]);
133: basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1];
134: basis[q][1] = xi[dim*q+0];
135: basis[q][2] = xi[dim*q+1];
136: }
138: /* 0->cell, 1->edge, 2->vert */
139: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
140: nel = pe - ps;
142: DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
143: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
144: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
146: DMGetCoordinatesLocal(dmc,&coorlocal);
147: DMGetCoordinateSection(dmc,&coordSection);
149: pcnt = 0;
150: for (e=0; e<nel; e++) {
151: DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);
153: for (q=0; q<npoints_q; q++) {
154: for (d=0; d<dim; d++) {
155: swarm_coor[dim*pcnt+d] = 0.0;
156: for (k=0; k<npe; k++) {
157: swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]);
158: }
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++) {
170: PetscFree(basis[q]);
171: }
172: PetscFree(basis);
173: return 0;
174: }
176: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(DM dm,DM dmc,PetscInt nsub)
177: {
178: PetscInt dim,nfaces,nbasis;
179: PetscInt q,npoints_q,e,nel,pcnt,ps,pe,d,k,r;
180: PetscTabulation T;
181: Vec coorlocal;
182: PetscSection coordSection;
183: PetscScalar *elcoor = NULL;
184: PetscReal *swarm_coor;
185: PetscInt *swarm_cellid;
186: const PetscReal *xiq;
187: PetscQuadrature quadrature;
188: PetscFE fe,feRef;
189: PetscBool is_simplex;
191: DMGetDimension(dmc,&dim);
192: is_simplex = PETSC_FALSE;
193: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
194: DMPlexGetConeSize(dmc, ps, &nfaces);
195: if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }
197: private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);
199: for (r=0; r<nsub; r++) {
200: PetscFERefine(fe,&feRef);
201: PetscFECopyQuadrature(feRef,fe);
202: PetscFEDestroy(&feRef);
203: }
205: PetscFEGetQuadrature(fe,&quadrature);
206: PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL);
207: PetscFEGetDimension(fe,&nbasis);
208: PetscFEGetCellTabulation(fe, 1, &T);
210: /* 0->cell, 1->edge, 2->vert */
211: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
212: nel = pe - ps;
214: DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
215: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
216: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
218: DMGetCoordinatesLocal(dmc,&coorlocal);
219: DMGetCoordinateSection(dmc,&coordSection);
221: pcnt = 0;
222: for (e=0; e<nel; e++) {
223: DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);
225: for (q=0; q<npoints_q; q++) {
226: for (d=0; d<dim; d++) {
227: swarm_coor[dim*pcnt+d] = 0.0;
228: for (k=0; k<nbasis; k++) {
229: swarm_coor[dim*pcnt+d] += T->T[0][q*nbasis + k] * PetscRealPart(elcoor[dim*k+d]);
230: }
231: }
232: swarm_cellid[pcnt] = e;
233: pcnt++;
234: }
235: DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);
236: }
237: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
238: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
240: PetscFEDestroy(&fe);
241: return 0;
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: DMGetDimension(dmc,&dim);
259: is_simplex = PETSC_FALSE;
260: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
261: DMPlexGetConeSize(dmc, ps, &nfaces);
262: if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }
265: PetscMalloc1(dim*npoints*npoints,&xi);
266: pcnt = 0;
267: ds = 1.0/((PetscReal)(npoints-1));
268: ds2 = 1.0/((PetscReal)(npoints));
269: for (jj = 0; jj<npoints; jj++) {
270: for (ii=0; ii<npoints-jj; ii++) {
271: xi[dim*pcnt+0] = ii * ds;
272: xi[dim*pcnt+1] = jj * ds;
274: xi[dim*pcnt+0] *= (1.0 - 1.2*ds2);
275: xi[dim*pcnt+1] *= (1.0 - 1.2*ds2);
277: xi[dim*pcnt+0] += 0.35*ds2;
278: xi[dim*pcnt+1] += 0.35*ds2;
279: pcnt++;
280: }
281: }
282: npoints_q = pcnt;
284: npe = 3; /* nodes per element (triangle) */
285: PetscMalloc1(npoints_q,&basis);
286: for (q=0; q<npoints_q; q++) {
287: PetscMalloc1(npe,&basis[q]);
289: basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1];
290: basis[q][1] = xi[dim*q+0];
291: basis[q][2] = xi[dim*q+1];
292: }
294: /* 0->cell, 1->edge, 2->vert */
295: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
296: nel = pe - ps;
298: DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
299: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
300: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
302: DMGetCoordinatesLocal(dmc,&coorlocal);
303: DMGetCoordinateSection(dmc,&coordSection);
305: pcnt = 0;
306: for (e=0; e<nel; e++) {
307: DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);
309: for (q=0; q<npoints_q; q++) {
310: for (d=0; d<dim; d++) {
311: swarm_coor[dim*pcnt+d] = 0.0;
312: for (k=0; k<npe; k++) {
313: swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]);
314: }
315: }
316: swarm_cellid[pcnt] = e;
317: pcnt++;
318: }
319: DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);
320: }
321: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
322: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
324: PetscFree(xi);
325: for (q=0; q<npoints_q; q++) {
326: PetscFree(basis[q]);
327: }
328: PetscFree(basis);
329: return 0;
330: }
332: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)
333: {
334: PetscInt dim;
336: DMGetDimension(celldm,&dim);
337: switch (layout) {
338: case DMSWARMPIC_LAYOUT_REGULAR:
340: private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm,celldm,layout_param);
341: break;
342: case DMSWARMPIC_LAYOUT_GAUSS:
343: {
344: PetscInt npoints,npoints1,ps,pe,nfaces;
345: const PetscReal *xi;
346: PetscBool is_simplex;
347: PetscQuadrature quadrature;
349: is_simplex = PETSC_FALSE;
350: DMPlexGetHeightStratum(celldm,0,&ps,&pe);
351: DMPlexGetConeSize(celldm, ps, &nfaces);
352: if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }
354: npoints1 = layout_param;
355: if (is_simplex) {
356: PetscDTStroudConicalQuadrature(dim,1,npoints1,-1.0,1.0,&quadrature);
357: } else {
358: PetscDTGaussTensorQuadrature(dim,1,npoints1,-1.0,1.0,&quadrature);
359: }
360: PetscQuadratureGetData(quadrature,NULL,NULL,&npoints,&xi,NULL);
361: private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm,celldm,npoints,(PetscReal*)xi);
362: PetscQuadratureDestroy(&quadrature);
363: }
364: break;
365: case DMSWARMPIC_LAYOUT_SUBDIVISION:
366: private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm,celldm,layout_param);
367: break;
368: }
369: return 0;
370: }
372: /*
373: typedef struct {
374: PetscReal x,y;
375: } Point2d;
377: static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s)
378: {
379: *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y);
380: return 0;
381: }
382: */
383: /*
384: static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v)
385: {
386: PetscReal s1,s2,s3;
387: PetscBool b1, b2, b3;
389: signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f;
390: signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f;
391: signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f;
393: *v = ((b1 == b2) && (b2 == b3));
394: return 0;
395: }
396: */
397: /*
398: static PetscErrorCode _ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ)
399: {
400: PetscReal x1,y1,x2,y2,x3,y3;
401: PetscReal c,b[2],A[2][2],inv[2][2],detJ,od;
403: x1 = coords[2*0+0];
404: x2 = coords[2*1+0];
405: x3 = coords[2*2+0];
407: y1 = coords[2*0+1];
408: y2 = coords[2*1+1];
409: y3 = coords[2*2+1];
411: c = x1 - 0.5*x1 - 0.5*x1 + 0.5*x2 + 0.5*x3;
412: b[0] = xp[0] - c;
413: c = y1 - 0.5*y1 - 0.5*y1 + 0.5*y2 + 0.5*y3;
414: b[1] = xp[1] - c;
416: A[0][0] = -0.5*x1 + 0.5*x2; A[0][1] = -0.5*x1 + 0.5*x3;
417: A[1][0] = -0.5*y1 + 0.5*y2; A[1][1] = -0.5*y1 + 0.5*y3;
419: detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0];
420: *dJ = PetscAbsReal(detJ);
421: od = 1.0/detJ;
423: inv[0][0] = A[1][1] * od;
424: inv[0][1] = -A[0][1] * od;
425: inv[1][0] = -A[1][0] * od;
426: inv[1][1] = A[0][0] * od;
428: xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1];
429: xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1];
430: return 0;
431: }
432: */
434: static PetscErrorCode ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscScalar coords[],PetscReal xip[],PetscReal *dJ)
435: {
436: PetscReal x1,y1,x2,y2,x3,y3;
437: PetscReal b[2],A[2][2],inv[2][2],detJ,od;
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; A[0][1] = x3-x1;
451: A[1][0] = y2-y1; A[1][1] = y3-y1;
453: detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0];
454: *dJ = PetscAbsReal(detJ);
455: od = 1.0/detJ;
457: inv[0][0] = A[1][1] * od;
458: inv[0][1] = -A[0][1] * od;
459: inv[1][0] = -A[1][0] * od;
460: inv[1][1] = A[0][0] * od;
462: xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1];
463: xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1];
464: return 0;
465: }
467: PetscErrorCode DMSwarmProjectField_ApproxP1_PLEX_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field)
468: {
469: const PetscReal PLEX_C_EPS = 1.0e-8;
470: Vec v_field_l,denom_l,coor_l,denom;
471: PetscInt k,p,e,npoints;
472: PetscInt *mpfield_cell;
473: PetscReal *mpfield_coor;
474: PetscReal xi_p[2];
475: PetscScalar Ni[3];
476: PetscSection coordSection;
477: PetscScalar *elcoor = NULL;
479: VecZeroEntries(v_field);
481: DMGetLocalVector(dm,&v_field_l);
482: DMGetGlobalVector(dm,&denom);
483: DMGetLocalVector(dm,&denom_l);
484: VecZeroEntries(v_field_l);
485: VecZeroEntries(denom);
486: VecZeroEntries(denom_l);
488: DMGetCoordinatesLocal(dm,&coor_l);
489: DMGetCoordinateSection(dm,&coordSection);
491: DMSwarmGetLocalSize(swarm,&npoints);
492: DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
493: DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
495: for (p=0; p<npoints; p++) {
496: PetscReal *coor_p,dJ;
497: PetscScalar elfield[3];
498: PetscBool point_located;
500: e = mpfield_cell[p];
501: coor_p = &mpfield_coor[2*p];
503: DMPlexVecGetClosure(dm,coordSection,coor_l,e,NULL,&elcoor);
505: /*
506: while (!point_located && (failed_counter < 25)) {
507: PointInTriangle(point, coords[0], coords[1], coords[2], &point_located);
508: point.x = coor_p[0];
509: point.y = coor_p[1];
510: point.x += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
511: point.y += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
512: failed_counter++;
513: }
515: if (!point_located) {
516: 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);
517: }
520: else {
521: _ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);
522: xi_p[0] = 0.5*(xi_p[0] + 1.0);
523: xi_p[1] = 0.5*(xi_p[1] + 1.0);
525: 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]);
527: }
528: */
530: ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);
531: /*
532: 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]);
533: */
534: /*
535: point_located = PETSC_TRUE;
536: if (xi_p[0] < 0.0) {
537: if (xi_p[0] > -PLEX_C_EPS) {
538: xi_p[0] = 0.0;
539: } else {
540: point_located = PETSC_FALSE;
541: }
542: }
543: if (xi_p[1] < 0.0) {
544: if (xi_p[1] > -PLEX_C_EPS) {
545: xi_p[1] = 0.0;
546: } else {
547: point_located = PETSC_FALSE;
548: }
549: }
550: if (xi_p[1] > (1.0-xi_p[0])) {
551: if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) {
552: xi_p[1] = 1.0 - xi_p[0];
553: } else {
554: point_located = PETSC_FALSE;
555: }
556: }
557: if (!point_located) {
558: PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]);
559: 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]);
560: }
562: */
564: Ni[0] = 1.0 - xi_p[0] - xi_p[1];
565: Ni[1] = xi_p[0];
566: Ni[2] = xi_p[1];
568: point_located = PETSC_TRUE;
569: for (k=0; k<3; k++) {
570: if (PetscRealPart(Ni[k]) < -PLEX_C_EPS) point_located = PETSC_FALSE;
571: if (PetscRealPart(Ni[k]) > (1.0+PLEX_C_EPS)) point_located = PETSC_FALSE;
572: }
573: if (!point_located) {
574: PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",(double)xi_p[0],(double)xi_p[1]);
575: 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]));
576: }
579: for (k=0; k<3; k++) {
580: Ni[k] = Ni[k] * dJ;
581: elfield[k] = Ni[k] * swarm_field[p];
582: }
583: DMPlexVecRestoreClosure(dm,coordSection,coor_l,e,NULL,&elcoor);
585: DMPlexVecSetClosure(dm, NULL,v_field_l, e, elfield, ADD_VALUES);
586: DMPlexVecSetClosure(dm, NULL,denom_l, e, Ni, ADD_VALUES);
587: }
589: DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
590: DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
592: DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);
593: DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);
594: DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);
595: DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);
597: VecPointwiseDivide(v_field,v_field,denom);
599: DMRestoreLocalVector(dm,&v_field_l);
600: DMRestoreLocalVector(dm,&denom_l);
601: DMRestoreGlobalVector(dm,&denom);
602: return 0;
603: }
605: PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[])
606: {
607: PetscInt f,dim;
609: DMGetDimension(swarm,&dim);
610: switch (dim) {
611: case 2:
612: for (f=0; f<nfields; f++) {
613: PetscReal *swarm_field;
615: DMSwarmDataFieldGetEntries(dfield[f],(void**)&swarm_field);
616: DMSwarmProjectField_ApproxP1_PLEX_2D(swarm,swarm_field,celldm,vecs[f]);
617: }
618: break;
619: case 3:
620: SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D");
621: default:
622: break;
623: }
624: return 0;
625: }
627: PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm,DM dmc,PetscInt npoints,PetscReal xi[])
628: {
629: PetscBool is_simplex,is_tensorcell;
630: PetscInt dim,nfaces,ps,pe,p,d,nbasis,pcnt,e,k,nel;
631: PetscFE fe;
632: PetscQuadrature quadrature;
633: PetscTabulation T;
634: PetscReal *xiq;
635: Vec coorlocal;
636: PetscSection coordSection;
637: PetscScalar *elcoor = NULL;
638: PetscReal *swarm_coor;
639: PetscInt *swarm_cellid;
641: DMGetDimension(dmc,&dim);
643: is_simplex = PETSC_FALSE;
644: is_tensorcell = PETSC_FALSE;
645: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
646: DMPlexGetConeSize(dmc, ps, &nfaces);
648: if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }
650: switch (dim) {
651: case 2:
652: if (nfaces == 4) { is_tensorcell = PETSC_TRUE; }
653: break;
654: case 3:
655: if (nfaces == 6) { is_tensorcell = PETSC_TRUE; }
656: break;
657: default:
658: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only support for 2D, 3D");
659: }
661: /* check points provided fail inside the reference cell */
662: if (is_simplex) {
663: for (p=0; p<npoints; p++) {
664: PetscReal sum;
665: for (d=0; d<dim; d++) {
667: }
668: sum = 0.0;
669: for (d=0; d<dim; d++) {
670: sum += xi[dim*p+d];
671: }
673: }
674: } else if (is_tensorcell) {
675: for (p=0; p<npoints; p++) {
676: for (d=0; d<dim; d++) {
678: }
679: }
680: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only support for d-simplex and d-tensorcell");
682: PetscQuadratureCreate(PetscObjectComm((PetscObject)dm),&quadrature);
683: PetscMalloc1(npoints*dim,&xiq);
684: PetscArraycpy(xiq,xi,npoints*dim);
685: PetscQuadratureSetData(quadrature,dim,1,npoints,(const PetscReal*)xiq,NULL);
686: private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);
687: PetscFESetQuadrature(fe,quadrature);
688: PetscFEGetDimension(fe,&nbasis);
689: PetscFEGetCellTabulation(fe, 1, &T);
691: /* for each cell, interpolate coordaintes and insert the interpolated points coordinates into swarm */
692: /* 0->cell, 1->edge, 2->vert */
693: DMPlexGetHeightStratum(dmc,0,&ps,&pe);
694: nel = pe - ps;
696: DMSwarmSetLocalSizes(dm,npoints*nel,-1);
697: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
698: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
700: DMGetCoordinatesLocal(dmc,&coorlocal);
701: DMGetCoordinateSection(dmc,&coordSection);
703: pcnt = 0;
704: for (e=0; e<nel; e++) {
705: DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);
707: for (p=0; p<npoints; p++) {
708: for (d=0; d<dim; d++) {
709: swarm_coor[dim*pcnt+d] = 0.0;
710: for (k=0; k<nbasis; k++) {
711: swarm_coor[dim*pcnt+d] += T->T[0][p*nbasis + k] * PetscRealPart(elcoor[dim*k+d]);
712: }
713: }
714: swarm_cellid[pcnt] = e;
715: pcnt++;
716: }
717: DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);
718: }
719: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
720: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
722: PetscQuadratureDestroy(&quadrature);
723: PetscFEDestroy(&fe);
724: return 0;
725: }