xref: /petsc/src/dm/impls/forest/tests/ex2.c (revision 2f613bf53f46f9356e00a2ca2bd69453be72fc31)
1c4762a1bSJed Brown static char help[] = "Create a mesh, refine and coarsen simultaneously, and transfer a field\n\n";
2c4762a1bSJed Brown 
3c4762a1bSJed Brown #include <petscds.h>
4c4762a1bSJed Brown #include <petscdmplex.h>
5c4762a1bSJed Brown #include <petscdmforest.h>
6c4762a1bSJed Brown #include <petscoptions.h>
7c4762a1bSJed Brown 
8c4762a1bSJed Brown static PetscErrorCode AddIdentityLabel(DM dm)
9c4762a1bSJed Brown {
10c4762a1bSJed Brown   PetscInt       pStart,pEnd,p;
11c4762a1bSJed Brown   PetscErrorCode ierr;
12c4762a1bSJed Brown 
13c4762a1bSJed Brown   PetscFunctionBegin;
14c4762a1bSJed Brown   ierr = DMCreateLabel(dm, "identity");CHKERRQ(ierr);
15c4762a1bSJed Brown   ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
16c4762a1bSJed Brown   for (p = pStart; p < pEnd; p++) {ierr = DMSetLabelValue(dm, "identity", p, p);CHKERRQ(ierr);}
17c4762a1bSJed Brown   PetscFunctionReturn(0);
18c4762a1bSJed Brown }
19c4762a1bSJed Brown 
20c4762a1bSJed Brown static PetscErrorCode CreateAdaptivityLabel(DM forest,DMLabel *adaptLabel)
21c4762a1bSJed Brown {
22c4762a1bSJed Brown   DMLabel        identLabel;
23c4762a1bSJed Brown   PetscInt       cStart, cEnd, c;
24c4762a1bSJed Brown   PetscErrorCode ierr;
25c4762a1bSJed Brown 
26c4762a1bSJed Brown   PetscFunctionBegin;
27c4762a1bSJed Brown   ierr = DMLabelCreate(PETSC_COMM_SELF,"adapt",adaptLabel);CHKERRQ(ierr);
28c4762a1bSJed Brown   ierr = DMLabelSetDefaultValue(*adaptLabel,DM_ADAPT_COARSEN);CHKERRQ(ierr);
29c4762a1bSJed Brown   ierr = DMGetLabel(forest,"identity",&identLabel);CHKERRQ(ierr);
30c4762a1bSJed Brown   ierr = DMForestGetCellChart(forest,&cStart,&cEnd);CHKERRQ(ierr);
31c4762a1bSJed Brown   for (c = cStart; c < cEnd; c++) {
32c4762a1bSJed Brown     PetscInt basePoint;
33c4762a1bSJed Brown 
34c4762a1bSJed Brown     ierr = DMLabelGetValue(identLabel,c,&basePoint);CHKERRQ(ierr);
35c4762a1bSJed Brown     if (!basePoint) {ierr = DMLabelSetValue(*adaptLabel,c,DM_ADAPT_REFINE);CHKERRQ(ierr);}
36c4762a1bSJed Brown   }
37c4762a1bSJed Brown   PetscFunctionReturn(0);
38c4762a1bSJed Brown }
39c4762a1bSJed Brown 
40c4762a1bSJed Brown static PetscErrorCode LinearFunction(PetscInt dim,PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
41c4762a1bSJed Brown {
42c4762a1bSJed Brown   PetscFunctionBeginUser;
43c4762a1bSJed Brown   u[0] = (x[0] * 2.0 + 1.) + (x[1] * 20.0 + 10.) + ((dim == 3) ? (x[2] * 200.0 + 100.) : 0.);
44c4762a1bSJed Brown   PetscFunctionReturn(0);
45c4762a1bSJed Brown }
46c4762a1bSJed Brown 
47c4762a1bSJed Brown static PetscErrorCode MultiaffineFunction(PetscInt dim,PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
48c4762a1bSJed Brown {
49c4762a1bSJed Brown   PetscFunctionBeginUser;
50c4762a1bSJed Brown   u[0] = (x[0] * 1.0 + 2.0) * (x[1] * 3.0 - 4.0) * ((dim == 3) ? (x[2] * 5.0 + 6.0) : 1.);
51c4762a1bSJed Brown   PetscFunctionReturn(0);
52c4762a1bSJed Brown }
53c4762a1bSJed Brown 
54c4762a1bSJed Brown static PetscErrorCode CoordsFunction(PetscInt dim,PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
55c4762a1bSJed Brown {
56c4762a1bSJed Brown   PetscInt f;
57c4762a1bSJed Brown 
58c4762a1bSJed Brown   PetscFunctionBeginUser;
59c4762a1bSJed Brown   for (f=0;f<Nf;f++) u[f] = x[f];
60c4762a1bSJed Brown   PetscFunctionReturn(0);
61c4762a1bSJed Brown }
62c4762a1bSJed Brown 
63c4762a1bSJed Brown typedef struct _bc_func_ctx
64c4762a1bSJed Brown {
65c4762a1bSJed Brown   PetscErrorCode (*func) (PetscInt,PetscReal,const PetscReal [], PetscInt, PetscScalar [], void *);
66c4762a1bSJed Brown   PetscInt dim;
67c4762a1bSJed Brown   PetscInt Nf;
68c4762a1bSJed Brown   void *ctx;
69c4762a1bSJed Brown }
70c4762a1bSJed Brown bc_func_ctx;
71c4762a1bSJed Brown 
72c4762a1bSJed Brown static PetscErrorCode bc_func_fv (PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
73c4762a1bSJed Brown {
74c4762a1bSJed Brown   bc_func_ctx    *bcCtx;
75c4762a1bSJed Brown   PetscErrorCode ierr;
76c4762a1bSJed Brown 
77c4762a1bSJed Brown   PetscFunctionBegin;
78c4762a1bSJed Brown   bcCtx = (bc_func_ctx *) ctx;
79c4762a1bSJed Brown   ierr = (bcCtx->func)(bcCtx->dim,time,c,bcCtx->Nf,xG,bcCtx->ctx);CHKERRQ(ierr);
80c4762a1bSJed Brown   PetscFunctionReturn(0);
81c4762a1bSJed Brown }
82c4762a1bSJed Brown 
83c4762a1bSJed Brown static PetscErrorCode IdentifyBadPoints (DM dm, Vec vec, PetscReal tol)
84c4762a1bSJed Brown {
85c4762a1bSJed Brown   DM             dmplex;
86c4762a1bSJed Brown   PetscInt       p, pStart, pEnd, maxDof;
87c4762a1bSJed Brown   Vec            vecLocal;
88c4762a1bSJed Brown   DMLabel        depthLabel;
89c4762a1bSJed Brown   PetscSection   section;
90c4762a1bSJed Brown   PetscErrorCode ierr;
91c4762a1bSJed Brown 
92c4762a1bSJed Brown   PetscFunctionBegin;
93c4762a1bSJed Brown   ierr = DMCreateLocalVector(dm, &vecLocal);CHKERRQ(ierr);
94c4762a1bSJed Brown   ierr = DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal);CHKERRQ(ierr);
95c4762a1bSJed Brown   ierr = DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal);CHKERRQ(ierr);
96c4762a1bSJed Brown   ierr = DMConvert(dm ,DMPLEX, &dmplex);CHKERRQ(ierr);
97c4762a1bSJed Brown   ierr = DMPlexGetChart(dmplex, &pStart, &pEnd);CHKERRQ(ierr);
98c4762a1bSJed Brown   ierr = DMPlexGetDepthLabel(dmplex, &depthLabel);CHKERRQ(ierr);
99c4762a1bSJed Brown   ierr = DMGetLocalSection(dmplex, &section);CHKERRQ(ierr);
100c4762a1bSJed Brown   ierr = PetscSectionGetMaxDof(section, &maxDof);CHKERRQ(ierr);
101c4762a1bSJed Brown   for (p = pStart; p < pEnd; p++) {
102c4762a1bSJed Brown     PetscInt     s, c, cSize, parent, childID, numChildren;
103c4762a1bSJed Brown     PetscInt     cl, closureSize, *closure = NULL;
104c4762a1bSJed Brown     PetscScalar *values = NULL;
105c4762a1bSJed Brown     PetscBool    bad = PETSC_FALSE;
106c4762a1bSJed Brown 
107c4762a1bSJed Brown     ierr = VecGetValuesSection(vecLocal, section, p, &values);CHKERRQ(ierr);
108c4762a1bSJed Brown     ierr = PetscSectionGetDof(section, p, &cSize);CHKERRQ(ierr);
109c4762a1bSJed Brown     for (c = 0; c < cSize; c++) {
110*2f613bf5SBarry Smith       PetscReal absDiff = PetscAbsScalar(values[c]);
111c4762a1bSJed Brown       if (absDiff > tol) {bad = PETSC_TRUE; break;}
112c4762a1bSJed Brown     }
113c4762a1bSJed Brown     if (!bad) continue;
114c4762a1bSJed Brown     ierr = PetscPrintf(PETSC_COMM_SELF, "Bad point %D\n", p);CHKERRQ(ierr);
115c4762a1bSJed Brown     ierr = DMLabelGetValue(depthLabel, p, &s);CHKERRQ(ierr);
116c4762a1bSJed Brown     ierr = PetscPrintf(PETSC_COMM_SELF, "  Depth %D\n", s);CHKERRQ(ierr);
117c4762a1bSJed Brown     ierr = DMPlexGetTransitiveClosure(dmplex, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
118c4762a1bSJed Brown     for (cl = 0; cl < closureSize; cl++) {
119c4762a1bSJed Brown       PetscInt cp = closure[2 * cl];
120c4762a1bSJed Brown       ierr = DMPlexGetTreeParent(dmplex, cp, &parent, &childID);CHKERRQ(ierr);
121c4762a1bSJed Brown       if (parent != cp) {
122c4762a1bSJed Brown         ierr = PetscPrintf(PETSC_COMM_SELF, "  Closure point %D (%D) child of %D (ID %D)\n", cl, cp, parent, childID);CHKERRQ(ierr);
123c4762a1bSJed Brown       }
124c4762a1bSJed Brown       ierr = DMPlexGetTreeChildren(dmplex, cp, &numChildren, NULL);CHKERRQ(ierr);
125c4762a1bSJed Brown       if (numChildren) {
126c4762a1bSJed Brown         ierr = PetscPrintf(PETSC_COMM_SELF, "  Closure point %D (%D) is parent\n", cl, cp);CHKERRQ(ierr);
127c4762a1bSJed Brown       }
128c4762a1bSJed Brown     }
129c4762a1bSJed Brown     ierr = DMPlexRestoreTransitiveClosure(dmplex, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
130c4762a1bSJed Brown     for (c = 0; c < cSize; c++) {
131*2f613bf5SBarry Smith       PetscReal absDiff = PetscAbsScalar(values[c]);
132c4762a1bSJed Brown       if (absDiff > tol) {
133c4762a1bSJed Brown         ierr = PetscPrintf(PETSC_COMM_SELF, "  Bad dof %D\n", c);CHKERRQ(ierr);
134c4762a1bSJed Brown       }
135c4762a1bSJed Brown     }
136c4762a1bSJed Brown   }
137c4762a1bSJed Brown   ierr = DMDestroy(&dmplex);CHKERRQ(ierr);
138c4762a1bSJed Brown   ierr = VecDestroy(&vecLocal);CHKERRQ(ierr);
139c4762a1bSJed Brown   PetscFunctionReturn(0);
140c4762a1bSJed Brown }
141c4762a1bSJed Brown 
142c4762a1bSJed Brown int main(int argc, char **argv)
143c4762a1bSJed Brown {
144c4762a1bSJed Brown   MPI_Comm       comm;
145c4762a1bSJed Brown   DM             base, preForest, postForest;
14630602db0SMatthew G. Knepley   PetscInt       dim, Nf = 1;
147c4762a1bSJed Brown   PetscInt       step, adaptSteps = 1;
148c4762a1bSJed Brown   PetscInt       preCount, postCount;
149c4762a1bSJed Brown   Vec            preVec, postVecTransfer, postVecExact;
150c4762a1bSJed Brown   PetscErrorCode (*funcs[1]) (PetscInt,PetscReal,const PetscReal [],PetscInt,PetscScalar [], void *) = {MultiaffineFunction};
151c4762a1bSJed Brown   void           *ctxs[1] = {NULL};
152c4762a1bSJed Brown   PetscReal      diff, tol = PETSC_SMALL;
153c4762a1bSJed Brown   PetscBool      linear = PETSC_FALSE;
154c4762a1bSJed Brown   PetscBool      coords = PETSC_FALSE;
155c4762a1bSJed Brown   PetscBool      useFV = PETSC_FALSE;
156c4762a1bSJed Brown   PetscBool      conv = PETSC_FALSE;
157c4762a1bSJed Brown   PetscBool      transfer_from_base[2] = {PETSC_TRUE,PETSC_FALSE};
158c4762a1bSJed Brown   PetscBool      use_bcs = PETSC_TRUE;
159c4762a1bSJed Brown   bc_func_ctx    bcCtx;
160c4762a1bSJed Brown   DMLabel        adaptLabel;
161c4762a1bSJed Brown   PetscErrorCode ierr;
162c4762a1bSJed Brown 
163c4762a1bSJed Brown   ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
164c4762a1bSJed Brown   comm = PETSC_COMM_WORLD;
165c4762a1bSJed Brown   ierr = PetscOptionsBegin(comm, "", "DMForestTransferVec() Test Options", "DMFOREST");CHKERRQ(ierr);
166c4762a1bSJed Brown   ierr = PetscOptionsBool("-linear","Transfer a simple linear function", "ex2.c", linear, &linear, NULL);CHKERRQ(ierr);
167c4762a1bSJed Brown   ierr = PetscOptionsBool("-coords","Transfer a simple coordinate function", "ex2.c", coords, &coords, NULL);CHKERRQ(ierr);
168c4762a1bSJed Brown   ierr = PetscOptionsBool("-use_fv","Use a finite volume approximation", "ex2.c", useFV, &useFV, NULL);CHKERRQ(ierr);
169c4762a1bSJed Brown   ierr = PetscOptionsBool("-test_convert","Test conversion to DMPLEX",NULL,conv,&conv,NULL);CHKERRQ(ierr);
170c4762a1bSJed Brown   ierr = PetscOptionsBool("-transfer_from_base","Transfer a vector from base DM to DMForest", "ex2.c", transfer_from_base[0], &transfer_from_base[0], NULL);CHKERRQ(ierr);
171c4762a1bSJed Brown   transfer_from_base[1] = transfer_from_base[0];
172c4762a1bSJed Brown   ierr = PetscOptionsBool("-transfer_from_base_steps","Transfer a vector from base DM to the latest DMForest after the adaptivity steps", "ex2.c", transfer_from_base[1], &transfer_from_base[1], NULL);CHKERRQ(ierr);
173c4762a1bSJed Brown   ierr = PetscOptionsBool("-use_bcs","Use dirichlet boundary conditions", "ex2.c", use_bcs, &use_bcs, NULL);CHKERRQ(ierr);
174c4762a1bSJed Brown   ierr = PetscOptionsBoundedInt("-adapt_steps","Number of adaptivity steps", "ex2.c", adaptSteps, &adaptSteps, NULL,0);CHKERRQ(ierr);
175c4762a1bSJed Brown   ierr = PetscOptionsEnd();CHKERRQ(ierr);
176c4762a1bSJed Brown 
177c4762a1bSJed Brown   tol = PetscMax(1.e-10,tol); /* XXX fix for quadruple precision -> why do I need to do this? */
178c4762a1bSJed Brown 
17930602db0SMatthew G. Knepley   /* the base mesh */
18030602db0SMatthew G. Knepley   ierr = DMCreate(comm, &base);CHKERRQ(ierr);
18130602db0SMatthew G. Knepley   ierr = DMSetType(base, DMPLEX);CHKERRQ(ierr);
18230602db0SMatthew G. Knepley   ierr = DMSetFromOptions(base);CHKERRQ(ierr);
18330602db0SMatthew G. Knepley 
18430602db0SMatthew G. Knepley   ierr = AddIdentityLabel(base);CHKERRQ(ierr);
18530602db0SMatthew G. Knepley   ierr = DMGetDimension(base, &dim);CHKERRQ(ierr);
186c4762a1bSJed Brown 
187c4762a1bSJed Brown   if (linear) {
188c4762a1bSJed Brown     funcs[0] = LinearFunction;
189c4762a1bSJed Brown   }
190c4762a1bSJed Brown   if (coords) {
191c4762a1bSJed Brown     funcs[0] = CoordsFunction;
192c4762a1bSJed Brown     Nf = dim;
193c4762a1bSJed Brown   }
194c4762a1bSJed Brown 
195c4762a1bSJed Brown   bcCtx.func = funcs[0];
196c4762a1bSJed Brown   bcCtx.dim  = dim;
197c4762a1bSJed Brown   bcCtx.Nf   = Nf;
198c4762a1bSJed Brown   bcCtx.ctx  = NULL;
199c4762a1bSJed Brown 
200c4762a1bSJed Brown   if (useFV) {
201c4762a1bSJed Brown     PetscFV      fv;
202c4762a1bSJed Brown     PetscLimiter limiter;
203c4762a1bSJed Brown     DM           baseFV;
204c4762a1bSJed Brown 
205c4762a1bSJed Brown     ierr = DMPlexConstructGhostCells(base,NULL,NULL,&baseFV);CHKERRQ(ierr);
20654fcfd0cSMatthew G. Knepley     ierr = DMViewFromOptions(baseFV, NULL, "-fv_dm_view");CHKERRQ(ierr);
207c4762a1bSJed Brown     ierr = DMDestroy(&base);CHKERRQ(ierr);
208c4762a1bSJed Brown     base = baseFV;
209c4762a1bSJed Brown     ierr = PetscFVCreate(comm, &fv);CHKERRQ(ierr);
210c4762a1bSJed Brown     ierr = PetscFVSetSpatialDimension(fv,dim);CHKERRQ(ierr);
211c4762a1bSJed Brown     ierr = PetscFVSetType(fv,PETSCFVLEASTSQUARES);CHKERRQ(ierr);
212c4762a1bSJed Brown     ierr = PetscFVSetNumComponents(fv,Nf);CHKERRQ(ierr);
213c4762a1bSJed Brown     ierr = PetscLimiterCreate(comm,&limiter);CHKERRQ(ierr);
214c4762a1bSJed Brown     ierr = PetscLimiterSetType(limiter,PETSCLIMITERNONE);CHKERRQ(ierr);
215c4762a1bSJed Brown     ierr = PetscFVSetLimiter(fv,limiter);CHKERRQ(ierr);
216c4762a1bSJed Brown     ierr = PetscLimiterDestroy(&limiter);CHKERRQ(ierr);
217c4762a1bSJed Brown     ierr = PetscFVSetFromOptions(fv);CHKERRQ(ierr);
218c4762a1bSJed Brown     ierr = DMSetField(base,0,NULL,(PetscObject)fv);CHKERRQ(ierr);
219c4762a1bSJed Brown     ierr = PetscFVDestroy(&fv);CHKERRQ(ierr);
220c4762a1bSJed Brown   } else {
221c4762a1bSJed Brown     PetscFE fe;
222c4762a1bSJed Brown 
223c4762a1bSJed Brown     ierr = PetscFECreateDefault(comm,dim,Nf,PETSC_FALSE,NULL,PETSC_DEFAULT,&fe);CHKERRQ(ierr);
224c4762a1bSJed Brown     ierr = DMSetField(base,0,NULL,(PetscObject)fe);CHKERRQ(ierr);
225c4762a1bSJed Brown     ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
226c4762a1bSJed Brown   }
227c4762a1bSJed Brown   ierr = DMCreateDS(base);CHKERRQ(ierr);
228c4762a1bSJed Brown 
229c4762a1bSJed Brown   if (use_bcs) {
230c4762a1bSJed Brown     PetscInt ids[] = {1, 2, 3, 4, 5, 6};
23145480ffeSMatthew G. Knepley     DMLabel  label;
232c4762a1bSJed Brown 
23345480ffeSMatthew G. Knepley     ierr = DMGetLabel(base, "marker", &label);CHKERRQ(ierr);
23445480ffeSMatthew G. Knepley     ierr = DMAddBoundary(base,DM_BC_ESSENTIAL, "bc", label, 2 * dim, ids, 0, 0, NULL, useFV ? (void(*)(void)) bc_func_fv : (void(*)(void)) funcs[0], NULL, useFV ? (void *) &bcCtx : NULL, NULL);CHKERRQ(ierr);
235c4762a1bSJed Brown   }
236c4762a1bSJed Brown   ierr = DMViewFromOptions(base,NULL,"-dm_base_view");CHKERRQ(ierr);
237c4762a1bSJed Brown 
238c4762a1bSJed Brown   /* the pre adaptivity forest */
239c4762a1bSJed Brown   ierr = DMCreate(comm,&preForest);CHKERRQ(ierr);
240c4762a1bSJed Brown   ierr = DMSetType(preForest,(dim == 2) ? DMP4EST : DMP8EST);CHKERRQ(ierr);
241c4762a1bSJed Brown   ierr = DMCopyDisc(base,preForest);CHKERRQ(ierr);
242c4762a1bSJed Brown   ierr = DMForestSetBaseDM(preForest,base);CHKERRQ(ierr);
243c4762a1bSJed Brown   ierr = DMForestSetMinimumRefinement(preForest,0);CHKERRQ(ierr);
244c4762a1bSJed Brown   ierr = DMForestSetInitialRefinement(preForest,1);CHKERRQ(ierr);
245c4762a1bSJed Brown   ierr = DMSetFromOptions(preForest);CHKERRQ(ierr);
246c4762a1bSJed Brown   ierr = DMSetUp(preForest);CHKERRQ(ierr);
247c4762a1bSJed Brown   ierr = DMViewFromOptions(preForest,NULL,"-dm_pre_view");CHKERRQ(ierr);
248c4762a1bSJed Brown 
249c4762a1bSJed Brown   /* the pre adaptivity field */
250c4762a1bSJed Brown   ierr = DMCreateGlobalVector(preForest,&preVec);CHKERRQ(ierr);
251c4762a1bSJed Brown   ierr = DMProjectFunction(preForest,0.,funcs,ctxs,INSERT_VALUES,preVec);CHKERRQ(ierr);
252c4762a1bSJed Brown   ierr = VecViewFromOptions(preVec,NULL,"-vec_pre_view");CHKERRQ(ierr);
253c4762a1bSJed Brown 
254c4762a1bSJed Brown   /* communicate between base and pre adaptivity forest */
255c4762a1bSJed Brown   if (transfer_from_base[0]) {
256c4762a1bSJed Brown     Vec baseVec, baseVecMapped;
257c4762a1bSJed Brown 
258c4762a1bSJed Brown     ierr = DMGetGlobalVector(base,&baseVec);CHKERRQ(ierr);
259c4762a1bSJed Brown     ierr = DMProjectFunction(base,0.,funcs,ctxs,INSERT_VALUES,baseVec);CHKERRQ(ierr);
260c4762a1bSJed Brown     ierr = PetscObjectSetName((PetscObject)baseVec,"Function Base");CHKERRQ(ierr);
261c4762a1bSJed Brown     ierr = VecViewFromOptions(baseVec,NULL,"-vec_base_view");CHKERRQ(ierr);
262c4762a1bSJed Brown 
263c4762a1bSJed Brown     ierr = DMGetGlobalVector(preForest,&baseVecMapped);CHKERRQ(ierr);
264c4762a1bSJed Brown     ierr = DMForestTransferVecFromBase(preForest,baseVec,baseVecMapped);CHKERRQ(ierr);
265c4762a1bSJed Brown     ierr = VecViewFromOptions(baseVecMapped,NULL,"-vec_map_base_view");CHKERRQ(ierr);
266c4762a1bSJed Brown 
267c4762a1bSJed Brown     /* compare */
268c4762a1bSJed Brown     ierr = VecAXPY(baseVecMapped,-1.,preVec);CHKERRQ(ierr);
269c4762a1bSJed Brown     ierr = VecViewFromOptions(baseVecMapped,NULL,"-vec_map_diff_view");CHKERRQ(ierr);
270c4762a1bSJed Brown     ierr = VecNorm(baseVecMapped,NORM_2,&diff);CHKERRQ(ierr);
271c4762a1bSJed Brown 
272c4762a1bSJed Brown     /* output */
273c4762a1bSJed Brown     if (diff < tol) {
274c4762a1bSJed Brown       ierr = PetscPrintf(comm,"DMForestTransferVecFromBase() passes.\n");CHKERRQ(ierr);
275c4762a1bSJed Brown     } else {
276c4762a1bSJed Brown       ierr = PetscPrintf(comm,"DMForestTransferVecFromBase() fails with error %g and tolerance %g\n",(double)diff,(double)tol);CHKERRQ(ierr);
277c4762a1bSJed Brown     }
278c4762a1bSJed Brown 
279c4762a1bSJed Brown     ierr = DMRestoreGlobalVector(base,&baseVec);CHKERRQ(ierr);
280c4762a1bSJed Brown     ierr = DMRestoreGlobalVector(preForest,&baseVecMapped);CHKERRQ(ierr);
281c4762a1bSJed Brown   }
282c4762a1bSJed Brown 
283c4762a1bSJed Brown   for (step = 0; step < adaptSteps; ++step) {
284c4762a1bSJed Brown 
285c4762a1bSJed Brown     if (!transfer_from_base[1]) {
286c4762a1bSJed Brown       ierr = PetscObjectGetReference((PetscObject)preForest,&preCount);CHKERRQ(ierr);
287c4762a1bSJed Brown     }
288c4762a1bSJed Brown 
289c4762a1bSJed Brown     /* adapt */
290c4762a1bSJed Brown     ierr = CreateAdaptivityLabel(preForest,&adaptLabel);CHKERRQ(ierr);
291c4762a1bSJed Brown     ierr = DMForestTemplate(preForest,comm,&postForest);CHKERRQ(ierr);
292c4762a1bSJed Brown     if (step) { ierr = DMForestSetAdaptivityLabel(postForest,adaptLabel);CHKERRQ(ierr); }
293c4762a1bSJed Brown     ierr = DMLabelDestroy(&adaptLabel);CHKERRQ(ierr);
294c4762a1bSJed Brown     ierr = DMSetUp(postForest);CHKERRQ(ierr);
295c4762a1bSJed Brown     ierr = DMViewFromOptions(postForest,NULL,"-dm_post_view");CHKERRQ(ierr);
296c4762a1bSJed Brown 
297c4762a1bSJed Brown     /* transfer */
298c4762a1bSJed Brown     ierr = DMCreateGlobalVector(postForest,&postVecTransfer);CHKERRQ(ierr);
299c4762a1bSJed Brown     ierr = DMForestTransferVec(preForest,preVec,postForest,postVecTransfer,PETSC_TRUE,0.0);CHKERRQ(ierr);
300c4762a1bSJed Brown     ierr = VecViewFromOptions(postVecTransfer,NULL,"-vec_post_transfer_view");CHKERRQ(ierr);
301c4762a1bSJed Brown 
302c4762a1bSJed Brown     /* the exact post adaptivity field */
303c4762a1bSJed Brown     ierr = DMCreateGlobalVector(postForest,&postVecExact);CHKERRQ(ierr);
304c4762a1bSJed Brown     ierr = DMProjectFunction(postForest,0.,funcs,ctxs,INSERT_VALUES,postVecExact);CHKERRQ(ierr);
305c4762a1bSJed Brown     ierr = VecViewFromOptions(postVecExact,NULL,"-vec_post_exact_view");CHKERRQ(ierr);
306c4762a1bSJed Brown 
307c4762a1bSJed Brown     /* compare */
308c4762a1bSJed Brown     ierr = VecAXPY(postVecExact,-1.,postVecTransfer);CHKERRQ(ierr);
309c4762a1bSJed Brown     ierr = VecViewFromOptions(postVecExact,NULL,"-vec_diff_view");CHKERRQ(ierr);
310c4762a1bSJed Brown     ierr = VecNorm(postVecExact,NORM_2,&diff);CHKERRQ(ierr);
311c4762a1bSJed Brown 
312c4762a1bSJed Brown     /* output */
313c4762a1bSJed Brown     if (diff < tol) {
314c4762a1bSJed Brown       ierr = PetscPrintf(comm,"DMForestTransferVec() passes.\n");CHKERRQ(ierr);
315c4762a1bSJed Brown     } else {
316c4762a1bSJed Brown       ierr = PetscPrintf(comm,"DMForestTransferVec() fails with error %g and tolerance %g\n",(double)diff,(double)tol);CHKERRQ(ierr);
317c4762a1bSJed Brown       ierr = IdentifyBadPoints(postForest, postVecExact, tol);CHKERRQ(ierr);
318c4762a1bSJed Brown     }
319c4762a1bSJed Brown     ierr = VecDestroy(&postVecExact);CHKERRQ(ierr);
320c4762a1bSJed Brown 
321c4762a1bSJed Brown     /* disconnect preForest from postForest if we don't test the transfer throughout the entire refinement process */
322c4762a1bSJed Brown     if (!transfer_from_base[1]) {
323c4762a1bSJed Brown       ierr = DMForestSetAdaptivityForest(postForest,NULL);CHKERRQ(ierr);
324c4762a1bSJed Brown       ierr = PetscObjectGetReference((PetscObject)preForest,&postCount);CHKERRQ(ierr);
325c4762a1bSJed Brown       if (postCount != preCount) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Adaptation not memory neutral: reference count increase from %d to %d\n",preCount,postCount);
326c4762a1bSJed Brown     }
327c4762a1bSJed Brown 
328c4762a1bSJed Brown     if (conv) {
329c4762a1bSJed Brown       DM dmConv;
330c4762a1bSJed Brown 
331c4762a1bSJed Brown       ierr = DMConvert(postForest,DMPLEX,&dmConv);CHKERRQ(ierr);
332c4762a1bSJed Brown       ierr = DMViewFromOptions(dmConv,NULL,"-dm_conv_view");CHKERRQ(ierr);
333c4762a1bSJed Brown       ierr = DMPlexCheckCellShape(dmConv,PETSC_TRUE,PETSC_DETERMINE);CHKERRQ(ierr);
334c4762a1bSJed Brown       ierr = DMDestroy(&dmConv);CHKERRQ(ierr);
335c4762a1bSJed Brown     }
336c4762a1bSJed Brown 
337c4762a1bSJed Brown     ierr = VecDestroy(&preVec);CHKERRQ(ierr);
338c4762a1bSJed Brown     ierr = DMDestroy(&preForest);CHKERRQ(ierr);
339c4762a1bSJed Brown 
340c4762a1bSJed Brown     preVec    = postVecTransfer;
341c4762a1bSJed Brown     preForest = postForest;
342c4762a1bSJed Brown   }
343c4762a1bSJed Brown 
344c4762a1bSJed Brown   if (transfer_from_base[1]) {
345c4762a1bSJed Brown     Vec baseVec, baseVecMapped;
346c4762a1bSJed Brown 
347c4762a1bSJed Brown     /* communicate between base and last adapted forest */
348c4762a1bSJed Brown     ierr = DMGetGlobalVector(base,&baseVec);CHKERRQ(ierr);
349c4762a1bSJed Brown     ierr = DMProjectFunction(base,0.,funcs,ctxs,INSERT_VALUES,baseVec);CHKERRQ(ierr);
350c4762a1bSJed Brown     ierr = PetscObjectSetName((PetscObject)baseVec,"Function Base");CHKERRQ(ierr);
351c4762a1bSJed Brown     ierr = VecViewFromOptions(baseVec,NULL,"-vec_base_view");CHKERRQ(ierr);
352c4762a1bSJed Brown 
353c4762a1bSJed Brown     ierr = DMGetGlobalVector(preForest,&baseVecMapped);CHKERRQ(ierr);
354c4762a1bSJed Brown     ierr = DMForestTransferVecFromBase(preForest,baseVec,baseVecMapped);CHKERRQ(ierr);
355c4762a1bSJed Brown     ierr = VecViewFromOptions(baseVecMapped,NULL,"-vec_map_base_view");CHKERRQ(ierr);
356c4762a1bSJed Brown 
357c4762a1bSJed Brown     /* compare */
358c4762a1bSJed Brown     ierr = VecAXPY(baseVecMapped,-1.,preVec);CHKERRQ(ierr);
359c4762a1bSJed Brown     ierr = VecViewFromOptions(baseVecMapped,NULL,"-vec_map_diff_view");CHKERRQ(ierr);
360c4762a1bSJed Brown     ierr = VecNorm(baseVecMapped,NORM_2,&diff);CHKERRQ(ierr);
361c4762a1bSJed Brown 
362c4762a1bSJed Brown     /* output */
363c4762a1bSJed Brown     if (diff < tol) {
364c4762a1bSJed Brown       ierr = PetscPrintf(comm,"DMForestTransferVecFromBase() passes.\n");CHKERRQ(ierr);
365c4762a1bSJed Brown     } else {
366c4762a1bSJed Brown       ierr = PetscPrintf(comm,"DMForestTransferVecFromBase() fails with error %g and tolerance %g\n",(double)diff,(double)tol);CHKERRQ(ierr);
367c4762a1bSJed Brown     }
368c4762a1bSJed Brown 
369c4762a1bSJed Brown     ierr = DMRestoreGlobalVector(base,&baseVec);CHKERRQ(ierr);
370c4762a1bSJed Brown     ierr = DMRestoreGlobalVector(preForest,&baseVecMapped);CHKERRQ(ierr);
371c4762a1bSJed Brown   }
372c4762a1bSJed Brown 
373c4762a1bSJed Brown   /* cleanup */
374c4762a1bSJed Brown   ierr = VecDestroy(&preVec);CHKERRQ(ierr);
375c4762a1bSJed Brown   ierr = DMDestroy(&preForest);CHKERRQ(ierr);
376c4762a1bSJed Brown   ierr = DMDestroy(&base);CHKERRQ(ierr);
377c4762a1bSJed Brown   ierr = PetscFinalize();
378c4762a1bSJed Brown   return ierr;
379c4762a1bSJed Brown }
380c4762a1bSJed Brown 
381c4762a1bSJed Brown /*TEST
38230602db0SMatthew G. Knepley   testset:
38330602db0SMatthew G. Knepley     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3 -petscspace_type tensor
384c4762a1bSJed Brown 
385c4762a1bSJed Brown     test:
386c4762a1bSJed Brown       output_file: output/ex2_2d.out
387c4762a1bSJed Brown       suffix: p4est_2d
38830602db0SMatthew G. Knepley       args: -petscspace_degree 2
389c4762a1bSJed Brown       nsize: 3
390c4762a1bSJed Brown       requires: p4est !single
391c4762a1bSJed Brown 
392c4762a1bSJed Brown     test:
393c4762a1bSJed Brown       output_file: output/ex2_2d.out
394c4762a1bSJed Brown       suffix: p4est_2d_deg4
39530602db0SMatthew G. Knepley       args: -petscspace_degree 4
396c4762a1bSJed Brown       requires: p4est !single
397c4762a1bSJed Brown 
398c4762a1bSJed Brown     test:
399c4762a1bSJed Brown       output_file: output/ex2_2d.out
400c4762a1bSJed Brown       suffix: p4est_2d_deg8
40130602db0SMatthew G. Knepley       args: -petscspace_degree 8
402c4762a1bSJed Brown       requires: p4est !single
403c4762a1bSJed Brown 
404c4762a1bSJed Brown     test:
405c4762a1bSJed Brown       output_file: output/ex2_steps2.out
406c4762a1bSJed Brown       suffix: p4est_2d_deg2_steps2
40730602db0SMatthew G. Knepley       args: -petscspace_degree 2 -coords -adapt_steps 2
408c4762a1bSJed Brown       nsize: 3
409c4762a1bSJed Brown       requires: p4est !single
410c4762a1bSJed Brown 
411c4762a1bSJed Brown     test:
412c4762a1bSJed Brown       output_file: output/ex2_steps3.out
413c4762a1bSJed Brown       suffix: p4est_2d_deg3_steps3
41430602db0SMatthew G. Knepley       args: -petscspace_degree 3 -coords -adapt_steps 3 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1
415c4762a1bSJed Brown       nsize: 3
416c4762a1bSJed Brown       requires: p4est !single
417c4762a1bSJed Brown 
418c4762a1bSJed Brown     test:
419c4762a1bSJed Brown       output_file: output/ex2_steps3.out
420c4762a1bSJed Brown       suffix: p4est_2d_deg3_steps3_L2_periodic
42130602db0SMatthew G. Knepley       args: -petscspace_degree 3 -petscdualspace_lagrange_continuity 0 -coords -adapt_steps 3 -dm_plex_box_bd periodic,periodic -use_bcs 0 -petscdualspace_lagrange_node_type equispaced
422c4762a1bSJed Brown       nsize: 3
423c4762a1bSJed Brown       requires: p4est !single
424c4762a1bSJed Brown 
425c4762a1bSJed Brown     test:
426c4762a1bSJed Brown       output_file: output/ex2_steps3.out
427c4762a1bSJed Brown       suffix: p4est_3d_deg2_steps3_L2_periodic
42830602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 2 -petscdualspace_lagrange_continuity 0 -coords -adapt_steps 3 -dm_plex_box_bd periodic,periodic,periodic -use_bcs 0
429c4762a1bSJed Brown       nsize: 3
430c4762a1bSJed Brown       requires: p4est !single
431c4762a1bSJed Brown 
432c4762a1bSJed Brown     test:
433c4762a1bSJed Brown       output_file: output/ex2_steps2.out
434c4762a1bSJed Brown       suffix: p4est_3d_deg2_steps2
43530602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 2 -coords -adapt_steps 2
436c4762a1bSJed Brown       nsize: 3
437c4762a1bSJed Brown       requires: p4est !single
438c4762a1bSJed Brown 
439c4762a1bSJed Brown     test:
440c4762a1bSJed Brown       output_file: output/ex2_steps3.out
441c4762a1bSJed Brown       suffix: p4est_3d_deg3_steps3
44230602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 3 -coords -adapt_steps 3 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1
443c4762a1bSJed Brown       nsize: 3
444c4762a1bSJed Brown       requires: p4est !single
445c4762a1bSJed Brown 
446c4762a1bSJed Brown     test:
447c4762a1bSJed Brown       output_file: output/ex2_3d.out
448c4762a1bSJed Brown       suffix: p4est_3d
44930602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 1
450c4762a1bSJed Brown       nsize: 3
451c4762a1bSJed Brown       requires: p4est !single
452c4762a1bSJed Brown 
453c4762a1bSJed Brown     test:
454c4762a1bSJed Brown       output_file: output/ex2_3d.out
455c4762a1bSJed Brown       suffix: p4est_3d_deg3
45630602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 3
457c4762a1bSJed Brown       nsize: 3
458c4762a1bSJed Brown       requires: p4est !single
459c4762a1bSJed Brown 
460c4762a1bSJed Brown     test:
461c4762a1bSJed Brown       output_file: output/ex2_2d.out
462c4762a1bSJed Brown       suffix: p4est_2d_deg2_coords
46330602db0SMatthew G. Knepley       args: -petscspace_degree 2 -coords
464c4762a1bSJed Brown       nsize: 3
465c4762a1bSJed Brown       requires: p4est !single
466c4762a1bSJed Brown 
467c4762a1bSJed Brown     test:
468c4762a1bSJed Brown       output_file: output/ex2_3d.out
469c4762a1bSJed Brown       suffix: p4est_3d_deg2_coords
47030602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 2 -coords
471c4762a1bSJed Brown       nsize: 3
472c4762a1bSJed Brown       requires: p4est !single
473c4762a1bSJed Brown 
474c4762a1bSJed Brown     test:
475c4762a1bSJed Brown       suffix: p4est_3d_nans
47630602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -dm_forest_partition_overlap 1 -test_convert -petscspace_degree 1
477c4762a1bSJed Brown       nsize: 2
478c4762a1bSJed Brown       requires: p4est !single
479c4762a1bSJed Brown 
480c4762a1bSJed Brown     test:
481c4762a1bSJed Brown       TODO: not broken, but the 3D case below is broken, so I do not trust this one
482c4762a1bSJed Brown       output_file: output/ex2_steps2.out
483c4762a1bSJed Brown       suffix: p4est_2d_tfb_distributed_nc
48430602db0SMatthew G. Knepley       args: -petscspace_degree 3 -dm_forest_maximum_refinement 2 -dm_p4est_refine_pattern hash -use_bcs 0 -coords -adapt_steps 2 -dm_distribute -petscpartitioner_type shell -petscpartitioner_shell_random
485c4762a1bSJed Brown       nsize: 3
486c4762a1bSJed Brown       requires: p4est !single
487c4762a1bSJed Brown 
488c4762a1bSJed Brown     test:
489c4762a1bSJed Brown       TODO: broken
490c4762a1bSJed Brown       output_file: output/ex2_steps2.out
491c4762a1bSJed Brown       suffix: p4est_3d_tfb_distributed_nc
49230602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 2 -dm_forest_maximum_refinement 2 -dm_p4est_refine_pattern hash -use_bcs 0 -coords -adapt_steps 2 -dm_distribute -petscpartitioner_type shell -petscpartitioner_shell_random
493c4762a1bSJed Brown       nsize: 3
494c4762a1bSJed Brown       requires: p4est !single
495c4762a1bSJed Brown 
49630602db0SMatthew G. Knepley   testset:
497012bc364SMatthew G. Knepley     args: -petscspace_type tensor -dm_coord_space 0 -dm_plex_transform_type refine_tobox
49830602db0SMatthew G. Knepley 
499c4762a1bSJed Brown     test:
500c4762a1bSJed Brown       TODO: broken
501c4762a1bSJed Brown       output_file: output/ex2_3d.out
502c4762a1bSJed Brown       suffix: p4est_3d_transfer_fails
503012bc364SMatthew G. Knepley       args: -petscspace_degree 1 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/doublet-tet.msh -adapt_steps 1 -dm_forest_initial_refinement 1 -use_bcs 0 -dm_refine
504c4762a1bSJed Brown       requires: p4est !single
505c4762a1bSJed Brown 
506c4762a1bSJed Brown     test:
507c4762a1bSJed Brown       TODO: broken
508c4762a1bSJed Brown       output_file: output/ex2_steps2_notfb.out
509c4762a1bSJed Brown       suffix: p4est_3d_transfer_fails_2
510012bc364SMatthew G. Knepley       args: -petscspace_degree 1 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/doublet-tet.msh -adapt_steps 2 -dm_forest_initial_refinement 0 -transfer_from_base 0 -use_bcs 0 -dm_refine
511c4762a1bSJed Brown       requires: p4est !single
512c4762a1bSJed Brown 
513c4762a1bSJed Brown     test:
514c4762a1bSJed Brown       output_file: output/ex2_steps2.out
515c4762a1bSJed Brown       suffix: p4est_3d_multi_transfer_s2t
516012bc364SMatthew G. Knepley       args: -petscspace_degree 3 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/doublet-tet.msh -adapt_steps 2 -dm_forest_initial_refinement 1 -petscdualspace_lagrange_continuity 0 -use_bcs 0 -dm_refine 1
517012bc364SMatthew G. Knepley       requires: p4est !single
518c4762a1bSJed Brown 
519c4762a1bSJed Brown     test:
520c4762a1bSJed Brown       output_file: output/ex2_steps2.out
521c4762a1bSJed Brown       suffix: p4est_3d_coords_transfer_s2t
522012bc364SMatthew G. Knepley       args: -petscspace_degree 3 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/doublet-tet.msh -adapt_steps 2 -dm_forest_initial_refinement 1 -petscdualspace_lagrange_continuity 0 -coords -use_bcs 0 -dm_refine 1
523012bc364SMatthew G. Knepley       requires: p4est !single
52430602db0SMatthew G. Knepley 
52530602db0SMatthew G. Knepley   testset:
52630602db0SMatthew G. Knepley     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3
52730602db0SMatthew G. Knepley 
52830602db0SMatthew G. Knepley     test:
52930602db0SMatthew G. Knepley       output_file: output/ex2_2d_fv.out
53030602db0SMatthew G. Knepley       suffix: p4est_2d_fv
53130602db0SMatthew G. Knepley       args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1
53230602db0SMatthew G. Knepley       nsize: 3
53330602db0SMatthew G. Knepley       requires: p4est !single
53430602db0SMatthew G. Knepley 
53530602db0SMatthew G. Knepley     test:
53630602db0SMatthew G. Knepley       TODO: broken (codimension adjacency)
53730602db0SMatthew G. Knepley       output_file: output/ex2_2d_fv.out
53830602db0SMatthew G. Knepley       suffix: p4est_2d_fv_adjcodim
53930602db0SMatthew G. Knepley       args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 -dm_forest_adjacency_codimension 1
54030602db0SMatthew G. Knepley       nsize: 2
54130602db0SMatthew G. Knepley       requires: p4est !single
54230602db0SMatthew G. Knepley 
54330602db0SMatthew G. Knepley     test:
54430602db0SMatthew G. Knepley       TODO: broken (dimension adjacency)
54530602db0SMatthew G. Knepley       output_file: output/ex2_2d_fv.out
54630602db0SMatthew G. Knepley       suffix: p4est_2d_fv_adjdim
54730602db0SMatthew G. Knepley       args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 -dm_forest_adjacency_dimension 1
54830602db0SMatthew G. Knepley       nsize: 2
54930602db0SMatthew G. Knepley       requires: p4est !single
55030602db0SMatthew G. Knepley 
55130602db0SMatthew G. Knepley     test:
55230602db0SMatthew G. Knepley       output_file: output/ex2_2d_fv.out
55330602db0SMatthew G. Knepley       suffix: p4est_2d_fv_zerocells
55430602db0SMatthew G. Knepley       args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1
55530602db0SMatthew G. Knepley       nsize: 10
55630602db0SMatthew G. Knepley       requires: p4est !single
55730602db0SMatthew G. Knepley 
55830602db0SMatthew G. Knepley     test:
55930602db0SMatthew G. Knepley       output_file: output/ex2_3d_fv.out
56030602db0SMatthew G. Knepley       suffix: p4est_3d_fv
56130602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1
56230602db0SMatthew G. Knepley       nsize: 3
563c4762a1bSJed Brown       requires: p4est !single
564c4762a1bSJed Brown 
565c4762a1bSJed Brown TEST*/
566