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