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