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