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