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

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

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

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

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

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

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

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

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

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

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

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

159:   PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1));
160:   PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
161:   PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));

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

169:     for (k = 0; k < npe; k++) {
170:       for (d = 0; d < dim; d++) elcoor[dim * k + d] = PetscRealPart(_coor[dim * element[k] + d]);
171:     }

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

187:   PetscCall(PetscFree(xi));
188:   PetscCall(PetscFree(elcoor));
189:   for (q = 0; q < npoints_q; q++) PetscCall(PetscFree(basis[q]));
190:   PetscCall(PetscFree(basis));
191:   PetscFunctionReturn(PETSC_SUCCESS);
192: }

194: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param)
195: {
196:   DMDAElementType etype;
197:   PetscInt        dim;

199:   PetscFunctionBegin;
200:   PetscCall(DMDAGetElementType(celldm, &etype));
201:   PetscCall(DMGetDimension(celldm, &dim));
202:   switch (etype) {
203:   case DMDA_ELEMENT_P1:
204:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DA support is not currently available for DMDA_ELEMENT_P1");
205:   case DMDA_ELEMENT_Q1:
206:     PetscCheck(dim != 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Support only available for dim = 2, 3");
207:     PetscCall(private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm, celldm, layout_param, layout));
208:     break;
209:   }
210:   PetscFunctionReturn(PETSC_SUCCESS);
211: }

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

225:   PetscFunctionBegin;
226:   PetscCall(VecZeroEntries(v_field));

228:   PetscCall(DMGetLocalVector(dm, &v_field_l));
229:   PetscCall(DMGetGlobalVector(dm, &denom));
230:   PetscCall(DMGetLocalVector(dm, &denom_l));
231:   PetscCall(VecZeroEntries(v_field_l));
232:   PetscCall(VecZeroEntries(denom));
233:   PetscCall(VecZeroEntries(denom_l));

235:   PetscCall(VecGetArray(v_field_l, &_field_l));
236:   PetscCall(VecGetArray(denom_l, &_denom_l));

238:   PetscCall(DMGetCoordinatesLocal(dm, &coor_l));
239:   PetscCall(VecGetArrayRead(coor_l, &_coor));

241:   PetscCall(DMDAGetElements(dm, &nel, &npe, &element_list));
242:   PetscCall(DMSwarmGetLocalSize(swarm, &npoints));
243:   PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
244:   PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));

246:   for (p = 0; p < npoints; p++) {
247:     PetscReal         *coor_p;
248:     const PetscScalar *x0;
249:     const PetscScalar *x2;
250:     PetscScalar        dx[2];

252:     e       = mpfield_cell[p];
253:     coor_p  = &mpfield_coor[2 * p];
254:     element = &element_list[npe * e];

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

260:     dx[0] = x2[0] - x0[0];
261:     dx[1] = x2[1] - x0[1];

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

266:     /* evaluate basis functions */
267:     Ni[0] = 0.25 * (1.0 - xi_p[0]) * (1.0 - xi_p[1]);
268:     Ni[1] = 0.25 * (1.0 + xi_p[0]) * (1.0 - xi_p[1]);
269:     Ni[2] = 0.25 * (1.0 + xi_p[0]) * (1.0 + xi_p[1]);
270:     Ni[3] = 0.25 * (1.0 - xi_p[0]) * (1.0 + xi_p[1]);

272:     for (k = 0; k < npe; k++) {
273:       _field_l[element[k]] += Ni[k] * swarm_field[p];
274:       _denom_l[element[k]] += Ni[k];
275:     }
276:   }

278:   PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));
279:   PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
280:   PetscCall(DMDARestoreElements(dm, &nel, &npe, &element_list));
281:   PetscCall(VecRestoreArrayRead(coor_l, &_coor));
282:   PetscCall(VecRestoreArray(v_field_l, &_field_l));
283:   PetscCall(VecRestoreArray(denom_l, &_denom_l));

285:   PetscCall(DMLocalToGlobalBegin(dm, v_field_l, ADD_VALUES, v_field));
286:   PetscCall(DMLocalToGlobalEnd(dm, v_field_l, ADD_VALUES, v_field));
287:   PetscCall(DMLocalToGlobalBegin(dm, denom_l, ADD_VALUES, denom));
288:   PetscCall(DMLocalToGlobalEnd(dm, denom_l, ADD_VALUES, denom));

290:   PetscCall(VecPointwiseDivide(v_field, v_field, denom));

292:   PetscCall(DMRestoreLocalVector(dm, &v_field_l));
293:   PetscCall(DMRestoreLocalVector(dm, &denom_l));
294:   PetscCall(DMRestoreGlobalVector(dm, &denom));
295:   PetscFunctionReturn(PETSC_SUCCESS);
296: }

298: PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm, DM celldm, PetscInt project_type, PetscInt nfields, DMSwarmDataField dfield[], Vec vecs[])
299: {
300:   PetscInt        f, dim;
301:   DMDAElementType etype;

303:   PetscFunctionBegin;
304:   PetscCall(DMDAGetElementType(celldm, &etype));
305:   PetscCheck(etype != DMDA_ELEMENT_P1, PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "Only Q1 DMDA supported");

307:   PetscCall(DMGetDimension(swarm, &dim));
308:   switch (dim) {
309:   case 2:
310:     for (f = 0; f < nfields; f++) {
311:       PetscReal *swarm_field;

313:       PetscCall(DMSwarmDataFieldGetEntries(dfield[f], (void **)&swarm_field));
314:       PetscCall(DMSwarmProjectField_ApproxQ1_DA_2D(swarm, swarm_field, celldm, vecs[f]));
315:     }
316:     break;
317:   case 3:
318:     SETERRQ(PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "No support for 3D");
319:   default:
320:     break;
321:   }
322:   PetscFunctionReturn(PETSC_SUCCESS);
323: }