Actual source code: swarmpic_da.c
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: {
9: PetscReal *xi;
10: PetscInt d,npoints=0,cnt;
11: PetscReal ds[] = {0.0,0.0,0.0};
12: PetscInt ii,jj,kk;
14: switch (dim) {
15: case 1:
16: npoints = np[0];
17: break;
18: case 2:
19: npoints = np[0]*np[1];
20: break;
21: case 3:
22: npoints = np[0]*np[1]*np[2];
23: break;
24: }
25: for (d=0; d<dim; d++) {
26: ds[d] = 2.0 / ((PetscReal)np[d]);
27: }
29: PetscMalloc1(dim*npoints,&xi);
30: switch (dim) {
31: case 1:
32: cnt = 0;
33: for (ii=0; ii<np[0]; ii++) {
34: xi[dim*cnt+0] = -1.0 + 0.5*ds[d] + ii*ds[0];
35: cnt++;
36: }
37: break;
39: case 2:
40: cnt = 0;
41: for (jj=0; jj<np[1]; jj++) {
42: for (ii=0; ii<np[0]; ii++) {
43: xi[dim*cnt+0] = -1.0 + 0.5*ds[0] + ii*ds[0];
44: xi[dim*cnt+1] = -1.0 + 0.5*ds[1] + jj*ds[1];
45: cnt++;
46: }
47: }
48: break;
50: case 3:
51: cnt = 0;
52: for (kk=0; kk<np[2]; kk++) {
53: for (jj=0; jj<np[1]; jj++) {
54: for (ii=0; ii<np[0]; ii++) {
55: xi[dim*cnt+0] = -1.0 + 0.5*ds[0] + ii*ds[0];
56: xi[dim*cnt+1] = -1.0 + 0.5*ds[1] + jj*ds[1];
57: xi[dim*cnt+2] = -1.0 + 0.5*ds[2] + kk*ds[2];
58: cnt++;
59: }
60: }
61: }
62: break;
63: }
64: *_npoints = npoints;
65: *_xi = xi;
66: return 0;
67: }
69: PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(PetscInt dim,PetscInt np_1d,PetscInt *_npoints,PetscReal **_xi)
70: {
71: PetscQuadrature quadrature;
72: const PetscReal *quadrature_xi;
73: PetscReal *xi;
74: PetscInt d,q,npoints_q;
76: PetscDTGaussTensorQuadrature(dim,1,np_1d,-1.0,1.0,&quadrature);
77: PetscQuadratureGetData(quadrature,NULL,NULL,&npoints_q,&quadrature_xi,NULL);
78: PetscMalloc1(dim*npoints_q,&xi);
79: for (q=0; q<npoints_q; q++) {
80: for (d=0; d<dim; d++) {
81: xi[dim*q+d] = quadrature_xi[dim*q+d];
82: }
83: }
84: PetscQuadratureDestroy(&quadrature);
85: *_npoints = npoints_q;
86: *_xi = xi;
87: return 0;
88: }
90: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA_Q1(DM dm,DM dmc,PetscInt npoints,DMSwarmPICLayoutType layout)
91: {
92: PetscInt dim,npoints_q;
93: PetscInt nel,npe,e,q,k,d;
94: const PetscInt *element_list;
95: PetscReal **basis;
96: PetscReal *xi;
97: Vec coor;
98: const PetscScalar *_coor;
99: PetscReal *elcoor;
100: PetscReal *swarm_coor;
101: PetscInt *swarm_cellid;
102: PetscInt pcnt;
104: DMGetDimension(dm,&dim);
105: switch (layout) {
106: case DMSWARMPIC_LAYOUT_REGULAR:
107: {
108: PetscInt np_dir[3];
109: np_dir[0] = np_dir[1] = np_dir[2] = npoints;
110: private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim,np_dir,&npoints_q,&xi);
111: }
112: break;
113: case DMSWARMPIC_LAYOUT_GAUSS:
114: private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(dim,npoints,&npoints_q,&xi);
115: break;
117: case DMSWARMPIC_LAYOUT_SUBDIVISION:
118: {
119: PetscInt s,nsub;
120: PetscInt np_dir[3];
121: nsub = npoints;
122: np_dir[0] = 1;
123: for (s=0; s<nsub; s++) {
124: np_dir[0] *= 2;
125: }
126: np_dir[1] = np_dir[0];
127: np_dir[2] = np_dir[0];
128: private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim,np_dir,&npoints_q,&xi);
129: }
130: break;
131: default:
132: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"A valid DMSwarmPIC layout must be provided");
133: }
135: DMDAGetElements(dmc,&nel,&npe,&element_list);
136: PetscMalloc1(dim*npe,&elcoor);
137: PetscMalloc1(npoints_q,&basis);
138: for (q=0; q<npoints_q; q++) {
139: PetscMalloc1(npe,&basis[q]);
141: switch (dim) {
142: case 1:
143: basis[q][0] = 0.5*(1.0 - xi[dim*q+0]);
144: basis[q][1] = 0.5*(1.0 + xi[dim*q+0]);
145: break;
146: case 2:
147: basis[q][0] = 0.25*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1]);
148: basis[q][1] = 0.25*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1]);
149: basis[q][2] = 0.25*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1]);
150: basis[q][3] = 0.25*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1]);
151: break;
153: case 3:
154: basis[q][0] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
155: basis[q][1] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
156: basis[q][2] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
157: basis[q][3] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
158: basis[q][4] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
159: basis[q][5] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
160: basis[q][6] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
161: basis[q][7] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
162: break;
163: }
164: }
166: DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
167: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
168: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
170: DMGetCoordinatesLocal(dmc,&coor);
171: VecGetArrayRead(coor,&_coor);
172: pcnt = 0;
173: for (e=0; e<nel; e++) {
174: const PetscInt *element = &element_list[npe*e];
176: for (k=0; k<npe; k++) {
177: for (d=0; d<dim; d++) {
178: elcoor[dim*k+d] = PetscRealPart(_coor[ dim*element[k] + d ]);
179: }
180: }
182: for (q=0; q<npoints_q; q++) {
183: for (d=0; d<dim; d++) {
184: swarm_coor[dim*pcnt+d] = 0.0;
185: }
186: for (k=0; k<npe; k++) {
187: for (d=0; d<dim; d++) {
188: swarm_coor[dim*pcnt+d] += basis[q][k] * elcoor[dim*k+d];
189: }
190: }
191: swarm_cellid[pcnt] = e;
192: pcnt++;
193: }
194: }
195: VecRestoreArrayRead(coor,&_coor);
196: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
197: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
198: DMDARestoreElements(dmc,&nel,&npe,&element_list);
200: PetscFree(xi);
201: PetscFree(elcoor);
202: for (q=0; q<npoints_q; q++) {
203: PetscFree(basis[q]);
204: }
205: PetscFree(basis);
206: return 0;
207: }
209: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)
210: {
211: DMDAElementType etype;
212: PetscInt dim;
214: DMDAGetElementType(celldm,&etype);
215: DMGetDimension(celldm,&dim);
216: switch (etype) {
217: case DMDA_ELEMENT_P1:
218: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DA support is not currently available for DMDA_ELEMENT_P1");
219: case DMDA_ELEMENT_Q1:
221: private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm,celldm,layout_param,layout);
222: break;
223: }
224: return 0;
225: }
227: PetscErrorCode DMSwarmProjectField_ApproxQ1_DA_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field)
228: {
229: Vec v_field_l,denom_l,coor_l,denom;
230: PetscScalar *_field_l,*_denom_l;
231: PetscInt k,p,e,npoints,nel,npe;
232: PetscInt *mpfield_cell;
233: PetscReal *mpfield_coor;
234: const PetscInt *element_list;
235: const PetscInt *element;
236: PetscScalar xi_p[2],Ni[4];
237: const PetscScalar *_coor;
239: VecZeroEntries(v_field);
241: DMGetLocalVector(dm,&v_field_l);
242: DMGetGlobalVector(dm,&denom);
243: DMGetLocalVector(dm,&denom_l);
244: VecZeroEntries(v_field_l);
245: VecZeroEntries(denom);
246: VecZeroEntries(denom_l);
248: VecGetArray(v_field_l,&_field_l);
249: VecGetArray(denom_l,&_denom_l);
251: DMGetCoordinatesLocal(dm,&coor_l);
252: VecGetArrayRead(coor_l,&_coor);
254: DMDAGetElements(dm,&nel,&npe,&element_list);
255: DMSwarmGetLocalSize(swarm,&npoints);
256: DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
257: DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
259: for (p=0; p<npoints; p++) {
260: PetscReal *coor_p;
261: const PetscScalar *x0;
262: const PetscScalar *x2;
263: PetscScalar dx[2];
265: e = mpfield_cell[p];
266: coor_p = &mpfield_coor[2*p];
267: element = &element_list[npe*e];
269: /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
270: x0 = &_coor[2*element[0]];
271: x2 = &_coor[2*element[2]];
273: dx[0] = x2[0] - x0[0];
274: dx[1] = x2[1] - x0[1];
276: xi_p[0] = 2.0 * (coor_p[0] - x0[0])/dx[0] - 1.0;
277: xi_p[1] = 2.0 * (coor_p[1] - x0[1])/dx[1] - 1.0;
279: /* evaluate basis functions */
280: Ni[0] = 0.25*(1.0 - xi_p[0])*(1.0 - xi_p[1]);
281: Ni[1] = 0.25*(1.0 + xi_p[0])*(1.0 - xi_p[1]);
282: Ni[2] = 0.25*(1.0 + xi_p[0])*(1.0 + xi_p[1]);
283: Ni[3] = 0.25*(1.0 - xi_p[0])*(1.0 + xi_p[1]);
285: for (k=0; k<npe; k++) {
286: _field_l[ element[k] ] += Ni[k] * swarm_field[p];
287: _denom_l[ element[k] ] += Ni[k];
288: }
289: }
291: DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
292: DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
293: DMDARestoreElements(dm,&nel,&npe,&element_list);
294: VecRestoreArrayRead(coor_l,&_coor);
295: VecRestoreArray(v_field_l,&_field_l);
296: VecRestoreArray(denom_l,&_denom_l);
298: DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);
299: DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);
300: DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);
301: DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);
303: VecPointwiseDivide(v_field,v_field,denom);
305: DMRestoreLocalVector(dm,&v_field_l);
306: DMRestoreLocalVector(dm,&denom_l);
307: DMRestoreGlobalVector(dm,&denom);
308: return 0;
309: }
311: PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[])
312: {
313: PetscInt f,dim;
314: DMDAElementType etype;
316: DMDAGetElementType(celldm,&etype);
319: DMGetDimension(swarm,&dim);
320: switch (dim) {
321: case 2:
322: for (f=0; f<nfields; f++) {
323: PetscReal *swarm_field;
325: DMSwarmDataFieldGetEntries(dfield[f],(void**)&swarm_field);
326: DMSwarmProjectField_ApproxQ1_DA_2D(swarm,swarm_field,celldm,vecs[f]);
327: }
328: break;
329: case 3:
330: SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D");
331: default:
332: break;
333: }
334: return 0;
335: }