1 static const char help[] = "Tests for injecting basis functions"; 2 3 #include <petscdmplex.h> 4 #include <petscfe.h> 5 #include <petscds.h> 6 7 typedef struct { 8 PetscInt its; /* Number of replications for timing */ 9 } AppCtx; 10 11 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 12 { 13 PetscFunctionBeginUser; 14 options->its = 1; 15 16 PetscOptionsBegin(comm, "", "FE Injection Options", "PETSCFE"); 17 PetscCall(PetscOptionsInt("-its", "The number of replications for timing", "ex1.c", options->its, &options->its, NULL)); 18 PetscOptionsEnd(); 19 PetscFunctionReturn(PETSC_SUCCESS); 20 } 21 22 static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 23 { 24 PetscInt d; 25 *u = 0.0; 26 for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0 * PETSC_PI * x[d]); 27 return PETSC_SUCCESS; 28 } 29 30 static void f0_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 31 { 32 PetscInt d; 33 for (d = 0; d < dim; ++d) f0[0] += -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]); 34 } 35 36 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 37 { 38 PetscInt d; 39 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 40 } 41 42 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 43 { 44 PetscInt d; 45 for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0; 46 } 47 48 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 49 { 50 PetscDS ds; 51 DMLabel label; 52 const PetscInt id = 1; 53 54 PetscFunctionBeginUser; 55 PetscCall(DMGetDS(dm, &ds)); 56 PetscCall(PetscDSSetResidual(ds, 0, f0_trig_u, f1_u)); 57 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu)); 58 PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user)); 59 PetscCall(DMGetLabel(dm, "marker", &label)); 60 if (label) PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL)); 61 PetscFunctionReturn(PETSC_SUCCESS); 62 } 63 64 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) 65 { 66 DM cdm = dm; 67 PetscFE fe; 68 char prefix[PETSC_MAX_PATH_LEN]; 69 PetscInt dim; 70 71 PetscFunctionBeginUser; 72 PetscCall(DMGetDimension(dm, &dim)); 73 PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name)); 74 PetscCall(DMCreateFEDefault(dm, dim, name ? prefix : NULL, -1, &fe)); 75 PetscCall(PetscObjectSetName((PetscObject)fe, name)); 76 /* Set discretization and boundary conditions for each mesh */ 77 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 78 PetscCall(DMCreateDS(dm)); 79 PetscCall((*setup)(dm, user)); 80 while (cdm) { 81 PetscCall(DMCopyDisc(dm, cdm)); 82 PetscCall(DMGetCoarseDM(cdm, &cdm)); 83 } 84 PetscCall(PetscFEDestroy(&fe)); 85 PetscFunctionReturn(PETSC_SUCCESS); 86 } 87 88 static PetscErrorCode PetscContainerUserDestroy_PetscFEGeom(void *ctx) 89 { 90 PetscFEGeom *geom = (PetscFEGeom *)ctx; 91 92 PetscFunctionBegin; 93 PetscCall(PetscFEGeomDestroy(&geom)); 94 PetscFunctionReturn(PETSC_SUCCESS); 95 } 96 97 PetscErrorCode CellRangeGetFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 98 { 99 char composeStr[33] = {0}; 100 PetscObjectId id; 101 PetscContainer container; 102 103 PetscFunctionBegin; 104 PetscCall(PetscObjectGetId((PetscObject)quad, &id)); 105 PetscCall(PetscSNPrintf(composeStr, 32, "CellRangeGetFEGeom_%" PetscInt64_FMT "\n", id)); 106 PetscCall(PetscObjectQuery((PetscObject)cellIS, composeStr, (PetscObject *)&container)); 107 if (container) { 108 PetscCall(PetscContainerGetPointer(container, (void **)geom)); 109 } else { 110 PetscCall(DMFieldCreateFEGeom(coordField, cellIS, quad, faceData, geom)); 111 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &container)); 112 PetscCall(PetscContainerSetPointer(container, (void *)*geom)); 113 PetscCall(PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom)); 114 PetscCall(PetscObjectCompose((PetscObject)cellIS, composeStr, (PetscObject)container)); 115 PetscCall(PetscContainerDestroy(&container)); 116 } 117 PetscFunctionReturn(PETSC_SUCCESS); 118 } 119 120 PetscErrorCode CellRangeRestoreFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 121 { 122 PetscFunctionBegin; 123 *geom = NULL; 124 PetscFunctionReturn(PETSC_SUCCESS); 125 } 126 127 static PetscErrorCode CreateFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms) 128 { 129 DMField coordField; 130 PetscInt Nf, f, maxDegree; 131 132 PetscFunctionBeginUser; 133 *affineQuad = NULL; 134 *affineGeom = NULL; 135 *quads = NULL; 136 *geoms = NULL; 137 PetscCall(PetscDSGetNumFields(ds, &Nf)); 138 PetscCall(DMGetCoordinateField(dm, &coordField)); 139 PetscCall(DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree)); 140 if (maxDegree <= 1) { 141 PetscCall(DMFieldCreateDefaultQuadrature(coordField, cellIS, affineQuad)); 142 if (*affineQuad) PetscCall(CellRangeGetFEGeom(cellIS, coordField, *affineQuad, PETSC_FALSE, affineGeom)); 143 } else { 144 PetscCall(PetscCalloc2(Nf, quads, Nf, geoms)); 145 for (f = 0; f < Nf; ++f) { 146 PetscFE fe; 147 148 PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe)); 149 PetscCall(PetscFEGetQuadrature(fe, &(*quads)[f])); 150 PetscCall(PetscObjectReference((PetscObject)(*quads)[f])); 151 PetscCall(CellRangeGetFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f])); 152 } 153 } 154 PetscFunctionReturn(PETSC_SUCCESS); 155 } 156 157 static PetscErrorCode DestroyFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms) 158 { 159 DMField coordField; 160 PetscInt Nf, f; 161 162 PetscFunctionBeginUser; 163 PetscCall(PetscDSGetNumFields(ds, &Nf)); 164 PetscCall(DMGetCoordinateField(dm, &coordField)); 165 if (*affineQuad) { 166 PetscCall(CellRangeRestoreFEGeom(cellIS, coordField, *affineQuad, PETSC_FALSE, affineGeom)); 167 PetscCall(PetscQuadratureDestroy(affineQuad)); 168 } else { 169 for (f = 0; f < Nf; ++f) { 170 PetscCall(CellRangeRestoreFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f])); 171 PetscCall(PetscQuadratureDestroy(&(*quads)[f])); 172 } 173 PetscCall(PetscFree2(*quads, *geoms)); 174 } 175 PetscFunctionReturn(PETSC_SUCCESS); 176 } 177 178 static PetscErrorCode TestEvaluation(DM dm) 179 { 180 PetscFE fe; 181 PetscSpace sp; 182 PetscReal *points; 183 PetscReal *B, *D, *H; 184 PetscInt dim, Nb, b, Nc, c, Np, p; 185 186 PetscFunctionBeginUser; 187 PetscCall(DMGetDimension(dm, &dim)); 188 PetscCall(DMGetField(dm, 0, NULL, (PetscObject *)&fe)); 189 Np = 6; 190 PetscCall(PetscMalloc1(Np * dim, &points)); 191 if (dim == 3) { 192 points[0] = -1.0; 193 points[1] = -1.0; 194 points[2] = -1.0; 195 points[3] = 1.0; 196 points[4] = -1.0; 197 points[5] = -1.0; 198 points[6] = -1.0; 199 points[7] = 1.0; 200 points[8] = -1.0; 201 points[9] = -1.0; 202 points[10] = -1.0; 203 points[11] = 1.0; 204 points[12] = 1.0; 205 points[13] = -1.0; 206 points[14] = 1.0; 207 points[15] = -1.0; 208 points[16] = 1.0; 209 points[17] = 1.0; 210 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only for 3D right now"); 211 PetscCall(PetscFEGetBasisSpace(fe, &sp)); 212 PetscCall(PetscSpaceGetDimension(sp, &Nb)); 213 PetscCall(PetscSpaceGetNumComponents(sp, &Nc)); 214 PetscCall(DMGetWorkArray(dm, Np * Nb * Nc, MPIU_REAL, &B)); 215 PetscCall(DMGetWorkArray(dm, Np * Nb * Nc * dim, MPIU_REAL, &D)); 216 PetscCall(DMGetWorkArray(dm, Np * Nb * Nc * dim * dim, MPIU_REAL, &H)); 217 PetscCall(PetscSpaceEvaluate(sp, Np, points, B, NULL, NULL /*D, H*/)); 218 for (p = 0; p < Np; ++p) { 219 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Point %" PetscInt_FMT "\n", p)); 220 for (b = 0; b < Nb; ++b) { 221 PetscCall(PetscPrintf(PETSC_COMM_SELF, "B[%" PetscInt_FMT "]:", b)); 222 for (c = 0; c < Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)B[(p * Nb + b) * Nc + c])); 223 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 224 #if 0 225 for (c = 0; c < Nc; ++c) { 226 PetscCall(PetscPrintf(PETSC_COMM_SELF, " D[%" PetscInt_FMT ",%" PetscInt_FMT "]:", b, c)); 227 for (d = 0; d < dim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double) B[((p*Nb+b)*Nc+c)*dim+d)])); 228 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 229 } 230 #endif 231 } 232 } 233 PetscCall(DMRestoreWorkArray(dm, Np * Nb, MPIU_REAL, &B)); 234 PetscCall(DMRestoreWorkArray(dm, Np * Nb * dim, MPIU_REAL, &D)); 235 PetscCall(DMRestoreWorkArray(dm, Np * Nb * dim * dim, MPIU_REAL, &H)); 236 PetscCall(PetscFree(points)); 237 PetscFunctionReturn(PETSC_SUCCESS); 238 } 239 240 static PetscErrorCode TestIntegration(DM dm, PetscInt cbs, PetscInt its) 241 { 242 PetscDS ds; 243 PetscFEGeom *chunkGeom = NULL; 244 PetscQuadrature affineQuad, *quads = NULL; 245 PetscFEGeom *affineGeom, **geoms = NULL; 246 PetscScalar *u, *elemVec; 247 IS cellIS; 248 PetscInt depth, cStart, cEnd, cell, chunkSize = cbs, Nch = 0, Nf, f, totDim, i, k; 249 250 PetscFunctionBeginUser; 251 PetscCall(DMPlexGetDepth(dm, &depth)); 252 PetscCall(DMGetStratumIS(dm, "depth", depth, &cellIS)); 253 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 254 PetscCall(DMGetCellDS(dm, cStart, &ds, NULL)); 255 PetscCall(PetscDSGetNumFields(ds, &Nf)); 256 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 257 PetscCall(CreateFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms)); 258 PetscCall(PetscMalloc2(chunkSize * totDim, &u, chunkSize * totDim, &elemVec)); 259 /* Assumptions: 260 - Single field 261 - No input data 262 - No auxiliary data 263 - No time-dependence 264 */ 265 for (i = 0; i < its; ++i) { 266 for (cell = cStart; cell < cEnd; cell += chunkSize, ++Nch) { 267 const PetscInt cS = cell, cE = PetscMin(cS + chunkSize, cEnd), Ne = cE - cS; 268 269 PetscCall(PetscArrayzero(elemVec, chunkSize * totDim)); 270 /* TODO Replace with DMPlexGetCellFields() */ 271 for (k = 0; k < chunkSize * totDim; ++k) u[k] = 1.0; 272 for (f = 0; f < Nf; ++f) { 273 PetscFormKey key; 274 PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f]; 275 /* PetscQuadrature quad = affineQuad ? affineQuad : quads[f]; */ 276 277 key.label = NULL; 278 key.value = 0; 279 key.field = f; 280 key.part = 0; 281 PetscCall(PetscFEGeomGetChunk(geom, cS, cE, &chunkGeom)); 282 PetscCall(PetscFEIntegrateResidual(ds, key, Ne, chunkGeom, u, NULL, NULL, NULL, 0.0, elemVec)); 283 } 284 } 285 } 286 PetscCall(PetscFEGeomRestoreChunk(affineGeom, cStart, cEnd, &chunkGeom)); 287 PetscCall(DestroyFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms)); 288 PetscCall(ISDestroy(&cellIS)); 289 PetscCall(PetscFree2(u, elemVec)); 290 PetscFunctionReturn(PETSC_SUCCESS); 291 } 292 293 static PetscErrorCode TestUnisolvence(DM dm) 294 { 295 Mat M; 296 Vec v; 297 298 PetscFunctionBeginUser; 299 PetscCall(DMGetLocalVector(dm, &v)); 300 PetscCall(DMRestoreLocalVector(dm, &v)); 301 PetscCall(DMCreateMassMatrix(dm, dm, &M)); 302 PetscCall(MatViewFromOptions(M, NULL, "-mass_view")); 303 PetscCall(MatDestroy(&M)); 304 PetscFunctionReturn(PETSC_SUCCESS); 305 } 306 307 int main(int argc, char **argv) 308 { 309 DM dm; 310 AppCtx ctx; 311 PetscMPIInt size; 312 313 PetscFunctionBeginUser; 314 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 315 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 316 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only."); 317 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx)); 318 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 319 PetscCall(DMSetType(dm, DMPLEX)); 320 PetscCall(DMSetFromOptions(dm)); 321 PetscCall(PetscObjectSetName((PetscObject)dm, "Mesh")); 322 PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_view")); 323 PetscCall(SetupDiscretization(dm, "field", SetupPrimalProblem, &ctx)); 324 PetscCall(TestEvaluation(dm)); 325 PetscCall(TestIntegration(dm, 1, ctx.its)); 326 PetscCall(TestUnisolvence(dm)); 327 PetscCall(DMDestroy(&dm)); 328 PetscCall(PetscFinalize()); 329 return 0; 330 } 331 332 /*TEST 333 test: 334 suffix: 0 335 args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism -field_petscspace_degree 0 336 337 test: 338 suffix: 2 339 args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism \ 340 -field_petscspace_type sum \ 341 -field_petscspace_variables 3 \ 342 -field_petscspace_components 3 \ 343 -field_petscspace_sum_spaces 2 \ 344 -field_petscspace_sum_concatenate false \ 345 -field_sumcomp_0_petscspace_variables 3 \ 346 -field_sumcomp_0_petscspace_components 3 \ 347 -field_sumcomp_0_petscspace_degree 1 \ 348 -field_sumcomp_1_petscspace_variables 3 \ 349 -field_sumcomp_1_petscspace_components 3 \ 350 -field_sumcomp_1_petscspace_type wxy \ 351 -field_petscdualspace_form_degree 0 \ 352 -field_petscdualspace_order 1 \ 353 -field_petscdualspace_components 3 354 355 TEST*/ 356