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" 6 7 PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(PetscInt dim, PetscInt np[], PetscInt *_npoints, PetscReal **_xi) { 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: npoints = np[0]; break; 16 case 2: npoints = np[0] * np[1]; break; 17 case 3: npoints = np[0] * np[1] * np[2]; break; 18 } 19 for (d = 0; d < dim; d++) { ds[d] = 2.0 / ((PetscReal)np[d]); } 20 21 PetscCall(PetscMalloc1(dim * npoints, &xi)); 22 switch (dim) { 23 case 1: 24 cnt = 0; 25 for (ii = 0; ii < np[0]; ii++) { 26 xi[dim * cnt + 0] = -1.0 + 0.5 * ds[d] + ii * ds[0]; 27 cnt++; 28 } 29 break; 30 31 case 2: 32 cnt = 0; 33 for (jj = 0; jj < np[1]; jj++) { 34 for (ii = 0; ii < np[0]; ii++) { 35 xi[dim * cnt + 0] = -1.0 + 0.5 * ds[0] + ii * ds[0]; 36 xi[dim * cnt + 1] = -1.0 + 0.5 * ds[1] + jj * ds[1]; 37 cnt++; 38 } 39 } 40 break; 41 42 case 3: 43 cnt = 0; 44 for (kk = 0; kk < np[2]; kk++) { 45 for (jj = 0; jj < np[1]; jj++) { 46 for (ii = 0; ii < np[0]; ii++) { 47 xi[dim * cnt + 0] = -1.0 + 0.5 * ds[0] + ii * ds[0]; 48 xi[dim * cnt + 1] = -1.0 + 0.5 * ds[1] + jj * ds[1]; 49 xi[dim * cnt + 2] = -1.0 + 0.5 * ds[2] + kk * ds[2]; 50 cnt++; 51 } 52 } 53 } 54 break; 55 } 56 *_npoints = npoints; 57 *_xi = xi; 58 PetscFunctionReturn(0); 59 } 60 61 PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(PetscInt dim, PetscInt np_1d, PetscInt *_npoints, PetscReal **_xi) { 62 PetscQuadrature quadrature; 63 const PetscReal *quadrature_xi; 64 PetscReal *xi; 65 PetscInt d, q, npoints_q; 66 67 PetscFunctionBegin; 68 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, np_1d, -1.0, 1.0, &quadrature)); 69 PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &quadrature_xi, NULL)); 70 PetscCall(PetscMalloc1(dim * npoints_q, &xi)); 71 for (q = 0; q < npoints_q; q++) { 72 for (d = 0; d < dim; d++) { xi[dim * q + d] = quadrature_xi[dim * q + d]; } 73 } 74 PetscCall(PetscQuadratureDestroy(&quadrature)); 75 *_npoints = npoints_q; 76 *_xi = xi; 77 PetscFunctionReturn(0); 78 } 79 80 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA_Q1(DM dm, DM dmc, PetscInt npoints, DMSwarmPICLayoutType layout) { 81 PetscInt dim, npoints_q; 82 PetscInt nel, npe, e, q, k, d; 83 const PetscInt *element_list; 84 PetscReal **basis; 85 PetscReal *xi; 86 Vec coor; 87 const PetscScalar *_coor; 88 PetscReal *elcoor; 89 PetscReal *swarm_coor; 90 PetscInt *swarm_cellid; 91 PetscInt pcnt; 92 93 PetscFunctionBegin; 94 PetscCall(DMGetDimension(dm, &dim)); 95 switch (layout) { 96 case DMSWARMPIC_LAYOUT_REGULAR: { 97 PetscInt np_dir[3]; 98 np_dir[0] = np_dir[1] = np_dir[2] = npoints; 99 PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim, np_dir, &npoints_q, &xi)); 100 } break; 101 case DMSWARMPIC_LAYOUT_GAUSS: PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(dim, npoints, &npoints_q, &xi)); break; 102 103 case DMSWARMPIC_LAYOUT_SUBDIVISION: { 104 PetscInt s, nsub; 105 PetscInt np_dir[3]; 106 nsub = npoints; 107 np_dir[0] = 1; 108 for (s = 0; s < nsub; s++) { np_dir[0] *= 2; } 109 np_dir[1] = np_dir[0]; 110 np_dir[2] = np_dir[0]; 111 PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim, np_dir, &npoints_q, &xi)); 112 } break; 113 default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "A valid DMSwarmPIC layout must be provided"); 114 } 115 116 PetscCall(DMDAGetElements(dmc, &nel, &npe, &element_list)); 117 PetscCall(PetscMalloc1(dim * npe, &elcoor)); 118 PetscCall(PetscMalloc1(npoints_q, &basis)); 119 for (q = 0; q < npoints_q; q++) { 120 PetscCall(PetscMalloc1(npe, &basis[q])); 121 122 switch (dim) { 123 case 1: 124 basis[q][0] = 0.5 * (1.0 - xi[dim * q + 0]); 125 basis[q][1] = 0.5 * (1.0 + xi[dim * q + 0]); 126 break; 127 case 2: 128 basis[q][0] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]); 129 basis[q][1] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]); 130 basis[q][2] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]); 131 basis[q][3] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]); 132 break; 133 134 case 3: 135 basis[q][0] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]); 136 basis[q][1] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]); 137 basis[q][2] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]); 138 basis[q][3] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]); 139 basis[q][4] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]); 140 basis[q][5] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]); 141 basis[q][6] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]); 142 basis[q][7] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]); 143 break; 144 } 145 } 146 147 PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1)); 148 PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor)); 149 PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid)); 150 151 PetscCall(DMGetCoordinatesLocal(dmc, &coor)); 152 PetscCall(VecGetArrayRead(coor, &_coor)); 153 pcnt = 0; 154 for (e = 0; e < nel; e++) { 155 const PetscInt *element = &element_list[npe * e]; 156 157 for (k = 0; k < npe; k++) { 158 for (d = 0; d < dim; d++) { elcoor[dim * k + d] = PetscRealPart(_coor[dim * element[k] + d]); } 159 } 160 161 for (q = 0; q < npoints_q; q++) { 162 for (d = 0; d < dim; d++) { swarm_coor[dim * pcnt + d] = 0.0; } 163 for (k = 0; k < npe; k++) { 164 for (d = 0; d < dim; d++) { swarm_coor[dim * pcnt + d] += basis[q][k] * elcoor[dim * k + d]; } 165 } 166 swarm_cellid[pcnt] = e; 167 pcnt++; 168 } 169 } 170 PetscCall(VecRestoreArrayRead(coor, &_coor)); 171 PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid)); 172 PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor)); 173 PetscCall(DMDARestoreElements(dmc, &nel, &npe, &element_list)); 174 175 PetscCall(PetscFree(xi)); 176 PetscCall(PetscFree(elcoor)); 177 for (q = 0; q < npoints_q; q++) { PetscCall(PetscFree(basis[q])); } 178 PetscCall(PetscFree(basis)); 179 PetscFunctionReturn(0); 180 } 181 182 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param) { 183 DMDAElementType etype; 184 PetscInt dim; 185 186 PetscFunctionBegin; 187 PetscCall(DMDAGetElementType(celldm, &etype)); 188 PetscCall(DMGetDimension(celldm, &dim)); 189 switch (etype) { 190 case DMDA_ELEMENT_P1: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DA support is not currently available for DMDA_ELEMENT_P1"); 191 case DMDA_ELEMENT_Q1: 192 PetscCheck(dim != 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Support only available for dim = 2, 3"); 193 PetscCall(private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm, celldm, layout_param, layout)); 194 break; 195 } 196 PetscFunctionReturn(0); 197 } 198 199 PetscErrorCode DMSwarmProjectField_ApproxQ1_DA_2D(DM swarm, PetscReal *swarm_field, DM dm, Vec v_field) { 200 Vec v_field_l, denom_l, coor_l, denom; 201 PetscScalar *_field_l, *_denom_l; 202 PetscInt k, p, e, npoints, nel, npe; 203 PetscInt *mpfield_cell; 204 PetscReal *mpfield_coor; 205 const PetscInt *element_list; 206 const PetscInt *element; 207 PetscScalar xi_p[2], Ni[4]; 208 const PetscScalar *_coor; 209 210 PetscFunctionBegin; 211 PetscCall(VecZeroEntries(v_field)); 212 213 PetscCall(DMGetLocalVector(dm, &v_field_l)); 214 PetscCall(DMGetGlobalVector(dm, &denom)); 215 PetscCall(DMGetLocalVector(dm, &denom_l)); 216 PetscCall(VecZeroEntries(v_field_l)); 217 PetscCall(VecZeroEntries(denom)); 218 PetscCall(VecZeroEntries(denom_l)); 219 220 PetscCall(VecGetArray(v_field_l, &_field_l)); 221 PetscCall(VecGetArray(denom_l, &_denom_l)); 222 223 PetscCall(DMGetCoordinatesLocal(dm, &coor_l)); 224 PetscCall(VecGetArrayRead(coor_l, &_coor)); 225 226 PetscCall(DMDAGetElements(dm, &nel, &npe, &element_list)); 227 PetscCall(DMSwarmGetLocalSize(swarm, &npoints)); 228 PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor)); 229 PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell)); 230 231 for (p = 0; p < npoints; p++) { 232 PetscReal *coor_p; 233 const PetscScalar *x0; 234 const PetscScalar *x2; 235 PetscScalar dx[2]; 236 237 e = mpfield_cell[p]; 238 coor_p = &mpfield_coor[2 * p]; 239 element = &element_list[npe * e]; 240 241 /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */ 242 x0 = &_coor[2 * element[0]]; 243 x2 = &_coor[2 * element[2]]; 244 245 dx[0] = x2[0] - x0[0]; 246 dx[1] = x2[1] - x0[1]; 247 248 xi_p[0] = 2.0 * (coor_p[0] - x0[0]) / dx[0] - 1.0; 249 xi_p[1] = 2.0 * (coor_p[1] - x0[1]) / dx[1] - 1.0; 250 251 /* evaluate basis functions */ 252 Ni[0] = 0.25 * (1.0 - xi_p[0]) * (1.0 - xi_p[1]); 253 Ni[1] = 0.25 * (1.0 + xi_p[0]) * (1.0 - xi_p[1]); 254 Ni[2] = 0.25 * (1.0 + xi_p[0]) * (1.0 + xi_p[1]); 255 Ni[3] = 0.25 * (1.0 - xi_p[0]) * (1.0 + xi_p[1]); 256 257 for (k = 0; k < npe; k++) { 258 _field_l[element[k]] += Ni[k] * swarm_field[p]; 259 _denom_l[element[k]] += Ni[k]; 260 } 261 } 262 263 PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell)); 264 PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor)); 265 PetscCall(DMDARestoreElements(dm, &nel, &npe, &element_list)); 266 PetscCall(VecRestoreArrayRead(coor_l, &_coor)); 267 PetscCall(VecRestoreArray(v_field_l, &_field_l)); 268 PetscCall(VecRestoreArray(denom_l, &_denom_l)); 269 270 PetscCall(DMLocalToGlobalBegin(dm, v_field_l, ADD_VALUES, v_field)); 271 PetscCall(DMLocalToGlobalEnd(dm, v_field_l, ADD_VALUES, v_field)); 272 PetscCall(DMLocalToGlobalBegin(dm, denom_l, ADD_VALUES, denom)); 273 PetscCall(DMLocalToGlobalEnd(dm, denom_l, ADD_VALUES, denom)); 274 275 PetscCall(VecPointwiseDivide(v_field, v_field, denom)); 276 277 PetscCall(DMRestoreLocalVector(dm, &v_field_l)); 278 PetscCall(DMRestoreLocalVector(dm, &denom_l)); 279 PetscCall(DMRestoreGlobalVector(dm, &denom)); 280 PetscFunctionReturn(0); 281 } 282 283 PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm, DM celldm, PetscInt project_type, PetscInt nfields, DMSwarmDataField dfield[], Vec vecs[]) { 284 PetscInt f, dim; 285 DMDAElementType etype; 286 287 PetscFunctionBegin; 288 PetscCall(DMDAGetElementType(celldm, &etype)); 289 PetscCheck(etype != DMDA_ELEMENT_P1, PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "Only Q1 DMDA supported"); 290 291 PetscCall(DMGetDimension(swarm, &dim)); 292 switch (dim) { 293 case 2: 294 for (f = 0; f < nfields; f++) { 295 PetscReal *swarm_field; 296 297 PetscCall(DMSwarmDataFieldGetEntries(dfield[f], (void **)&swarm_field)); 298 PetscCall(DMSwarmProjectField_ApproxQ1_DA_2D(swarm, swarm_field, celldm, vecs[f])); 299 } 300 break; 301 case 3: SETERRQ(PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "No support for 3D"); 302 default: break; 303 } 304 PetscFunctionReturn(0); 305 } 306