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