xref: /petsc/src/dm/impls/forest/tests/ex2.c (revision 40cbb1a031ea8f2be4fe2b92dc842b003ad37be3)
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 
12c4762a1bSJed Brown   PetscFunctionBegin;
139566063dSJacob Faibussowitsch   PetscCall(DMCreateLabel(dm, "identity"));
149566063dSJacob Faibussowitsch   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
159566063dSJacob Faibussowitsch   for (p = pStart; p < pEnd; p++) PetscCall(DMSetLabelValue(dm, "identity", p, p));
16c4762a1bSJed Brown   PetscFunctionReturn(0);
17c4762a1bSJed Brown }
18c4762a1bSJed Brown 
19c4762a1bSJed Brown static PetscErrorCode CreateAdaptivityLabel(DM forest,DMLabel *adaptLabel)
20c4762a1bSJed Brown {
21c4762a1bSJed Brown   DMLabel        identLabel;
22c4762a1bSJed Brown   PetscInt       cStart, cEnd, c;
23c4762a1bSJed Brown 
24c4762a1bSJed Brown   PetscFunctionBegin;
259566063dSJacob Faibussowitsch   PetscCall(DMLabelCreate(PETSC_COMM_SELF,"adapt",adaptLabel));
269566063dSJacob Faibussowitsch   PetscCall(DMLabelSetDefaultValue(*adaptLabel,DM_ADAPT_COARSEN));
279566063dSJacob Faibussowitsch   PetscCall(DMGetLabel(forest,"identity",&identLabel));
289566063dSJacob Faibussowitsch   PetscCall(DMForestGetCellChart(forest,&cStart,&cEnd));
29c4762a1bSJed Brown   for (c = cStart; c < cEnd; c++) {
30c4762a1bSJed Brown     PetscInt basePoint;
31c4762a1bSJed Brown 
329566063dSJacob Faibussowitsch     PetscCall(DMLabelGetValue(identLabel,c,&basePoint));
339566063dSJacob Faibussowitsch     if (!basePoint) PetscCall(DMLabelSetValue(*adaptLabel,c,DM_ADAPT_REFINE));
34c4762a1bSJed Brown   }
35c4762a1bSJed Brown   PetscFunctionReturn(0);
36c4762a1bSJed Brown }
37c4762a1bSJed Brown 
38c4762a1bSJed Brown static PetscErrorCode LinearFunction(PetscInt dim,PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
39c4762a1bSJed Brown {
40c4762a1bSJed Brown   PetscFunctionBeginUser;
41c4762a1bSJed Brown   u[0] = (x[0] * 2.0 + 1.) + (x[1] * 20.0 + 10.) + ((dim == 3) ? (x[2] * 200.0 + 100.) : 0.);
42c4762a1bSJed Brown   PetscFunctionReturn(0);
43c4762a1bSJed Brown }
44c4762a1bSJed Brown 
45c4762a1bSJed Brown static PetscErrorCode MultiaffineFunction(PetscInt dim,PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
46c4762a1bSJed Brown {
47c4762a1bSJed Brown   PetscFunctionBeginUser;
48c4762a1bSJed Brown   u[0] = (x[0] * 1.0 + 2.0) * (x[1] * 3.0 - 4.0) * ((dim == 3) ? (x[2] * 5.0 + 6.0) : 1.);
49c4762a1bSJed Brown   PetscFunctionReturn(0);
50c4762a1bSJed Brown }
51c4762a1bSJed Brown 
52c4762a1bSJed Brown static PetscErrorCode CoordsFunction(PetscInt dim,PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
53c4762a1bSJed Brown {
54c4762a1bSJed Brown   PetscInt f;
55c4762a1bSJed Brown 
56c4762a1bSJed Brown   PetscFunctionBeginUser;
57c4762a1bSJed Brown   for (f=0;f<Nf;f++) u[f] = x[f];
58c4762a1bSJed Brown   PetscFunctionReturn(0);
59c4762a1bSJed Brown }
60c4762a1bSJed Brown 
61c4762a1bSJed Brown typedef struct _bc_func_ctx
62c4762a1bSJed Brown {
63c4762a1bSJed Brown   PetscErrorCode (*func) (PetscInt,PetscReal,const PetscReal [], PetscInt, PetscScalar [], void *);
64c4762a1bSJed Brown   PetscInt dim;
65c4762a1bSJed Brown   PetscInt Nf;
66c4762a1bSJed Brown   void *ctx;
67c4762a1bSJed Brown }
68c4762a1bSJed Brown bc_func_ctx;
69c4762a1bSJed Brown 
70c4762a1bSJed Brown static PetscErrorCode bc_func_fv (PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
71c4762a1bSJed Brown {
72c4762a1bSJed Brown   bc_func_ctx    *bcCtx;
73c4762a1bSJed Brown 
74c4762a1bSJed Brown   PetscFunctionBegin;
75c4762a1bSJed Brown   bcCtx = (bc_func_ctx *) ctx;
769566063dSJacob Faibussowitsch   PetscCall((bcCtx->func)(bcCtx->dim,time,c,bcCtx->Nf,xG,bcCtx->ctx));
77c4762a1bSJed Brown   PetscFunctionReturn(0);
78c4762a1bSJed Brown }
79c4762a1bSJed Brown 
80c4762a1bSJed Brown static PetscErrorCode IdentifyBadPoints (DM dm, Vec vec, PetscReal tol)
81c4762a1bSJed Brown {
82c4762a1bSJed Brown   DM             dmplex;
83c4762a1bSJed Brown   PetscInt       p, pStart, pEnd, maxDof;
84c4762a1bSJed Brown   Vec            vecLocal;
85c4762a1bSJed Brown   DMLabel        depthLabel;
86c4762a1bSJed Brown   PetscSection   section;
87c4762a1bSJed Brown 
88c4762a1bSJed Brown   PetscFunctionBegin;
899566063dSJacob Faibussowitsch   PetscCall(DMCreateLocalVector(dm, &vecLocal));
909566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal));
919566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal));
929566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm ,DMPLEX, &dmplex));
939566063dSJacob Faibussowitsch   PetscCall(DMPlexGetChart(dmplex, &pStart, &pEnd));
949566063dSJacob Faibussowitsch   PetscCall(DMPlexGetDepthLabel(dmplex, &depthLabel));
959566063dSJacob Faibussowitsch   PetscCall(DMGetLocalSection(dmplex, &section));
969566063dSJacob Faibussowitsch   PetscCall(PetscSectionGetMaxDof(section, &maxDof));
97c4762a1bSJed Brown   for (p = pStart; p < pEnd; p++) {
98c4762a1bSJed Brown     PetscInt     s, c, cSize, parent, childID, numChildren;
99c4762a1bSJed Brown     PetscInt     cl, closureSize, *closure = NULL;
100c4762a1bSJed Brown     PetscScalar *values = NULL;
101c4762a1bSJed Brown     PetscBool    bad = PETSC_FALSE;
102c4762a1bSJed Brown 
1039566063dSJacob Faibussowitsch     PetscCall(VecGetValuesSection(vecLocal, section, p, &values));
1049566063dSJacob Faibussowitsch     PetscCall(PetscSectionGetDof(section, p, &cSize));
105c4762a1bSJed Brown     for (c = 0; c < cSize; c++) {
1062f613bf5SBarry Smith       PetscReal absDiff = PetscAbsScalar(values[c]);
107c4762a1bSJed Brown       if (absDiff > tol) {bad = PETSC_TRUE; break;}
108c4762a1bSJed Brown     }
109c4762a1bSJed Brown     if (!bad) continue;
11063a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Bad point %" PetscInt_FMT "\n", p));
1119566063dSJacob Faibussowitsch     PetscCall(DMLabelGetValue(depthLabel, p, &s));
11263a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Depth %" PetscInt_FMT "\n", s));
1139566063dSJacob Faibussowitsch     PetscCall(DMPlexGetTransitiveClosure(dmplex, p, PETSC_TRUE, &closureSize, &closure));
114c4762a1bSJed Brown     for (cl = 0; cl < closureSize; cl++) {
115c4762a1bSJed Brown       PetscInt cp = closure[2 * cl];
1169566063dSJacob Faibussowitsch       PetscCall(DMPlexGetTreeParent(dmplex, cp, &parent, &childID));
117c4762a1bSJed Brown       if (parent != cp) {
11863a3b9bcSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Closure point %" PetscInt_FMT " (%" PetscInt_FMT ") child of %" PetscInt_FMT " (ID %" PetscInt_FMT ")\n", cl, cp, parent, childID));
119c4762a1bSJed Brown       }
1209566063dSJacob Faibussowitsch       PetscCall(DMPlexGetTreeChildren(dmplex, cp, &numChildren, NULL));
121c4762a1bSJed Brown       if (numChildren) {
12263a3b9bcSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Closure point %" PetscInt_FMT " (%" PetscInt_FMT ") is parent\n", cl, cp));
123c4762a1bSJed Brown       }
124c4762a1bSJed Brown     }
1259566063dSJacob Faibussowitsch     PetscCall(DMPlexRestoreTransitiveClosure(dmplex, p, PETSC_TRUE, &closureSize, &closure));
126c4762a1bSJed Brown     for (c = 0; c < cSize; c++) {
1272f613bf5SBarry Smith       PetscReal absDiff = PetscAbsScalar(values[c]);
128c4762a1bSJed Brown       if (absDiff > tol) {
12963a3b9bcSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Bad dof %" PetscInt_FMT "\n", c));
130c4762a1bSJed Brown       }
131c4762a1bSJed Brown     }
132c4762a1bSJed Brown   }
1339566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dmplex));
1349566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&vecLocal));
135c4762a1bSJed Brown   PetscFunctionReturn(0);
136c4762a1bSJed Brown }
137c4762a1bSJed Brown 
138c4762a1bSJed Brown int main(int argc, char **argv)
139c4762a1bSJed Brown {
140c4762a1bSJed Brown   MPI_Comm       comm;
141c4762a1bSJed Brown   DM             base, preForest, postForest;
14230602db0SMatthew G. Knepley   PetscInt       dim, Nf = 1;
143c4762a1bSJed Brown   PetscInt       step, adaptSteps = 1;
144c4762a1bSJed Brown   PetscInt       preCount, postCount;
145c4762a1bSJed Brown   Vec            preVec, postVecTransfer, postVecExact;
146c4762a1bSJed Brown   PetscErrorCode (*funcs[1]) (PetscInt,PetscReal,const PetscReal [],PetscInt,PetscScalar [], void *) = {MultiaffineFunction};
147c4762a1bSJed Brown   void           *ctxs[1] = {NULL};
148c4762a1bSJed Brown   PetscReal      diff, tol = PETSC_SMALL;
149c4762a1bSJed Brown   PetscBool      linear = PETSC_FALSE;
150c4762a1bSJed Brown   PetscBool      coords = PETSC_FALSE;
151c4762a1bSJed Brown   PetscBool      useFV = PETSC_FALSE;
152c4762a1bSJed Brown   PetscBool      conv = PETSC_FALSE;
153c4762a1bSJed Brown   PetscBool      transfer_from_base[2] = {PETSC_TRUE,PETSC_FALSE};
154c4762a1bSJed Brown   PetscBool      use_bcs = PETSC_TRUE;
155c4762a1bSJed Brown   bc_func_ctx    bcCtx;
156c4762a1bSJed Brown   DMLabel        adaptLabel;
157c4762a1bSJed Brown 
1589566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &argv, NULL,help));
159c4762a1bSJed Brown   comm = PETSC_COMM_WORLD;
160d0609cedSBarry Smith   PetscOptionsBegin(comm, "", "DMForestTransferVec() Test Options", "DMFOREST");
1619566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-linear","Transfer a simple linear function", "ex2.c", linear, &linear, NULL));
1629566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-coords","Transfer a simple coordinate function", "ex2.c", coords, &coords, NULL));
1639566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-use_fv","Use a finite volume approximation", "ex2.c", useFV, &useFV, NULL));
1649566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-test_convert","Test conversion to DMPLEX",NULL,conv,&conv,NULL));
1659566063dSJacob 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));
166c4762a1bSJed Brown   transfer_from_base[1] = transfer_from_base[0];
1679566063dSJacob 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));
1689566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBool("-use_bcs","Use dirichlet boundary conditions", "ex2.c", use_bcs, &use_bcs, NULL));
1699566063dSJacob Faibussowitsch   PetscCall(PetscOptionsBoundedInt("-adapt_steps","Number of adaptivity steps", "ex2.c", adaptSteps, &adaptSteps, NULL,0));
170d0609cedSBarry Smith   PetscOptionsEnd();
171c4762a1bSJed Brown 
172c4762a1bSJed Brown   tol = PetscMax(1.e-10,tol); /* XXX fix for quadruple precision -> why do I need to do this? */
173c4762a1bSJed Brown 
17430602db0SMatthew G. Knepley   /* the base mesh */
1759566063dSJacob Faibussowitsch   PetscCall(DMCreate(comm, &base));
1769566063dSJacob Faibussowitsch   PetscCall(DMSetType(base, DMPLEX));
1779566063dSJacob Faibussowitsch   PetscCall(DMSetFromOptions(base));
17830602db0SMatthew G. Knepley 
1799566063dSJacob Faibussowitsch   PetscCall(AddIdentityLabel(base));
1809566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(base, &dim));
181c4762a1bSJed Brown 
182c4762a1bSJed Brown   if (linear) {
183c4762a1bSJed Brown     funcs[0] = LinearFunction;
184c4762a1bSJed Brown   }
185c4762a1bSJed Brown   if (coords) {
186c4762a1bSJed Brown     funcs[0] = CoordsFunction;
187c4762a1bSJed Brown     Nf = dim;
188c4762a1bSJed Brown   }
189c4762a1bSJed Brown 
190c4762a1bSJed Brown   bcCtx.func = funcs[0];
191c4762a1bSJed Brown   bcCtx.dim  = dim;
192c4762a1bSJed Brown   bcCtx.Nf   = Nf;
193c4762a1bSJed Brown   bcCtx.ctx  = NULL;
194c4762a1bSJed Brown 
195c4762a1bSJed Brown   if (useFV) {
196c4762a1bSJed Brown     PetscFV      fv;
197c4762a1bSJed Brown     PetscLimiter limiter;
198c4762a1bSJed Brown     DM           baseFV;
199c4762a1bSJed Brown 
2009566063dSJacob Faibussowitsch     PetscCall(DMPlexConstructGhostCells(base,NULL,NULL,&baseFV));
2019566063dSJacob Faibussowitsch     PetscCall(DMViewFromOptions(baseFV, NULL, "-fv_dm_view"));
2029566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&base));
203c4762a1bSJed Brown     base = baseFV;
204*40cbb1a0SMatthew G. Knepley     PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv));
2059566063dSJacob Faibussowitsch     PetscCall(PetscFVSetSpatialDimension(fv,dim));
2069566063dSJacob Faibussowitsch     PetscCall(PetscFVSetType(fv,PETSCFVLEASTSQUARES));
2079566063dSJacob Faibussowitsch     PetscCall(PetscFVSetNumComponents(fv,Nf));
2089566063dSJacob Faibussowitsch     PetscCall(PetscLimiterCreate(comm,&limiter));
2099566063dSJacob Faibussowitsch     PetscCall(PetscLimiterSetType(limiter,PETSCLIMITERNONE));
2109566063dSJacob Faibussowitsch     PetscCall(PetscFVSetLimiter(fv,limiter));
2119566063dSJacob Faibussowitsch     PetscCall(PetscLimiterDestroy(&limiter));
2129566063dSJacob Faibussowitsch     PetscCall(PetscFVSetFromOptions(fv));
2139566063dSJacob Faibussowitsch     PetscCall(DMSetField(base,0,NULL,(PetscObject)fv));
2149566063dSJacob Faibussowitsch     PetscCall(PetscFVDestroy(&fv));
215c4762a1bSJed Brown   } else {
216c4762a1bSJed Brown     PetscFE fe;
217c4762a1bSJed Brown 
2189566063dSJacob Faibussowitsch     PetscCall(PetscFECreateDefault(comm,dim,Nf,PETSC_FALSE,NULL,PETSC_DEFAULT,&fe));
2199566063dSJacob Faibussowitsch     PetscCall(DMSetField(base,0,NULL,(PetscObject)fe));
2209566063dSJacob Faibussowitsch     PetscCall(PetscFEDestroy(&fe));
221c4762a1bSJed Brown   }
2229566063dSJacob Faibussowitsch   PetscCall(DMCreateDS(base));
223c4762a1bSJed Brown 
224c4762a1bSJed Brown   if (use_bcs) {
225c4762a1bSJed Brown     PetscInt ids[] = {1, 2, 3, 4, 5, 6};
22645480ffeSMatthew G. Knepley     DMLabel  label;
227c4762a1bSJed Brown 
2289566063dSJacob Faibussowitsch     PetscCall(DMGetLabel(base, "marker", &label));
2299566063dSJacob 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));
230c4762a1bSJed Brown   }
2319566063dSJacob Faibussowitsch   PetscCall(DMViewFromOptions(base,NULL,"-dm_base_view"));
232c4762a1bSJed Brown 
233c4762a1bSJed Brown   /* the pre adaptivity forest */
2349566063dSJacob Faibussowitsch   PetscCall(DMCreate(comm,&preForest));
2359566063dSJacob Faibussowitsch   PetscCall(DMSetType(preForest,(dim == 2) ? DMP4EST : DMP8EST));
2369566063dSJacob Faibussowitsch   PetscCall(DMCopyDisc(base,preForest));
2379566063dSJacob Faibussowitsch   PetscCall(DMForestSetBaseDM(preForest,base));
2389566063dSJacob Faibussowitsch   PetscCall(DMForestSetMinimumRefinement(preForest,0));
2399566063dSJacob Faibussowitsch   PetscCall(DMForestSetInitialRefinement(preForest,1));
2409566063dSJacob Faibussowitsch   PetscCall(DMSetFromOptions(preForest));
2419566063dSJacob Faibussowitsch   PetscCall(DMSetUp(preForest));
2429566063dSJacob Faibussowitsch   PetscCall(DMViewFromOptions(preForest,NULL,"-dm_pre_view"));
243c4762a1bSJed Brown 
244c4762a1bSJed Brown   /* the pre adaptivity field */
2459566063dSJacob Faibussowitsch   PetscCall(DMCreateGlobalVector(preForest,&preVec));
2469566063dSJacob Faibussowitsch   PetscCall(DMProjectFunction(preForest,0.,funcs,ctxs,INSERT_VALUES,preVec));
2479566063dSJacob Faibussowitsch   PetscCall(VecViewFromOptions(preVec,NULL,"-vec_pre_view"));
248c4762a1bSJed Brown 
249c4762a1bSJed Brown   /* communicate between base and pre adaptivity forest */
250c4762a1bSJed Brown   if (transfer_from_base[0]) {
251c4762a1bSJed Brown     Vec baseVec, baseVecMapped;
252c4762a1bSJed Brown 
2539566063dSJacob Faibussowitsch     PetscCall(DMGetGlobalVector(base,&baseVec));
2549566063dSJacob Faibussowitsch     PetscCall(DMProjectFunction(base,0.,funcs,ctxs,INSERT_VALUES,baseVec));
2559566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetName((PetscObject)baseVec,"Function Base"));
2569566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(baseVec,NULL,"-vec_base_view"));
257c4762a1bSJed Brown 
2589566063dSJacob Faibussowitsch     PetscCall(DMGetGlobalVector(preForest,&baseVecMapped));
2599566063dSJacob Faibussowitsch     PetscCall(DMForestTransferVecFromBase(preForest,baseVec,baseVecMapped));
2609566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(baseVecMapped,NULL,"-vec_map_base_view"));
261c4762a1bSJed Brown 
262c4762a1bSJed Brown     /* compare */
2639566063dSJacob Faibussowitsch     PetscCall(VecAXPY(baseVecMapped,-1.,preVec));
2649566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(baseVecMapped,NULL,"-vec_map_diff_view"));
2659566063dSJacob Faibussowitsch     PetscCall(VecNorm(baseVecMapped,NORM_2,&diff));
266c4762a1bSJed Brown 
267c4762a1bSJed Brown     /* output */
268c4762a1bSJed Brown     if (diff < tol) {
2699566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(comm,"DMForestTransferVecFromBase() passes.\n"));
270c4762a1bSJed Brown     } else {
2719566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(comm,"DMForestTransferVecFromBase() fails with error %g and tolerance %g\n",(double)diff,(double)tol));
272c4762a1bSJed Brown     }
273c4762a1bSJed Brown 
2749566063dSJacob Faibussowitsch     PetscCall(DMRestoreGlobalVector(base,&baseVec));
2759566063dSJacob Faibussowitsch     PetscCall(DMRestoreGlobalVector(preForest,&baseVecMapped));
276c4762a1bSJed Brown   }
277c4762a1bSJed Brown 
278c4762a1bSJed Brown   for (step = 0; step < adaptSteps; ++step) {
279c4762a1bSJed Brown 
280c4762a1bSJed Brown     if (!transfer_from_base[1]) {
2819566063dSJacob Faibussowitsch       PetscCall(PetscObjectGetReference((PetscObject)preForest,&preCount));
282c4762a1bSJed Brown     }
283c4762a1bSJed Brown 
284c4762a1bSJed Brown     /* adapt */
2859566063dSJacob Faibussowitsch     PetscCall(CreateAdaptivityLabel(preForest,&adaptLabel));
2869566063dSJacob Faibussowitsch     PetscCall(DMForestTemplate(preForest,comm,&postForest));
2879566063dSJacob Faibussowitsch     if (step) PetscCall(DMForestSetAdaptivityLabel(postForest,adaptLabel));
2889566063dSJacob Faibussowitsch     PetscCall(DMLabelDestroy(&adaptLabel));
2899566063dSJacob Faibussowitsch     PetscCall(DMSetUp(postForest));
2909566063dSJacob Faibussowitsch     PetscCall(DMViewFromOptions(postForest,NULL,"-dm_post_view"));
291c4762a1bSJed Brown 
292c4762a1bSJed Brown     /* transfer */
2939566063dSJacob Faibussowitsch     PetscCall(DMCreateGlobalVector(postForest,&postVecTransfer));
2949566063dSJacob Faibussowitsch     PetscCall(DMForestTransferVec(preForest,preVec,postForest,postVecTransfer,PETSC_TRUE,0.0));
2959566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(postVecTransfer,NULL,"-vec_post_transfer_view"));
296c4762a1bSJed Brown 
297c4762a1bSJed Brown     /* the exact post adaptivity field */
2989566063dSJacob Faibussowitsch     PetscCall(DMCreateGlobalVector(postForest,&postVecExact));
2999566063dSJacob Faibussowitsch     PetscCall(DMProjectFunction(postForest,0.,funcs,ctxs,INSERT_VALUES,postVecExact));
3009566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(postVecExact,NULL,"-vec_post_exact_view"));
301c4762a1bSJed Brown 
302c4762a1bSJed Brown     /* compare */
3039566063dSJacob Faibussowitsch     PetscCall(VecAXPY(postVecExact,-1.,postVecTransfer));
3049566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(postVecExact,NULL,"-vec_diff_view"));
3059566063dSJacob Faibussowitsch     PetscCall(VecNorm(postVecExact,NORM_2,&diff));
306c4762a1bSJed Brown 
307c4762a1bSJed Brown     /* output */
308c4762a1bSJed Brown     if (diff < tol) {
3099566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(comm,"DMForestTransferVec() passes.\n"));
310c4762a1bSJed Brown     } else {
3119566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(comm,"DMForestTransferVec() fails with error %g and tolerance %g\n",(double)diff,(double)tol));
3129566063dSJacob Faibussowitsch       PetscCall(IdentifyBadPoints(postForest, postVecExact, tol));
313c4762a1bSJed Brown     }
3149566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&postVecExact));
315c4762a1bSJed Brown 
316c4762a1bSJed Brown     /* disconnect preForest from postForest if we don't test the transfer throughout the entire refinement process */
317c4762a1bSJed Brown     if (!transfer_from_base[1]) {
3189566063dSJacob Faibussowitsch       PetscCall(DMForestSetAdaptivityForest(postForest,NULL));
3199566063dSJacob Faibussowitsch       PetscCall(PetscObjectGetReference((PetscObject)preForest,&postCount));
32063a3b9bcSJacob 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);
321c4762a1bSJed Brown     }
322c4762a1bSJed Brown 
323c4762a1bSJed Brown     if (conv) {
324c4762a1bSJed Brown       DM dmConv;
325c4762a1bSJed Brown 
3269566063dSJacob Faibussowitsch       PetscCall(DMConvert(postForest,DMPLEX,&dmConv));
3279566063dSJacob Faibussowitsch       PetscCall(DMViewFromOptions(dmConv,NULL,"-dm_conv_view"));
3289566063dSJacob Faibussowitsch       PetscCall(DMPlexCheckCellShape(dmConv,PETSC_TRUE,PETSC_DETERMINE));
3299566063dSJacob Faibussowitsch       PetscCall(DMDestroy(&dmConv));
330c4762a1bSJed Brown     }
331c4762a1bSJed Brown 
3329566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&preVec));
3339566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&preForest));
334c4762a1bSJed Brown 
335c4762a1bSJed Brown     preVec    = postVecTransfer;
336c4762a1bSJed Brown     preForest = postForest;
337c4762a1bSJed Brown   }
338c4762a1bSJed Brown 
339c4762a1bSJed Brown   if (transfer_from_base[1]) {
340c4762a1bSJed Brown     Vec baseVec, baseVecMapped;
341c4762a1bSJed Brown 
342c4762a1bSJed Brown     /* communicate between base and last adapted forest */
3439566063dSJacob Faibussowitsch     PetscCall(DMGetGlobalVector(base,&baseVec));
3449566063dSJacob Faibussowitsch     PetscCall(DMProjectFunction(base,0.,funcs,ctxs,INSERT_VALUES,baseVec));
3459566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetName((PetscObject)baseVec,"Function Base"));
3469566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(baseVec,NULL,"-vec_base_view"));
347c4762a1bSJed Brown 
3489566063dSJacob Faibussowitsch     PetscCall(DMGetGlobalVector(preForest,&baseVecMapped));
3499566063dSJacob Faibussowitsch     PetscCall(DMForestTransferVecFromBase(preForest,baseVec,baseVecMapped));
3509566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(baseVecMapped,NULL,"-vec_map_base_view"));
351c4762a1bSJed Brown 
352c4762a1bSJed Brown     /* compare */
3539566063dSJacob Faibussowitsch     PetscCall(VecAXPY(baseVecMapped,-1.,preVec));
3549566063dSJacob Faibussowitsch     PetscCall(VecViewFromOptions(baseVecMapped,NULL,"-vec_map_diff_view"));
3559566063dSJacob Faibussowitsch     PetscCall(VecNorm(baseVecMapped,NORM_2,&diff));
356c4762a1bSJed Brown 
357c4762a1bSJed Brown     /* output */
358c4762a1bSJed Brown     if (diff < tol) {
3599566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(comm,"DMForestTransferVecFromBase() passes.\n"));
360c4762a1bSJed Brown     } else {
3619566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(comm,"DMForestTransferVecFromBase() fails with error %g and tolerance %g\n",(double)diff,(double)tol));
362c4762a1bSJed Brown     }
363c4762a1bSJed Brown 
3649566063dSJacob Faibussowitsch     PetscCall(DMRestoreGlobalVector(base,&baseVec));
3659566063dSJacob Faibussowitsch     PetscCall(DMRestoreGlobalVector(preForest,&baseVecMapped));
366c4762a1bSJed Brown   }
367c4762a1bSJed Brown 
368c4762a1bSJed Brown   /* cleanup */
3699566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&preVec));
3709566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&preForest));
3719566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&base));
3729566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
373b122ec5aSJacob Faibussowitsch   return 0;
374c4762a1bSJed Brown }
375c4762a1bSJed Brown 
376c4762a1bSJed Brown /*TEST
37730602db0SMatthew G. Knepley   testset:
37830602db0SMatthew G. Knepley     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3 -petscspace_type tensor
379c4762a1bSJed Brown 
380c4762a1bSJed Brown     test:
381c4762a1bSJed Brown       output_file: output/ex2_2d.out
382c4762a1bSJed Brown       suffix: p4est_2d
38330602db0SMatthew G. Knepley       args: -petscspace_degree 2
384c4762a1bSJed Brown       nsize: 3
385c4762a1bSJed Brown       requires: p4est !single
386c4762a1bSJed Brown 
387c4762a1bSJed Brown     test:
388c4762a1bSJed Brown       output_file: output/ex2_2d.out
389c4762a1bSJed Brown       suffix: p4est_2d_deg4
39030602db0SMatthew G. Knepley       args: -petscspace_degree 4
391c4762a1bSJed Brown       requires: p4est !single
392c4762a1bSJed Brown 
393c4762a1bSJed Brown     test:
394c4762a1bSJed Brown       output_file: output/ex2_2d.out
395c4762a1bSJed Brown       suffix: p4est_2d_deg8
39630602db0SMatthew G. Knepley       args: -petscspace_degree 8
397c4762a1bSJed Brown       requires: p4est !single
398c4762a1bSJed Brown 
399c4762a1bSJed Brown     test:
400c4762a1bSJed Brown       output_file: output/ex2_steps2.out
401c4762a1bSJed Brown       suffix: p4est_2d_deg2_steps2
40230602db0SMatthew G. Knepley       args: -petscspace_degree 2 -coords -adapt_steps 2
403c4762a1bSJed Brown       nsize: 3
404c4762a1bSJed Brown       requires: p4est !single
405c4762a1bSJed Brown 
406c4762a1bSJed Brown     test:
407c4762a1bSJed Brown       output_file: output/ex2_steps3.out
408c4762a1bSJed Brown       suffix: p4est_2d_deg3_steps3
40930602db0SMatthew G. Knepley       args: -petscspace_degree 3 -coords -adapt_steps 3 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1
410c4762a1bSJed Brown       nsize: 3
411c4762a1bSJed Brown       requires: p4est !single
412c4762a1bSJed Brown 
413c4762a1bSJed Brown     test:
414c4762a1bSJed Brown       output_file: output/ex2_steps3.out
415c4762a1bSJed Brown       suffix: p4est_2d_deg3_steps3_L2_periodic
41630602db0SMatthew 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
417c4762a1bSJed Brown       nsize: 3
418c4762a1bSJed Brown       requires: p4est !single
419c4762a1bSJed Brown 
420c4762a1bSJed Brown     test:
421c4762a1bSJed Brown       output_file: output/ex2_steps3.out
422c4762a1bSJed Brown       suffix: p4est_3d_deg2_steps3_L2_periodic
42330602db0SMatthew 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
424c4762a1bSJed Brown       nsize: 3
425c4762a1bSJed Brown       requires: p4est !single
426c4762a1bSJed Brown 
427c4762a1bSJed Brown     test:
428c4762a1bSJed Brown       output_file: output/ex2_steps2.out
429c4762a1bSJed Brown       suffix: p4est_3d_deg2_steps2
43030602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 2 -coords -adapt_steps 2
431c4762a1bSJed Brown       nsize: 3
432c4762a1bSJed Brown       requires: p4est !single
433c4762a1bSJed Brown 
434c4762a1bSJed Brown     test:
435c4762a1bSJed Brown       output_file: output/ex2_steps3.out
436c4762a1bSJed Brown       suffix: p4est_3d_deg3_steps3
43730602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 3 -coords -adapt_steps 3 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1
438c4762a1bSJed Brown       nsize: 3
439c4762a1bSJed Brown       requires: p4est !single
440c4762a1bSJed Brown 
441c4762a1bSJed Brown     test:
442c4762a1bSJed Brown       output_file: output/ex2_3d.out
443c4762a1bSJed Brown       suffix: p4est_3d
44430602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 1
445c4762a1bSJed Brown       nsize: 3
446c4762a1bSJed Brown       requires: p4est !single
447c4762a1bSJed Brown 
448c4762a1bSJed Brown     test:
449c4762a1bSJed Brown       output_file: output/ex2_3d.out
450c4762a1bSJed Brown       suffix: p4est_3d_deg3
45130602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 3
452c4762a1bSJed Brown       nsize: 3
453c4762a1bSJed Brown       requires: p4est !single
454c4762a1bSJed Brown 
455c4762a1bSJed Brown     test:
456c4762a1bSJed Brown       output_file: output/ex2_2d.out
457c4762a1bSJed Brown       suffix: p4est_2d_deg2_coords
45830602db0SMatthew G. Knepley       args: -petscspace_degree 2 -coords
459c4762a1bSJed Brown       nsize: 3
460c4762a1bSJed Brown       requires: p4est !single
461c4762a1bSJed Brown 
462c4762a1bSJed Brown     test:
463c4762a1bSJed Brown       output_file: output/ex2_3d.out
464c4762a1bSJed Brown       suffix: p4est_3d_deg2_coords
46530602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -petscspace_degree 2 -coords
466c4762a1bSJed Brown       nsize: 3
467c4762a1bSJed Brown       requires: p4est !single
468c4762a1bSJed Brown 
469c4762a1bSJed Brown     test:
470c4762a1bSJed Brown       suffix: p4est_3d_nans
47130602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -dm_forest_partition_overlap 1 -test_convert -petscspace_degree 1
472c4762a1bSJed Brown       nsize: 2
473c4762a1bSJed Brown       requires: p4est !single
474c4762a1bSJed Brown 
475c4762a1bSJed Brown     test:
476c4762a1bSJed Brown       TODO: not broken, but the 3D case below is broken, so I do not trust this one
477c4762a1bSJed Brown       output_file: output/ex2_steps2.out
478c4762a1bSJed Brown       suffix: p4est_2d_tfb_distributed_nc
479e600fa54SMatthew 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
480c4762a1bSJed Brown       nsize: 3
481c4762a1bSJed Brown       requires: p4est !single
482c4762a1bSJed Brown 
483c4762a1bSJed Brown     test:
484c4762a1bSJed Brown       TODO: broken
485c4762a1bSJed Brown       output_file: output/ex2_steps2.out
486c4762a1bSJed Brown       suffix: p4est_3d_tfb_distributed_nc
487e600fa54SMatthew 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
488c4762a1bSJed Brown       nsize: 3
489c4762a1bSJed Brown       requires: p4est !single
490c4762a1bSJed Brown 
49130602db0SMatthew G. Knepley   testset:
492012bc364SMatthew G. Knepley     args: -petscspace_type tensor -dm_coord_space 0 -dm_plex_transform_type refine_tobox
49330602db0SMatthew G. Knepley 
494c4762a1bSJed Brown     test:
495c4762a1bSJed Brown       TODO: broken
496c4762a1bSJed Brown       output_file: output/ex2_3d.out
497c4762a1bSJed Brown       suffix: p4est_3d_transfer_fails
498012bc364SMatthew 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
499c4762a1bSJed Brown       requires: p4est !single
500c4762a1bSJed Brown 
501c4762a1bSJed Brown     test:
502c4762a1bSJed Brown       TODO: broken
503c4762a1bSJed Brown       output_file: output/ex2_steps2_notfb.out
504c4762a1bSJed Brown       suffix: p4est_3d_transfer_fails_2
505012bc364SMatthew 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
506c4762a1bSJed Brown       requires: p4est !single
507c4762a1bSJed Brown 
508c4762a1bSJed Brown     test:
509c4762a1bSJed Brown       output_file: output/ex2_steps2.out
510c4762a1bSJed Brown       suffix: p4est_3d_multi_transfer_s2t
511012bc364SMatthew 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
512012bc364SMatthew G. Knepley       requires: p4est !single
513c4762a1bSJed Brown 
514c4762a1bSJed Brown     test:
515c4762a1bSJed Brown       output_file: output/ex2_steps2.out
516c4762a1bSJed Brown       suffix: p4est_3d_coords_transfer_s2t
517012bc364SMatthew 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
518012bc364SMatthew G. Knepley       requires: p4est !single
51930602db0SMatthew G. Knepley 
52030602db0SMatthew G. Knepley   testset:
52130602db0SMatthew G. Knepley     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3
52230602db0SMatthew G. Knepley 
52330602db0SMatthew G. Knepley     test:
52430602db0SMatthew G. Knepley       output_file: output/ex2_2d_fv.out
52530602db0SMatthew G. Knepley       suffix: p4est_2d_fv
52630602db0SMatthew G. Knepley       args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1
52730602db0SMatthew G. Knepley       nsize: 3
52830602db0SMatthew G. Knepley       requires: p4est !single
52930602db0SMatthew G. Knepley 
53030602db0SMatthew G. Knepley     test:
53130602db0SMatthew G. Knepley       TODO: broken (codimension adjacency)
53230602db0SMatthew G. Knepley       output_file: output/ex2_2d_fv.out
53330602db0SMatthew G. Knepley       suffix: p4est_2d_fv_adjcodim
53430602db0SMatthew G. Knepley       args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 -dm_forest_adjacency_codimension 1
53530602db0SMatthew G. Knepley       nsize: 2
53630602db0SMatthew G. Knepley       requires: p4est !single
53730602db0SMatthew G. Knepley 
53830602db0SMatthew G. Knepley     test:
53930602db0SMatthew G. Knepley       TODO: broken (dimension adjacency)
54030602db0SMatthew G. Knepley       output_file: output/ex2_2d_fv.out
54130602db0SMatthew G. Knepley       suffix: p4est_2d_fv_adjdim
54230602db0SMatthew G. Knepley       args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 -dm_forest_adjacency_dimension 1
54330602db0SMatthew G. Knepley       nsize: 2
54430602db0SMatthew G. Knepley       requires: p4est !single
54530602db0SMatthew G. Knepley 
54630602db0SMatthew G. Knepley     test:
54730602db0SMatthew G. Knepley       output_file: output/ex2_2d_fv.out
54830602db0SMatthew G. Knepley       suffix: p4est_2d_fv_zerocells
54930602db0SMatthew G. Knepley       args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1
55030602db0SMatthew G. Knepley       nsize: 10
55130602db0SMatthew G. Knepley       requires: p4est !single
55230602db0SMatthew G. Knepley 
55330602db0SMatthew G. Knepley     test:
55430602db0SMatthew G. Knepley       output_file: output/ex2_3d_fv.out
55530602db0SMatthew G. Knepley       suffix: p4est_3d_fv
55630602db0SMatthew G. Knepley       args: -dm_plex_dim 3 -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1
55730602db0SMatthew G. Knepley       nsize: 3
558c4762a1bSJed Brown       requires: p4est !single
559c4762a1bSJed Brown 
560c4762a1bSJed Brown TEST*/
561