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
private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(PetscInt dim,PetscInt np[],PetscInt * _npoints,PetscReal ** _xi)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
private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(PetscInt dim,PetscInt np_1d,PetscInt * _npoints,PetscReal ** _xi)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
private_DMSwarmInsertPointsUsingCellDM_DA_Q1(DM dm,DM dmc,PetscInt npoints,DMSwarmPICLayoutType layout)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
private_DMSwarmInsertPointsUsingCellDM_DA(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)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
DMSwarmProjectField_ApproxQ1_DA_2D(DM swarm,PetscReal * swarm_field,DM dm,Vec v_field)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
private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[])331 PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField 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 = DMSwarmDataFieldGetEntries(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