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: }