1 #include <petscdm.h> 2 #include <petscdmplex.h> 3 #include <petscdmswarm.h> 4 #include "../src/dm/impls/swarm/data_bucket.h" 5 6 PetscBool SwarmProjcite = PETSC_FALSE; 7 const char SwarmProjCitation[] = "@article{PusztayKnepleyAdams2022,\n" 8 "title = {Conservative Projection Between FEM and Particle Bases},\n" 9 "author = {Joseph V. Pusztay and Matthew G. Knepley and Mark F. Adams},\n" 10 "journal = {SIAM Journal on Scientific Computing},\n" 11 "volume = {44},\n" 12 "number = {4},\n" 13 "pages = {C310--C319},\n" 14 "doi = {10.1137/21M145407},\n" 15 "year = {2022}\n}\n"; 16 17 PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM, DM, PetscInt, PetscReal *xi); 18 19 static PetscErrorCode private_PetscFECreateDefault_scalar_pk1(DM dm, PetscInt dim, PetscBool isSimplex, PetscInt qorder, PetscFE *fem) 20 { 21 const PetscInt Nc = 1; 22 PetscQuadrature q, fq; 23 DM K; 24 PetscSpace P; 25 PetscDualSpace Q; 26 PetscInt order, quadPointsPerEdge; 27 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 28 29 PetscFunctionBegin; 30 /* Create space */ 31 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)dm), &P)); 32 /* PetscCall(PetscObjectSetOptionsPrefix((PetscObject) P, prefix)); */ 33 PetscCall(PetscSpacePolynomialSetTensor(P, tensor)); 34 /* PetscCall(PetscSpaceSetFromOptions(P)); */ 35 PetscCall(PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL)); 36 PetscCall(PetscSpaceSetDegree(P, 1, PETSC_DETERMINE)); 37 PetscCall(PetscSpaceSetNumComponents(P, Nc)); 38 PetscCall(PetscSpaceSetNumVariables(P, dim)); 39 PetscCall(PetscSpaceSetUp(P)); 40 PetscCall(PetscSpaceGetDegree(P, &order, NULL)); 41 PetscCall(PetscSpacePolynomialGetTensor(P, &tensor)); 42 /* Create dual space */ 43 PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)dm), &Q)); 44 PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE)); 45 /* PetscCall(PetscObjectSetOptionsPrefix((PetscObject) Q, prefix)); */ 46 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K)); 47 PetscCall(PetscDualSpaceSetDM(Q, K)); 48 PetscCall(DMDestroy(&K)); 49 PetscCall(PetscDualSpaceSetNumComponents(Q, Nc)); 50 PetscCall(PetscDualSpaceSetOrder(Q, order)); 51 PetscCall(PetscDualSpaceLagrangeSetTensor(Q, tensor)); 52 /* PetscCall(PetscDualSpaceSetFromOptions(Q)); */ 53 PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE)); 54 PetscCall(PetscDualSpaceSetUp(Q)); 55 /* Create element */ 56 PetscCall(PetscFECreate(PetscObjectComm((PetscObject)dm), fem)); 57 /* PetscCall(PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix)); */ 58 /* PetscCall(PetscFESetFromOptions(*fem)); */ 59 PetscCall(PetscFESetType(*fem, PETSCFEBASIC)); 60 PetscCall(PetscFESetBasisSpace(*fem, P)); 61 PetscCall(PetscFESetDualSpace(*fem, Q)); 62 PetscCall(PetscFESetNumComponents(*fem, Nc)); 63 PetscCall(PetscFESetUp(*fem)); 64 PetscCall(PetscSpaceDestroy(&P)); 65 PetscCall(PetscDualSpaceDestroy(&Q)); 66 /* Create quadrature (with specified order if given) */ 67 qorder = qorder >= 0 ? qorder : order; 68 quadPointsPerEdge = PetscMax(qorder + 1, 1); 69 if (isSimplex) { 70 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q)); 71 PetscCall(PetscDTStroudConicalQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, &fq)); 72 } else { 73 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q)); 74 PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, &fq)); 75 } 76 PetscCall(PetscFESetQuadrature(*fem, q)); 77 PetscCall(PetscFESetFaceQuadrature(*fem, fq)); 78 PetscCall(PetscQuadratureDestroy(&q)); 79 PetscCall(PetscQuadratureDestroy(&fq)); 80 PetscFunctionReturn(PETSC_SUCCESS); 81 } 82 83 static PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(DM dm, DM dmc, PetscInt nsub) 84 { 85 PetscInt dim, nfaces, nbasis; 86 PetscInt q, npoints_q, e, nel, pcnt, ps, pe, d, k, r, Nfc; 87 DMSwarmCellDM celldm; 88 PetscTabulation T; 89 Vec coorlocal; 90 PetscSection coordSection; 91 PetscScalar *elcoor = NULL; 92 PetscReal *swarm_coor; 93 PetscInt *swarm_cellid; 94 const PetscReal *xiq; 95 PetscQuadrature quadrature; 96 PetscFE fe, feRef; 97 PetscBool is_simplex; 98 const char **coordFields, *cellid; 99 100 PetscFunctionBegin; 101 PetscCall(DMGetDimension(dmc, &dim)); 102 is_simplex = PETSC_FALSE; 103 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe)); 104 PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces)); 105 if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE; 106 107 PetscCall(private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe)); 108 109 for (r = 0; r < nsub; r++) { 110 PetscCall(PetscFERefine(fe, &feRef)); 111 PetscCall(PetscFECopyQuadrature(feRef, fe)); 112 PetscCall(PetscFEDestroy(&feRef)); 113 } 114 115 PetscCall(PetscFEGetQuadrature(fe, &quadrature)); 116 PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL)); 117 PetscCall(PetscFEGetDimension(fe, &nbasis)); 118 PetscCall(PetscFEGetCellTabulation(fe, 1, &T)); 119 120 /* 0->cell, 1->edge, 2->vert */ 121 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe)); 122 nel = pe - ps; 123 124 PetscCall(DMSwarmGetCellDMActive(dmc, &celldm)); 125 PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields)); 126 PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc); 127 PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid)); 128 129 PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1)); 130 PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor)); 131 PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid)); 132 133 PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal)); 134 PetscCall(DMGetCoordinateSection(dmc, &coordSection)); 135 136 pcnt = 0; 137 for (e = 0; e < nel; e++) { 138 PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor)); 139 140 for (q = 0; q < npoints_q; q++) { 141 for (d = 0; d < dim; d++) { 142 swarm_coor[dim * pcnt + d] = 0.0; 143 for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][q * nbasis + k] * PetscRealPart(elcoor[dim * k + d]); 144 } 145 swarm_cellid[pcnt] = e; 146 pcnt++; 147 } 148 PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor)); 149 } 150 PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid)); 151 PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor)); 152 153 PetscCall(PetscFEDestroy(&fe)); 154 PetscFunctionReturn(PETSC_SUCCESS); 155 } 156 157 static PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm, DM dmc, PetscInt npoints) 158 { 159 PetscInt dim; 160 PetscInt ii, jj, q, npoints_q, e, nel, npe, pcnt, ps, pe, d, k, nfaces, Nfc; 161 PetscReal *xi, ds, ds2; 162 PetscReal **basis; 163 DMSwarmCellDM celldm; 164 Vec coorlocal; 165 PetscSection coordSection; 166 PetscScalar *elcoor = NULL; 167 PetscReal *swarm_coor; 168 PetscInt *swarm_cellid; 169 PetscBool is_simplex; 170 const char **coordFields, *cellid; 171 172 PetscFunctionBegin; 173 PetscCall(DMGetDimension(dmc, &dim)); 174 PetscCheck(dim == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only 2D is supported"); 175 is_simplex = PETSC_FALSE; 176 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe)); 177 PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces)); 178 if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE; 179 PetscCheck(is_simplex, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only the simplex is supported"); 180 181 PetscCall(PetscMalloc1(dim * npoints * npoints, &xi)); 182 pcnt = 0; 183 ds = 1.0 / (PetscReal)(npoints - 1); 184 ds2 = 1.0 / (PetscReal)npoints; 185 for (jj = 0; jj < npoints; jj++) { 186 for (ii = 0; ii < npoints - jj; ii++) { 187 xi[dim * pcnt + 0] = ii * ds; 188 xi[dim * pcnt + 1] = jj * ds; 189 190 xi[dim * pcnt + 0] *= (1.0 - 1.2 * ds2); 191 xi[dim * pcnt + 1] *= (1.0 - 1.2 * ds2); 192 193 xi[dim * pcnt + 0] += 0.35 * ds2; 194 xi[dim * pcnt + 1] += 0.35 * ds2; 195 pcnt++; 196 } 197 } 198 npoints_q = pcnt; 199 200 npe = 3; /* nodes per element (triangle) */ 201 PetscCall(PetscMalloc1(npoints_q, &basis)); 202 for (q = 0; q < npoints_q; q++) { 203 PetscCall(PetscMalloc1(npe, &basis[q])); 204 205 basis[q][0] = 1.0 - xi[dim * q + 0] - xi[dim * q + 1]; 206 basis[q][1] = xi[dim * q + 0]; 207 basis[q][2] = xi[dim * q + 1]; 208 } 209 210 /* 0->cell, 1->edge, 2->vert */ 211 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe)); 212 nel = pe - ps; 213 214 PetscCall(DMSwarmGetCellDMActive(dmc, &celldm)); 215 PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields)); 216 PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc); 217 PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid)); 218 219 PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1)); 220 PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor)); 221 PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid)); 222 223 PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal)); 224 PetscCall(DMGetCoordinateSection(dmc, &coordSection)); 225 226 pcnt = 0; 227 for (e = 0; e < nel; e++) { 228 PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor)); 229 230 for (q = 0; q < npoints_q; q++) { 231 for (d = 0; d < dim; d++) { 232 swarm_coor[dim * pcnt + d] = 0.0; 233 for (k = 0; k < npe; k++) swarm_coor[dim * pcnt + d] += basis[q][k] * PetscRealPart(elcoor[dim * k + d]); 234 } 235 swarm_cellid[pcnt] = e; 236 pcnt++; 237 } 238 PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor)); 239 } 240 PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid)); 241 PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor)); 242 243 PetscCall(PetscFree(xi)); 244 for (q = 0; q < npoints_q; q++) PetscCall(PetscFree(basis[q])); 245 PetscCall(PetscFree(basis)); 246 PetscFunctionReturn(PETSC_SUCCESS); 247 } 248 249 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param) 250 { 251 PetscInt dim; 252 253 PetscFunctionBegin; 254 PetscCall(DMGetDimension(celldm, &dim)); 255 switch (layout) { 256 case DMSWARMPIC_LAYOUT_REGULAR: 257 PetscCheck(dim != 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No 3D support for REGULAR+PLEX"); 258 PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm, celldm, layout_param)); 259 break; 260 case DMSWARMPIC_LAYOUT_GAUSS: { 261 PetscQuadrature quad, facequad; 262 const PetscReal *xi; 263 DMPolytopeType ct; 264 PetscInt cStart, Nq; 265 266 PetscCall(DMPlexGetHeightStratum(celldm, 0, &cStart, NULL)); 267 PetscCall(DMPlexGetCellType(celldm, cStart, &ct)); 268 PetscCall(PetscDTCreateDefaultQuadrature(ct, layout_param, &quad, &facequad)); 269 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &xi, NULL)); 270 PetscCall(private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm, celldm, Nq, (PetscReal *)xi)); 271 PetscCall(PetscQuadratureDestroy(&quad)); 272 PetscCall(PetscQuadratureDestroy(&facequad)); 273 } break; 274 case DMSWARMPIC_LAYOUT_SUBDIVISION: 275 PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm, celldm, layout_param)); 276 break; 277 } 278 PetscFunctionReturn(PETSC_SUCCESS); 279 } 280 281 PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm, DM dmc, PetscInt npoints, PetscReal xi[]) 282 { 283 PetscBool is_simplex, is_tensorcell; 284 PetscInt dim, nfaces, ps, pe, p, d, nbasis, pcnt, e, k, nel, Nfc; 285 DMSwarmCellDM celldm; 286 PetscFE fe; 287 PetscQuadrature quadrature; 288 PetscTabulation T; 289 PetscReal *xiq; 290 Vec coorlocal; 291 PetscSection coordSection; 292 PetscScalar *elcoor = NULL; 293 PetscReal *swarm_coor; 294 PetscInt *swarm_cellid; 295 const char **coordFields, *cellid; 296 297 PetscFunctionBegin; 298 PetscCall(DMGetDimension(dmc, &dim)); 299 300 is_simplex = PETSC_FALSE; 301 is_tensorcell = PETSC_FALSE; 302 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe)); 303 PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces)); 304 305 if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE; 306 307 switch (dim) { 308 case 2: 309 if (nfaces == 4) is_tensorcell = PETSC_TRUE; 310 break; 311 case 3: 312 if (nfaces == 6) is_tensorcell = PETSC_TRUE; 313 break; 314 default: 315 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for 2D, 3D"); 316 } 317 318 /* check points provided fail inside the reference cell */ 319 if (is_simplex) { 320 for (p = 0; p < npoints; p++) { 321 PetscReal sum; 322 for (d = 0; d < dim; d++) PetscCheck(xi[dim * p + d] >= -1.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the simplex domain"); 323 sum = 0.0; 324 for (d = 0; d < dim; d++) sum += xi[dim * p + d]; 325 PetscCheck(sum <= 0.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the simplex domain"); 326 } 327 } else if (is_tensorcell) { 328 for (p = 0; p < npoints; p++) { 329 for (d = 0; d < dim; d++) PetscCheck(PetscAbsReal(xi[dim * p + d]) <= 1.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the tensor domain [-1,1]^d"); 330 } 331 } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for d-simplex and d-tensorcell"); 332 333 PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)dm), &quadrature)); 334 PetscCall(PetscMalloc1(npoints * dim, &xiq)); 335 PetscCall(PetscArraycpy(xiq, xi, npoints * dim)); 336 PetscCall(PetscQuadratureSetData(quadrature, dim, 1, npoints, (const PetscReal *)xiq, NULL)); 337 PetscCall(private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe)); 338 PetscCall(PetscFESetQuadrature(fe, quadrature)); 339 PetscCall(PetscFEGetDimension(fe, &nbasis)); 340 PetscCall(PetscFEGetCellTabulation(fe, 1, &T)); 341 342 /* for each cell, interpolate coordaintes and insert the interpolated points coordinates into swarm */ 343 /* 0->cell, 1->edge, 2->vert */ 344 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe)); 345 nel = pe - ps; 346 347 PetscCall(DMSwarmGetCellDMActive(dm, &celldm)); 348 PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields)); 349 PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dmc), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc); 350 PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid)); 351 352 PetscCall(DMSwarmSetLocalSizes(dm, npoints * nel, -1)); 353 PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor)); 354 PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid)); 355 356 PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal)); 357 PetscCall(DMGetCoordinateSection(dmc, &coordSection)); 358 359 pcnt = 0; 360 for (e = 0; e < nel; e++) { 361 PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor)); 362 363 for (p = 0; p < npoints; p++) { 364 for (d = 0; d < dim; d++) { 365 swarm_coor[dim * pcnt + d] = 0.0; 366 for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][p * nbasis + k] * PetscRealPart(elcoor[dim * k + d]); 367 } 368 swarm_cellid[pcnt] = e; 369 pcnt++; 370 } 371 PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor)); 372 } 373 PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid)); 374 PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor)); 375 376 PetscCall(PetscQuadratureDestroy(&quadrature)); 377 PetscCall(PetscFEDestroy(&fe)); 378 PetscFunctionReturn(PETSC_SUCCESS); 379 } 380