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++) ds[d] = 2.0 / ((PetscReal)np[d]);

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

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

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

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

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

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

100:   DMGetDimension(dm, &dim);
101:   switch (layout) {
102:   case DMSWARMPIC_LAYOUT_REGULAR: {
103:     PetscInt np_dir[3];
104:     np_dir[0] = np_dir[1] = np_dir[2] = npoints;
105:     private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim, np_dir, &npoints_q, &xi);
106:   } break;
107:   case DMSWARMPIC_LAYOUT_GAUSS:
108:     private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(dim, npoints, &npoints_q, &xi);
109:     break;

111:   case DMSWARMPIC_LAYOUT_SUBDIVISION: {
112:     PetscInt s, nsub;
113:     PetscInt np_dir[3];
114:     nsub      = npoints;
115:     np_dir[0] = 1;
116:     for (s = 0; s < nsub; s++) np_dir[0] *= 2;
117:     np_dir[1] = np_dir[0];
118:     np_dir[2] = np_dir[0];
119:     private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim, np_dir, &npoints_q, &xi);
120:   } break;
121:   default:
122:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "A valid DMSwarmPIC layout must be provided");
123:   }

125:   DMDAGetElements(dmc, &nel, &npe, &element_list);
126:   PetscMalloc1(dim * npe, &elcoor);
127:   PetscMalloc1(npoints_q, &basis);
128:   for (q = 0; q < npoints_q; q++) {
129:     PetscMalloc1(npe, &basis[q]);

131:     switch (dim) {
132:     case 1:
133:       basis[q][0] = 0.5 * (1.0 - xi[dim * q + 0]);
134:       basis[q][1] = 0.5 * (1.0 + xi[dim * q + 0]);
135:       break;
136:     case 2:
137:       basis[q][0] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
138:       basis[q][1] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
139:       basis[q][2] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
140:       basis[q][3] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
141:       break;

143:     case 3:
144:       basis[q][0] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
145:       basis[q][1] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
146:       basis[q][2] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
147:       basis[q][3] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
148:       basis[q][4] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
149:       basis[q][5] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
150:       basis[q][6] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
151:       basis[q][7] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
152:       break;
153:     }
154:   }

156:   DMSwarmSetLocalSizes(dm, npoints_q * nel, -1);
157:   DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
158:   DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);

160:   DMGetCoordinatesLocal(dmc, &coor);
161:   VecGetArrayRead(coor, &_coor);
162:   pcnt = 0;
163:   for (e = 0; e < nel; e++) {
164:     const PetscInt *element = &element_list[npe * e];

166:     for (k = 0; k < npe; k++) {
167:       for (d = 0; d < dim; d++) elcoor[dim * k + d] = PetscRealPart(_coor[dim * element[k] + d]);
168:     }

170:     for (q = 0; q < npoints_q; q++) {
171:       for (d = 0; d < dim; d++) swarm_coor[dim * pcnt + d] = 0.0;
172:       for (k = 0; k < npe; k++) {
173:         for (d = 0; d < dim; d++) swarm_coor[dim * pcnt + d] += basis[q][k] * elcoor[dim * k + d];
174:       }
175:       swarm_cellid[pcnt] = e;
176:       pcnt++;
177:     }
178:   }
179:   VecRestoreArrayRead(coor, &_coor);
180:   DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
181:   DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
182:   DMDARestoreElements(dmc, &nel, &npe, &element_list);

184:   PetscFree(xi);
185:   PetscFree(elcoor);
186:   for (q = 0; q < npoints_q; q++) PetscFree(basis[q]);
187:   PetscFree(basis);
188:   return 0;
189: }

191: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param)
192: {
193:   DMDAElementType etype;
194:   PetscInt        dim;

196:   DMDAGetElementType(celldm, &etype);
197:   DMGetDimension(celldm, &dim);
198:   switch (etype) {
199:   case DMDA_ELEMENT_P1:
200:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DA support is not currently available for DMDA_ELEMENT_P1");
201:   case DMDA_ELEMENT_Q1:
203:     private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm, celldm, layout_param, layout);
204:     break;
205:   }
206:   return 0;
207: }

209: PetscErrorCode DMSwarmProjectField_ApproxQ1_DA_2D(DM swarm, PetscReal *swarm_field, DM dm, Vec v_field)
210: {
211:   Vec                v_field_l, denom_l, coor_l, denom;
212:   PetscScalar       *_field_l, *_denom_l;
213:   PetscInt           k, p, e, npoints, nel, npe;
214:   PetscInt          *mpfield_cell;
215:   PetscReal         *mpfield_coor;
216:   const PetscInt    *element_list;
217:   const PetscInt    *element;
218:   PetscScalar        xi_p[2], Ni[4];
219:   const PetscScalar *_coor;

221:   VecZeroEntries(v_field);

223:   DMGetLocalVector(dm, &v_field_l);
224:   DMGetGlobalVector(dm, &denom);
225:   DMGetLocalVector(dm, &denom_l);
226:   VecZeroEntries(v_field_l);
227:   VecZeroEntries(denom);
228:   VecZeroEntries(denom_l);

230:   VecGetArray(v_field_l, &_field_l);
231:   VecGetArray(denom_l, &_denom_l);

233:   DMGetCoordinatesLocal(dm, &coor_l);
234:   VecGetArrayRead(coor_l, &_coor);

236:   DMDAGetElements(dm, &nel, &npe, &element_list);
237:   DMSwarmGetLocalSize(swarm, &npoints);
238:   DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor);
239:   DMSwarmGetField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell);

241:   for (p = 0; p < npoints; p++) {
242:     PetscReal         *coor_p;
243:     const PetscScalar *x0;
244:     const PetscScalar *x2;
245:     PetscScalar        dx[2];

247:     e       = mpfield_cell[p];
248:     coor_p  = &mpfield_coor[2 * p];
249:     element = &element_list[npe * e];

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

255:     dx[0] = x2[0] - x0[0];
256:     dx[1] = x2[1] - x0[1];

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

261:     /* evaluate basis functions */
262:     Ni[0] = 0.25 * (1.0 - xi_p[0]) * (1.0 - xi_p[1]);
263:     Ni[1] = 0.25 * (1.0 + xi_p[0]) * (1.0 - xi_p[1]);
264:     Ni[2] = 0.25 * (1.0 + xi_p[0]) * (1.0 + xi_p[1]);
265:     Ni[3] = 0.25 * (1.0 - xi_p[0]) * (1.0 + xi_p[1]);

267:     for (k = 0; k < npe; k++) {
268:       _field_l[element[k]] += Ni[k] * swarm_field[p];
269:       _denom_l[element[k]] += Ni[k];
270:     }
271:   }

273:   DMSwarmRestoreField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell);
274:   DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor);
275:   DMDARestoreElements(dm, &nel, &npe, &element_list);
276:   VecRestoreArrayRead(coor_l, &_coor);
277:   VecRestoreArray(v_field_l, &_field_l);
278:   VecRestoreArray(denom_l, &_denom_l);

280:   DMLocalToGlobalBegin(dm, v_field_l, ADD_VALUES, v_field);
281:   DMLocalToGlobalEnd(dm, v_field_l, ADD_VALUES, v_field);
282:   DMLocalToGlobalBegin(dm, denom_l, ADD_VALUES, denom);
283:   DMLocalToGlobalEnd(dm, denom_l, ADD_VALUES, denom);

285:   VecPointwiseDivide(v_field, v_field, denom);

287:   DMRestoreLocalVector(dm, &v_field_l);
288:   DMRestoreLocalVector(dm, &denom_l);
289:   DMRestoreGlobalVector(dm, &denom);
290:   return 0;
291: }

293: PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm, DM celldm, PetscInt project_type, PetscInt nfields, DMSwarmDataField dfield[], Vec vecs[])
294: {
295:   PetscInt        f, dim;
296:   DMDAElementType etype;

298:   DMDAGetElementType(celldm, &etype);

301:   DMGetDimension(swarm, &dim);
302:   switch (dim) {
303:   case 2:
304:     for (f = 0; f < nfields; f++) {
305:       PetscReal *swarm_field;

307:       DMSwarmDataFieldGetEntries(dfield[f], (void **)&swarm_field);
308:       DMSwarmProjectField_ApproxQ1_DA_2D(swarm, swarm_field, celldm, vecs[f]);
309:     }
310:     break;
311:   case 3:
312:     SETERRQ(PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "No support for 3D");
313:   default:
314:     break;
315:   }
316:   return 0;
317: }