Actual source code: swarmpic_da.c
petsc-3.14.6 2021-03-30
1: #include <petscdm.h>
2: #include <petscdmda.h>
3: #include <petscdmswarm.h>
4: #include <petsc/private/dmswarmimpl.h>
5: #include "../src/dm/impls/swarm/data_bucket.h"
7: PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(PetscInt dim,PetscInt np[],PetscInt *_npoints,PetscReal **_xi)
8: {
10: PetscReal *xi;
11: PetscInt d,npoints=0,cnt;
12: PetscReal ds[] = {0.0,0.0,0.0};
13: PetscInt ii,jj,kk;
16: switch (dim) {
17: case 1:
18: npoints = np[0];
19: break;
20: case 2:
21: npoints = np[0]*np[1];
22: break;
23: case 3:
24: npoints = np[0]*np[1]*np[2];
25: break;
26: }
27: for (d=0; d<dim; d++) {
28: ds[d] = 2.0 / ((PetscReal)np[d]);
29: }
31: PetscMalloc1(dim*npoints,&xi);
33: switch (dim) {
34: case 1:
35: cnt = 0;
36: for (ii=0; ii<np[0]; ii++) {
37: xi[dim*cnt+0] = -1.0 + 0.5*ds[d] + ii*ds[0];
38: cnt++;
39: }
40: break;
42: case 2:
43: cnt = 0;
44: for (jj=0; jj<np[1]; jj++) {
45: for (ii=0; ii<np[0]; ii++) {
46: xi[dim*cnt+0] = -1.0 + 0.5*ds[0] + ii*ds[0];
47: xi[dim*cnt+1] = -1.0 + 0.5*ds[1] + jj*ds[1];
48: cnt++;
49: }
50: }
51: break;
53: case 3:
54: cnt = 0;
55: for (kk=0; kk<np[2]; kk++) {
56: for (jj=0; jj<np[1]; jj++) {
57: for (ii=0; ii<np[0]; ii++) {
58: xi[dim*cnt+0] = -1.0 + 0.5*ds[0] + ii*ds[0];
59: xi[dim*cnt+1] = -1.0 + 0.5*ds[1] + jj*ds[1];
60: xi[dim*cnt+2] = -1.0 + 0.5*ds[2] + kk*ds[2];
61: cnt++;
62: }
63: }
64: }
65: break;
66: }
68: *_npoints = npoints;
69: *_xi = xi;
71: return(0);
72: }
74: PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(PetscInt dim,PetscInt np_1d,PetscInt *_npoints,PetscReal **_xi)
75: {
77: PetscQuadrature quadrature;
78: const PetscReal *quadrature_xi;
79: PetscReal *xi;
80: PetscInt d,q,npoints_q;
83: PetscDTGaussTensorQuadrature(dim,1,np_1d,-1.0,1.0,&quadrature);
84: PetscQuadratureGetData(quadrature,NULL,NULL,&npoints_q,&quadrature_xi,NULL);
85: PetscMalloc1(dim*npoints_q,&xi);
86: for (q=0; q<npoints_q; q++) {
87: for (d=0; d<dim; d++) {
88: xi[dim*q+d] = quadrature_xi[dim*q+d];
89: }
90: }
91: PetscQuadratureDestroy(&quadrature);
93: *_npoints = npoints_q;
94: *_xi = xi;
96: return(0);
97: }
99: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA_Q1(DM dm,DM dmc,PetscInt npoints,DMSwarmPICLayoutType layout)
100: {
102: PetscInt dim,npoints_q;
103: PetscInt nel,npe,e,q,k,d;
104: const PetscInt *element_list;
105: PetscReal **basis;
106: PetscReal *xi;
107: Vec coor;
108: const PetscScalar *_coor;
109: PetscReal *elcoor;
110: PetscReal *swarm_coor;
111: PetscInt *swarm_cellid;
112: PetscInt pcnt;
115: DMGetDimension(dm,&dim);
116: switch (layout) {
117: case DMSWARMPIC_LAYOUT_REGULAR:
118: {
119: PetscInt np_dir[3];
120: np_dir[0] = np_dir[1] = np_dir[2] = npoints;
121: private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim,np_dir,&npoints_q,&xi);
122: }
123: break;
124: case DMSWARMPIC_LAYOUT_GAUSS:
125: private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(dim,npoints,&npoints_q,&xi);
126: break;
128: case DMSWARMPIC_LAYOUT_SUBDIVISION:
129: {
130: PetscInt s,nsub;
131: PetscInt np_dir[3];
132: nsub = npoints;
133: np_dir[0] = 1;
134: for (s=0; s<nsub; s++) {
135: np_dir[0] *= 2;
136: }
137: np_dir[1] = np_dir[0];
138: np_dir[2] = np_dir[0];
139: private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim,np_dir,&npoints_q,&xi);
140: }
141: break;
142: default:
143: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"A valid DMSwarmPIC layout must be provided");
144: }
146: DMDAGetElements(dmc,&nel,&npe,&element_list);
148: PetscMalloc1(dim*npe,&elcoor);
149: PetscMalloc1(npoints_q,&basis);
150: for (q=0; q<npoints_q; q++) {
151: PetscMalloc1(npe,&basis[q]);
153: switch (dim) {
154: case 1:
155: basis[q][0] = 0.5*(1.0 - xi[dim*q+0]);
156: basis[q][1] = 0.5*(1.0 + xi[dim*q+0]);
157: break;
158: case 2:
159: basis[q][0] = 0.25*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1]);
160: basis[q][1] = 0.25*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1]);
161: basis[q][2] = 0.25*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1]);
162: basis[q][3] = 0.25*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1]);
163: break;
165: case 3:
166: basis[q][0] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
167: basis[q][1] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
168: basis[q][2] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
169: basis[q][3] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
170: basis[q][4] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
171: basis[q][5] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
172: basis[q][6] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
173: basis[q][7] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
174: break;
175: }
176: }
178: DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
179: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
180: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
182: DMGetCoordinatesLocal(dmc,&coor);
183: VecGetArrayRead(coor,&_coor);
184: pcnt = 0;
185: for (e=0; e<nel; e++) {
186: const PetscInt *element = &element_list[npe*e];
188: for (k=0; k<npe; k++) {
189: for (d=0; d<dim; d++) {
190: elcoor[dim*k+d] = PetscRealPart(_coor[ dim*element[k] + d ]);
191: }
192: }
194: for (q=0; q<npoints_q; q++) {
195: for (d=0; d<dim; d++) {
196: swarm_coor[dim*pcnt+d] = 0.0;
197: }
198: for (k=0; k<npe; k++) {
199: for (d=0; d<dim; d++) {
200: swarm_coor[dim*pcnt+d] += basis[q][k] * elcoor[dim*k+d];
201: }
202: }
203: swarm_cellid[pcnt] = e;
204: pcnt++;
205: }
206: }
207: VecRestoreArrayRead(coor,&_coor);
208: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
209: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
210: DMDARestoreElements(dmc,&nel,&npe,&element_list);
212: PetscFree(xi);
213: PetscFree(elcoor);
214: for (q=0; q<npoints_q; q++) {
215: PetscFree(basis[q]);
216: }
217: PetscFree(basis);
219: return(0);
220: }
222: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)
223: {
225: DMDAElementType etype;
226: PetscInt dim;
229: DMDAGetElementType(celldm,&etype);
230: DMGetDimension(celldm,&dim);
231: switch (etype) {
232: case DMDA_ELEMENT_P1:
233: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DA support is not currently available for DMDA_ELEMENT_P1");
234: case DMDA_ELEMENT_Q1:
235: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Support only available for dim = 2, 3");
236: private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm,celldm,layout_param,layout);
237: break;
238: }
239: return(0);
240: }
242: PetscErrorCode DMSwarmProjectField_ApproxQ1_DA_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field)
243: {
245: Vec v_field_l,denom_l,coor_l,denom;
246: PetscScalar *_field_l,*_denom_l;
247: PetscInt k,p,e,npoints,nel,npe;
248: PetscInt *mpfield_cell;
249: PetscReal *mpfield_coor;
250: const PetscInt *element_list;
251: const PetscInt *element;
252: PetscScalar xi_p[2],Ni[4];
253: const PetscScalar *_coor;
256: VecZeroEntries(v_field);
258: DMGetLocalVector(dm,&v_field_l);
259: DMGetGlobalVector(dm,&denom);
260: DMGetLocalVector(dm,&denom_l);
261: VecZeroEntries(v_field_l);
262: VecZeroEntries(denom);
263: VecZeroEntries(denom_l);
265: VecGetArray(v_field_l,&_field_l);
266: VecGetArray(denom_l,&_denom_l);
268: DMGetCoordinatesLocal(dm,&coor_l);
269: VecGetArrayRead(coor_l,&_coor);
271: DMDAGetElements(dm,&nel,&npe,&element_list);
272: DMSwarmGetLocalSize(swarm,&npoints);
273: DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
274: DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
276: for (p=0; p<npoints; p++) {
277: PetscReal *coor_p;
278: const PetscScalar *x0;
279: const PetscScalar *x2;
280: PetscScalar dx[2];
282: e = mpfield_cell[p];
283: coor_p = &mpfield_coor[2*p];
284: element = &element_list[npe*e];
286: /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
287: x0 = &_coor[2*element[0]];
288: x2 = &_coor[2*element[2]];
290: dx[0] = x2[0] - x0[0];
291: dx[1] = x2[1] - x0[1];
293: xi_p[0] = 2.0 * (coor_p[0] - x0[0])/dx[0] - 1.0;
294: xi_p[1] = 2.0 * (coor_p[1] - x0[1])/dx[1] - 1.0;
296: /* evaluate basis functions */
297: Ni[0] = 0.25*(1.0 - xi_p[0])*(1.0 - xi_p[1]);
298: Ni[1] = 0.25*(1.0 + xi_p[0])*(1.0 - xi_p[1]);
299: Ni[2] = 0.25*(1.0 + xi_p[0])*(1.0 + xi_p[1]);
300: Ni[3] = 0.25*(1.0 - xi_p[0])*(1.0 + xi_p[1]);
302: for (k=0; k<npe; k++) {
303: _field_l[ element[k] ] += Ni[k] * swarm_field[p];
304: _denom_l[ element[k] ] += Ni[k];
305: }
306: }
308: DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
309: DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
310: DMDARestoreElements(dm,&nel,&npe,&element_list);
311: VecRestoreArrayRead(coor_l,&_coor);
312: VecRestoreArray(v_field_l,&_field_l);
313: VecRestoreArray(denom_l,&_denom_l);
315: DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);
316: DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);
317: DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);
318: DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);
320: VecPointwiseDivide(v_field,v_field,denom);
322: DMRestoreLocalVector(dm,&v_field_l);
323: DMRestoreLocalVector(dm,&denom_l);
324: DMRestoreGlobalVector(dm,&denom);
326: return(0);
327: }
329: PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[])
330: {
332: PetscInt f,dim;
333: DMDAElementType etype;
336: DMDAGetElementType(celldm,&etype);
337: if (etype == DMDA_ELEMENT_P1) SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"Only Q1 DMDA supported");
339: DMGetDimension(swarm,&dim);
340: switch (dim) {
341: case 2:
342: for (f=0; f<nfields; f++) {
343: PetscReal *swarm_field;
345: DMSwarmDataFieldGetEntries(dfield[f],(void**)&swarm_field);
346: DMSwarmProjectField_ApproxQ1_DA_2D(swarm,swarm_field,celldm,vecs[f]);
347: }
348: break;
349: case 3:
350: SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D");
351: default:
352: break;
353: }
354: return(0);
355: }