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 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 236 comm = PETSC_COMM_WORLD; 237 PetscOptionsBegin(comm, "", "DMField Tutorial Options", "DM"); 238 PetscCall(PetscOptionsFList("-dm_type","DM implementation on which to define field","ex1.c",DMList,type,type,256,NULL)); 239 PetscCall(PetscOptionsRangeInt("-dim","DM intrinsic dimension", "ex1.c", dim, &dim, NULL,1,3)); 240 PetscCall(PetscOptionsBoundedInt("-num_components","Number of components in field", "ex1.c", nc, &nc, NULL,1)); 241 PetscCall(PetscOptionsBoundedInt("-num_quad_points","Number of quadrature points per dimension", "ex1.c", pointsPerEdge, &pointsPerEdge, NULL,1)); 242 PetscCall(PetscOptionsBoundedInt("-num_point_tests", "Number of test points for DMFieldEvaluate()", "ex1.c", numPoint, &numPoint, NULL,0)); 243 PetscCall(PetscOptionsBoundedInt("-num_fe_tests", "Number of test cells for DMFieldEvaluateFE()", "ex1.c", numFE, &numFE, NULL,0)); 244 PetscCall(PetscOptionsBoundedInt("-num_fv_tests", "Number of test cells for DMFieldEvaluateFV()", "ex1.c", numFV, &numFV, NULL,0)); 245 PetscCall(PetscOptionsBool("-test_shell", "Test the DMFIELDSHELL implementation of DMField", "ex1.c", testShell, &testShell, NULL)); 246 PetscOptionsEnd(); 247 248 PetscCheck(dim <= 3,comm,PETSC_ERR_ARG_OUTOFRANGE,"This examples works for dim <= 3, not %" PetscInt_FMT,dim); 249 PetscCall(PetscStrncmp(type,DMPLEX,256,&isplex)); 250 PetscCall(PetscStrncmp(type,DMDA,256,&isda)); 251 252 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rand)); 253 PetscCall(PetscRandomSetFromOptions(rand)); 254 if (isplex) { 255 PetscInt overlap = 0; 256 Vec fieldvec; 257 PetscInt cells[3] = {3,3,3}; 258 PetscBool simplex; 259 PetscFE fe; 260 261 PetscOptionsBegin(comm, "", "DMField DMPlex Options", "DM"); 262 PetscCall(PetscOptionsBoundedInt("-overlap","DMPlex parallel overlap","ex1.c",overlap,&overlap,NULL,0)); 263 PetscOptionsEnd(); 264 if (0) { 265 PetscCall(DMPlexCreateBoxMesh(comm,2,PETSC_TRUE,cells,NULL,NULL,NULL,PETSC_TRUE,&dm)); 266 } else { 267 PetscCall(DMCreate(comm, &dm)); 268 PetscCall(DMSetType(dm, DMPLEX)); 269 PetscCall(DMSetFromOptions(dm)); 270 CHKMEMQ; 271 } 272 PetscCall(DMGetDimension(dm, &dim)); 273 PetscCall(DMPlexIsSimplex(dm, &simplex)); 274 if (simplex) PetscCall(PetscDTStroudConicalQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad)); 275 else PetscCall(PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad)); 276 PetscCall(DMPlexGetHeightStratum(dm,0,&cStart,&cEnd)); 277 if (testShell) { 278 Vec ctxVec; 279 PetscInt i; 280 PetscScalar *array; 281 282 PetscCall(VecCreateSeq(PETSC_COMM_SELF, nc, &ctxVec)); 283 PetscCall(VecSetUp(ctxVec)); 284 PetscCall(VecGetArray(ctxVec,&array)); 285 for (i = 0; i < nc; i++) array[i] = i + 1.; 286 PetscCall(VecRestoreArray(ctxVec,&array)); 287 PetscCall(DMFieldCreateShell(dm, nc, DMFIELD_VERTEX, (void *) ctxVec, &field)); 288 PetscCall(DMFieldShellSetEvaluate(field, TestShellEvaluate)); 289 PetscCall(DMFieldShellSetDestroy(field, TestShellDestroy)); 290 } else { 291 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF,dim,nc,simplex,NULL,PETSC_DEFAULT,&fe)); 292 PetscCall(PetscFESetName(fe,"MyPetscFE")); 293 PetscCall(DMSetField(dm,0,NULL,(PetscObject)fe)); 294 PetscCall(PetscFEDestroy(&fe)); 295 PetscCall(DMCreateDS(dm)); 296 PetscCall(DMCreateLocalVector(dm,&fieldvec)); 297 { 298 PetscErrorCode (*func[1]) (PetscInt,PetscReal,const PetscReal [],PetscInt, PetscScalar *,void *); 299 void *ctxs[1]; 300 301 func[0] = radiusSquared; 302 ctxs[0] = NULL; 303 304 PetscCall(DMProjectFunctionLocal(dm,0.0,func,ctxs,INSERT_ALL_VALUES,fieldvec)); 305 } 306 PetscCall(DMFieldCreateDS(dm,0,fieldvec,&field)); 307 PetscCall(VecDestroy(&fieldvec)); 308 } 309 } else if (isda) { 310 PetscInt i; 311 PetscScalar *cv; 312 313 switch (dim) { 314 case 1: 315 PetscCall(DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm)); 316 break; 317 case 2: 318 PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 3, 3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dm)); 319 break; 320 default: 321 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)); 322 break; 323 } 324 PetscCall(DMSetUp(dm)); 325 PetscCall(DMDAGetHeightStratum(dm,0,&cStart,&cEnd)); 326 PetscCall(PetscMalloc1(nc * (1 << dim),&cv)); 327 for (i = 0; i < nc * (1 << dim); i++) { 328 PetscReal rv; 329 330 PetscCall(PetscRandomGetValueReal(rand,&rv)); 331 cv[i] = rv; 332 } 333 PetscCall(DMFieldCreateDA(dm,nc,cv,&field)); 334 PetscCall(PetscFree(cv)); 335 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad)); 336 } else SETERRQ(comm,PETSC_ERR_SUP,"This test does not run for DM type %s",type); 337 338 PetscCall(PetscObjectSetName((PetscObject)dm,"mesh")); 339 PetscCall(DMViewFromOptions(dm,NULL,"-dm_view")); 340 PetscCall(DMDestroy(&dm)); 341 PetscCall(PetscObjectSetName((PetscObject)field,"field")); 342 PetscCall(PetscObjectViewFromOptions((PetscObject)field,NULL,"-dmfield_view")); 343 if (numPoint) PetscCall(TestEvaluate(field,numPoint,rand)); 344 if (numFE) PetscCall(TestEvaluateFE(field,numFE,cStart,cEnd,quad,rand)); 345 if (numFV) PetscCall(TestEvaluateFV(field,numFV,cStart,cEnd,rand)); 346 PetscCall(DMFieldDestroy(&field)); 347 PetscCall(PetscQuadratureDestroy(&quad)); 348 PetscCall(PetscRandomDestroy(&rand)); 349 PetscCall(PetscFinalize()); 350 return 0; 351 } 352 353 /*TEST 354 355 test: 356 suffix: da 357 requires: !complex 358 args: -dm_type da -dim 2 -num_components 2 -num_point_tests 2 -num_fe_tests 2 -num_fv_tests 2 -dmfield_view 359 360 test: 361 suffix: da_1 362 requires: !complex 363 args: -dm_type da -dim 1 -num_fe_tests 2 364 365 test: 366 suffix: da_2 367 requires: !complex 368 args: -dm_type da -dim 2 -num_fe_tests 2 369 370 test: 371 suffix: da_3 372 requires: !complex 373 args: -dm_type da -dim 3 -num_fe_tests 2 374 375 test: 376 suffix: ds 377 requires: !complex triangle 378 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 379 380 test: 381 suffix: ds_simplex_0 382 requires: !complex triangle 383 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 0 384 385 test: 386 suffix: ds_simplex_1 387 requires: !complex triangle 388 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 1 389 390 test: 391 suffix: ds_simplex_2 392 requires: !complex triangle 393 args: -dm_type plex -dm_plex_dim 2 -dm_plex_box_faces 3,3 -num_fe_tests 2 -petscspace_degree 2 394 395 test: 396 suffix: ds_tensor_2_0 397 requires: !complex 398 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 399 400 test: 401 suffix: ds_tensor_2_1 402 requires: !complex 403 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 404 405 test: 406 suffix: ds_tensor_2_2 407 requires: !complex 408 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 409 410 test: 411 suffix: ds_tensor_3_0 412 requires: !complex 413 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 414 415 test: 416 suffix: ds_tensor_3_1 417 requires: !complex 418 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 419 420 test: 421 suffix: ds_tensor_3_2 422 requires: !complex 423 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 424 425 test: 426 suffix: shell 427 requires: !complex triangle 428 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 429 430 TEST*/ 431