xref: /petsc/src/dm/impls/swarm/swarmpic_da.c (revision ccb4e88a40f0b86eaeca07ff64c64e4de2fae686) 
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"
6 
7 PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(PetscInt dim,PetscInt np[],PetscInt *_npoints,PetscReal **_xi)
8 {
9   PetscErrorCode ierr;
10   PetscReal *xi;
11   PetscInt d,npoints=0,cnt;
12   PetscReal ds[] = {0.0,0.0,0.0};
13   PetscInt ii,jj,kk;
14 
15   PetscFunctionBegin;
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   }
30 
31   ierr = PetscMalloc1(dim*npoints,&xi);CHKERRQ(ierr);
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   }
67 
68   *_npoints = npoints;
69   *_xi = xi;
70 
71   PetscFunctionReturn(0);
72 }
73 
74 PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(PetscInt dim,PetscInt np_1d,PetscInt *_npoints,PetscReal **_xi)
75 {
76   PetscErrorCode ierr;
77   PetscQuadrature quadrature;
78   const PetscReal *quadrature_xi;
79   PetscReal *xi;
80   PetscInt d,q,npoints_q;
81 
82   PetscFunctionBegin;
83   ierr = PetscDTGaussTensorQuadrature(dim,1,np_1d,-1.0,1.0,&quadrature);CHKERRQ(ierr);
84   ierr = PetscQuadratureGetData(quadrature,NULL,NULL,&npoints_q,&quadrature_xi,NULL);CHKERRQ(ierr);
85   ierr = PetscMalloc1(dim*npoints_q,&xi);CHKERRQ(ierr);
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   ierr = PetscQuadratureDestroy(&quadrature);CHKERRQ(ierr);
92 
93   *_npoints = npoints_q;
94   *_xi = xi;
95 
96   PetscFunctionReturn(0);
97 }
98 
99 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA_Q1(DM dm,DM dmc,PetscInt npoints,DMSwarmPICLayoutType layout)
100 {
101   PetscErrorCode ierr;
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 
114   PetscFunctionBegin;
115   ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
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       ierr = private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim,np_dir,&npoints_q,&xi);CHKERRQ(ierr);
122     }
123       break;
124     case DMSWARMPIC_LAYOUT_GAUSS:
125       ierr = private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(dim,npoints,&npoints_q,&xi);CHKERRQ(ierr);
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       ierr = private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim,np_dir,&npoints_q,&xi);CHKERRQ(ierr);
140     }
141       break;
142     default:
143       SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"A valid DMSwarmPIC layout must be provided");
144   }
145 
146   ierr = DMDAGetElements(dmc,&nel,&npe,&element_list);CHKERRQ(ierr);
147 
148   ierr = PetscMalloc1(dim*npe,&elcoor);CHKERRQ(ierr);
149   ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr);
150   for (q=0; q<npoints_q; q++) {
151     ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr);
152 
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;
164 
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   }
177 
178   ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr);
179   ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
180   ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
181 
182   ierr = DMGetCoordinatesLocal(dmc,&coor);CHKERRQ(ierr);
183   ierr = VecGetArrayRead(coor,&_coor);CHKERRQ(ierr);
184   pcnt = 0;
185   for (e=0; e<nel; e++) {
186     const PetscInt *element = &element_list[npe*e];
187 
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     }
193 
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   ierr = VecRestoreArrayRead(coor,&_coor);CHKERRQ(ierr);
208   ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
209   ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
210   ierr = DMDARestoreElements(dmc,&nel,&npe,&element_list);CHKERRQ(ierr);
211 
212   ierr = PetscFree(xi);CHKERRQ(ierr);
213   ierr = PetscFree(elcoor);CHKERRQ(ierr);
214   for (q=0; q<npoints_q; q++) {
215     ierr = PetscFree(basis[q]);CHKERRQ(ierr);
216   }
217   ierr = PetscFree(basis);CHKERRQ(ierr);
218 
219   PetscFunctionReturn(0);
220 }
221 
222 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)
223 {
224   PetscErrorCode ierr;
225   DMDAElementType etype;
226   PetscInt dim;
227 
228   PetscFunctionBegin;
229   ierr = DMDAGetElementType(celldm,&etype);CHKERRQ(ierr);
230   ierr = DMGetDimension(celldm,&dim);CHKERRQ(ierr);
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       ierr = private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm,celldm,layout_param,layout);CHKERRQ(ierr);
237       break;
238   }
239   PetscFunctionReturn(0);
240 }
241 
242 PetscErrorCode DMSwarmProjectField_ApproxQ1_DA_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field)
243 {
244   PetscErrorCode ierr;
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;
254 
255   PetscFunctionBegin;
256   ierr = VecZeroEntries(v_field);CHKERRQ(ierr);
257 
258   ierr = DMGetLocalVector(dm,&v_field_l);CHKERRQ(ierr);
259   ierr = DMGetGlobalVector(dm,&denom);CHKERRQ(ierr);
260   ierr = DMGetLocalVector(dm,&denom_l);CHKERRQ(ierr);
261   ierr = VecZeroEntries(v_field_l);CHKERRQ(ierr);
262   ierr = VecZeroEntries(denom);CHKERRQ(ierr);
263   ierr = VecZeroEntries(denom_l);CHKERRQ(ierr);
264 
265   ierr = VecGetArray(v_field_l,&_field_l);CHKERRQ(ierr);
266   ierr = VecGetArray(denom_l,&_denom_l);CHKERRQ(ierr);
267 
268   ierr = DMGetCoordinatesLocal(dm,&coor_l);CHKERRQ(ierr);
269   ierr = VecGetArrayRead(coor_l,&_coor);CHKERRQ(ierr);
270 
271   ierr = DMDAGetElements(dm,&nel,&npe,&element_list);CHKERRQ(ierr);
272   ierr = DMSwarmGetLocalSize(swarm,&npoints);CHKERRQ(ierr);
273   ierr = DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr);
274   ierr = DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr);
275 
276   for (p=0; p<npoints; p++) {
277     PetscReal *coor_p;
278     const PetscScalar *x0;
279     const PetscScalar *x2;
280     PetscScalar dx[2];
281 
282     e = mpfield_cell[p];
283     coor_p = &mpfield_coor[2*p];
284     element = &element_list[npe*e];
285 
286     /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
287     x0 = &_coor[2*element[0]];
288     x2 = &_coor[2*element[2]];
289 
290     dx[0] = x2[0] - x0[0];
291     dx[1] = x2[1] - x0[1];
292 
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;
295 
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]);
301 
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   }
307 
308   ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr);
309   ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr);
310   ierr = DMDARestoreElements(dm,&nel,&npe,&element_list);CHKERRQ(ierr);
311   ierr = VecRestoreArrayRead(coor_l,&_coor);CHKERRQ(ierr);
312   ierr = VecRestoreArray(v_field_l,&_field_l);CHKERRQ(ierr);
313   ierr = VecRestoreArray(denom_l,&_denom_l);CHKERRQ(ierr);
314 
315   ierr = DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr);
316   ierr = DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr);
317   ierr = DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr);
318   ierr = DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr);
319 
320   ierr = VecPointwiseDivide(v_field,v_field,denom);CHKERRQ(ierr);
321 
322   ierr = DMRestoreLocalVector(dm,&v_field_l);CHKERRQ(ierr);
323   ierr = DMRestoreLocalVector(dm,&denom_l);CHKERRQ(ierr);
324   ierr = DMRestoreGlobalVector(dm,&denom);CHKERRQ(ierr);
325 
326   PetscFunctionReturn(0);
327 }
328 
329 PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[])
330 {
331   PetscErrorCode ierr;
332   PetscInt f,dim;
333   DMDAElementType etype;
334 
335   PetscFunctionBegin;
336   ierr = DMDAGetElementType(celldm,&etype);CHKERRQ(ierr);
337   if (etype == DMDA_ELEMENT_P1) SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"Only Q1 DMDA supported");
338 
339   ierr = DMGetDimension(swarm,&dim);CHKERRQ(ierr);
340   switch (dim) {
341     case 2:
342       for (f=0; f<nfields; f++) {
343         PetscReal *swarm_field;
344 
345         ierr = DMSwarmDataFieldGetEntries(dfield[f],(void**)&swarm_field);CHKERRQ(ierr);
346         ierr = DMSwarmProjectField_ApproxQ1_DA_2D(swarm,swarm_field,celldm,vecs[f]);CHKERRQ(ierr);
347       }
348       break;
349     case 3:
350       SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D");
351     default:
352       break;
353   }
354   PetscFunctionReturn(0);
355 }
356