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