xref: /petsc/src/dm/impls/swarm/swarmpic_da.c (revision 5e71baeff2f3138f93cd4f5927dfd596eb8325cc) 
1 #include <petscdm.h>
2 #include <petscdmda.h>
3 #include <petscdmswarm.h>
4 #include <petsc/private/dmswarmimpl.h>
5 #include "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       break;
145   }
146 
147   ierr = DMDAGetElements(dmc,&nel,&npe,&element_list);CHKERRQ(ierr);
148 
149   ierr = PetscMalloc1(dim*npe,&elcoor);CHKERRQ(ierr);
150   ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr);
151   for (q=0; q<npoints_q; q++) {
152     ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr);
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   ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr);
180   ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
181   ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
182 
183   ierr = DMGetCoordinatesLocal(dmc,&coor);CHKERRQ(ierr);
184   ierr = VecGetArrayRead(coor,&_coor);CHKERRQ(ierr);
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   ierr = VecRestoreArrayRead(coor,&_coor);CHKERRQ(ierr);
209   ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
210   ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
211   ierr = DMDARestoreElements(dmc,&nel,&npe,&element_list);CHKERRQ(ierr);
212 
213   ierr = PetscFree(xi);CHKERRQ(ierr);
214   ierr = PetscFree(elcoor);CHKERRQ(ierr);
215   for (q=0; q<npoints_q; q++) {
216     ierr = PetscFree(basis[q]);CHKERRQ(ierr);
217   }
218   ierr = PetscFree(basis);CHKERRQ(ierr);
219 
220   PetscFunctionReturn(0);
221 }
222 
223 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)
224 {
225   PetscErrorCode ierr;
226   DMDAElementType etype;
227   PetscInt dim;
228 
229   PetscFunctionBegin;
230   ierr = DMDAGetElementType(celldm,&etype);CHKERRQ(ierr);
231   ierr = DMGetDimension(celldm,&dim);CHKERRQ(ierr);
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       ierr = private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm,celldm,layout_param,layout);CHKERRQ(ierr);
239       break;
240   }
241   PetscFunctionReturn(0);
242 }
243 
244 PetscErrorCode DMSwarmProjectField_ApproxQ1_DA_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field)
245 {
246   PetscErrorCode ierr;
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 
257   PetscFunctionBegin;
258   ierr = VecZeroEntries(v_field);CHKERRQ(ierr);
259 
260   ierr = DMGetLocalVector(dm,&v_field_l);CHKERRQ(ierr);
261   ierr = DMGetGlobalVector(dm,&denom);CHKERRQ(ierr);
262   ierr = DMGetLocalVector(dm,&denom_l);CHKERRQ(ierr);
263   ierr = VecZeroEntries(v_field_l);CHKERRQ(ierr);
264   ierr = VecZeroEntries(denom);CHKERRQ(ierr);
265   ierr = VecZeroEntries(denom_l);CHKERRQ(ierr);
266 
267   ierr = VecGetArray(v_field_l,&_field_l);CHKERRQ(ierr);
268   ierr = VecGetArray(denom_l,&_denom_l);CHKERRQ(ierr);
269 
270   ierr = DMGetCoordinatesLocal(dm,&coor_l);CHKERRQ(ierr);
271   ierr = VecGetArrayRead(coor_l,&_coor);CHKERRQ(ierr);
272 
273   ierr = DMDAGetElements(dm,&nel,&npe,&element_list);CHKERRQ(ierr);
274   ierr = DMSwarmGetLocalSize(swarm,&npoints);CHKERRQ(ierr);
275   ierr = DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr);
276   ierr = DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr);
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   ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr);
311   ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr);
312   ierr = DMDARestoreElements(dm,&nel,&npe,&element_list);CHKERRQ(ierr);
313   ierr = VecRestoreArrayRead(coor_l,&_coor);CHKERRQ(ierr);
314   ierr = VecRestoreArray(v_field_l,&_field_l);CHKERRQ(ierr);
315   ierr = VecRestoreArray(denom_l,&_denom_l);CHKERRQ(ierr);
316 
317   ierr = DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr);
318   ierr = DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr);
319   ierr = DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr);
320   ierr = DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr);
321 
322   ierr = VecPointwiseDivide(v_field,v_field,denom);CHKERRQ(ierr);
323 
324   ierr = DMRestoreLocalVector(dm,&v_field_l);CHKERRQ(ierr);
325   ierr = DMRestoreLocalVector(dm,&denom_l);CHKERRQ(ierr);
326   ierr = DMRestoreGlobalVector(dm,&denom);CHKERRQ(ierr);
327 
328   PetscFunctionReturn(0);
329 }
330 
331 PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DataField dfield[],Vec vecs[])
332 {
333   PetscErrorCode ierr;
334   PetscInt f,dim;
335   DMDAElementType etype;
336 
337   PetscFunctionBegin;
338   ierr = DMDAGetElementType(celldm,&etype);CHKERRQ(ierr);
339   if (etype == DMDA_ELEMENT_P1) SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"Only Q1 DMDA supported");
340 
341   ierr = DMGetDimension(swarm,&dim);CHKERRQ(ierr);
342   switch (dim) {
343     case 2:
344       for (f=0; f<nfields; f++) {
345         PetscReal *swarm_field;
346 
347         ierr = DataFieldGetEntries(dfield[f],(void**)&swarm_field);CHKERRQ(ierr);
348         ierr = DMSwarmProjectField_ApproxQ1_DA_2D(swarm,swarm_field,celldm,vecs[f]);CHKERRQ(ierr);
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   PetscFunctionReturn(0);
358 }
359