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