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(0); 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(VecCreateMPI(PetscObjectComm((PetscObject)field),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(0); 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(0); 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(0); 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++) { 154 u[i] = (i + 1) * r2; 155 } 156 PetscFunctionReturn(0); 157 } 158 159 static PetscErrorCode TestShellEvaluate(DMField field, Vec points, PetscDataType type, void *B, void *D, void *H) 160 { 161 Vec ctxVec = NULL; 162 const PetscScalar *mult; 163 PetscInt dim; 164 const PetscScalar *x; 165 PetscInt Nc, n, i, j, k, l; 166 167 PetscFunctionBegin; 168 PetscCall(DMFieldGetNumComponents(field, &Nc)); 169 PetscCall(DMFieldShellGetContext(field, &ctxVec)); 170 PetscCall(VecGetBlockSize(points, &dim)); 171 PetscCall(VecGetLocalSize(points, &n)); 172 n /= Nc; 173 PetscCall(VecGetArrayRead(ctxVec, &mult)); 174 PetscCall(VecGetArrayRead(points, &x)); 175 for (i = 0; i < n; i++) { 176 PetscReal r2 = 0.; 177 178 for (j = 0; j < dim; j++) {r2 += PetscSqr(PetscRealPart(x[i * dim + j]));} 179 for (j = 0; j < Nc; j++) { 180 PetscReal m = PetscRealPart(mult[j]); 181 if (B) { 182 if (type == PETSC_SCALAR) { 183 ((PetscScalar *)B)[i * Nc + j] = m * r2; 184 } else { 185 ((PetscReal *)B)[i * Nc + j] = m * r2; 186 } 187 } 188 if (D) { 189 if (type == PETSC_SCALAR) { 190 for (k = 0; k < dim; k++) ((PetscScalar *)D)[(i * Nc + j) * dim + k] = 2. * m * x[i * dim + k]; 191 } else { 192 for (k = 0; k < dim; k++) ((PetscReal *)D)[(i * Nc + j) * dim + k] = 2. * m * PetscRealPart(x[i * dim + k]); 193 } 194 } 195 if (H) { 196 if (type == PETSC_SCALAR) { 197 for (k = 0; k < dim; k++) for (l = 0; l < dim; l++) ((PetscScalar *)H)[((i * Nc + j) * dim + k) * dim + l] = (k == l) ? 2. * m : 0.; 198 } else { 199 for (k = 0; k < dim; k++) 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(0); 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(0); 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 PetscErrorCode ierr; 235 236 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 237 comm = PETSC_COMM_WORLD; 238 ierr = PetscOptionsBegin(comm, "", "DMField Tutorial Options", "DM");PetscCall(ierr); 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 ierr = PetscOptionsEnd();PetscCall(ierr); 248 249 PetscCheckFalse(dim > 3,comm,PETSC_ERR_ARG_OUTOFRANGE,"This examples works for dim <= 3, not %D",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 ierr = PetscOptionsBegin(comm, "", "DMField DMPlex Options", "DM");PetscCall(ierr); 263 PetscCall(PetscOptionsBoundedInt("-overlap","DMPlex parallel overlap","ex1.c",overlap,&overlap,NULL,0)); 264 ierr = PetscOptionsEnd();PetscCall(ierr); 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) { 276 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad)); 277 } else { 278 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad)); 279 } 280 PetscCall(DMPlexGetHeightStratum(dm,0,&cStart,&cEnd)); 281 if (testShell) { 282 Vec ctxVec; 283 PetscInt i; 284 PetscScalar *array; 285 286 PetscCall(VecCreateSeq(PETSC_COMM_SELF, nc, &ctxVec)); 287 PetscCall(VecSetUp(ctxVec)); 288 PetscCall(VecGetArray(ctxVec,&array)); 289 for (i = 0; i < nc; i++) array[i] = i + 1.; 290 PetscCall(VecRestoreArray(ctxVec,&array)); 291 PetscCall(DMFieldCreateShell(dm, nc, DMFIELD_VERTEX, (void *) ctxVec, &field)); 292 PetscCall(DMFieldShellSetEvaluate(field, TestShellEvaluate)); 293 PetscCall(DMFieldShellSetDestroy(field, TestShellDestroy)); 294 } else { 295 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF,dim,nc,simplex,NULL,PETSC_DEFAULT,&fe)); 296 PetscCall(PetscFESetName(fe,"MyPetscFE")); 297 PetscCall(DMSetField(dm,0,NULL,(PetscObject)fe)); 298 PetscCall(PetscFEDestroy(&fe)); 299 PetscCall(DMCreateDS(dm)); 300 PetscCall(DMCreateLocalVector(dm,&fieldvec)); 301 { 302 PetscErrorCode (*func[1]) (PetscInt,PetscReal,const PetscReal [],PetscInt, PetscScalar *,void *); 303 void *ctxs[1]; 304 305 func[0] = radiusSquared; 306 ctxs[0] = NULL; 307 308 PetscCall(DMProjectFunctionLocal(dm,0.0,func,ctxs,INSERT_ALL_VALUES,fieldvec)); 309 } 310 PetscCall(DMFieldCreateDS(dm,0,fieldvec,&field)); 311 PetscCall(VecDestroy(&fieldvec)); 312 } 313 } else if (isda) { 314 PetscInt i; 315 PetscScalar *cv; 316 317 switch (dim) { 318 case 1: 319 PetscCall(DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm)); 320 break; 321 case 2: 322 PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm)); 323 break; 324 default: 325 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)); 326 break; 327 } 328 PetscCall(DMSetUp(dm)); 329 PetscCall(DMDAGetHeightStratum(dm,0,&cStart,&cEnd)); 330 PetscCall(PetscMalloc1(nc * (1 << dim),&cv)); 331 for (i = 0; i < nc * (1 << dim); i++) { 332 PetscReal rv; 333 334 PetscCall(PetscRandomGetValueReal(rand,&rv)); 335 cv[i] = rv; 336 } 337 PetscCall(DMFieldCreateDA(dm,nc,cv,&field)); 338 PetscCall(PetscFree(cv)); 339 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad)); 340 } else SETERRQ(comm,PETSC_ERR_SUP,"This test does not run for DM type %s",type); 341 342 PetscCall(PetscObjectSetName((PetscObject)dm,"mesh")); 343 PetscCall(DMViewFromOptions(dm,NULL,"-dm_view")); 344 PetscCall(DMDestroy(&dm)); 345 PetscCall(PetscObjectSetName((PetscObject)field,"field")); 346 PetscCall(PetscObjectViewFromOptions((PetscObject)field,NULL,"-dmfield_view")); 347 if (numPoint) PetscCall(TestEvaluate(field,numPoint,rand)); 348 if (numFE) PetscCall(TestEvaluateFE(field,numFE,cStart,cEnd,quad,rand)); 349 if (numFV) PetscCall(TestEvaluateFV(field,numFV,cStart,cEnd,rand)); 350 PetscCall(DMFieldDestroy(&field)); 351 PetscCall(PetscQuadratureDestroy(&quad)); 352 PetscCall(PetscRandomDestroy(&rand)); 353 PetscCall(PetscFinalize()); 354 return 0; 355 } 356 357 /*TEST 358 359 test: 360 suffix: da 361 requires: !complex 362 args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view 363 364 test: 365 suffix: da_1 366 requires: !complex 367 args: -dm_type da -dim 1 -num_fe_tests 2 368 369 test: 370 suffix: da_2 371 requires: !complex 372 args: -dm_type da -dim 2 -num_fe_tests 2 373 374 test: 375 suffix: da_3 376 requires: !complex 377 args: -dm_type da -dim 3 -num_fe_tests 2 378 379 test: 380 suffix: ds 381 requires: !complex triangle 382 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 383 384 test: 385 suffix: ds_simplex_0 386 requires: !complex triangle 387 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 0 388 389 test: 390 suffix: ds_simplex_1 391 requires: !complex triangle 392 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 1 393 394 test: 395 suffix: ds_simplex_2 396 requires: !complex triangle 397 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 2 398 399 test: 400 suffix: ds_tensor_2_0 401 requires: !complex 402 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 403 404 test: 405 suffix: ds_tensor_2_1 406 requires: !complex 407 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 408 409 test: 410 suffix: ds_tensor_2_2 411 requires: !complex 412 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 413 414 test: 415 suffix: ds_tensor_3_0 416 requires: !complex 417 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 418 419 test: 420 suffix: ds_tensor_3_1 421 requires: !complex 422 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 423 424 test: 425 suffix: ds_tensor_3_2 426 requires: !complex 427 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 428 429 test: 430 suffix: shell 431 requires: !complex triangle 432 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 433 434 TEST*/ 435