1 #include <petscdm.h>
2 #include <petscdmda.h> /*I "petscdmda.h" I*/
3 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/
4 #include "../src/dm/impls/swarm/data_bucket.h"
5
private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(PetscInt dim,PetscInt np[],PetscInt * _npoints,PetscReal ** _xi)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;
12
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]);
26
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;
36
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;
47
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 }
66
private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(PetscInt dim,PetscInt np_1d,PetscInt * _npoints,PetscReal ** _xi)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;
73
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 }
86
private_DMSwarmInsertPointsUsingCellDM_DA_Q1(DM dm,DM dmc,PetscInt npoints,DMSwarmPICLayoutType layout)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;
102
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;
114
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 }
128
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]));
134
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;
146
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 }
159
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));
164
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));
168
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];
174
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 }
178
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));
192
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 }
199
private_DMSwarmInsertPointsUsingCellDM_DA(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)200 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param)
201 {
202 DMDAElementType etype;
203 PetscInt dim;
204
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 }
218