Actual source code: swarmpic_da.c

  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"

  7: PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(PetscInt dim,PetscInt np[],PetscInt *_npoints,PetscReal **_xi)
  8: {
  9:   PetscReal      *xi;
 10:   PetscInt       d,npoints=0,cnt;
 11:   PetscReal      ds[] = {0.0,0.0,0.0};
 12:   PetscInt       ii,jj,kk;

 14:   switch (dim) {
 15:     case 1:
 16:       npoints = np[0];
 17:       break;
 18:     case 2:
 19:       npoints = np[0]*np[1];
 20:       break;
 21:     case 3:
 22:       npoints = np[0]*np[1]*np[2];
 23:       break;
 24:   }
 25:   for (d=0; d<dim; d++) {
 26:     ds[d] = 2.0 / ((PetscReal)np[d]);
 27:   }

 29:   PetscMalloc1(dim*npoints,&xi);
 30:   switch (dim) {
 31:     case 1:
 32:       cnt = 0;
 33:       for (ii=0; ii<np[0]; ii++) {
 34:         xi[dim*cnt+0] = -1.0 + 0.5*ds[d] + ii*ds[0];
 35:         cnt++;
 36:       }
 37:       break;

 39:     case 2:
 40:       cnt = 0;
 41:       for (jj=0; jj<np[1]; jj++) {
 42:         for (ii=0; ii<np[0]; ii++) {
 43:           xi[dim*cnt+0] = -1.0 + 0.5*ds[0] + ii*ds[0];
 44:           xi[dim*cnt+1] = -1.0 + 0.5*ds[1] + jj*ds[1];
 45:           cnt++;
 46:         }
 47:       }
 48:       break;

 50:     case 3:
 51:       cnt = 0;
 52:       for (kk=0; kk<np[2]; kk++) {
 53:         for (jj=0; jj<np[1]; jj++) {
 54:           for (ii=0; ii<np[0]; ii++) {
 55:             xi[dim*cnt+0] = -1.0 + 0.5*ds[0] + ii*ds[0];
 56:             xi[dim*cnt+1] = -1.0 + 0.5*ds[1] + jj*ds[1];
 57:             xi[dim*cnt+2] = -1.0 + 0.5*ds[2] + kk*ds[2];
 58:             cnt++;
 59:           }
 60:         }
 61:       }
 62:       break;
 63:   }
 64:   *_npoints = npoints;
 65:   *_xi = xi;
 66:   return 0;
 67: }

 69: PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(PetscInt dim,PetscInt np_1d,PetscInt *_npoints,PetscReal **_xi)
 70: {
 71:   PetscQuadrature quadrature;
 72:   const PetscReal *quadrature_xi;
 73:   PetscReal       *xi;
 74:   PetscInt        d,q,npoints_q;

 76:   PetscDTGaussTensorQuadrature(dim,1,np_1d,-1.0,1.0,&quadrature);
 77:   PetscQuadratureGetData(quadrature,NULL,NULL,&npoints_q,&quadrature_xi,NULL);
 78:   PetscMalloc1(dim*npoints_q,&xi);
 79:   for (q=0; q<npoints_q; q++) {
 80:     for (d=0; d<dim; d++) {
 81:       xi[dim*q+d] = quadrature_xi[dim*q+d];
 82:     }
 83:   }
 84:   PetscQuadratureDestroy(&quadrature);
 85:   *_npoints = npoints_q;
 86:   *_xi = xi;
 87:   return 0;
 88: }

 90: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA_Q1(DM dm,DM dmc,PetscInt npoints,DMSwarmPICLayoutType layout)
 91: {
 92:   PetscInt          dim,npoints_q;
 93:   PetscInt          nel,npe,e,q,k,d;
 94:   const PetscInt    *element_list;
 95:   PetscReal         **basis;
 96:   PetscReal         *xi;
 97:   Vec               coor;
 98:   const PetscScalar *_coor;
 99:   PetscReal         *elcoor;
100:   PetscReal         *swarm_coor;
101:   PetscInt          *swarm_cellid;
102:   PetscInt          pcnt;

104:   DMGetDimension(dm,&dim);
105:   switch (layout) {
106:     case DMSWARMPIC_LAYOUT_REGULAR:
107:     {
108:       PetscInt np_dir[3];
109:       np_dir[0] = np_dir[1] = np_dir[2] = npoints;
110:       private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim,np_dir,&npoints_q,&xi);
111:     }
112:       break;
113:     case DMSWARMPIC_LAYOUT_GAUSS:
114:       private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(dim,npoints,&npoints_q,&xi);
115:       break;

117:     case DMSWARMPIC_LAYOUT_SUBDIVISION:
118:     {
119:       PetscInt s,nsub;
120:       PetscInt np_dir[3];
121:       nsub = npoints;
122:       np_dir[0] = 1;
123:       for (s=0; s<nsub; s++) {
124:         np_dir[0] *= 2;
125:       }
126:       np_dir[1] = np_dir[0];
127:       np_dir[2] = np_dir[0];
128:       private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim,np_dir,&npoints_q,&xi);
129:     }
130:       break;
131:     default:
132:       SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"A valid DMSwarmPIC layout must be provided");
133:   }

135:   DMDAGetElements(dmc,&nel,&npe,&element_list);
136:   PetscMalloc1(dim*npe,&elcoor);
137:   PetscMalloc1(npoints_q,&basis);
138:   for (q=0; q<npoints_q; q++) {
139:     PetscMalloc1(npe,&basis[q]);

141:     switch (dim) {
142:       case 1:
143:         basis[q][0] = 0.5*(1.0 - xi[dim*q+0]);
144:         basis[q][1] = 0.5*(1.0 + xi[dim*q+0]);
145:         break;
146:       case 2:
147:         basis[q][0] = 0.25*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1]);
148:         basis[q][1] = 0.25*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1]);
149:         basis[q][2] = 0.25*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1]);
150:         basis[q][3] = 0.25*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1]);
151:         break;

153:       case 3:
154:         basis[q][0] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
155:         basis[q][1] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
156:         basis[q][2] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
157:         basis[q][3] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
158:         basis[q][4] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
159:         basis[q][5] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
160:         basis[q][6] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
161:         basis[q][7] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
162:         break;
163:     }
164:   }

166:   DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
167:   DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
168:   DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);

170:   DMGetCoordinatesLocal(dmc,&coor);
171:   VecGetArrayRead(coor,&_coor);
172:   pcnt = 0;
173:   for (e=0; e<nel; e++) {
174:     const PetscInt *element = &element_list[npe*e];

176:     for (k=0; k<npe; k++) {
177:       for (d=0; d<dim; d++) {
178:         elcoor[dim*k+d] = PetscRealPart(_coor[ dim*element[k] + d ]);
179:       }
180:     }

182:     for (q=0; q<npoints_q; q++) {
183:       for (d=0; d<dim; d++) {
184:         swarm_coor[dim*pcnt+d] = 0.0;
185:       }
186:       for (k=0; k<npe; k++) {
187:         for (d=0; d<dim; d++) {
188:           swarm_coor[dim*pcnt+d] += basis[q][k] * elcoor[dim*k+d];
189:         }
190:       }
191:       swarm_cellid[pcnt] = e;
192:       pcnt++;
193:     }
194:   }
195:   VecRestoreArrayRead(coor,&_coor);
196:   DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
197:   DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
198:   DMDARestoreElements(dmc,&nel,&npe,&element_list);

200:   PetscFree(xi);
201:   PetscFree(elcoor);
202:   for (q=0; q<npoints_q; q++) {
203:     PetscFree(basis[q]);
204:   }
205:   PetscFree(basis);
206:   return 0;
207: }

209: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)
210: {
211:   DMDAElementType etype;
212:   PetscInt        dim;

214:   DMDAGetElementType(celldm,&etype);
215:   DMGetDimension(celldm,&dim);
216:   switch (etype) {
217:     case DMDA_ELEMENT_P1:
218:       SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DA support is not currently available for DMDA_ELEMENT_P1");
219:     case DMDA_ELEMENT_Q1:
221:       private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm,celldm,layout_param,layout);
222:       break;
223:   }
224:   return 0;
225: }

227: PetscErrorCode DMSwarmProjectField_ApproxQ1_DA_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field)
228: {
229:   Vec               v_field_l,denom_l,coor_l,denom;
230:   PetscScalar       *_field_l,*_denom_l;
231:   PetscInt          k,p,e,npoints,nel,npe;
232:   PetscInt          *mpfield_cell;
233:   PetscReal         *mpfield_coor;
234:   const PetscInt    *element_list;
235:   const PetscInt    *element;
236:   PetscScalar       xi_p[2],Ni[4];
237:   const PetscScalar *_coor;

239:   VecZeroEntries(v_field);

241:   DMGetLocalVector(dm,&v_field_l);
242:   DMGetGlobalVector(dm,&denom);
243:   DMGetLocalVector(dm,&denom_l);
244:   VecZeroEntries(v_field_l);
245:   VecZeroEntries(denom);
246:   VecZeroEntries(denom_l);

248:   VecGetArray(v_field_l,&_field_l);
249:   VecGetArray(denom_l,&_denom_l);

251:   DMGetCoordinatesLocal(dm,&coor_l);
252:   VecGetArrayRead(coor_l,&_coor);

254:   DMDAGetElements(dm,&nel,&npe,&element_list);
255:   DMSwarmGetLocalSize(swarm,&npoints);
256:   DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
257:   DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);

259:   for (p=0; p<npoints; p++) {
260:     PetscReal         *coor_p;
261:     const PetscScalar *x0;
262:     const PetscScalar *x2;
263:     PetscScalar       dx[2];

265:     e = mpfield_cell[p];
266:     coor_p = &mpfield_coor[2*p];
267:     element = &element_list[npe*e];

269:     /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
270:     x0 = &_coor[2*element[0]];
271:     x2 = &_coor[2*element[2]];

273:     dx[0] = x2[0] - x0[0];
274:     dx[1] = x2[1] - x0[1];

276:     xi_p[0] = 2.0 * (coor_p[0] - x0[0])/dx[0] - 1.0;
277:     xi_p[1] = 2.0 * (coor_p[1] - x0[1])/dx[1] - 1.0;

279:     /* evaluate basis functions */
280:     Ni[0] = 0.25*(1.0 - xi_p[0])*(1.0 - xi_p[1]);
281:     Ni[1] = 0.25*(1.0 + xi_p[0])*(1.0 - xi_p[1]);
282:     Ni[2] = 0.25*(1.0 + xi_p[0])*(1.0 + xi_p[1]);
283:     Ni[3] = 0.25*(1.0 - xi_p[0])*(1.0 + xi_p[1]);

285:     for (k=0; k<npe; k++) {
286:       _field_l[ element[k] ] += Ni[k] * swarm_field[p];
287:       _denom_l[ element[k] ] += Ni[k];
288:     }
289:   }

291:   DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
292:   DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
293:   DMDARestoreElements(dm,&nel,&npe,&element_list);
294:   VecRestoreArrayRead(coor_l,&_coor);
295:   VecRestoreArray(v_field_l,&_field_l);
296:   VecRestoreArray(denom_l,&_denom_l);

298:   DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);
299:   DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);
300:   DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);
301:   DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);

303:   VecPointwiseDivide(v_field,v_field,denom);

305:   DMRestoreLocalVector(dm,&v_field_l);
306:   DMRestoreLocalVector(dm,&denom_l);
307:   DMRestoreGlobalVector(dm,&denom);
308:   return 0;
309: }

311: PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[])
312: {
313:   PetscInt        f,dim;
314:   DMDAElementType etype;

316:   DMDAGetElementType(celldm,&etype);

319:   DMGetDimension(swarm,&dim);
320:   switch (dim) {
321:     case 2:
322:       for (f=0; f<nfields; f++) {
323:         PetscReal *swarm_field;

325:         DMSwarmDataFieldGetEntries(dfield[f],(void**)&swarm_field);
326:         DMSwarmProjectField_ApproxQ1_DA_2D(swarm,swarm_field,celldm,vecs[f]);
327:       }
328:       break;
329:     case 3:
330:       SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D");
331:     default:
332:       break;
333:   }
334:   return 0;
335: }