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