1 static char help[] = "Demonstration of creating and viewing DMFields objects.\n\n"; 2 3 #include <petscdmfield.h> 4 #include <petscdmplex.h> 5 #include <petscdmda.h> 6 7 static PetscErrorCode ViewResults(PetscViewer viewer, PetscInt N, PetscInt dim, PetscScalar *B, PetscScalar *D, PetscScalar *H, PetscReal *rB, PetscReal *rD, PetscReal *rH) { 8 PetscFunctionBegin; 9 PetscCall(PetscViewerASCIIPrintf(viewer, "B:\n")); 10 PetscCall(PetscScalarView(N, B, viewer)); 11 PetscCall(PetscViewerASCIIPrintf(viewer, "D:\n")); 12 PetscCall(PetscScalarView(N * dim, D, viewer)); 13 PetscCall(PetscViewerASCIIPrintf(viewer, "H:\n")); 14 PetscCall(PetscScalarView(N * dim * dim, H, viewer)); 15 16 PetscCall(PetscViewerASCIIPrintf(viewer, "rB:\n")); 17 PetscCall(PetscRealView(N, rB, viewer)); 18 PetscCall(PetscViewerASCIIPrintf(viewer, "rD:\n")); 19 PetscCall(PetscRealView(N * dim, rD, viewer)); 20 PetscCall(PetscViewerASCIIPrintf(viewer, "rH:\n")); 21 PetscCall(PetscRealView(N * dim * dim, rH, viewer)); 22 PetscFunctionReturn(0); 23 } 24 25 static PetscErrorCode TestEvaluate(DMField field, PetscInt n, PetscRandom rand) { 26 DM dm; 27 PetscInt dim, i, nc; 28 PetscScalar *B, *D, *H; 29 PetscReal *rB, *rD, *rH; 30 Vec points; 31 PetscScalar *array; 32 PetscViewer viewer; 33 MPI_Comm comm; 34 35 PetscFunctionBegin; 36 comm = PetscObjectComm((PetscObject)field); 37 PetscCall(DMFieldGetNumComponents(field, &nc)); 38 PetscCall(DMFieldGetDM(field, &dm)); 39 PetscCall(DMGetDimension(dm, &dim)); 40 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)field), n * dim, PETSC_DETERMINE, &points)); 41 PetscCall(VecSetBlockSize(points, dim)); 42 PetscCall(VecGetArray(points, &array)); 43 for (i = 0; i < n * dim; i++) PetscCall(PetscRandomGetValue(rand, &array[i])); 44 PetscCall(VecRestoreArray(points, &array)); 45 PetscCall(PetscMalloc6(n * nc, &B, n * nc, &rB, n * nc * dim, &D, n * nc * dim, &rD, n * nc * dim * dim, &H, n * nc * dim * dim, &rH)); 46 PetscCall(DMFieldEvaluate(field, points, PETSC_SCALAR, B, D, H)); 47 PetscCall(DMFieldEvaluate(field, points, PETSC_REAL, rB, rD, rH)); 48 viewer = PETSC_VIEWER_STDOUT_(comm); 49 50 PetscCall(PetscObjectSetName((PetscObject)points, "Test Points")); 51 PetscCall(VecView(points, viewer)); 52 PetscCall(ViewResults(viewer, n * nc, dim, B, D, H, rB, rD, rH)); 53 54 PetscCall(PetscFree6(B, rB, D, rD, H, rH)); 55 PetscCall(VecDestroy(&points)); 56 PetscFunctionReturn(0); 57 } 58 59 static PetscErrorCode TestEvaluateFE(DMField field, PetscInt n, PetscInt cStart, PetscInt cEnd, PetscQuadrature quad, PetscRandom rand) { 60 DM dm; 61 PetscInt dim, i, nc, nq; 62 PetscInt N; 63 PetscScalar *B, *D, *H; 64 PetscReal *rB, *rD, *rH; 65 PetscInt *cells; 66 IS cellIS; 67 PetscViewer viewer; 68 MPI_Comm comm; 69 70 PetscFunctionBegin; 71 comm = PetscObjectComm((PetscObject)field); 72 PetscCall(DMFieldGetNumComponents(field, &nc)); 73 PetscCall(DMFieldGetDM(field, &dm)); 74 PetscCall(DMGetDimension(dm, &dim)); 75 PetscCall(PetscRandomSetInterval(rand, (PetscScalar)cStart, (PetscScalar)cEnd)); 76 PetscCall(PetscMalloc1(n, &cells)); 77 for (i = 0; i < n; i++) { 78 PetscReal rc; 79 80 PetscCall(PetscRandomGetValueReal(rand, &rc)); 81 cells[i] = PetscFloorReal(rc); 82 } 83 PetscCall(ISCreateGeneral(comm, n, cells, PETSC_OWN_POINTER, &cellIS)); 84 PetscCall(PetscObjectSetName((PetscObject)cellIS, "FE Test Cells")); 85 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &nq, NULL, NULL)); 86 N = n * nq * nc; 87 PetscCall(PetscMalloc6(N, &B, N, &rB, N * dim, &D, N * dim, &rD, N * dim * dim, &H, N * dim * dim, &rH)); 88 PetscCall(DMFieldEvaluateFE(field, cellIS, quad, PETSC_SCALAR, B, D, H)); 89 PetscCall(DMFieldEvaluateFE(field, cellIS, quad, PETSC_REAL, rB, rD, rH)); 90 viewer = PETSC_VIEWER_STDOUT_(comm); 91 92 PetscCall(PetscObjectSetName((PetscObject)quad, "Test quadrature")); 93 PetscCall(PetscQuadratureView(quad, viewer)); 94 PetscCall(ISView(cellIS, viewer)); 95 PetscCall(ViewResults(viewer, N, dim, B, D, H, rB, rD, rH)); 96 97 PetscCall(PetscFree6(B, rB, D, rD, H, rH)); 98 PetscCall(ISDestroy(&cellIS)); 99 PetscFunctionReturn(0); 100 } 101 102 static PetscErrorCode TestEvaluateFV(DMField field, PetscInt n, PetscInt cStart, PetscInt cEnd, PetscRandom rand) { 103 DM dm; 104 PetscInt dim, i, nc; 105 PetscInt N; 106 PetscScalar *B, *D, *H; 107 PetscReal *rB, *rD, *rH; 108 PetscInt *cells; 109 IS cellIS; 110 PetscViewer viewer; 111 MPI_Comm comm; 112 113 PetscFunctionBegin; 114 comm = PetscObjectComm((PetscObject)field); 115 PetscCall(DMFieldGetNumComponents(field, &nc)); 116 PetscCall(DMFieldGetDM(field, &dm)); 117 PetscCall(DMGetDimension(dm, &dim)); 118 PetscCall(PetscRandomSetInterval(rand, (PetscScalar)cStart, (PetscScalar)cEnd)); 119 PetscCall(PetscMalloc1(n, &cells)); 120 for (i = 0; i < n; i++) { 121 PetscReal rc; 122 123 PetscCall(PetscRandomGetValueReal(rand, &rc)); 124 cells[i] = PetscFloorReal(rc); 125 } 126 PetscCall(ISCreateGeneral(comm, n, cells, PETSC_OWN_POINTER, &cellIS)); 127 PetscCall(PetscObjectSetName((PetscObject)cellIS, "FV Test Cells")); 128 N = n * nc; 129 PetscCall(PetscMalloc6(N, &B, N, &rB, N * dim, &D, N * dim, &rD, N * dim * dim, &H, N * dim * dim, &rH)); 130 PetscCall(DMFieldEvaluateFV(field, cellIS, PETSC_SCALAR, B, D, H)); 131 PetscCall(DMFieldEvaluateFV(field, cellIS, PETSC_REAL, rB, rD, rH)); 132 viewer = PETSC_VIEWER_STDOUT_(comm); 133 134 PetscCall(ISView(cellIS, viewer)); 135 PetscCall(ViewResults(viewer, N, dim, B, D, H, rB, rD, rH)); 136 137 PetscCall(PetscFree6(B, rB, D, rD, H, rH)); 138 PetscCall(ISDestroy(&cellIS)); 139 PetscFunctionReturn(0); 140 } 141 142 static PetscErrorCode radiusSquared(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) { 143 PetscInt i; 144 PetscReal r2 = 0.; 145 146 PetscFunctionBegin; 147 for (i = 0; i < dim; i++) r2 += PetscSqr(x[i]); 148 for (i = 0; i < Nf; i++) u[i] = (i + 1) * r2; 149 PetscFunctionReturn(0); 150 } 151 152 static PetscErrorCode TestShellEvaluate(DMField field, Vec points, PetscDataType type, void *B, void *D, void *H) { 153 Vec ctxVec = NULL; 154 const PetscScalar *mult; 155 PetscInt dim; 156 const PetscScalar *x; 157 PetscInt Nc, n, i, j, k, l; 158 159 PetscFunctionBegin; 160 PetscCall(DMFieldGetNumComponents(field, &Nc)); 161 PetscCall(DMFieldShellGetContext(field, &ctxVec)); 162 PetscCall(VecGetBlockSize(points, &dim)); 163 PetscCall(VecGetLocalSize(points, &n)); 164 n /= Nc; 165 PetscCall(VecGetArrayRead(ctxVec, &mult)); 166 PetscCall(VecGetArrayRead(points, &x)); 167 for (i = 0; i < n; i++) { 168 PetscReal r2 = 0.; 169 170 for (j = 0; j < dim; j++) r2 += PetscSqr(PetscRealPart(x[i * dim + j])); 171 for (j = 0; j < Nc; j++) { 172 PetscReal m = PetscRealPart(mult[j]); 173 if (B) { 174 if (type == PETSC_SCALAR) { 175 ((PetscScalar *)B)[i * Nc + j] = m * r2; 176 } else { 177 ((PetscReal *)B)[i * Nc + j] = m * r2; 178 } 179 } 180 if (D) { 181 if (type == PETSC_SCALAR) { 182 for (k = 0; k < dim; k++) ((PetscScalar *)D)[(i * Nc + j) * dim + k] = 2. * m * x[i * dim + k]; 183 } else { 184 for (k = 0; k < dim; k++) ((PetscReal *)D)[(i * Nc + j) * dim + k] = 2. * m * PetscRealPart(x[i * dim + k]); 185 } 186 } 187 if (H) { 188 if (type == PETSC_SCALAR) { 189 for (k = 0; k < dim; k++) 190 for (l = 0; l < dim; l++) ((PetscScalar *)H)[((i * Nc + j) * dim + k) * dim + l] = (k == l) ? 2. * m : 0.; 191 } else { 192 for (k = 0; k < dim; k++) 193 for (l = 0; l < dim; l++) ((PetscReal *)H)[((i * Nc + j) * dim + k) * dim + l] = (k == l) ? 2. * m : 0.; 194 } 195 } 196 } 197 } 198 PetscCall(VecRestoreArrayRead(points, &x)); 199 PetscCall(VecRestoreArrayRead(ctxVec, &mult)); 200 PetscFunctionReturn(0); 201 } 202 203 static PetscErrorCode TestShellDestroy(DMField field) { 204 Vec ctxVec = NULL; 205 206 PetscFunctionBegin; 207 PetscCall(DMFieldShellGetContext(field, &ctxVec)); 208 PetscCall(VecDestroy(&ctxVec)); 209 PetscFunctionReturn(0); 210 } 211 212 int main(int argc, char **argv) { 213 DM dm = NULL; 214 MPI_Comm comm; 215 char type[256] = DMPLEX; 216 PetscBool isda, isplex; 217 PetscInt dim = 2; 218 DMField field = NULL; 219 PetscInt nc = 1; 220 PetscInt cStart = -1, cEnd = -1; 221 PetscRandom rand; 222 PetscQuadrature quad = NULL; 223 PetscInt pointsPerEdge = 2; 224 PetscInt numPoint = 0, numFE = 0, numFV = 0; 225 PetscBool testShell = PETSC_FALSE; 226 227 PetscFunctionBeginUser; 228 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 229 comm = PETSC_COMM_WORLD; 230 PetscOptionsBegin(comm, "", "DMField Tutorial Options", "DM"); 231 PetscCall(PetscOptionsFList("-dm_type", "DM implementation on which to define field", "ex1.c", DMList, type, type, 256, NULL)); 232 PetscCall(PetscOptionsRangeInt("-dim", "DM intrinsic dimension", "ex1.c", dim, &dim, NULL, 1, 3)); 233 PetscCall(PetscOptionsBoundedInt("-num_components", "Number of components in field", "ex1.c", nc, &nc, NULL, 1)); 234 PetscCall(PetscOptionsBoundedInt("-num_quad_points", "Number of quadrature points per dimension", "ex1.c", pointsPerEdge, &pointsPerEdge, NULL, 1)); 235 PetscCall(PetscOptionsBoundedInt("-num_point_tests", "Number of test points for DMFieldEvaluate()", "ex1.c", numPoint, &numPoint, NULL, 0)); 236 PetscCall(PetscOptionsBoundedInt("-num_fe_tests", "Number of test cells for DMFieldEvaluateFE()", "ex1.c", numFE, &numFE, NULL, 0)); 237 PetscCall(PetscOptionsBoundedInt("-num_fv_tests", "Number of test cells for DMFieldEvaluateFV()", "ex1.c", numFV, &numFV, NULL, 0)); 238 PetscCall(PetscOptionsBool("-test_shell", "Test the DMFIELDSHELL implementation of DMField", "ex1.c", testShell, &testShell, NULL)); 239 PetscOptionsEnd(); 240 241 PetscCheck(dim <= 3, comm, PETSC_ERR_ARG_OUTOFRANGE, "This examples works for dim <= 3, not %" PetscInt_FMT, dim); 242 PetscCall(PetscStrncmp(type, DMPLEX, 256, &isplex)); 243 PetscCall(PetscStrncmp(type, DMDA, 256, &isda)); 244 245 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand)); 246 PetscCall(PetscRandomSetFromOptions(rand)); 247 if (isplex) { 248 PetscInt overlap = 0; 249 Vec fieldvec; 250 PetscInt cells[3] = {3, 3, 3}; 251 PetscBool simplex; 252 PetscFE fe; 253 254 PetscOptionsBegin(comm, "", "DMField DMPlex Options", "DM"); 255 PetscCall(PetscOptionsBoundedInt("-overlap", "DMPlex parallel overlap", "ex1.c", overlap, &overlap, NULL, 0)); 256 PetscOptionsEnd(); 257 if (0) { 258 PetscCall(DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, cells, NULL, NULL, NULL, PETSC_TRUE, &dm)); 259 } else { 260 PetscCall(DMCreate(comm, &dm)); 261 PetscCall(DMSetType(dm, DMPLEX)); 262 PetscCall(DMSetFromOptions(dm)); 263 CHKMEMQ; 264 } 265 PetscCall(DMGetDimension(dm, &dim)); 266 PetscCall(DMPlexIsSimplex(dm, &simplex)); 267 if (simplex) PetscCall(PetscDTStroudConicalQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad)); 268 else PetscCall(PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad)); 269 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 270 if (testShell) { 271 Vec ctxVec; 272 PetscInt i; 273 PetscScalar *array; 274 275 PetscCall(VecCreateSeq(PETSC_COMM_SELF, nc, &ctxVec)); 276 PetscCall(VecSetUp(ctxVec)); 277 PetscCall(VecGetArray(ctxVec, &array)); 278 for (i = 0; i < nc; i++) array[i] = i + 1.; 279 PetscCall(VecRestoreArray(ctxVec, &array)); 280 PetscCall(DMFieldCreateShell(dm, nc, DMFIELD_VERTEX, (void *)ctxVec, &field)); 281 PetscCall(DMFieldShellSetEvaluate(field, TestShellEvaluate)); 282 PetscCall(DMFieldShellSetDestroy(field, TestShellDestroy)); 283 } else { 284 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, nc, simplex, NULL, PETSC_DEFAULT, &fe)); 285 PetscCall(PetscFESetName(fe, "MyPetscFE")); 286 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 287 PetscCall(PetscFEDestroy(&fe)); 288 PetscCall(DMCreateDS(dm)); 289 PetscCall(DMCreateLocalVector(dm, &fieldvec)); 290 { 291 PetscErrorCode (*func[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 292 void *ctxs[1]; 293 294 func[0] = radiusSquared; 295 ctxs[0] = NULL; 296 297 PetscCall(DMProjectFunctionLocal(dm, 0.0, func, ctxs, INSERT_ALL_VALUES, fieldvec)); 298 } 299 PetscCall(DMFieldCreateDS(dm, 0, fieldvec, &field)); 300 PetscCall(VecDestroy(&fieldvec)); 301 } 302 } else if (isda) { 303 PetscInt i; 304 PetscScalar *cv; 305 306 switch (dim) { 307 case 1: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm)); break; 308 case 2: PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm)); break; 309 default: PetscCall(DMDACreate3d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, NULL, &dm)); break; 310 } 311 PetscCall(DMSetUp(dm)); 312 PetscCall(DMDAGetHeightStratum(dm, 0, &cStart, &cEnd)); 313 PetscCall(PetscMalloc1(nc * (1 << dim), &cv)); 314 for (i = 0; i < nc * (1 << dim); i++) { 315 PetscReal rv; 316 317 PetscCall(PetscRandomGetValueReal(rand, &rv)); 318 cv[i] = rv; 319 } 320 PetscCall(DMFieldCreateDA(dm, nc, cv, &field)); 321 PetscCall(PetscFree(cv)); 322 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad)); 323 } else SETERRQ(comm, PETSC_ERR_SUP, "This test does not run for DM type %s", type); 324 325 PetscCall(PetscObjectSetName((PetscObject)dm, "mesh")); 326 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 327 PetscCall(DMDestroy(&dm)); 328 PetscCall(PetscObjectSetName((PetscObject)field, "field")); 329 PetscCall(PetscObjectViewFromOptions((PetscObject)field, NULL, "-dmfield_view")); 330 if (numPoint) PetscCall(TestEvaluate(field, numPoint, rand)); 331 if (numFE) PetscCall(TestEvaluateFE(field, numFE, cStart, cEnd, quad, rand)); 332 if (numFV) PetscCall(TestEvaluateFV(field, numFV, cStart, cEnd, rand)); 333 PetscCall(DMFieldDestroy(&field)); 334 PetscCall(PetscQuadratureDestroy(&quad)); 335 PetscCall(PetscRandomDestroy(&rand)); 336 PetscCall(PetscFinalize()); 337 return 0; 338 } 339 340 /*TEST 341 342 test: 343 suffix: da 344 requires: !complex 345 args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view 346 347 test: 348 suffix: da_1 349 requires: !complex 350 args: -dm_type da -dim 1 -num_fe_tests 2 351 352 test: 353 suffix: da_2 354 requires: !complex 355 args: -dm_type da -dim 2 -num_fe_tests 2 356 357 test: 358 suffix: da_3 359 requires: !complex 360 args: -dm_type da -dim 3 -num_fe_tests 2 361 362 test: 363 suffix: ds 364 requires: !complex triangle 365 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view -petscspace_degree 2 -num_quad_points 1 366 367 test: 368 suffix: ds_simplex_0 369 requires: !complex triangle 370 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 0 371 372 test: 373 suffix: ds_simplex_1 374 requires: !complex triangle 375 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 1 376 377 test: 378 suffix: ds_simplex_2 379 requires: !complex triangle 380 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 2 381 382 test: 383 suffix: ds_tensor_2_0 384 requires: !complex 385 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_poly_tensor 1 -petscspace_degree 0 -dm_plex_simplex 0 386 387 test: 388 suffix: ds_tensor_2_1 389 requires: !complex 390 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_poly_tensor 1 -petscspace_degree 1 -dm_plex_simplex 0 391 392 test: 393 suffix: ds_tensor_2_2 394 requires: !complex 395 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_poly_tensor 1 -petscspace_degree 2 -dm_plex_simplex 0 396 397 test: 398 suffix: ds_tensor_3_0 399 requires: !complex 400 args: -dm_type plex -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -num_fe_tests 2 -petscspace_poly_tensor 1 -petscspace_degree 0 -dm_plex_simplex 0 401 402 test: 403 suffix: ds_tensor_3_1 404 requires: !complex 405 args: -dm_type plex -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -num_fe_tests 2 -petscspace_poly_tensor 1 -petscspace_degree 1 -dm_plex_simplex 0 406 407 test: 408 suffix: ds_tensor_3_2 409 requires: !complex 410 args: -dm_type plex -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -num_fe_tests 2 -petscspace_poly_tensor 1 -petscspace_degree 2 -dm_plex_simplex 0 411 412 test: 413 suffix: shell 414 requires: !complex triangle 415 args: -dm_coord_space 0 -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view -num_quad_points 1 -test_shell 416 417 TEST*/ 418