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