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