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