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 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