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 CHKERRQ(PetscViewerASCIIPrintf(viewer,"B:\n")); 11 CHKERRQ(PetscScalarView(N,B,viewer)); 12 CHKERRQ(PetscViewerASCIIPrintf(viewer,"D:\n")); 13 CHKERRQ(PetscScalarView(N*dim,D,viewer)); 14 CHKERRQ(PetscViewerASCIIPrintf(viewer,"H:\n")); 15 CHKERRQ(PetscScalarView(N*dim*dim,H,viewer)); 16 17 CHKERRQ(PetscViewerASCIIPrintf(viewer,"rB:\n")); 18 CHKERRQ(PetscRealView(N,rB,viewer)); 19 CHKERRQ(PetscViewerASCIIPrintf(viewer,"rD:\n")); 20 CHKERRQ(PetscRealView(N*dim,rD,viewer)); 21 CHKERRQ(PetscViewerASCIIPrintf(viewer,"rH:\n")); 22 CHKERRQ(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 CHKERRQ(DMFieldGetNumComponents(field,&nc)); 40 CHKERRQ(DMFieldGetDM(field,&dm)); 41 CHKERRQ(DMGetDimension(dm,&dim)); 42 CHKERRQ(VecCreateMPI(PetscObjectComm((PetscObject)field),n * dim,PETSC_DETERMINE,&points)); 43 CHKERRQ(VecSetBlockSize(points,dim)); 44 CHKERRQ(VecGetArray(points,&array)); 45 for (i = 0; i < n * dim; i++) CHKERRQ(PetscRandomGetValue(rand,&array[i])); 46 CHKERRQ(VecRestoreArray(points,&array)); 47 CHKERRQ(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 CHKERRQ(DMFieldEvaluate(field,points,PETSC_SCALAR,B,D,H)); 49 CHKERRQ(DMFieldEvaluate(field,points,PETSC_REAL,rB,rD,rH)); 50 viewer = PETSC_VIEWER_STDOUT_(comm); 51 52 CHKERRQ(PetscObjectSetName((PetscObject)points,"Test Points")); 53 CHKERRQ(VecView(points,viewer)); 54 CHKERRQ(ViewResults(viewer,n*nc,dim,B,D,H,rB,rD,rH)); 55 56 CHKERRQ(PetscFree6(B,rB,D,rD,H,rH)); 57 CHKERRQ(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 CHKERRQ(DMFieldGetNumComponents(field,&nc)); 76 CHKERRQ(DMFieldGetDM(field,&dm)); 77 CHKERRQ(DMGetDimension(dm,&dim)); 78 CHKERRQ(PetscRandomSetInterval(rand,(PetscScalar) cStart, (PetscScalar) cEnd)); 79 CHKERRQ(PetscMalloc1(n,&cells)); 80 for (i = 0; i < n; i++) { 81 PetscReal rc; 82 83 CHKERRQ(PetscRandomGetValueReal(rand,&rc)); 84 cells[i] = PetscFloorReal(rc); 85 } 86 CHKERRQ(ISCreateGeneral(comm,n,cells,PETSC_OWN_POINTER,&cellIS)); 87 CHKERRQ(PetscObjectSetName((PetscObject)cellIS,"FE Test Cells")); 88 CHKERRQ(PetscQuadratureGetData(quad,NULL,NULL,&nq,NULL,NULL)); 89 N = n * nq * nc; 90 CHKERRQ(PetscMalloc6(N,&B,N,&rB,N*dim,&D,N*dim,&rD,N*dim*dim,&H,N*dim*dim,&rH)); 91 CHKERRQ(DMFieldEvaluateFE(field,cellIS,quad,PETSC_SCALAR,B,D,H)); 92 CHKERRQ(DMFieldEvaluateFE(field,cellIS,quad,PETSC_REAL,rB,rD,rH)); 93 viewer = PETSC_VIEWER_STDOUT_(comm); 94 95 CHKERRQ(PetscObjectSetName((PetscObject)quad,"Test quadrature")); 96 CHKERRQ(PetscQuadratureView(quad,viewer)); 97 CHKERRQ(ISView(cellIS,viewer)); 98 CHKERRQ(ViewResults(viewer,N,dim,B,D,H,rB,rD,rH)); 99 100 CHKERRQ(PetscFree6(B,rB,D,rD,H,rH)); 101 CHKERRQ(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 CHKERRQ(DMFieldGetNumComponents(field,&nc)); 120 CHKERRQ(DMFieldGetDM(field,&dm)); 121 CHKERRQ(DMGetDimension(dm,&dim)); 122 CHKERRQ(PetscRandomSetInterval(rand,(PetscScalar) cStart, (PetscScalar) cEnd)); 123 CHKERRQ(PetscMalloc1(n,&cells)); 124 for (i = 0; i < n; i++) { 125 PetscReal rc; 126 127 CHKERRQ(PetscRandomGetValueReal(rand,&rc)); 128 cells[i] = PetscFloorReal(rc); 129 } 130 CHKERRQ(ISCreateGeneral(comm,n,cells,PETSC_OWN_POINTER,&cellIS)); 131 CHKERRQ(PetscObjectSetName((PetscObject)cellIS,"FV Test Cells")); 132 N = n * nc; 133 CHKERRQ(PetscMalloc6(N,&B,N,&rB,N*dim,&D,N*dim,&rD,N*dim*dim,&H,N*dim*dim,&rH)); 134 CHKERRQ(DMFieldEvaluateFV(field,cellIS,PETSC_SCALAR,B,D,H)); 135 CHKERRQ(DMFieldEvaluateFV(field,cellIS,PETSC_REAL,rB,rD,rH)); 136 viewer = PETSC_VIEWER_STDOUT_(comm); 137 138 CHKERRQ(ISView(cellIS,viewer)); 139 CHKERRQ(ViewResults(viewer,N,dim,B,D,H,rB,rD,rH)); 140 141 CHKERRQ(PetscFree6(B,rB,D,rD,H,rH)); 142 CHKERRQ(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 CHKERRQ(DMFieldGetNumComponents(field, &Nc)); 169 CHKERRQ(DMFieldShellGetContext(field, &ctxVec)); 170 CHKERRQ(VecGetBlockSize(points, &dim)); 171 CHKERRQ(VecGetLocalSize(points, &n)); 172 n /= Nc; 173 CHKERRQ(VecGetArrayRead(ctxVec, &mult)); 174 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(points, &x)); 205 CHKERRQ(VecRestoreArrayRead(ctxVec, &mult)); 206 PetscFunctionReturn(0); 207 } 208 209 static PetscErrorCode TestShellDestroy(DMField field) 210 { 211 Vec ctxVec = NULL; 212 213 PetscFunctionBegin; 214 CHKERRQ(DMFieldShellGetContext(field, &ctxVec)); 215 CHKERRQ(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 CHKERRQ(PetscInitialize(&argc, &argv, NULL, help)); 237 comm = PETSC_COMM_WORLD; 238 ierr = PetscOptionsBegin(comm, "", "DMField Tutorial Options", "DM");CHKERRQ(ierr); 239 CHKERRQ(PetscOptionsFList("-dm_type","DM implementation on which to define field","ex1.c",DMList,type,type,256,NULL)); 240 CHKERRQ(PetscOptionsRangeInt("-dim","DM intrinsic dimension", "ex1.c", dim, &dim, NULL,1,3)); 241 CHKERRQ(PetscOptionsBoundedInt("-num_components","Number of components in field", "ex1.c", nc, &nc, NULL,1)); 242 CHKERRQ(PetscOptionsBoundedInt("-num_quad_points","Number of quadrature points per dimension", "ex1.c", pointsPerEdge, &pointsPerEdge, NULL,1)); 243 CHKERRQ(PetscOptionsBoundedInt("-num_point_tests", "Number of test points for DMFieldEvaluate()", "ex1.c", numPoint, &numPoint, NULL,0)); 244 CHKERRQ(PetscOptionsBoundedInt("-num_fe_tests", "Number of test cells for DMFieldEvaluateFE()", "ex1.c", numFE, &numFE, NULL,0)); 245 CHKERRQ(PetscOptionsBoundedInt("-num_fv_tests", "Number of test cells for DMFieldEvaluateFV()", "ex1.c", numFV, &numFV, NULL,0)); 246 CHKERRQ(PetscOptionsBool("-test_shell", "Test the DMFIELDSHELL implementation of DMField", "ex1.c", testShell, &testShell, NULL)); 247 ierr = PetscOptionsEnd();CHKERRQ(ierr); 248 249 PetscCheckFalse(dim > 3,comm,PETSC_ERR_ARG_OUTOFRANGE,"This examples works for dim <= 3, not %D",dim); 250 CHKERRQ(PetscStrncmp(type,DMPLEX,256,&isplex)); 251 CHKERRQ(PetscStrncmp(type,DMDA,256,&isda)); 252 253 CHKERRQ(PetscRandomCreate(PETSC_COMM_SELF,&rand)); 254 CHKERRQ(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");CHKERRQ(ierr); 263 CHKERRQ(PetscOptionsBoundedInt("-overlap","DMPlex parallel overlap","ex1.c",overlap,&overlap,NULL,0)); 264 ierr = PetscOptionsEnd();CHKERRQ(ierr); 265 if (0) { 266 CHKERRQ(DMPlexCreateBoxMesh(comm,2,PETSC_TRUE,cells,NULL,NULL,NULL,PETSC_TRUE,&dm)); 267 } else { 268 CHKERRQ(DMCreate(comm, &dm)); 269 CHKERRQ(DMSetType(dm, DMPLEX)); 270 CHKERRQ(DMSetFromOptions(dm)); 271 CHKMEMQ; 272 } 273 CHKERRQ(DMGetDimension(dm, &dim)); 274 CHKERRQ(DMPlexIsSimplex(dm, &simplex)); 275 if (simplex) { 276 CHKERRQ(PetscDTStroudConicalQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad)); 277 } else { 278 CHKERRQ(PetscDTGaussTensorQuadrature(dim, 1, pointsPerEdge, -1.0, 1.0, &quad)); 279 } 280 CHKERRQ(DMPlexGetHeightStratum(dm,0,&cStart,&cEnd)); 281 if (testShell) { 282 Vec ctxVec; 283 PetscInt i; 284 PetscScalar *array; 285 286 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF, nc, &ctxVec)); 287 CHKERRQ(VecSetUp(ctxVec)); 288 CHKERRQ(VecGetArray(ctxVec,&array)); 289 for (i = 0; i < nc; i++) array[i] = i + 1.; 290 CHKERRQ(VecRestoreArray(ctxVec,&array)); 291 CHKERRQ(DMFieldCreateShell(dm, nc, DMFIELD_VERTEX, (void *) ctxVec, &field)); 292 CHKERRQ(DMFieldShellSetEvaluate(field, TestShellEvaluate)); 293 CHKERRQ(DMFieldShellSetDestroy(field, TestShellDestroy)); 294 } else { 295 CHKERRQ(PetscFECreateDefault(PETSC_COMM_SELF,dim,nc,simplex,NULL,PETSC_DEFAULT,&fe)); 296 CHKERRQ(PetscFESetName(fe,"MyPetscFE")); 297 CHKERRQ(DMSetField(dm,0,NULL,(PetscObject)fe)); 298 CHKERRQ(PetscFEDestroy(&fe)); 299 CHKERRQ(DMCreateDS(dm)); 300 CHKERRQ(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 CHKERRQ(DMProjectFunctionLocal(dm,0.0,func,ctxs,INSERT_ALL_VALUES,fieldvec)); 309 } 310 CHKERRQ(DMFieldCreateDS(dm,0,fieldvec,&field)); 311 CHKERRQ(VecDestroy(&fieldvec)); 312 } 313 } else if (isda) { 314 PetscInt i; 315 PetscScalar *cv; 316 317 switch (dim) { 318 case 1: 319 CHKERRQ(DMDACreate1d(comm, DM_BOUNDARY_NONE, 3, 1, 1, NULL, &dm)); 320 break; 321 case 2: 322 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMSetUp(dm)); 329 CHKERRQ(DMDAGetHeightStratum(dm,0,&cStart,&cEnd)); 330 CHKERRQ(PetscMalloc1(nc * (1 << dim),&cv)); 331 for (i = 0; i < nc * (1 << dim); i++) { 332 PetscReal rv; 333 334 CHKERRQ(PetscRandomGetValueReal(rand,&rv)); 335 cv[i] = rv; 336 } 337 CHKERRQ(DMFieldCreateDA(dm,nc,cv,&field)); 338 CHKERRQ(PetscFree(cv)); 339 CHKERRQ(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 CHKERRQ(PetscObjectSetName((PetscObject)dm,"mesh")); 343 CHKERRQ(DMViewFromOptions(dm,NULL,"-dm_view")); 344 CHKERRQ(DMDestroy(&dm)); 345 CHKERRQ(PetscObjectSetName((PetscObject)field,"field")); 346 CHKERRQ(PetscObjectViewFromOptions((PetscObject)field,NULL,"-dmfield_view")); 347 if (numPoint) CHKERRQ(TestEvaluate(field,numPoint,rand)); 348 if (numFE) CHKERRQ(TestEvaluateFE(field,numFE,cStart,cEnd,quad,rand)); 349 if (numFV) CHKERRQ(TestEvaluateFV(field,numFV,cStart,cEnd,rand)); 350 CHKERRQ(DMFieldDestroy(&field)); 351 CHKERRQ(PetscQuadratureDestroy(&quad)); 352 CHKERRQ(PetscRandomDestroy(&rand)); 353 CHKERRQ(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