Actual source code: swarmpic_da.c
1: #include <petscdm.h>
2: #include <petscdmda.h>
3: #include <petsc/private/dmswarmimpl.h>
4: #include "../src/dm/impls/swarm/data_bucket.h"
6: static PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(PetscInt dim, PetscInt np[], PetscInt *_npoints, PetscReal **_xi)
7: {
8: PetscReal *xi;
9: PetscInt d, npoints = 0, cnt;
10: PetscReal ds[] = {0.0, 0.0, 0.0};
11: PetscInt ii, jj, kk;
13: PetscFunctionBegin;
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: PetscCall(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: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: static 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: PetscFunctionBegin;
75: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, np_1d, -1.0, 1.0, &quadrature));
76: PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &quadrature_xi, NULL));
77: PetscCall(PetscMalloc1(dim * npoints_q, &xi));
78: for (q = 0; q < npoints_q; q++) {
79: for (d = 0; d < dim; d++) xi[dim * q + d] = quadrature_xi[dim * q + d];
80: }
81: PetscCall(PetscQuadratureDestroy(&quadrature));
82: *_npoints = npoints_q;
83: *_xi = xi;
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: static PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA_Q1(DM dm, DM dmc, PetscInt npoints, DMSwarmPICLayoutType layout)
88: {
89: DMSwarmCellDM celldm;
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, Nfc;
101: const char **coordFields, *cellid;
103: PetscFunctionBegin;
104: PetscCall(DMGetDimension(dm, &dim));
105: switch (layout) {
106: case DMSWARMPIC_LAYOUT_REGULAR: {
107: PetscInt np_dir[3];
108: np_dir[0] = np_dir[1] = np_dir[2] = npoints;
109: PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim, np_dir, &npoints_q, &xi));
110: } break;
111: case DMSWARMPIC_LAYOUT_GAUSS:
112: PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(dim, npoints, &npoints_q, &xi));
113: break;
115: case DMSWARMPIC_LAYOUT_SUBDIVISION: {
116: PetscInt s, nsub;
117: PetscInt np_dir[3];
118: nsub = npoints;
119: np_dir[0] = 1;
120: for (s = 0; s < nsub; s++) np_dir[0] *= 2;
121: np_dir[1] = np_dir[0];
122: np_dir[2] = np_dir[0];
123: PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim, np_dir, &npoints_q, &xi));
124: } break;
125: default:
126: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "A valid DMSwarmPIC layout must be provided");
127: }
129: PetscCall(DMDAGetElements(dmc, &nel, &npe, &element_list));
130: PetscCall(PetscMalloc1(dim * npe, &elcoor));
131: PetscCall(PetscMalloc1(npoints_q, &basis));
132: for (q = 0; q < npoints_q; q++) {
133: PetscCall(PetscMalloc1(npe, &basis[q]));
135: switch (dim) {
136: case 1:
137: basis[q][0] = 0.5 * (1.0 - xi[dim * q + 0]);
138: basis[q][1] = 0.5 * (1.0 + xi[dim * q + 0]);
139: break;
140: case 2:
141: basis[q][0] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
142: basis[q][1] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
143: basis[q][2] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
144: basis[q][3] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
145: break;
147: case 3:
148: basis[q][0] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
149: basis[q][1] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
150: basis[q][2] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
151: basis[q][3] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
152: basis[q][4] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
153: basis[q][5] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
154: basis[q][6] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
155: basis[q][7] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
156: break;
157: }
158: }
160: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
161: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
162: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
163: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
165: PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1));
166: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
167: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
169: PetscCall(DMGetCoordinatesLocal(dmc, &coor));
170: PetscCall(VecGetArrayRead(coor, &_coor));
171: pcnt = 0;
172: for (e = 0; e < nel; e++) {
173: const PetscInt *element = &element_list[npe * e];
175: for (k = 0; k < npe; k++) {
176: for (d = 0; d < dim; d++) elcoor[dim * k + d] = PetscRealPart(_coor[dim * element[k] + d]);
177: }
179: for (q = 0; q < npoints_q; q++) {
180: for (d = 0; d < dim; d++) swarm_coor[dim * pcnt + d] = 0.0;
181: for (k = 0; k < npe; k++) {
182: for (d = 0; d < dim; d++) swarm_coor[dim * pcnt + d] += basis[q][k] * elcoor[dim * k + d];
183: }
184: swarm_cellid[pcnt] = e;
185: pcnt++;
186: }
187: }
188: PetscCall(VecRestoreArrayRead(coor, &_coor));
189: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
190: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
191: PetscCall(DMDARestoreElements(dmc, &nel, &npe, &element_list));
193: PetscCall(PetscFree(xi));
194: PetscCall(PetscFree(elcoor));
195: for (q = 0; q < npoints_q; q++) PetscCall(PetscFree(basis[q]));
196: PetscCall(PetscFree(basis));
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param)
201: {
202: DMDAElementType etype;
203: PetscInt dim;
205: PetscFunctionBegin;
206: PetscCall(DMDAGetElementType(celldm, &etype));
207: PetscCall(DMGetDimension(celldm, &dim));
208: switch (etype) {
209: case DMDA_ELEMENT_P1:
210: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DA support is not currently available for DMDA_ELEMENT_P1");
211: case DMDA_ELEMENT_Q1:
212: PetscCheck(dim != 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Support only available for dim = 2, 3");
213: PetscCall(private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm, celldm, layout_param, layout));
214: break;
215: }
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }