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