xref: /petsc/src/dm/field/tutorials/ex1.c (revision c4762a1b19cd2af06abeed90e8f9d34fb975dd94)
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