Actual source code: swarmpic_da.c
petsc-3.8.4 2018-03-24
1: #include <petscdm.h>
2: #include <petscdmda.h>
3: #include <petscdmswarm.h>
4: #include <petsc/private/dmswarmimpl.h>
5: #include "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);
32:
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;
41:
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;
52:
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;
81:
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;
95:
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;
113:
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;
127:
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: break;
145: }
146:
147: DMDAGetElements(dmc,&nel,&npe,&element_list);
148:
149: PetscMalloc1(dim*npe,&elcoor);
150: PetscMalloc1(npoints_q,&basis);
151: for (q=0; q<npoints_q; q++) {
152: PetscMalloc1(npe,&basis[q]);
153:
154: switch (dim) {
155: case 1:
156: basis[q][0] = 0.5*(1.0 - xi[dim*q+0]);
157: basis[q][1] = 0.5*(1.0 + xi[dim*q+0]);
158: break;
159: case 2:
160: basis[q][0] = 0.25*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1]);
161: basis[q][1] = 0.25*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1]);
162: basis[q][2] = 0.25*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1]);
163: basis[q][3] = 0.25*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1]);
164: break;
165:
166: case 3:
167: basis[q][0] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
168: basis[q][1] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
169: basis[q][2] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
170: basis[q][3] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
171: basis[q][4] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
172: basis[q][5] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
173: basis[q][6] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
174: basis[q][7] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
175: break;
176: }
177: }
178:
179: DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
180: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
181: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
182:
183: DMGetCoordinatesLocal(dmc,&coor);
184: VecGetArrayRead(coor,&_coor);
185: pcnt = 0;
186: for (e=0; e<nel; e++) {
187: const PetscInt *element = &element_list[npe*e];
188:
189: for (k=0; k<npe; k++) {
190: for (d=0; d<dim; d++) {
191: elcoor[dim*k+d] = PetscRealPart(_coor[ dim*element[k] + d ]);
192: }
193: }
194:
195: for (q=0; q<npoints_q; q++) {
196: for (d=0; d<dim; d++) {
197: swarm_coor[dim*pcnt+d] = 0.0;
198: }
199: for (k=0; k<npe; k++) {
200: for (d=0; d<dim; d++) {
201: swarm_coor[dim*pcnt+d] += basis[q][k] * elcoor[dim*k+d];
202: }
203: }
204: swarm_cellid[pcnt] = e;
205: pcnt++;
206: }
207: }
208: VecRestoreArrayRead(coor,&_coor);
209: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
210: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
211: DMDARestoreElements(dmc,&nel,&npe,&element_list);
212:
213: PetscFree(xi);
214: PetscFree(elcoor);
215: for (q=0; q<npoints_q; q++) {
216: PetscFree(basis[q]);
217: }
218: PetscFree(basis);
219:
220: return(0);
221: }
223: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)
224: {
226: DMDAElementType etype;
227: PetscInt dim;
228:
230: DMDAGetElementType(celldm,&etype);
231: DMGetDimension(celldm,&dim);
232: switch (etype) {
233: case DMDA_ELEMENT_P1:
234: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DA support is not currently available for DMDA_ELEMENT_P1");
235: break;
236: case DMDA_ELEMENT_Q1:
237: if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Support only available for dim = 2, 3");
238: private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm,celldm,layout_param,layout);
239: break;
240: }
241: return(0);
242: }
244: PetscErrorCode DMSwarmProjectField_ApproxQ1_DA_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field)
245: {
247: Vec v_field_l,denom_l,coor_l,denom;
248: PetscScalar *_field_l,*_denom_l;
249: PetscInt k,p,e,npoints,nel,npe;
250: PetscInt *mpfield_cell;
251: PetscReal *mpfield_coor;
252: const PetscInt *element_list;
253: const PetscInt *element;
254: PetscScalar xi_p[2],Ni[4];
255: const PetscScalar *_coor;
256:
258: VecZeroEntries(v_field);
259:
260: DMGetLocalVector(dm,&v_field_l);
261: DMGetGlobalVector(dm,&denom);
262: DMGetLocalVector(dm,&denom_l);
263: VecZeroEntries(v_field_l);
264: VecZeroEntries(denom);
265: VecZeroEntries(denom_l);
266:
267: VecGetArray(v_field_l,&_field_l);
268: VecGetArray(denom_l,&_denom_l);
269:
270: DMGetCoordinatesLocal(dm,&coor_l);
271: VecGetArrayRead(coor_l,&_coor);
272:
273: DMDAGetElements(dm,&nel,&npe,&element_list);
274: DMSwarmGetLocalSize(swarm,&npoints);
275: DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
276: DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
277:
278: for (p=0; p<npoints; p++) {
279: PetscReal *coor_p;
280: const PetscScalar *x0;
281: const PetscScalar *x2;
282: PetscScalar dx[2];
283:
284: e = mpfield_cell[p];
285: coor_p = &mpfield_coor[2*p];
286: element = &element_list[npe*e];
287:
288: /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
289: x0 = &_coor[2*element[0]];
290: x2 = &_coor[2*element[2]];
291:
292: dx[0] = x2[0] - x0[0];
293: dx[1] = x2[1] - x0[1];
294:
295: xi_p[0] = 2.0 * (coor_p[0] - x0[0])/dx[0] - 1.0;
296: xi_p[1] = 2.0 * (coor_p[1] - x0[1])/dx[1] - 1.0;
297:
298: /* evaluate basis functions */
299: Ni[0] = 0.25*(1.0 - xi_p[0])*(1.0 - xi_p[1]);
300: Ni[1] = 0.25*(1.0 + xi_p[0])*(1.0 - xi_p[1]);
301: Ni[2] = 0.25*(1.0 + xi_p[0])*(1.0 + xi_p[1]);
302: Ni[3] = 0.25*(1.0 - xi_p[0])*(1.0 + xi_p[1]);
303:
304: for (k=0; k<npe; k++) {
305: _field_l[ element[k] ] += Ni[k] * swarm_field[p];
306: _denom_l[ element[k] ] += Ni[k];
307: }
308: }
309:
310: DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
311: DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
312: DMDARestoreElements(dm,&nel,&npe,&element_list);
313: VecRestoreArrayRead(coor_l,&_coor);
314: VecRestoreArray(v_field_l,&_field_l);
315: VecRestoreArray(denom_l,&_denom_l);
316:
317: DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);
318: DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);
319: DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);
320: DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);
321:
322: VecPointwiseDivide(v_field,v_field,denom);
323:
324: DMRestoreLocalVector(dm,&v_field_l);
325: DMRestoreLocalVector(dm,&denom_l);
326: DMRestoreGlobalVector(dm,&denom);
327:
328: return(0);
329: }
331: PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DataField dfield[],Vec vecs[])
332: {
334: PetscInt f,dim;
335: DMDAElementType etype;
336:
338: DMDAGetElementType(celldm,&etype);
339: if (etype == DMDA_ELEMENT_P1) SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"Only Q1 DMDA supported");
340:
341: DMGetDimension(swarm,&dim);
342: switch (dim) {
343: case 2:
344: for (f=0; f<nfields; f++) {
345: PetscReal *swarm_field;
346:
347: DataFieldGetEntries(dfield[f],(void**)&swarm_field);
348: DMSwarmProjectField_ApproxQ1_DA_2D(swarm,swarm_field,celldm,vecs[f]);
349: }
350: break;
351: case 3:
352: SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D");
353: break;
354: default:
355: break;
356: }
357: return(0);
358: }