xref: /petsc/src/dm/impls/forest/tests/ex2.c (revision 48a46eb9bd028bec07ec0f396b1a3abb43f14558)
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 
89371c9d4SSatish Balay static PetscErrorCode AddIdentityLabel(DM dm) {
9c4762a1bSJed Brown   PetscInt pStart, pEnd, p;
10c4762a1bSJed Brown 
11c4762a1bSJed Brown   PetscFunctionBegin;
129566063dSJacob Faibussowitsch   PetscCall(DMCreateLabel(dm, "identity"));
139566063dSJacob Faibussowitsch   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
149566063dSJacob Faibussowitsch   for (p = pStart; p < pEnd; p++) PetscCall(DMSetLabelValue(dm, "identity", p, p));
15c4762a1bSJed Brown   PetscFunctionReturn(0);
16c4762a1bSJed Brown }
17c4762a1bSJed Brown 
189371c9d4SSatish Balay static PetscErrorCode CreateAdaptivityLabel(DM forest, DMLabel *adaptLabel) {
19c4762a1bSJed Brown   DMLabel  identLabel;
20c4762a1bSJed Brown   PetscInt cStart, cEnd, c;
21c4762a1bSJed Brown 
22c4762a1bSJed Brown   PetscFunctionBegin;
239566063dSJacob Faibussowitsch   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", adaptLabel));
249566063dSJacob Faibussowitsch   PetscCall(DMLabelSetDefaultValue(*adaptLabel, DM_ADAPT_COARSEN));
259566063dSJacob Faibussowitsch   PetscCall(DMGetLabel(forest, "identity", &identLabel));
269566063dSJacob Faibussowitsch   PetscCall(DMForestGetCellChart(forest, &cStart, &cEnd));
27c4762a1bSJed Brown   for (c = cStart; c < cEnd; c++) {
28c4762a1bSJed Brown     PetscInt basePoint;
29c4762a1bSJed Brown 
309566063dSJacob Faibussowitsch     PetscCall(DMLabelGetValue(identLabel, c, &basePoint));
319566063dSJacob Faibussowitsch     if (!basePoint) PetscCall(DMLabelSetValue(*adaptLabel, c, DM_ADAPT_REFINE));
32c4762a1bSJed Brown   }
33c4762a1bSJed Brown   PetscFunctionReturn(0);
34c4762a1bSJed Brown }
35c4762a1bSJed Brown 
369371c9d4SSatish Balay static PetscErrorCode LinearFunction(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) {
37c4762a1bSJed Brown   PetscFunctionBeginUser;
38c4762a1bSJed Brown   u[0] = (x[0] * 2.0 + 1.) + (x[1] * 20.0 + 10.) + ((dim == 3) ? (x[2] * 200.0 + 100.) : 0.);
39c4762a1bSJed Brown   PetscFunctionReturn(0);
40c4762a1bSJed Brown }
41c4762a1bSJed Brown 
429371c9d4SSatish Balay static PetscErrorCode MultiaffineFunction(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) {
43c4762a1bSJed Brown   PetscFunctionBeginUser;
44c4762a1bSJed Brown   u[0] = (x[0] * 1.0 + 2.0) * (x[1] * 3.0 - 4.0) * ((dim == 3) ? (x[2] * 5.0 + 6.0) : 1.);
45c4762a1bSJed Brown   PetscFunctionReturn(0);
46c4762a1bSJed Brown }
47c4762a1bSJed Brown 
489371c9d4SSatish Balay static PetscErrorCode CoordsFunction(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) {
49c4762a1bSJed Brown   PetscInt f;
50c4762a1bSJed Brown 
51c4762a1bSJed Brown   PetscFunctionBeginUser;
52c4762a1bSJed Brown   for (f = 0; f < Nf; f++) u[f] = x[f];
53c4762a1bSJed Brown   PetscFunctionReturn(0);
54c4762a1bSJed Brown }
55c4762a1bSJed Brown 
569371c9d4SSatish Balay typedef struct _bc_func_ctx {
57c4762a1bSJed Brown   PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
58c4762a1bSJed Brown   PetscInt dim;
59c4762a1bSJed Brown   PetscInt Nf;
60c4762a1bSJed Brown   void    *ctx;
619371c9d4SSatish Balay } bc_func_ctx;
62c4762a1bSJed Brown 
639371c9d4SSatish Balay static PetscErrorCode bc_func_fv(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx) {
64c4762a1bSJed Brown   bc_func_ctx *bcCtx;
65c4762a1bSJed Brown 
66c4762a1bSJed Brown   PetscFunctionBegin;
67c4762a1bSJed Brown   bcCtx = (bc_func_ctx *)ctx;
689566063dSJacob Faibussowitsch   PetscCall((bcCtx->func)(bcCtx->dim, time, c, bcCtx->Nf, xG, bcCtx->ctx));
69c4762a1bSJed Brown   PetscFunctionReturn(0);
70c4762a1bSJed Brown }
71c4762a1bSJed Brown 
729371c9d4SSatish Balay static PetscErrorCode IdentifyBadPoints(DM dm, Vec vec, PetscReal tol) {
73c4762a1bSJed Brown   DM           dmplex;
74c4762a1bSJed Brown   PetscInt     p, pStart, pEnd, maxDof;
75c4762a1bSJed Brown   Vec          vecLocal;
76c4762a1bSJed Brown   DMLabel      depthLabel;
77c4762a1bSJed Brown   PetscSection section;
78c4762a1bSJed Brown 
79c4762a1bSJed Brown   PetscFunctionBegin;
809566063dSJacob Faibussowitsch   PetscCall(DMCreateLocalVector(dm, &vecLocal));
819566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal));
829566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal));
839566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm, DMPLEX, &dmplex));
849566063dSJacob Faibussowitsch   PetscCall(DMPlexGetChart(dmplex, &pStart, &pEnd));
859566063dSJacob Faibussowitsch   PetscCall(DMPlexGetDepthLabel(dmplex, &depthLabel));
869566063dSJacob Faibussowitsch   PetscCall(DMGetLocalSection(dmplex, &section));
879566063dSJacob Faibussowitsch   PetscCall(PetscSectionGetMaxDof(section, &maxDof));
88c4762a1bSJed Brown   for (p = pStart; p < pEnd; p++) {
89c4762a1bSJed Brown     PetscInt     s, c, cSize, parent, childID, numChildren;
90c4762a1bSJed Brown     PetscInt     cl, closureSize, *closure = NULL;
91c4762a1bSJed Brown     PetscScalar *values = NULL;
92c4762a1bSJed Brown     PetscBool    bad    = PETSC_FALSE;
93c4762a1bSJed Brown 
949566063dSJacob Faibussowitsch     PetscCall(VecGetValuesSection(vecLocal, section, p, &values));
959566063dSJacob Faibussowitsch     PetscCall(PetscSectionGetDof(section, p, &cSize));
96c4762a1bSJed Brown     for (c = 0; c < cSize; c++) {
972f613bf5SBarry Smith       PetscReal absDiff = PetscAbsScalar(values[c]);
989371c9d4SSatish Balay       if (absDiff > tol) {
999371c9d4SSatish Balay         bad = PETSC_TRUE;
1009371c9d4SSatish Balay         break;
1019371c9d4SSatish Balay       }
102c4762a1bSJed Brown     }
103c4762a1bSJed Brown     if (!bad) continue;
10463a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Bad point %" PetscInt_FMT "\n", p));
1059566063dSJacob Faibussowitsch     PetscCall(DMLabelGetValue(depthLabel, p, &s));
10663a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Depth %" PetscInt_FMT "\n", s));
1079566063dSJacob Faibussowitsch     PetscCall(DMPlexGetTransitiveClosure(dmplex, p, PETSC_TRUE, &closureSize, &closure));
108c4762a1bSJed Brown     for (cl = 0; cl < closureSize; cl++) {
109c4762a1bSJed Brown       PetscInt cp = closure[2 * cl];
1109566063dSJacob Faibussowitsch       PetscCall(DMPlexGetTreeParent(dmplex, cp, &parent, &childID));
111*48a46eb9SPierre Jolivet       if (parent != cp) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Closure point %" PetscInt_FMT " (%" PetscInt_FMT ") child of %" PetscInt_FMT " (ID %" PetscInt_FMT ")\n", cl, cp, parent, childID));
1129566063dSJacob Faibussowitsch       PetscCall(DMPlexGetTreeChildren(dmplex, cp, &numChildren, NULL));
113*48a46eb9SPierre Jolivet       if (numChildren) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Closure point %" PetscInt_FMT " (%" PetscInt_FMT ") is parent\n", cl, cp));
114c4762a1bSJed Brown     }
1159566063dSJacob Faibussowitsch     PetscCall(DMPlexRestoreTransitiveClosure(dmplex, p, PETSC_TRUE, &closureSize, &closure));
116c4762a1bSJed Brown     for (c = 0; c < cSize; c++) {
1172f613bf5SBarry Smith       PetscReal absDiff = PetscAbsScalar(values[c]);
118*48a46eb9SPierre Jolivet       if (absDiff > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Bad dof %" PetscInt_FMT "\n", c));
119c4762a1bSJed Brown     }
120c4762a1bSJed Brown   }
1219566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dmplex));
1229566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&vecLocal));
123c4762a1bSJed Brown   PetscFunctionReturn(0);
124c4762a1bSJed Brown }
125c4762a1bSJed Brown 
1269371c9d4SSatish Balay int main(int argc, char **argv) {
127c4762a1bSJed Brown   MPI_Comm comm;
128c4762a1bSJed Brown   DM       base, preForest, postForest;
12930602db0SMatthew G. Knepley   PetscInt dim, Nf          = 1;
130c4762a1bSJed Brown   PetscInt step, adaptSteps = 1;
131c4762a1bSJed Brown   PetscInt preCount, postCount;
132c4762a1bSJed Brown   Vec      preVec, postVecTransfer, postVecExact;
133c4762a1bSJed Brown   PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *) = {MultiaffineFunction};
134c4762a1bSJed Brown   void       *ctxs[1]                                                                                 = {NULL};
135c4762a1bSJed Brown   PetscReal   diff, tol = PETSC_SMALL;
136c4762a1bSJed Brown   PetscBool   linear                = PETSC_FALSE;
137c4762a1bSJed Brown   PetscBool   coords                = PETSC_FALSE;
138c4762a1bSJed Brown   PetscBool   useFV                 = PETSC_FALSE;
139c4762a1bSJed Brown   PetscBool   conv                  = PETSC_FALSE;
140c4762a1bSJed Brown   PetscBool   transfer_from_base[2] = {PETSC_TRUE, PETSC_FALSE};
141c4762a1bSJed Brown   PetscBool   use_bcs               = PETSC_TRUE;
142c4762a1bSJed Brown   bc_func_ctx bcCtx;
143c4762a1bSJed Brown   DMLabel     adaptLabel;
144c4762a1bSJed Brown 
145327415f7SBarry Smith   PetscFunctionBeginUser;
1469566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
147c4762a1bSJed Brown   comm = PETSC_COMM_WORLD;
148d0609cedSBarry Smith   PetscOptionsBegin(comm, "", "DMForestTransferVec() Test Options", "DMFOREST");
1499566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-linear", "Transfer a simple linear function", "ex2.c", linear, &linear, NULL));
1509566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-coords", "Transfer a simple coordinate function", "ex2.c", coords, &coords, NULL));
1519566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-use_fv", "Use a finite volume approximation", "ex2.c", useFV, &useFV, NULL));
1529566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-test_convert", "Test conversion to DMPLEX", NULL, conv, &conv, NULL));
1539566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-transfer_from_base", "Transfer a vector from base DM to DMForest", "ex2.c", transfer_from_base[0], &transfer_from_base[0], NULL));
154c4762a1bSJed Brown   transfer_from_base[1] = transfer_from_base[0];
1559566063dSJacob Faibussowitsch   PetscCall(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));
1569566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-use_bcs", "Use dirichlet boundary conditions", "ex2.c", use_bcs, &use_bcs, NULL));
1579566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBoundedInt("-adapt_steps", "Number of adaptivity steps", "ex2.c", adaptSteps, &adaptSteps, NULL, 0));
158d0609cedSBarry Smith   PetscOptionsEnd();
159c4762a1bSJed Brown 
160c4762a1bSJed Brown   tol = PetscMax(1.e-10, tol); /* XXX fix for quadruple precision -> why do I need to do this? */
161c4762a1bSJed Brown 
16230602db0SMatthew G. Knepley   /* the base mesh */
1639566063dSJacob Faibussowitsch   PetscCall(DMCreate(comm, &base));
1649566063dSJacob Faibussowitsch   PetscCall(DMSetType(base, DMPLEX));
1659566063dSJacob Faibussowitsch   PetscCall(DMSetFromOptions(base));
16630602db0SMatthew G. Knepley 
1679566063dSJacob Faibussowitsch   PetscCall(AddIdentityLabel(base));
1689566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(base, &dim));
169c4762a1bSJed Brown 
1709371c9d4SSatish Balay   if (linear) { funcs[0] = LinearFunction; }
171c4762a1bSJed Brown   if (coords) {
172c4762a1bSJed Brown     funcs[0] = CoordsFunction;
173c4762a1bSJed Brown     Nf       = dim;
174c4762a1bSJed Brown   }
175c4762a1bSJed Brown 
176c4762a1bSJed Brown   bcCtx.func = funcs[0];
177c4762a1bSJed Brown   bcCtx.dim  = dim;
178c4762a1bSJed Brown   bcCtx.Nf   = Nf;
179c4762a1bSJed Brown   bcCtx.ctx  = NULL;
180c4762a1bSJed Brown 
181c4762a1bSJed Brown   if (useFV) {
182c4762a1bSJed Brown     PetscFV      fv;
183c4762a1bSJed Brown     PetscLimiter limiter;
184c4762a1bSJed Brown     DM           baseFV;
185c4762a1bSJed Brown 
1869566063dSJacob Faibussowitsch     PetscCall(DMPlexConstructGhostCells(base, NULL, NULL, &baseFV));
1879566063dSJacob Faibussowitsch     PetscCall(DMViewFromOptions(baseFV, NULL, "-fv_dm_view"));
1889566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&base));
189c4762a1bSJed Brown     base = baseFV;
19040cbb1a0SMatthew G. Knepley     PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv));
1919566063dSJacob Faibussowitsch     PetscCall(PetscFVSetSpatialDimension(fv, dim));
1929566063dSJacob Faibussowitsch     PetscCall(PetscFVSetType(fv, PETSCFVLEASTSQUARES));
1939566063dSJacob Faibussowitsch     PetscCall(PetscFVSetNumComponents(fv, Nf));
1949566063dSJacob Faibussowitsch     PetscCall(PetscLimiterCreate(comm, &limiter));
1959566063dSJacob Faibussowitsch     PetscCall(PetscLimiterSetType(limiter, PETSCLIMITERNONE));
1969566063dSJacob Faibussowitsch     PetscCall(PetscFVSetLimiter(fv, limiter));
1979566063dSJacob Faibussowitsch     PetscCall(PetscLimiterDestroy(&limiter));
1989566063dSJacob Faibussowitsch     PetscCall(PetscFVSetFromOptions(fv));
1999566063dSJacob Faibussowitsch     PetscCall(DMSetField(base, 0, NULL, (PetscObject)fv));
2009566063dSJacob Faibussowitsch     PetscCall(PetscFVDestroy(&fv));
201c4762a1bSJed Brown   } else {
202c4762a1bSJed Brown     PetscFE fe;
203c4762a1bSJed Brown 
2049566063dSJacob Faibussowitsch     PetscCall(PetscFECreateDefault(comm, dim, Nf, PETSC_FALSE, NULL, PETSC_DEFAULT, &fe));
2059566063dSJacob Faibussowitsch     PetscCall(DMSetField(base, 0, NULL, (PetscObject)fe));
2069566063dSJacob Faibussowitsch     PetscCall(PetscFEDestroy(&fe));
207c4762a1bSJed Brown   }
2089566063dSJacob Faibussowitsch   PetscCall(DMCreateDS(base));
209c4762a1bSJed Brown 
210c4762a1bSJed Brown   if (use_bcs) {
211c4762a1bSJed Brown     PetscInt ids[] = {1, 2, 3, 4, 5, 6};
21245480ffeSMatthew G. Knepley     DMLabel  label;
213c4762a1bSJed Brown 
2149566063dSJacob Faibussowitsch     PetscCall(DMGetLabel(base, "marker", &label));
2159566063dSJacob Faibussowitsch     PetscCall(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));
216c4762a1bSJed Brown   }
2179566063dSJacob Faibussowitsch   PetscCall(DMViewFromOptions(base, NULL, "-dm_base_view"));
218c4762a1bSJed Brown 
219c4762a1bSJed Brown   /* the pre adaptivity forest */
2209566063dSJacob Faibussowitsch   PetscCall(DMCreate(comm, &preForest));
2219566063dSJacob Faibussowitsch   PetscCall(DMSetType(preForest, (dim == 2) ? DMP4EST : DMP8EST));
2229566063dSJacob Faibussowitsch   PetscCall(DMCopyDisc(base, preForest));
2239566063dSJacob Faibussowitsch   PetscCall(DMForestSetBaseDM(preForest, base));
2249566063dSJacob Faibussowitsch   PetscCall(DMForestSetMinimumRefinement(preForest, 0));
2259566063dSJacob Faibussowitsch   PetscCall(DMForestSetInitialRefinement(preForest, 1));
2269566063dSJacob Faibussowitsch   PetscCall(DMSetFromOptions(preForest));
2279566063dSJacob Faibussowitsch   PetscCall(DMSetUp(preForest));
2289566063dSJacob Faibussowitsch   PetscCall(DMViewFromOptions(preForest, NULL, "-dm_pre_view"));
229c4762a1bSJed Brown 
230c4762a1bSJed Brown   /* the pre adaptivity field */
2319566063dSJacob Faibussowitsch   PetscCall(DMCreateGlobalVector(preForest, &preVec));
2329566063dSJacob Faibussowitsch   PetscCall(DMProjectFunction(preForest, 0., funcs, ctxs, INSERT_VALUES, preVec));
2339566063dSJacob Faibussowitsch   PetscCall(VecViewFromOptions(preVec, NULL, "-vec_pre_view"));
234c4762a1bSJed Brown 
235c4762a1bSJed Brown   /* communicate between base and pre adaptivity forest */
236c4762a1bSJed Brown   if (transfer_from_base[0]) {
237c4762a1bSJed Brown     Vec baseVec, baseVecMapped;
238c4762a1bSJed Brown 
2399566063dSJacob Faibussowitsch     PetscCall(DMGetGlobalVector(base, &baseVec));
2409566063dSJacob Faibussowitsch     PetscCall(DMProjectFunction(base, 0., funcs, ctxs, INSERT_VALUES, baseVec));
2419566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetName((PetscObject)baseVec, "Function Base"));
2429566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(baseVec, NULL, "-vec_base_view"));
243c4762a1bSJed Brown 
2449566063dSJacob Faibussowitsch     PetscCall(DMGetGlobalVector(preForest, &baseVecMapped));
2459566063dSJacob Faibussowitsch     PetscCall(DMForestTransferVecFromBase(preForest, baseVec, baseVecMapped));
2469566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(baseVecMapped, NULL, "-vec_map_base_view"));
247c4762a1bSJed Brown 
248c4762a1bSJed Brown     /* compare */
2499566063dSJacob Faibussowitsch     PetscCall(VecAXPY(baseVecMapped, -1., preVec));
2509566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(baseVecMapped, NULL, "-vec_map_diff_view"));
2519566063dSJacob Faibussowitsch     PetscCall(VecNorm(baseVecMapped, NORM_2, &diff));
252c4762a1bSJed Brown 
253c4762a1bSJed Brown     /* output */
254c4762a1bSJed Brown     if (diff < tol) {
2559566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(comm, "DMForestTransferVecFromBase() passes.\n"));
256c4762a1bSJed Brown     } else {
2579566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(comm, "DMForestTransferVecFromBase() fails with error %g and tolerance %g\n", (double)diff, (double)tol));
258c4762a1bSJed Brown     }
259c4762a1bSJed Brown 
2609566063dSJacob Faibussowitsch     PetscCall(DMRestoreGlobalVector(base, &baseVec));
2619566063dSJacob Faibussowitsch     PetscCall(DMRestoreGlobalVector(preForest, &baseVecMapped));
262c4762a1bSJed Brown   }
263c4762a1bSJed Brown 
264c4762a1bSJed Brown   for (step = 0; step < adaptSteps; ++step) {
265*48a46eb9SPierre Jolivet     if (!transfer_from_base[1]) PetscCall(PetscObjectGetReference((PetscObject)preForest, &preCount));
266c4762a1bSJed Brown 
267c4762a1bSJed Brown     /* adapt */
2689566063dSJacob Faibussowitsch     PetscCall(CreateAdaptivityLabel(preForest, &adaptLabel));
2699566063dSJacob Faibussowitsch     PetscCall(DMForestTemplate(preForest, comm, &postForest));
2709566063dSJacob Faibussowitsch     if (step) PetscCall(DMForestSetAdaptivityLabel(postForest, adaptLabel));
2719566063dSJacob Faibussowitsch     PetscCall(DMLabelDestroy(&adaptLabel));
2729566063dSJacob Faibussowitsch     PetscCall(DMSetUp(postForest));
2739566063dSJacob Faibussowitsch     PetscCall(DMViewFromOptions(postForest, NULL, "-dm_post_view"));
274c4762a1bSJed Brown 
275c4762a1bSJed Brown     /* transfer */
2769566063dSJacob Faibussowitsch     PetscCall(DMCreateGlobalVector(postForest, &postVecTransfer));
2779566063dSJacob Faibussowitsch     PetscCall(DMForestTransferVec(preForest, preVec, postForest, postVecTransfer, PETSC_TRUE, 0.0));
2789566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(postVecTransfer, NULL, "-vec_post_transfer_view"));
279c4762a1bSJed Brown 
280c4762a1bSJed Brown     /* the exact post adaptivity field */
2819566063dSJacob Faibussowitsch     PetscCall(DMCreateGlobalVector(postForest, &postVecExact));
2829566063dSJacob Faibussowitsch     PetscCall(DMProjectFunction(postForest, 0., funcs, ctxs, INSERT_VALUES, postVecExact));
2839566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(postVecExact, NULL, "-vec_post_exact_view"));
284c4762a1bSJed Brown 
285c4762a1bSJed Brown     /* compare */
2869566063dSJacob Faibussowitsch     PetscCall(VecAXPY(postVecExact, -1., postVecTransfer));
2879566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(postVecExact, NULL, "-vec_diff_view"));
2889566063dSJacob Faibussowitsch     PetscCall(VecNorm(postVecExact, NORM_2, &diff));
289c4762a1bSJed Brown 
290c4762a1bSJed Brown     /* output */
291c4762a1bSJed Brown     if (diff < tol) {
2929566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(comm, "DMForestTransferVec() passes.\n"));
293c4762a1bSJed Brown     } else {
2949566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(comm, "DMForestTransferVec() fails with error %g and tolerance %g\n", (double)diff, (double)tol));
2959566063dSJacob Faibussowitsch       PetscCall(IdentifyBadPoints(postForest, postVecExact, tol));
296c4762a1bSJed Brown     }
2979566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&postVecExact));
298c4762a1bSJed Brown 
299c4762a1bSJed Brown     /* disconnect preForest from postForest if we don't test the transfer throughout the entire refinement process */
300c4762a1bSJed Brown     if (!transfer_from_base[1]) {
3019566063dSJacob Faibussowitsch       PetscCall(DMForestSetAdaptivityForest(postForest, NULL));
3029566063dSJacob Faibussowitsch       PetscCall(PetscObjectGetReference((PetscObject)preForest, &postCount));
30363a3b9bcSJacob Faibussowitsch       PetscCheck(postCount == preCount, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Adaptation not memory neutral: reference count increase from %" PetscInt_FMT " to %" PetscInt_FMT, preCount, postCount);
304c4762a1bSJed Brown     }
305c4762a1bSJed Brown 
306c4762a1bSJed Brown     if (conv) {
307c4762a1bSJed Brown       DM dmConv;
308c4762a1bSJed Brown 
3099566063dSJacob Faibussowitsch       PetscCall(DMConvert(postForest, DMPLEX, &dmConv));
3109566063dSJacob Faibussowitsch       PetscCall(DMViewFromOptions(dmConv, NULL, "-dm_conv_view"));
3119566063dSJacob Faibussowitsch       PetscCall(DMPlexCheckCellShape(dmConv, PETSC_TRUE, PETSC_DETERMINE));
3129566063dSJacob Faibussowitsch       PetscCall(DMDestroy(&dmConv));
313c4762a1bSJed Brown     }
314c4762a1bSJed Brown 
3159566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&preVec));
3169566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&preForest));
317c4762a1bSJed Brown 
318c4762a1bSJed Brown     preVec    = postVecTransfer;
319c4762a1bSJed Brown     preForest = postForest;
320c4762a1bSJed Brown   }
321c4762a1bSJed Brown 
322c4762a1bSJed Brown   if (transfer_from_base[1]) {
323c4762a1bSJed Brown     Vec baseVec, baseVecMapped;
324c4762a1bSJed Brown 
325c4762a1bSJed Brown     /* communicate between base and last adapted forest */
3269566063dSJacob Faibussowitsch     PetscCall(DMGetGlobalVector(base, &baseVec));
3279566063dSJacob Faibussowitsch     PetscCall(DMProjectFunction(base, 0., funcs, ctxs, INSERT_VALUES, baseVec));
3289566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetName((PetscObject)baseVec, "Function Base"));
3299566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(baseVec, NULL, "-vec_base_view"));
330c4762a1bSJed Brown 
3319566063dSJacob Faibussowitsch     PetscCall(DMGetGlobalVector(preForest, &baseVecMapped));
3329566063dSJacob Faibussowitsch     PetscCall(DMForestTransferVecFromBase(preForest, baseVec, baseVecMapped));
3339566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(baseVecMapped, NULL, "-vec_map_base_view"));
334c4762a1bSJed Brown 
335c4762a1bSJed Brown     /* compare */
3369566063dSJacob Faibussowitsch     PetscCall(VecAXPY(baseVecMapped, -1., preVec));
3379566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(baseVecMapped, NULL, "-vec_map_diff_view"));
3389566063dSJacob Faibussowitsch     PetscCall(VecNorm(baseVecMapped, NORM_2, &diff));
339c4762a1bSJed Brown 
340c4762a1bSJed Brown     /* output */
341c4762a1bSJed Brown     if (diff < tol) {
3429566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(comm, "DMForestTransferVecFromBase() passes.\n"));
343c4762a1bSJed Brown     } else {
3449566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(comm, "DMForestTransferVecFromBase() fails with error %g and tolerance %g\n", (double)diff, (double)tol));
345c4762a1bSJed Brown     }
346c4762a1bSJed Brown 
3479566063dSJacob Faibussowitsch     PetscCall(DMRestoreGlobalVector(base, &baseVec));
3489566063dSJacob Faibussowitsch     PetscCall(DMRestoreGlobalVector(preForest, &baseVecMapped));
349c4762a1bSJed Brown   }
350c4762a1bSJed Brown 
351c4762a1bSJed Brown   /* cleanup */
3529566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&preVec));
3539566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&preForest));
3549566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&base));
3559566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
356b122ec5aSJacob Faibussowitsch   return 0;
357c4762a1bSJed Brown }
358c4762a1bSJed Brown 
359c4762a1bSJed Brown /*TEST
36030602db0SMatthew G. Knepley   testset:
36130602db0SMatthew G. Knepley     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3 -petscspace_type tensor
362c4762a1bSJed Brown 
363c4762a1bSJed Brown     test:
364c4762a1bSJed Brown       output_file: output/ex2_2d.out
365c4762a1bSJed Brown       suffix: p4est_2d
36630602db0SMatthew G. Knepley       args: -petscspace_degree 2
367c4762a1bSJed Brown       nsize: 3
368c4762a1bSJed Brown       requires: p4est !single
369c4762a1bSJed Brown 
370c4762a1bSJed Brown     test:
371c4762a1bSJed Brown       output_file: output/ex2_2d.out
372c4762a1bSJed Brown       suffix: p4est_2d_deg4
37330602db0SMatthew G. Knepley       args: -petscspace_degree 4
374c4762a1bSJed Brown       requires: p4est !single
375c4762a1bSJed Brown 
376c4762a1bSJed Brown     test:
377c4762a1bSJed Brown       output_file: output/ex2_2d.out
378c4762a1bSJed Brown       suffix: p4est_2d_deg8
37930602db0SMatthew G. Knepley       args: -petscspace_degree 8
380c4762a1bSJed Brown       requires: p4est !single
381c4762a1bSJed Brown 
382c4762a1bSJed Brown     test:
383c4762a1bSJed Brown       output_file: output/ex2_steps2.out
384c4762a1bSJed Brown       suffix: p4est_2d_deg2_steps2
38530602db0SMatthew G. Knepley       args: -petscspace_degree 2 -coords -adapt_steps 2
386c4762a1bSJed Brown       nsize: 3
387c4762a1bSJed Brown       requires: p4est !single
388c4762a1bSJed Brown 
389c4762a1bSJed Brown     test:
390c4762a1bSJed Brown       output_file: output/ex2_steps3.out
391c4762a1bSJed Brown       suffix: p4est_2d_deg3_steps3
39230602db0SMatthew G. Knepley       args: -petscspace_degree 3 -coords -adapt_steps 3 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1
393c4762a1bSJed Brown       nsize: 3
394c4762a1bSJed Brown       requires: p4est !single
395c4762a1bSJed Brown 
396c4762a1bSJed Brown     test:
397c4762a1bSJed Brown       output_file: output/ex2_steps3.out
398c4762a1bSJed Brown       suffix: p4est_2d_deg3_steps3_L2_periodic
39930602db0SMatthew 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
400c4762a1bSJed Brown       nsize: 3
401c4762a1bSJed Brown       requires: p4est !single
402c4762a1bSJed Brown 
403c4762a1bSJed Brown     test:
404c4762a1bSJed Brown       output_file: output/ex2_steps3.out
405c4762a1bSJed Brown       suffix: p4est_3d_deg2_steps3_L2_periodic
40630602db0SMatthew 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
407c4762a1bSJed Brown       nsize: 3
408c4762a1bSJed Brown       requires: p4est !single
409c4762a1bSJed Brown 
410c4762a1bSJed Brown     test:
411c4762a1bSJed Brown       output_file: output/ex2_steps2.out
412c4762a1bSJed Brown       suffix: p4est_3d_deg2_steps2
41330602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 2 -coords -adapt_steps 2
414c4762a1bSJed Brown       nsize: 3
415c4762a1bSJed Brown       requires: p4est !single
416c4762a1bSJed Brown 
417c4762a1bSJed Brown     test:
418c4762a1bSJed Brown       output_file: output/ex2_steps3.out
419c4762a1bSJed Brown       suffix: p4est_3d_deg3_steps3
42030602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 3 -coords -adapt_steps 3 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1
421c4762a1bSJed Brown       nsize: 3
422c4762a1bSJed Brown       requires: p4est !single
423c4762a1bSJed Brown 
424c4762a1bSJed Brown     test:
425c4762a1bSJed Brown       output_file: output/ex2_3d.out
426c4762a1bSJed Brown       suffix: p4est_3d
42730602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 1
428c4762a1bSJed Brown       nsize: 3
429c4762a1bSJed Brown       requires: p4est !single
430c4762a1bSJed Brown 
431c4762a1bSJed Brown     test:
432c4762a1bSJed Brown       output_file: output/ex2_3d.out
433c4762a1bSJed Brown       suffix: p4est_3d_deg3
43430602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 3
435c4762a1bSJed Brown       nsize: 3
436c4762a1bSJed Brown       requires: p4est !single
437c4762a1bSJed Brown 
438c4762a1bSJed Brown     test:
439c4762a1bSJed Brown       output_file: output/ex2_2d.out
440c4762a1bSJed Brown       suffix: p4est_2d_deg2_coords
44130602db0SMatthew G. Knepley       args: -petscspace_degree 2 -coords
442c4762a1bSJed Brown       nsize: 3
443c4762a1bSJed Brown       requires: p4est !single
444c4762a1bSJed Brown 
445c4762a1bSJed Brown     test:
446c4762a1bSJed Brown       output_file: output/ex2_3d.out
447c4762a1bSJed Brown       suffix: p4est_3d_deg2_coords
44830602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 2 -coords
449c4762a1bSJed Brown       nsize: 3
450c4762a1bSJed Brown       requires: p4est !single
451c4762a1bSJed Brown 
452c4762a1bSJed Brown     test:
453c4762a1bSJed Brown       suffix: p4est_3d_nans
45430602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -dm_forest_partition_overlap 1 -test_convert -petscspace_degree 1
455c4762a1bSJed Brown       nsize: 2
456c4762a1bSJed Brown       requires: p4est !single
457c4762a1bSJed Brown 
458c4762a1bSJed Brown     test:
459c4762a1bSJed Brown       TODO: not broken, but the 3D case below is broken, so I do not trust this one
460c4762a1bSJed Brown       output_file: output/ex2_steps2.out
461c4762a1bSJed Brown       suffix: p4est_2d_tfb_distributed_nc
462e600fa54SMatthew G. Knepley       args: -petscspace_degree 3 -dm_forest_maximum_refinement 2 -dm_p4est_refine_pattern hash -use_bcs 0 -coords -adapt_steps 2 -petscpartitioner_type shell -petscpartitioner_shell_random
463c4762a1bSJed Brown       nsize: 3
464c4762a1bSJed Brown       requires: p4est !single
465c4762a1bSJed Brown 
466c4762a1bSJed Brown     test:
467c4762a1bSJed Brown       TODO: broken
468c4762a1bSJed Brown       output_file: output/ex2_steps2.out
469c4762a1bSJed Brown       suffix: p4est_3d_tfb_distributed_nc
470e600fa54SMatthew 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 -petscpartitioner_type shell -petscpartitioner_shell_random
471c4762a1bSJed Brown       nsize: 3
472c4762a1bSJed Brown       requires: p4est !single
473c4762a1bSJed Brown 
47430602db0SMatthew G. Knepley   testset:
475012bc364SMatthew G. Knepley     args: -petscspace_type tensor -dm_coord_space 0 -dm_plex_transform_type refine_tobox
47630602db0SMatthew G. Knepley 
477c4762a1bSJed Brown     test:
478c4762a1bSJed Brown       TODO: broken
479c4762a1bSJed Brown       output_file: output/ex2_3d.out
480c4762a1bSJed Brown       suffix: p4est_3d_transfer_fails
481012bc364SMatthew 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
482c4762a1bSJed Brown       requires: p4est !single
483c4762a1bSJed Brown 
484c4762a1bSJed Brown     test:
485c4762a1bSJed Brown       TODO: broken
486c4762a1bSJed Brown       output_file: output/ex2_steps2_notfb.out
487c4762a1bSJed Brown       suffix: p4est_3d_transfer_fails_2
488012bc364SMatthew 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
489c4762a1bSJed Brown       requires: p4est !single
490c4762a1bSJed Brown 
491c4762a1bSJed Brown     test:
492c4762a1bSJed Brown       output_file: output/ex2_steps2.out
493c4762a1bSJed Brown       suffix: p4est_3d_multi_transfer_s2t
494012bc364SMatthew 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
495012bc364SMatthew G. Knepley       requires: p4est !single
496c4762a1bSJed Brown 
497c4762a1bSJed Brown     test:
498c4762a1bSJed Brown       output_file: output/ex2_steps2.out
499c4762a1bSJed Brown       suffix: p4est_3d_coords_transfer_s2t
500012bc364SMatthew 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
501012bc364SMatthew G. Knepley       requires: p4est !single
50230602db0SMatthew G. Knepley 
50330602db0SMatthew G. Knepley   testset:
50430602db0SMatthew G. Knepley     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3
50530602db0SMatthew G. Knepley 
50630602db0SMatthew G. Knepley     test:
50730602db0SMatthew G. Knepley       output_file: output/ex2_2d_fv.out
50830602db0SMatthew G. Knepley       suffix: p4est_2d_fv
50930602db0SMatthew G. Knepley       args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1
51030602db0SMatthew G. Knepley       nsize: 3
51130602db0SMatthew G. Knepley       requires: p4est !single
51230602db0SMatthew G. Knepley 
51330602db0SMatthew G. Knepley     test:
51430602db0SMatthew G. Knepley       TODO: broken (codimension adjacency)
51530602db0SMatthew G. Knepley       output_file: output/ex2_2d_fv.out
51630602db0SMatthew G. Knepley       suffix: p4est_2d_fv_adjcodim
51730602db0SMatthew G. Knepley       args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 -dm_forest_adjacency_codimension 1
51830602db0SMatthew G. Knepley       nsize: 2
51930602db0SMatthew G. Knepley       requires: p4est !single
52030602db0SMatthew G. Knepley 
52130602db0SMatthew G. Knepley     test:
52230602db0SMatthew G. Knepley       TODO: broken (dimension adjacency)
52330602db0SMatthew G. Knepley       output_file: output/ex2_2d_fv.out
52430602db0SMatthew G. Knepley       suffix: p4est_2d_fv_adjdim
52530602db0SMatthew G. Knepley       args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 -dm_forest_adjacency_dimension 1
52630602db0SMatthew G. Knepley       nsize: 2
52730602db0SMatthew G. Knepley       requires: p4est !single
52830602db0SMatthew G. Knepley 
52930602db0SMatthew G. Knepley     test:
53030602db0SMatthew G. Knepley       output_file: output/ex2_2d_fv.out
53130602db0SMatthew G. Knepley       suffix: p4est_2d_fv_zerocells
53230602db0SMatthew G. Knepley       args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1
53330602db0SMatthew G. Knepley       nsize: 10
53430602db0SMatthew G. Knepley       requires: p4est !single
53530602db0SMatthew G. Knepley 
53630602db0SMatthew G. Knepley     test:
53730602db0SMatthew G. Knepley       output_file: output/ex2_3d_fv.out
53830602db0SMatthew G. Knepley       suffix: p4est_3d_fv
53930602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1
54030602db0SMatthew G. Knepley       nsize: 3
541c4762a1bSJed Brown       requires: p4est !single
542c4762a1bSJed Brown 
543c4762a1bSJed Brown TEST*/
544