xref: /petsc/src/dm/impls/forest/tests/ex2.c (revision bef158480efac06de457f7a665168877ab3c2fd7)
1 static char help[] = "Create a mesh, refine and coarsen simultaneously, and transfer a field\n\n";
2 
3 #include <petscds.h>
4 #include <petscdmplex.h>
5 #include <petscdmforest.h>
6 #include <petscoptions.h>
7 
8 static PetscErrorCode AddIdentityLabel(DM dm)
9 {
10   PetscInt       pStart,pEnd,p;
11   PetscErrorCode ierr;
12 
13   PetscFunctionBegin;
14   ierr = DMCreateLabel(dm, "identity");CHKERRQ(ierr);
15   ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
16   for (p = pStart; p < pEnd; p++) {ierr = DMSetLabelValue(dm, "identity", p, p);CHKERRQ(ierr);}
17   PetscFunctionReturn(0);
18 }
19 
20 static PetscErrorCode CreateAdaptivityLabel(DM forest,DMLabel *adaptLabel)
21 {
22   DMLabel        identLabel;
23   PetscInt       cStart, cEnd, c;
24   PetscErrorCode ierr;
25 
26   PetscFunctionBegin;
27   ierr = DMLabelCreate(PETSC_COMM_SELF,"adapt",adaptLabel);CHKERRQ(ierr);
28   ierr = DMLabelSetDefaultValue(*adaptLabel,DM_ADAPT_COARSEN);CHKERRQ(ierr);
29   ierr = DMGetLabel(forest,"identity",&identLabel);CHKERRQ(ierr);
30   ierr = DMForestGetCellChart(forest,&cStart,&cEnd);CHKERRQ(ierr);
31   for (c = cStart; c < cEnd; c++) {
32     PetscInt basePoint;
33 
34     ierr = DMLabelGetValue(identLabel,c,&basePoint);CHKERRQ(ierr);
35     if (!basePoint) {ierr = DMLabelSetValue(*adaptLabel,c,DM_ADAPT_REFINE);CHKERRQ(ierr);}
36   }
37   PetscFunctionReturn(0);
38 }
39 
40 static PetscErrorCode LinearFunction(PetscInt dim,PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
41 {
42   PetscFunctionBeginUser;
43   u[0] = (x[0] * 2.0 + 1.) + (x[1] * 20.0 + 10.) + ((dim == 3) ? (x[2] * 200.0 + 100.) : 0.);
44   PetscFunctionReturn(0);
45 }
46 
47 static PetscErrorCode MultiaffineFunction(PetscInt dim,PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
48 {
49   PetscFunctionBeginUser;
50   u[0] = (x[0] * 1.0 + 2.0) * (x[1] * 3.0 - 4.0) * ((dim == 3) ? (x[2] * 5.0 + 6.0) : 1.);
51   PetscFunctionReturn(0);
52 }
53 
54 static PetscErrorCode CoordsFunction(PetscInt dim,PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
55 {
56   PetscInt f;
57 
58   PetscFunctionBeginUser;
59   for (f=0;f<Nf;f++) u[f] = x[f];
60   PetscFunctionReturn(0);
61 }
62 
63 typedef struct _bc_func_ctx
64 {
65   PetscErrorCode (*func) (PetscInt,PetscReal,const PetscReal [], PetscInt, PetscScalar [], void *);
66   PetscInt dim;
67   PetscInt Nf;
68   void *ctx;
69 }
70 bc_func_ctx;
71 
72 static PetscErrorCode bc_func_fv (PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
73 {
74   bc_func_ctx    *bcCtx;
75   PetscErrorCode ierr;
76 
77   PetscFunctionBegin;
78   bcCtx = (bc_func_ctx *) ctx;
79   ierr = (bcCtx->func)(bcCtx->dim,time,c,bcCtx->Nf,xG,bcCtx->ctx);CHKERRQ(ierr);
80   PetscFunctionReturn(0);
81 }
82 
83 static PetscErrorCode IdentifyBadPoints (DM dm, Vec vec, PetscReal tol)
84 {
85   DM             dmplex;
86   PetscInt       p, pStart, pEnd, maxDof;
87   Vec            vecLocal;
88   DMLabel        depthLabel;
89   PetscSection   section;
90   PetscErrorCode ierr;
91 
92   PetscFunctionBegin;
93   ierr = DMCreateLocalVector(dm, &vecLocal);CHKERRQ(ierr);
94   ierr = DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal);CHKERRQ(ierr);
95   ierr = DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal);CHKERRQ(ierr);
96   ierr = DMConvert(dm ,DMPLEX, &dmplex);CHKERRQ(ierr);
97   ierr = DMPlexGetChart(dmplex, &pStart, &pEnd);CHKERRQ(ierr);
98   ierr = DMPlexGetDepthLabel(dmplex, &depthLabel);CHKERRQ(ierr);
99   ierr = DMGetLocalSection(dmplex, &section);CHKERRQ(ierr);
100   ierr = PetscSectionGetMaxDof(section, &maxDof);CHKERRQ(ierr);
101   for (p = pStart; p < pEnd; p++) {
102     PetscInt     s, c, cSize, parent, childID, numChildren;
103     PetscInt     cl, closureSize, *closure = NULL;
104     PetscScalar *values = NULL;
105     PetscBool    bad = PETSC_FALSE;
106 
107     ierr = VecGetValuesSection(vecLocal, section, p, &values);CHKERRQ(ierr);
108     ierr = PetscSectionGetDof(section, p, &cSize);CHKERRQ(ierr);
109     for (c = 0; c < cSize; c++) {
110       PetscReal absDiff = PetscAbsScalar(values[c]);CHKERRQ(ierr);
111       if (absDiff > tol) {bad = PETSC_TRUE; break;}
112     }
113     if (!bad) continue;
114     ierr = PetscPrintf(PETSC_COMM_SELF, "Bad point %D\n", p);CHKERRQ(ierr);
115     ierr = DMLabelGetValue(depthLabel, p, &s);CHKERRQ(ierr);
116     ierr = PetscPrintf(PETSC_COMM_SELF, "  Depth %D\n", s);CHKERRQ(ierr);
117     ierr = DMPlexGetTransitiveClosure(dmplex, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
118     for (cl = 0; cl < closureSize; cl++) {
119       PetscInt cp = closure[2 * cl];
120       ierr = DMPlexGetTreeParent(dmplex, cp, &parent, &childID);CHKERRQ(ierr);
121       if (parent != cp) {
122         ierr = PetscPrintf(PETSC_COMM_SELF, "  Closure point %D (%D) child of %D (ID %D)\n", cl, cp, parent, childID);CHKERRQ(ierr);
123       }
124       ierr = DMPlexGetTreeChildren(dmplex, cp, &numChildren, NULL);CHKERRQ(ierr);
125       if (numChildren) {
126         ierr = PetscPrintf(PETSC_COMM_SELF, "  Closure point %D (%D) is parent\n", cl, cp);CHKERRQ(ierr);
127       }
128     }
129     ierr = DMPlexRestoreTransitiveClosure(dmplex, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
130     for (c = 0; c < cSize; c++) {
131       PetscReal absDiff = PetscAbsScalar(values[c]);CHKERRQ(ierr);
132       if (absDiff > tol) {
133         ierr = PetscPrintf(PETSC_COMM_SELF, "  Bad dof %D\n", c);CHKERRQ(ierr);
134       }
135     }
136   }
137   ierr = DMDestroy(&dmplex);CHKERRQ(ierr);
138   ierr = VecDestroy(&vecLocal);CHKERRQ(ierr);
139   PetscFunctionReturn(0);
140 }
141 
142 int main(int argc, char **argv)
143 {
144   MPI_Comm       comm;
145   DM             base, preForest, postForest;
146   PetscInt       dim = 2, Nf = 1;
147   PetscInt       step, adaptSteps = 1;
148   PetscInt       preCount, postCount;
149   Vec            preVec, postVecTransfer, postVecExact;
150   PetscErrorCode (*funcs[1]) (PetscInt,PetscReal,const PetscReal [],PetscInt,PetscScalar [], void *) = {MultiaffineFunction};
151   char           filename[PETSC_MAX_PATH_LEN];
152   void           *ctxs[1] = {NULL};
153   const PetscInt cells[] = {3, 3, 3};
154   PetscReal      diff, tol = PETSC_SMALL;
155   DMBoundaryType periodicity[] = {DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE};
156   PetscBool      linear = PETSC_FALSE;
157   PetscBool      coords = PETSC_FALSE;
158   PetscBool      useFV = PETSC_FALSE;
159   PetscBool      conv = PETSC_FALSE;
160   PetscBool      distribute_base = PETSC_FALSE;
161   PetscBool      transfer_from_base[2] = {PETSC_TRUE,PETSC_FALSE};
162   PetscBool      use_bcs = PETSC_TRUE;
163   PetscBool      periodic = PETSC_FALSE;
164   bc_func_ctx    bcCtx;
165   DMLabel        adaptLabel;
166   size_t         len;
167   PetscErrorCode ierr;
168 
169   filename[0] = '\0';
170   ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
171   comm = PETSC_COMM_WORLD;
172   ierr = PetscOptionsBegin(comm, "", "DMForestTransferVec() Test Options", "DMFOREST");CHKERRQ(ierr);
173   ierr = PetscOptionsRangeInt("-dim", "The dimension (2 or 3)", "ex2.c", dim, &dim, NULL,2,3);CHKERRQ(ierr);
174   ierr = PetscOptionsBool("-linear","Transfer a simple linear function", "ex2.c", linear, &linear, NULL);CHKERRQ(ierr);
175   ierr = PetscOptionsBool("-coords","Transfer a simple coordinate function", "ex2.c", coords, &coords, NULL);CHKERRQ(ierr);
176   ierr = PetscOptionsBool("-use_fv","Use a finite volume approximation", "ex2.c", useFV, &useFV, NULL);CHKERRQ(ierr);
177   ierr = PetscOptionsBool("-test_convert","Test conversion to DMPLEX",NULL,conv,&conv,NULL);CHKERRQ(ierr);
178   ierr = PetscOptionsString("-filename", "Read the base mesh from file", "ex2.c", filename, filename, sizeof(filename), NULL);CHKERRQ(ierr);
179   ierr = PetscOptionsBool("-distribute_base","Distribute base DM", "ex2.c", distribute_base, &distribute_base, NULL);CHKERRQ(ierr);
180   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);
181   transfer_from_base[1] = transfer_from_base[0];
182   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);
183   ierr = PetscOptionsBool("-use_bcs","Use dirichlet boundary conditions", "ex2.c", use_bcs, &use_bcs, NULL);CHKERRQ(ierr);
184   ierr = PetscOptionsBoundedInt("-adapt_steps","Number of adaptivity steps", "ex2.c", adaptSteps, &adaptSteps, NULL,0);CHKERRQ(ierr);
185   ierr = PetscOptionsBool("-periodic","Use periodic box mesh", "ex2.c", periodic, &periodic, NULL);CHKERRQ(ierr);
186   ierr = PetscOptionsEnd();CHKERRQ(ierr);
187 
188   tol = PetscMax(1.e-10,tol); /* XXX fix for quadruple precision -> why do I need to do this? */
189 
190   if (periodic) periodicity[0] = periodicity[1] = periodicity[2] = DM_BOUNDARY_PERIODIC;
191 
192   if (linear) {
193     funcs[0] = LinearFunction;
194   }
195   if (coords) {
196     funcs[0] = CoordsFunction;
197     Nf = dim;
198   }
199 
200   bcCtx.func = funcs[0];
201   bcCtx.dim  = dim;
202   bcCtx.Nf   = Nf;
203   bcCtx.ctx  = NULL;
204 
205   /* the base mesh */
206   ierr = PetscStrlen(filename, &len);CHKERRQ(ierr);
207   if (!len) {
208     ierr = DMPlexCreateBoxMesh(comm,dim,PETSC_FALSE,cells,NULL,NULL,periodicity,PETSC_TRUE,&base);CHKERRQ(ierr);
209   } else {
210     DM tdm = NULL;
211 
212     ierr = DMPlexCreateFromFile(comm,filename,PETSC_TRUE,&base);CHKERRQ(ierr);
213     ierr = DMGetDimension(base,&dim);CHKERRQ(ierr);
214     ierr = DMPlexSetCellRefinerType(base, DM_REFINER_TO_BOX);CHKERRQ(ierr);
215     ierr = DMRefine(base, PETSC_COMM_WORLD, &tdm);CHKERRQ(ierr);
216     if (tdm) {
217       ierr = DMDestroy(&base);CHKERRQ(ierr);
218       base = tdm;
219     }
220     use_bcs = PETSC_FALSE;
221   }
222   ierr = AddIdentityLabel(base);CHKERRQ(ierr);
223   if (distribute_base) {
224     DM               baseParallel;
225     PetscPartitioner part;
226 
227     ierr = DMPlexGetPartitioner(base,&part);CHKERRQ(ierr);
228     ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);
229     ierr = DMPlexDistribute(base,0,NULL,&baseParallel);CHKERRQ(ierr);
230     if (baseParallel) {
231       ierr = DMDestroy(&base);CHKERRQ(ierr);
232       base = baseParallel;
233     }
234   }
235   if (useFV) {
236     PetscFV      fv;
237     PetscLimiter limiter;
238     DM           baseFV;
239 
240     ierr = DMPlexConstructGhostCells(base,NULL,NULL,&baseFV);CHKERRQ(ierr);
241     ierr = DMViewFromOptions(baseFV, NULL, "-fv_dm_view");CHKERRQ(ierr);
242     ierr = DMDestroy(&base);CHKERRQ(ierr);
243     base = baseFV;
244     ierr = PetscFVCreate(comm, &fv);CHKERRQ(ierr);
245     ierr = PetscFVSetSpatialDimension(fv,dim);CHKERRQ(ierr);
246     ierr = PetscFVSetType(fv,PETSCFVLEASTSQUARES);CHKERRQ(ierr);
247     ierr = PetscFVSetNumComponents(fv,Nf);CHKERRQ(ierr);
248     ierr = PetscLimiterCreate(comm,&limiter);CHKERRQ(ierr);
249     ierr = PetscLimiterSetType(limiter,PETSCLIMITERNONE);CHKERRQ(ierr);
250     ierr = PetscFVSetLimiter(fv,limiter);CHKERRQ(ierr);
251     ierr = PetscLimiterDestroy(&limiter);CHKERRQ(ierr);
252     ierr = PetscFVSetFromOptions(fv);CHKERRQ(ierr);
253     ierr = DMSetField(base,0,NULL,(PetscObject)fv);CHKERRQ(ierr);
254     ierr = PetscFVDestroy(&fv);CHKERRQ(ierr);
255   } else {
256     PetscFE fe;
257 
258     ierr = PetscFECreateDefault(comm,dim,Nf,PETSC_FALSE,NULL,PETSC_DEFAULT,&fe);CHKERRQ(ierr);
259     ierr = DMSetField(base,0,NULL,(PetscObject)fe);CHKERRQ(ierr);
260     ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
261   }
262   ierr = DMCreateDS(base);CHKERRQ(ierr);
263 
264   if (use_bcs) {
265     PetscInt ids[] = {1, 2, 3, 4, 5, 6};
266 
267     ierr = DMAddBoundary(base,DM_BC_ESSENTIAL, "bc", "marker", 0, 0, NULL, useFV ? (void(*)(void)) bc_func_fv : (void(*)(void)) funcs[0], NULL, 2 * dim, ids, useFV ? (void *) &bcCtx : NULL);CHKERRQ(ierr);
268   }
269   ierr = DMViewFromOptions(base,NULL,"-dm_base_view");CHKERRQ(ierr);
270 
271   /* the pre adaptivity forest */
272   ierr = DMCreate(comm,&preForest);CHKERRQ(ierr);
273   ierr = DMSetType(preForest,(dim == 2) ? DMP4EST : DMP8EST);CHKERRQ(ierr);
274   ierr = DMCopyDisc(base,preForest);CHKERRQ(ierr);
275   ierr = DMForestSetBaseDM(preForest,base);CHKERRQ(ierr);
276   ierr = DMForestSetMinimumRefinement(preForest,0);CHKERRQ(ierr);
277   ierr = DMForestSetInitialRefinement(preForest,1);CHKERRQ(ierr);
278   ierr = DMSetFromOptions(preForest);CHKERRQ(ierr);
279   ierr = DMSetUp(preForest);CHKERRQ(ierr);
280   ierr = DMViewFromOptions(preForest,NULL,"-dm_pre_view");CHKERRQ(ierr);
281 
282   /* the pre adaptivity field */
283   ierr = DMCreateGlobalVector(preForest,&preVec);CHKERRQ(ierr);
284   ierr = DMProjectFunction(preForest,0.,funcs,ctxs,INSERT_VALUES,preVec);CHKERRQ(ierr);
285   ierr = VecViewFromOptions(preVec,NULL,"-vec_pre_view");CHKERRQ(ierr);
286 
287   /* communicate between base and pre adaptivity forest */
288   if (transfer_from_base[0]) {
289     Vec baseVec, baseVecMapped;
290 
291     ierr = DMGetGlobalVector(base,&baseVec);CHKERRQ(ierr);
292     ierr = DMProjectFunction(base,0.,funcs,ctxs,INSERT_VALUES,baseVec);CHKERRQ(ierr);
293     ierr = PetscObjectSetName((PetscObject)baseVec,"Function Base");CHKERRQ(ierr);
294     ierr = VecViewFromOptions(baseVec,NULL,"-vec_base_view");CHKERRQ(ierr);
295 
296     ierr = DMGetGlobalVector(preForest,&baseVecMapped);CHKERRQ(ierr);
297     ierr = DMForestTransferVecFromBase(preForest,baseVec,baseVecMapped);CHKERRQ(ierr);
298     ierr = VecViewFromOptions(baseVecMapped,NULL,"-vec_map_base_view");CHKERRQ(ierr);
299 
300     /* compare */
301     ierr = VecAXPY(baseVecMapped,-1.,preVec);CHKERRQ(ierr);
302     ierr = VecViewFromOptions(baseVecMapped,NULL,"-vec_map_diff_view");CHKERRQ(ierr);
303     ierr = VecNorm(baseVecMapped,NORM_2,&diff);CHKERRQ(ierr);
304 
305     /* output */
306     if (diff < tol) {
307       ierr = PetscPrintf(comm,"DMForestTransferVecFromBase() passes.\n");CHKERRQ(ierr);
308     } else {
309       ierr = PetscPrintf(comm,"DMForestTransferVecFromBase() fails with error %g and tolerance %g\n",(double)diff,(double)tol);CHKERRQ(ierr);
310     }
311 
312     ierr = DMRestoreGlobalVector(base,&baseVec);CHKERRQ(ierr);
313     ierr = DMRestoreGlobalVector(preForest,&baseVecMapped);CHKERRQ(ierr);
314   }
315 
316   for (step = 0; step < adaptSteps; ++step) {
317 
318     if (!transfer_from_base[1]) {
319       ierr = PetscObjectGetReference((PetscObject)preForest,&preCount);CHKERRQ(ierr);
320     }
321 
322     /* adapt */
323     ierr = CreateAdaptivityLabel(preForest,&adaptLabel);CHKERRQ(ierr);
324     ierr = DMForestTemplate(preForest,comm,&postForest);CHKERRQ(ierr);
325     if (step) { ierr = DMForestSetAdaptivityLabel(postForest,adaptLabel);CHKERRQ(ierr); }
326     ierr = DMLabelDestroy(&adaptLabel);CHKERRQ(ierr);
327     ierr = DMSetUp(postForest);CHKERRQ(ierr);
328     ierr = DMViewFromOptions(postForest,NULL,"-dm_post_view");CHKERRQ(ierr);
329 
330     /* transfer */
331     ierr = DMCreateGlobalVector(postForest,&postVecTransfer);CHKERRQ(ierr);
332     ierr = DMForestTransferVec(preForest,preVec,postForest,postVecTransfer,PETSC_TRUE,0.0);CHKERRQ(ierr);
333     ierr = VecViewFromOptions(postVecTransfer,NULL,"-vec_post_transfer_view");CHKERRQ(ierr);
334 
335     /* the exact post adaptivity field */
336     ierr = DMCreateGlobalVector(postForest,&postVecExact);CHKERRQ(ierr);
337     ierr = DMProjectFunction(postForest,0.,funcs,ctxs,INSERT_VALUES,postVecExact);CHKERRQ(ierr);
338     ierr = VecViewFromOptions(postVecExact,NULL,"-vec_post_exact_view");CHKERRQ(ierr);
339 
340     /* compare */
341     ierr = VecAXPY(postVecExact,-1.,postVecTransfer);CHKERRQ(ierr);
342     ierr = VecViewFromOptions(postVecExact,NULL,"-vec_diff_view");CHKERRQ(ierr);
343     ierr = VecNorm(postVecExact,NORM_2,&diff);CHKERRQ(ierr);
344 
345     /* output */
346     if (diff < tol) {
347       ierr = PetscPrintf(comm,"DMForestTransferVec() passes.\n");CHKERRQ(ierr);
348     } else {
349       ierr = PetscPrintf(comm,"DMForestTransferVec() fails with error %g and tolerance %g\n",(double)diff,(double)tol);CHKERRQ(ierr);
350       ierr = IdentifyBadPoints(postForest, postVecExact, tol);CHKERRQ(ierr);
351     }
352     ierr = VecDestroy(&postVecExact);CHKERRQ(ierr);
353 
354     /* disconnect preForest from postForest if we don't test the transfer throughout the entire refinement process */
355     if (!transfer_from_base[1]) {
356       ierr = DMForestSetAdaptivityForest(postForest,NULL);CHKERRQ(ierr);
357       ierr = PetscObjectGetReference((PetscObject)preForest,&postCount);CHKERRQ(ierr);
358       if (postCount != preCount) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Adaptation not memory neutral: reference count increase from %d to %d\n",preCount,postCount);
359     }
360 
361     if (conv) {
362       DM dmConv;
363 
364       ierr = DMConvert(postForest,DMPLEX,&dmConv);CHKERRQ(ierr);
365       ierr = DMViewFromOptions(dmConv,NULL,"-dm_conv_view");CHKERRQ(ierr);
366       ierr = DMPlexCheckCellShape(dmConv,PETSC_TRUE,PETSC_DETERMINE);CHKERRQ(ierr);
367       ierr = DMDestroy(&dmConv);CHKERRQ(ierr);
368     }
369 
370     ierr = VecDestroy(&preVec);CHKERRQ(ierr);
371     ierr = DMDestroy(&preForest);CHKERRQ(ierr);
372 
373     preVec    = postVecTransfer;
374     preForest = postForest;
375   }
376 
377   if (transfer_from_base[1]) {
378     Vec baseVec, baseVecMapped;
379 
380     /* communicate between base and last adapted forest */
381     ierr = DMGetGlobalVector(base,&baseVec);CHKERRQ(ierr);
382     ierr = DMProjectFunction(base,0.,funcs,ctxs,INSERT_VALUES,baseVec);CHKERRQ(ierr);
383     ierr = PetscObjectSetName((PetscObject)baseVec,"Function Base");CHKERRQ(ierr);
384     ierr = VecViewFromOptions(baseVec,NULL,"-vec_base_view");CHKERRQ(ierr);
385 
386     ierr = DMGetGlobalVector(preForest,&baseVecMapped);CHKERRQ(ierr);
387     ierr = DMForestTransferVecFromBase(preForest,baseVec,baseVecMapped);CHKERRQ(ierr);
388     ierr = VecViewFromOptions(baseVecMapped,NULL,"-vec_map_base_view");CHKERRQ(ierr);
389 
390     /* compare */
391     ierr = VecAXPY(baseVecMapped,-1.,preVec);CHKERRQ(ierr);
392     ierr = VecViewFromOptions(baseVecMapped,NULL,"-vec_map_diff_view");CHKERRQ(ierr);
393     ierr = VecNorm(baseVecMapped,NORM_2,&diff);CHKERRQ(ierr);
394 
395     /* output */
396     if (diff < tol) {
397       ierr = PetscPrintf(comm,"DMForestTransferVecFromBase() passes.\n");CHKERRQ(ierr);
398     } else {
399       ierr = PetscPrintf(comm,"DMForestTransferVecFromBase() fails with error %g and tolerance %g\n",(double)diff,(double)tol);CHKERRQ(ierr);
400     }
401 
402     ierr = DMRestoreGlobalVector(base,&baseVec);CHKERRQ(ierr);
403     ierr = DMRestoreGlobalVector(preForest,&baseVecMapped);CHKERRQ(ierr);
404   }
405 
406   /* cleanup */
407   ierr = VecDestroy(&preVec);CHKERRQ(ierr);
408   ierr = DMDestroy(&preForest);CHKERRQ(ierr);
409   ierr = DMDestroy(&base);CHKERRQ(ierr);
410   ierr = PetscFinalize();
411   return ierr;
412 }
413 
414 /*TEST
415 
416      test:
417        output_file: output/ex2_2d.out
418        suffix: p4est_2d
419        args: -petscspace_type tensor -petscspace_degree 2 -dim 2
420        nsize: 3
421        requires: p4est !single
422 
423      test:
424        output_file: output/ex2_2d.out
425        suffix: p4est_2d_deg4
426        args: -petscspace_type tensor -petscspace_degree 4 -dim 2
427        requires: p4est !single
428 
429      test:
430        output_file: output/ex2_2d.out
431        suffix: p4est_2d_deg8
432        args: -petscspace_type tensor -petscspace_degree 8 -dim 2
433        requires: p4est !single
434 
435      test:
436        output_file: output/ex2_steps2.out
437        suffix: p4est_2d_deg2_steps2
438        args: -petscspace_type tensor -petscspace_degree 2 -dim 2 -coords -adapt_steps 2
439        nsize: 3
440        requires: p4est !single
441 
442      test:
443        output_file: output/ex2_steps3.out
444        suffix: p4est_2d_deg3_steps3
445        args: -petscspace_type tensor -petscspace_degree 3 -dim 2 -coords -adapt_steps 3 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1
446        nsize: 3
447        requires: p4est !single
448 
449      test:
450        output_file: output/ex2_steps3.out
451        suffix: p4est_2d_deg3_steps3_L2_periodic
452        args: -petscspace_type tensor -petscspace_degree 3 -petscdualspace_lagrange_continuity 0 -dim 2 -coords -adapt_steps 3 -periodic -use_bcs 0 -petscdualspace_lagrange_node_type equispaced
453        nsize: 3
454        requires: p4est !single
455 
456      test:
457        output_file: output/ex2_steps3.out
458        suffix: p4est_3d_deg2_steps3_L2_periodic
459        args: -petscspace_type tensor -petscspace_degree 2 -petscdualspace_lagrange_continuity 0 -dim 3 -coords -adapt_steps 3 -periodic -use_bcs 0
460        nsize: 3
461        requires: p4est !single
462 
463      test:
464        output_file: output/ex2_steps2.out
465        suffix: p4est_3d_deg2_steps2
466        args: -petscspace_type tensor -petscspace_degree 2 -dim 3 -coords -adapt_steps 2
467        nsize: 3
468        requires: p4est !single
469 
470      test:
471        output_file: output/ex2_steps3.out
472        suffix: p4est_3d_deg3_steps3
473        args: -petscspace_type tensor -petscspace_degree 3 -dim 3 -coords -adapt_steps 3 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1
474        nsize: 3
475        requires: p4est !single
476 
477      test:
478        output_file: output/ex2_2d_fv.out
479        suffix: p4est_2d_fv
480        args: -transfer_from_base 0 -use_fv -linear -dim 2 -dm_forest_partition_overlap 1
481        nsize: 3
482        requires: p4est !single
483 
484      test:
485        TODO: broken (codimension adjacency)
486        output_file: output/ex2_2d_fv.out
487        suffix: p4est_2d_fv_adjcodim
488        args: -transfer_from_base 0 -use_fv -linear -dim 2 -dm_forest_partition_overlap 1 -dm_forest_adjacency_codimension 1
489        nsize: 2
490        requires: p4est !single
491 
492      test:
493        TODO: broken (dimension adjacency)
494        output_file: output/ex2_2d_fv.out
495        suffix: p4est_2d_fv_adjdim
496        args: -transfer_from_base 0 -use_fv -linear -dim 2 -dm_forest_partition_overlap 1 -dm_forest_adjacency_dimension 1
497        nsize: 2
498        requires: p4est !single
499 
500      test:
501        output_file: output/ex2_2d_fv.out
502        suffix: p4est_2d_fv_zerocells
503        args: -transfer_from_base 0 -use_fv -linear -dim 2 -dm_forest_partition_overlap 1
504        nsize: 10
505        requires: p4est !single
506 
507      test:
508        output_file: output/ex2_3d.out
509        suffix: p4est_3d
510        args: -petscspace_type tensor -petscspace_degree 1 -dim 3
511        nsize: 3
512        requires: p4est !single
513 
514      test:
515        output_file: output/ex2_3d.out
516        suffix: p4est_3d_deg3
517        args: -petscspace_type tensor -petscspace_degree 3 -dim 3
518        nsize: 3
519        requires: p4est !single
520 
521      test:
522        output_file: output/ex2_2d.out
523        suffix: p4est_2d_deg2_coords
524        args: -petscspace_type tensor -petscspace_degree 2 -dim 2 -coords
525        nsize: 3
526        requires: p4est !single
527 
528      test:
529        output_file: output/ex2_3d.out
530        suffix: p4est_3d_deg2_coords
531        args: -petscspace_type tensor -petscspace_degree 2 -dim 3 -coords
532        nsize: 3
533        requires: p4est !single
534 
535      test:
536        output_file: output/ex2_3d_fv.out
537        suffix: p4est_3d_fv
538        args: -transfer_from_base 0 -use_fv -linear -dim 3 -dm_forest_partition_overlap 1
539        nsize: 3
540        requires: p4est !single
541 
542      test:
543        suffix: p4est_3d_nans
544        args: -dim 3 -dm_forest_partition_overlap 1 -test_convert -petscspace_type tensor -petscspace_degree 1
545        nsize: 2
546        requires: p4est !single
547 
548      test:
549        TODO: not broken, but the 3D case below is broken, so I do not trust this one
550        output_file: output/ex2_steps2.out
551        suffix: p4est_2d_tfb_distributed_nc
552        args: -petscspace_type tensor -petscspace_degree 3 -dim 2 -dm_forest_maximum_refinement 2 -dm_p4est_refine_pattern hash -use_bcs 0 -coords -adapt_steps 2 -distribute_base -petscpartitioner_type shell -petscpartitioner_shell_random
553        nsize: 3
554        requires: p4est !single
555 
556      test:
557        TODO: broken
558        output_file: output/ex2_steps2.out
559        suffix: p4est_3d_tfb_distributed_nc
560        args: -petscspace_type tensor -petscspace_degree 2 -dim 3 -dm_forest_maximum_refinement 2 -dm_p4est_refine_pattern hash -use_bcs 0 -coords -adapt_steps 2 -distribute_base -petscpartitioner_type shell -petscpartitioner_shell_random
561        nsize: 3
562        requires: p4est !single
563 
564      test:
565        TODO: broken
566        output_file: output/ex2_3d.out
567        suffix: p4est_3d_transfer_fails
568        args: -petscspace_type tensor -petscspace_degree 1 -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/doublet-tet.msh -adapt_steps 1 -dm_forest_initial_refinement 1
569        requires: p4est !single
570 
571      test:
572        TODO: broken
573        output_file: output/ex2_steps2_notfb.out
574        suffix: p4est_3d_transfer_fails_2
575        args: -petscspace_type tensor -petscspace_degree 1 -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/doublet-tet.msh -adapt_steps 2 -dm_forest_initial_refinement 0 -transfer_from_base 0
576        requires: p4est !single
577 
578      test:
579        output_file: output/ex2_steps2.out
580        suffix: p4est_3d_multi_transfer_s2t
581        args: -petscspace_type tensor -petscspace_degree 3 -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
582        requires: p4est !single
583 
584      test:
585        output_file: output/ex2_steps2.out
586        suffix: p4est_3d_coords_transfer_s2t
587        args: -petscspace_type tensor -petscspace_degree 3 -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
588        requires: p4est !single
589 
590 TEST*/
591