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 { 9 PetscReal *xi; 10 PetscInt d, npoints = 0, cnt; 11 PetscReal ds[] = {0.0, 0.0, 0.0}; 12 PetscInt ii, jj, kk; 13 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]); 27 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; 37 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; 48 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(0); 66 } 67 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; 74 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(0); 86 } 87 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; 101 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; 113 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 } 127 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])); 133 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; 145 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 } 158 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)); 162 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]; 168 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 } 172 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)); 186 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(0); 192 } 193 194 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param) 195 { 196 DMDAElementType etype; 197 PetscInt dim; 198 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(0); 211 } 212 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; 224 225 PetscFunctionBegin; 226 PetscCall(VecZeroEntries(v_field)); 227 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)); 234 235 PetscCall(VecGetArray(v_field_l, &_field_l)); 236 PetscCall(VecGetArray(denom_l, &_denom_l)); 237 238 PetscCall(DMGetCoordinatesLocal(dm, &coor_l)); 239 PetscCall(VecGetArrayRead(coor_l, &_coor)); 240 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)); 245 246 for (p = 0; p < npoints; p++) { 247 PetscReal *coor_p; 248 const PetscScalar *x0; 249 const PetscScalar *x2; 250 PetscScalar dx[2]; 251 252 e = mpfield_cell[p]; 253 coor_p = &mpfield_coor[2 * p]; 254 element = &element_list[npe * e]; 255 256 /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */ 257 x0 = &_coor[2 * element[0]]; 258 x2 = &_coor[2 * element[2]]; 259 260 dx[0] = x2[0] - x0[0]; 261 dx[1] = x2[1] - x0[1]; 262 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; 265 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]); 271 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 } 277 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)); 284 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)); 289 290 PetscCall(VecPointwiseDivide(v_field, v_field, denom)); 291 292 PetscCall(DMRestoreLocalVector(dm, &v_field_l)); 293 PetscCall(DMRestoreLocalVector(dm, &denom_l)); 294 PetscCall(DMRestoreGlobalVector(dm, &denom)); 295 PetscFunctionReturn(0); 296 } 297 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; 302 303 PetscFunctionBegin; 304 PetscCall(DMDAGetElementType(celldm, &etype)); 305 PetscCheck(etype != DMDA_ELEMENT_P1, PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "Only Q1 DMDA supported"); 306 307 PetscCall(DMGetDimension(swarm, &dim)); 308 switch (dim) { 309 case 2: 310 for (f = 0; f < nfields; f++) { 311 PetscReal *swarm_field; 312 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(0); 323 } 324