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
AddIdentityLabel(DM dm)8d71ae5a4SJacob Faibussowitsch static PetscErrorCode AddIdentityLabel(DM dm)
9d71ae5a4SJacob Faibussowitsch {
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));
163ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
17c4762a1bSJed Brown }
18c4762a1bSJed Brown
CreateAdaptivityLabel(DM forest,DMLabel * adaptLabel)19d71ae5a4SJacob Faibussowitsch static PetscErrorCode CreateAdaptivityLabel(DM forest, DMLabel *adaptLabel)
20d71ae5a4SJacob Faibussowitsch {
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 }
353ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
36c4762a1bSJed Brown }
37c4762a1bSJed Brown
LinearFunction(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar u[],PetscCtx ctx)38*2a8381b2SBarry Smith static PetscErrorCode LinearFunction(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], PetscCtx ctx)
39d71ae5a4SJacob Faibussowitsch {
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.);
423ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
43c4762a1bSJed Brown }
44c4762a1bSJed Brown
MultiaffineFunction(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar u[],PetscCtx ctx)45*2a8381b2SBarry Smith static PetscErrorCode MultiaffineFunction(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], PetscCtx ctx)
46d71ae5a4SJacob Faibussowitsch {
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.);
493ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
50c4762a1bSJed Brown }
51c4762a1bSJed Brown
CoordsFunction(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar u[],PetscCtx ctx)52*2a8381b2SBarry Smith static PetscErrorCode CoordsFunction(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], PetscCtx ctx)
53d71ae5a4SJacob Faibussowitsch {
54c4762a1bSJed Brown PetscInt f;
55c4762a1bSJed Brown
56c4762a1bSJed Brown PetscFunctionBeginUser;
57c4762a1bSJed Brown for (f = 0; f < Nf; f++) u[f] = x[f];
583ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
59c4762a1bSJed Brown }
60c4762a1bSJed Brown
619371c9d4SSatish Balay typedef struct _bc_func_ctx {
62c4762a1bSJed Brown PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
63c4762a1bSJed Brown PetscInt dim;
64c4762a1bSJed Brown PetscInt Nf;
65c4762a1bSJed Brown void *ctx;
669371c9d4SSatish Balay } bc_func_ctx;
67c4762a1bSJed Brown
bc_func_fv(PetscReal time,const PetscReal * c,const PetscReal * n,const PetscScalar * xI,PetscScalar * xG,PetscCtx ctx)68*2a8381b2SBarry Smith static PetscErrorCode bc_func_fv(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, PetscCtx ctx)
69d71ae5a4SJacob Faibussowitsch {
70c4762a1bSJed Brown bc_func_ctx *bcCtx;
71c4762a1bSJed Brown
72c4762a1bSJed Brown PetscFunctionBegin;
73c4762a1bSJed Brown bcCtx = (bc_func_ctx *)ctx;
74f4f49eeaSPierre Jolivet PetscCall(bcCtx->func(bcCtx->dim, time, c, bcCtx->Nf, xG, bcCtx->ctx));
753ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
76c4762a1bSJed Brown }
77c4762a1bSJed Brown
IdentifyBadPoints(DM dm,Vec vec,PetscReal tol)78d71ae5a4SJacob Faibussowitsch static PetscErrorCode IdentifyBadPoints(DM dm, Vec vec, PetscReal tol)
79d71ae5a4SJacob Faibussowitsch {
80c4762a1bSJed Brown DM dmplex;
81c4762a1bSJed Brown PetscInt p, pStart, pEnd, maxDof;
82c4762a1bSJed Brown Vec vecLocal;
83c4762a1bSJed Brown DMLabel depthLabel;
84c4762a1bSJed Brown PetscSection section;
85c4762a1bSJed Brown
86c4762a1bSJed Brown PetscFunctionBegin;
879566063dSJacob Faibussowitsch PetscCall(DMCreateLocalVector(dm, &vecLocal));
889566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal));
899566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal));
909566063dSJacob Faibussowitsch PetscCall(DMConvert(dm, DMPLEX, &dmplex));
919566063dSJacob Faibussowitsch PetscCall(DMPlexGetChart(dmplex, &pStart, &pEnd));
929566063dSJacob Faibussowitsch PetscCall(DMPlexGetDepthLabel(dmplex, &depthLabel));
939566063dSJacob Faibussowitsch PetscCall(DMGetLocalSection(dmplex, §ion));
949566063dSJacob Faibussowitsch PetscCall(PetscSectionGetMaxDof(section, &maxDof));
95c4762a1bSJed Brown for (p = pStart; p < pEnd; p++) {
96c4762a1bSJed Brown PetscInt s, c, cSize, parent, childID, numChildren;
97c4762a1bSJed Brown PetscInt cl, closureSize, *closure = NULL;
98c4762a1bSJed Brown PetscScalar *values = NULL;
99c4762a1bSJed Brown PetscBool bad = PETSC_FALSE;
100c4762a1bSJed Brown
1019566063dSJacob Faibussowitsch PetscCall(VecGetValuesSection(vecLocal, section, p, &values));
1029566063dSJacob Faibussowitsch PetscCall(PetscSectionGetDof(section, p, &cSize));
103c4762a1bSJed Brown for (c = 0; c < cSize; c++) {
1042f613bf5SBarry Smith PetscReal absDiff = PetscAbsScalar(values[c]);
1059371c9d4SSatish Balay if (absDiff > tol) {
1069371c9d4SSatish Balay bad = PETSC_TRUE;
1079371c9d4SSatish Balay break;
1089371c9d4SSatish Balay }
109c4762a1bSJed Brown }
110c4762a1bSJed Brown if (!bad) continue;
11163a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_SELF, "Bad point %" PetscInt_FMT "\n", p));
1129566063dSJacob Faibussowitsch PetscCall(DMLabelGetValue(depthLabel, p, &s));
11363a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_SELF, " Depth %" PetscInt_FMT "\n", s));
1149566063dSJacob Faibussowitsch PetscCall(DMPlexGetTransitiveClosure(dmplex, p, PETSC_TRUE, &closureSize, &closure));
115c4762a1bSJed Brown for (cl = 0; cl < closureSize; cl++) {
116c4762a1bSJed Brown PetscInt cp = closure[2 * cl];
1179566063dSJacob Faibussowitsch PetscCall(DMPlexGetTreeParent(dmplex, cp, &parent, &childID));
11848a46eb9SPierre 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));
1199566063dSJacob Faibussowitsch PetscCall(DMPlexGetTreeChildren(dmplex, cp, &numChildren, NULL));
12048a46eb9SPierre Jolivet if (numChildren) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Closure point %" PetscInt_FMT " (%" PetscInt_FMT ") is parent\n", cl, cp));
121c4762a1bSJed Brown }
1229566063dSJacob Faibussowitsch PetscCall(DMPlexRestoreTransitiveClosure(dmplex, p, PETSC_TRUE, &closureSize, &closure));
123c4762a1bSJed Brown for (c = 0; c < cSize; c++) {
1242f613bf5SBarry Smith PetscReal absDiff = PetscAbsScalar(values[c]);
12548a46eb9SPierre Jolivet if (absDiff > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Bad dof %" PetscInt_FMT "\n", c));
126c4762a1bSJed Brown }
127c4762a1bSJed Brown }
1289566063dSJacob Faibussowitsch PetscCall(DMDestroy(&dmplex));
1299566063dSJacob Faibussowitsch PetscCall(VecDestroy(&vecLocal));
1303ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
131c4762a1bSJed Brown }
132c4762a1bSJed Brown
main(int argc,char ** argv)133d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv)
134d71ae5a4SJacob Faibussowitsch {
135c4762a1bSJed Brown MPI_Comm comm;
136c4762a1bSJed Brown DM base, preForest, postForest;
13730602db0SMatthew G. Knepley PetscInt dim, Nf = 1;
138c4762a1bSJed Brown PetscInt step, adaptSteps = 1;
139c4762a1bSJed Brown PetscInt preCount, postCount;
140c4762a1bSJed Brown Vec preVec, postVecTransfer, postVecExact;
141c4762a1bSJed Brown PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *) = {MultiaffineFunction};
142c4762a1bSJed Brown void *ctxs[1] = {NULL};
143c4762a1bSJed Brown PetscReal diff, tol = PETSC_SMALL;
144c4762a1bSJed Brown PetscBool linear = PETSC_FALSE;
145c4762a1bSJed Brown PetscBool coords = PETSC_FALSE;
146c4762a1bSJed Brown PetscBool useFV = PETSC_FALSE;
147c4762a1bSJed Brown PetscBool conv = PETSC_FALSE;
148c4762a1bSJed Brown PetscBool transfer_from_base[2] = {PETSC_TRUE, PETSC_FALSE};
149c4762a1bSJed Brown PetscBool use_bcs = PETSC_TRUE;
150c4762a1bSJed Brown bc_func_ctx bcCtx;
151c4762a1bSJed Brown DMLabel adaptLabel;
152c4762a1bSJed Brown
153327415f7SBarry Smith PetscFunctionBeginUser;
1549566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, NULL, help));
155c4762a1bSJed Brown comm = PETSC_COMM_WORLD;
156d0609cedSBarry Smith PetscOptionsBegin(comm, "", "DMForestTransferVec() Test Options", "DMFOREST");
1579566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-linear", "Transfer a simple linear function", "ex2.c", linear, &linear, NULL));
1589566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-coords", "Transfer a simple coordinate function", "ex2.c", coords, &coords, NULL));
1599566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-use_fv", "Use a finite volume approximation", "ex2.c", useFV, &useFV, NULL));
1609566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-test_convert", "Test conversion to DMPLEX", NULL, conv, &conv, NULL));
1619566063dSJacob 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));
162c4762a1bSJed Brown transfer_from_base[1] = transfer_from_base[0];
1639566063dSJacob 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));
1649566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-use_bcs", "Use dirichlet boundary conditions", "ex2.c", use_bcs, &use_bcs, NULL));
1659566063dSJacob Faibussowitsch PetscCall(PetscOptionsBoundedInt("-adapt_steps", "Number of adaptivity steps", "ex2.c", adaptSteps, &adaptSteps, NULL, 0));
166d0609cedSBarry Smith PetscOptionsEnd();
167c4762a1bSJed Brown
168c4762a1bSJed Brown tol = PetscMax(1.e-10, tol); /* XXX fix for quadruple precision -> why do I need to do this? */
169c4762a1bSJed Brown
17030602db0SMatthew G. Knepley /* the base mesh */
1719566063dSJacob Faibussowitsch PetscCall(DMCreate(comm, &base));
1729566063dSJacob Faibussowitsch PetscCall(DMSetType(base, DMPLEX));
1739566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(base));
17430602db0SMatthew G. Knepley
1759566063dSJacob Faibussowitsch PetscCall(AddIdentityLabel(base));
1769566063dSJacob Faibussowitsch PetscCall(DMGetDimension(base, &dim));
177c4762a1bSJed Brown
178ad540459SPierre Jolivet if (linear) funcs[0] = LinearFunction;
179c4762a1bSJed Brown if (coords) {
180c4762a1bSJed Brown funcs[0] = CoordsFunction;
181c4762a1bSJed Brown Nf = dim;
182c4762a1bSJed Brown }
183c4762a1bSJed Brown
184c4762a1bSJed Brown bcCtx.func = funcs[0];
185c4762a1bSJed Brown bcCtx.dim = dim;
186c4762a1bSJed Brown bcCtx.Nf = Nf;
187c4762a1bSJed Brown bcCtx.ctx = NULL;
188c4762a1bSJed Brown
189c4762a1bSJed Brown if (useFV) {
190c4762a1bSJed Brown PetscFV fv;
191c4762a1bSJed Brown PetscLimiter limiter;
192c4762a1bSJed Brown DM baseFV;
193c4762a1bSJed Brown
1949566063dSJacob Faibussowitsch PetscCall(DMPlexConstructGhostCells(base, NULL, NULL, &baseFV));
1959566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(baseFV, NULL, "-fv_dm_view"));
1969566063dSJacob Faibussowitsch PetscCall(DMDestroy(&base));
197c4762a1bSJed Brown base = baseFV;
19840cbb1a0SMatthew G. Knepley PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv));
1999566063dSJacob Faibussowitsch PetscCall(PetscFVSetSpatialDimension(fv, dim));
2009566063dSJacob Faibussowitsch PetscCall(PetscFVSetType(fv, PETSCFVLEASTSQUARES));
2019566063dSJacob Faibussowitsch PetscCall(PetscFVSetNumComponents(fv, Nf));
2029566063dSJacob Faibussowitsch PetscCall(PetscLimiterCreate(comm, &limiter));
2039566063dSJacob Faibussowitsch PetscCall(PetscLimiterSetType(limiter, PETSCLIMITERNONE));
2049566063dSJacob Faibussowitsch PetscCall(PetscFVSetLimiter(fv, limiter));
2059566063dSJacob Faibussowitsch PetscCall(PetscLimiterDestroy(&limiter));
2069566063dSJacob Faibussowitsch PetscCall(PetscFVSetFromOptions(fv));
2079566063dSJacob Faibussowitsch PetscCall(DMSetField(base, 0, NULL, (PetscObject)fv));
2089566063dSJacob Faibussowitsch PetscCall(PetscFVDestroy(&fv));
209c4762a1bSJed Brown } else {
210c4762a1bSJed Brown PetscFE fe;
211c4762a1bSJed Brown
2129566063dSJacob Faibussowitsch PetscCall(PetscFECreateDefault(comm, dim, Nf, PETSC_FALSE, NULL, PETSC_DEFAULT, &fe));
2139566063dSJacob Faibussowitsch PetscCall(DMSetField(base, 0, NULL, (PetscObject)fe));
2149566063dSJacob Faibussowitsch PetscCall(PetscFEDestroy(&fe));
215c4762a1bSJed Brown }
2169566063dSJacob Faibussowitsch PetscCall(DMCreateDS(base));
217c4762a1bSJed Brown
218c4762a1bSJed Brown if (use_bcs) {
219c4762a1bSJed Brown PetscInt ids[] = {1, 2, 3, 4, 5, 6};
22045480ffeSMatthew G. Knepley DMLabel label;
221c4762a1bSJed Brown
2229566063dSJacob Faibussowitsch PetscCall(DMGetLabel(base, "marker", &label));
22357d50842SBarry Smith PetscCall(DMAddBoundary(base, DM_BC_ESSENTIAL, "bc", label, 2 * dim, ids, 0, 0, NULL, useFV ? (PetscVoidFn *)bc_func_fv : (PetscVoidFn *)funcs[0], NULL, useFV ? (void *)&bcCtx : NULL, NULL));
224c4762a1bSJed Brown }
2259566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(base, NULL, "-dm_base_view"));
226c4762a1bSJed Brown
227c4762a1bSJed Brown /* the pre adaptivity forest */
2289566063dSJacob Faibussowitsch PetscCall(DMCreate(comm, &preForest));
2299566063dSJacob Faibussowitsch PetscCall(DMSetType(preForest, (dim == 2) ? DMP4EST : DMP8EST));
2309566063dSJacob Faibussowitsch PetscCall(DMCopyDisc(base, preForest));
2319566063dSJacob Faibussowitsch PetscCall(DMForestSetBaseDM(preForest, base));
2329566063dSJacob Faibussowitsch PetscCall(DMForestSetMinimumRefinement(preForest, 0));
2339566063dSJacob Faibussowitsch PetscCall(DMForestSetInitialRefinement(preForest, 1));
2349566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(preForest));
2359566063dSJacob Faibussowitsch PetscCall(DMSetUp(preForest));
2369566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(preForest, NULL, "-dm_pre_view"));
237c4762a1bSJed Brown
238c4762a1bSJed Brown /* the pre adaptivity field */
2399566063dSJacob Faibussowitsch PetscCall(DMCreateGlobalVector(preForest, &preVec));
2409566063dSJacob Faibussowitsch PetscCall(DMProjectFunction(preForest, 0., funcs, ctxs, INSERT_VALUES, preVec));
2419566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(preVec, NULL, "-vec_pre_view"));
242c4762a1bSJed Brown
243c4762a1bSJed Brown /* communicate between base and pre adaptivity forest */
244c4762a1bSJed Brown if (transfer_from_base[0]) {
245c4762a1bSJed Brown Vec baseVec, baseVecMapped;
246c4762a1bSJed Brown
2479566063dSJacob Faibussowitsch PetscCall(DMGetGlobalVector(base, &baseVec));
2489566063dSJacob Faibussowitsch PetscCall(DMProjectFunction(base, 0., funcs, ctxs, INSERT_VALUES, baseVec));
2499566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)baseVec, "Function Base"));
2509566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(baseVec, NULL, "-vec_base_view"));
251c4762a1bSJed Brown
2529566063dSJacob Faibussowitsch PetscCall(DMGetGlobalVector(preForest, &baseVecMapped));
2539566063dSJacob Faibussowitsch PetscCall(DMForestTransferVecFromBase(preForest, baseVec, baseVecMapped));
2549566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(baseVecMapped, NULL, "-vec_map_base_view"));
255c4762a1bSJed Brown
256c4762a1bSJed Brown /* compare */
2579566063dSJacob Faibussowitsch PetscCall(VecAXPY(baseVecMapped, -1., preVec));
2589566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(baseVecMapped, NULL, "-vec_map_diff_view"));
2599566063dSJacob Faibussowitsch PetscCall(VecNorm(baseVecMapped, NORM_2, &diff));
260c4762a1bSJed Brown
261c4762a1bSJed Brown /* output */
262c4762a1bSJed Brown if (diff < tol) {
2639566063dSJacob Faibussowitsch PetscCall(PetscPrintf(comm, "DMForestTransferVecFromBase() passes.\n"));
264c4762a1bSJed Brown } else {
2659566063dSJacob Faibussowitsch PetscCall(PetscPrintf(comm, "DMForestTransferVecFromBase() fails with error %g and tolerance %g\n", (double)diff, (double)tol));
266c4762a1bSJed Brown }
267c4762a1bSJed Brown
2689566063dSJacob Faibussowitsch PetscCall(DMRestoreGlobalVector(base, &baseVec));
2699566063dSJacob Faibussowitsch PetscCall(DMRestoreGlobalVector(preForest, &baseVecMapped));
270c4762a1bSJed Brown }
271c4762a1bSJed Brown
272c4762a1bSJed Brown for (step = 0; step < adaptSteps; ++step) {
27348a46eb9SPierre Jolivet if (!transfer_from_base[1]) PetscCall(PetscObjectGetReference((PetscObject)preForest, &preCount));
274c4762a1bSJed Brown
275c4762a1bSJed Brown /* adapt */
2769566063dSJacob Faibussowitsch PetscCall(CreateAdaptivityLabel(preForest, &adaptLabel));
2779566063dSJacob Faibussowitsch PetscCall(DMForestTemplate(preForest, comm, &postForest));
2789566063dSJacob Faibussowitsch if (step) PetscCall(DMForestSetAdaptivityLabel(postForest, adaptLabel));
2799566063dSJacob Faibussowitsch PetscCall(DMLabelDestroy(&adaptLabel));
2809566063dSJacob Faibussowitsch PetscCall(DMSetUp(postForest));
2819566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(postForest, NULL, "-dm_post_view"));
282c4762a1bSJed Brown
283c4762a1bSJed Brown /* transfer */
2849566063dSJacob Faibussowitsch PetscCall(DMCreateGlobalVector(postForest, &postVecTransfer));
2859566063dSJacob Faibussowitsch PetscCall(DMForestTransferVec(preForest, preVec, postForest, postVecTransfer, PETSC_TRUE, 0.0));
2869566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(postVecTransfer, NULL, "-vec_post_transfer_view"));
287c4762a1bSJed Brown
288c4762a1bSJed Brown /* the exact post adaptivity field */
2899566063dSJacob Faibussowitsch PetscCall(DMCreateGlobalVector(postForest, &postVecExact));
2909566063dSJacob Faibussowitsch PetscCall(DMProjectFunction(postForest, 0., funcs, ctxs, INSERT_VALUES, postVecExact));
2919566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(postVecExact, NULL, "-vec_post_exact_view"));
292c4762a1bSJed Brown
293c4762a1bSJed Brown /* compare */
2949566063dSJacob Faibussowitsch PetscCall(VecAXPY(postVecExact, -1., postVecTransfer));
2959566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(postVecExact, NULL, "-vec_diff_view"));
2969566063dSJacob Faibussowitsch PetscCall(VecNorm(postVecExact, NORM_2, &diff));
297c4762a1bSJed Brown
298c4762a1bSJed Brown /* output */
299c4762a1bSJed Brown if (diff < tol) {
3009566063dSJacob Faibussowitsch PetscCall(PetscPrintf(comm, "DMForestTransferVec() passes.\n"));
301c4762a1bSJed Brown } else {
3029566063dSJacob Faibussowitsch PetscCall(PetscPrintf(comm, "DMForestTransferVec() fails with error %g and tolerance %g\n", (double)diff, (double)tol));
3039566063dSJacob Faibussowitsch PetscCall(IdentifyBadPoints(postForest, postVecExact, tol));
304c4762a1bSJed Brown }
3059566063dSJacob Faibussowitsch PetscCall(VecDestroy(&postVecExact));
306c4762a1bSJed Brown
307c4762a1bSJed Brown /* disconnect preForest from postForest if we don't test the transfer throughout the entire refinement process */
308c4762a1bSJed Brown if (!transfer_from_base[1]) {
3099566063dSJacob Faibussowitsch PetscCall(DMForestSetAdaptivityForest(postForest, NULL));
3109566063dSJacob Faibussowitsch PetscCall(PetscObjectGetReference((PetscObject)preForest, &postCount));
31163a3b9bcSJacob 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);
312c4762a1bSJed Brown }
313c4762a1bSJed Brown
314c4762a1bSJed Brown if (conv) {
315c4762a1bSJed Brown DM dmConv;
316c4762a1bSJed Brown
3179566063dSJacob Faibussowitsch PetscCall(DMConvert(postForest, DMPLEX, &dmConv));
3189566063dSJacob Faibussowitsch PetscCall(DMViewFromOptions(dmConv, NULL, "-dm_conv_view"));
3199566063dSJacob Faibussowitsch PetscCall(DMPlexCheckCellShape(dmConv, PETSC_TRUE, PETSC_DETERMINE));
3209566063dSJacob Faibussowitsch PetscCall(DMDestroy(&dmConv));
321c4762a1bSJed Brown }
322c4762a1bSJed Brown
3239566063dSJacob Faibussowitsch PetscCall(VecDestroy(&preVec));
3249566063dSJacob Faibussowitsch PetscCall(DMDestroy(&preForest));
325c4762a1bSJed Brown
326c4762a1bSJed Brown preVec = postVecTransfer;
327c4762a1bSJed Brown preForest = postForest;
328c4762a1bSJed Brown }
329c4762a1bSJed Brown
330c4762a1bSJed Brown if (transfer_from_base[1]) {
331c4762a1bSJed Brown Vec baseVec, baseVecMapped;
332c4762a1bSJed Brown
333c4762a1bSJed Brown /* communicate between base and last adapted forest */
3349566063dSJacob Faibussowitsch PetscCall(DMGetGlobalVector(base, &baseVec));
3359566063dSJacob Faibussowitsch PetscCall(DMProjectFunction(base, 0., funcs, ctxs, INSERT_VALUES, baseVec));
3369566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)baseVec, "Function Base"));
3379566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(baseVec, NULL, "-vec_base_view"));
338c4762a1bSJed Brown
3399566063dSJacob Faibussowitsch PetscCall(DMGetGlobalVector(preForest, &baseVecMapped));
3409566063dSJacob Faibussowitsch PetscCall(DMForestTransferVecFromBase(preForest, baseVec, baseVecMapped));
3419566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(baseVecMapped, NULL, "-vec_map_base_view"));
342c4762a1bSJed Brown
343c4762a1bSJed Brown /* compare */
3449566063dSJacob Faibussowitsch PetscCall(VecAXPY(baseVecMapped, -1., preVec));
3459566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(baseVecMapped, NULL, "-vec_map_diff_view"));
3469566063dSJacob Faibussowitsch PetscCall(VecNorm(baseVecMapped, NORM_2, &diff));
347c4762a1bSJed Brown
348c4762a1bSJed Brown /* output */
349c4762a1bSJed Brown if (diff < tol) {
3509566063dSJacob Faibussowitsch PetscCall(PetscPrintf(comm, "DMForestTransferVecFromBase() passes.\n"));
351c4762a1bSJed Brown } else {
3529566063dSJacob Faibussowitsch PetscCall(PetscPrintf(comm, "DMForestTransferVecFromBase() fails with error %g and tolerance %g\n", (double)diff, (double)tol));
353c4762a1bSJed Brown }
354c4762a1bSJed Brown
3559566063dSJacob Faibussowitsch PetscCall(DMRestoreGlobalVector(base, &baseVec));
3569566063dSJacob Faibussowitsch PetscCall(DMRestoreGlobalVector(preForest, &baseVecMapped));
357c4762a1bSJed Brown }
358c4762a1bSJed Brown
359c4762a1bSJed Brown /* cleanup */
3609566063dSJacob Faibussowitsch PetscCall(VecDestroy(&preVec));
3619566063dSJacob Faibussowitsch PetscCall(DMDestroy(&preForest));
3629566063dSJacob Faibussowitsch PetscCall(DMDestroy(&base));
3639566063dSJacob Faibussowitsch PetscCall(PetscFinalize());
364b122ec5aSJacob Faibussowitsch return 0;
365c4762a1bSJed Brown }
366c4762a1bSJed Brown
367c4762a1bSJed Brown /*TEST
36830602db0SMatthew G. Knepley testset:
36930602db0SMatthew G. Knepley args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3 -petscspace_type tensor
370c4762a1bSJed Brown
371c4762a1bSJed Brown test:
372c4762a1bSJed Brown output_file: output/ex2_2d.out
373c4762a1bSJed Brown suffix: p4est_2d
37430602db0SMatthew G. Knepley args: -petscspace_degree 2
375c4762a1bSJed Brown nsize: 3
376c4762a1bSJed Brown requires: p4est !single
377c4762a1bSJed Brown
378c4762a1bSJed Brown test:
379c4762a1bSJed Brown output_file: output/ex2_2d.out
380c4762a1bSJed Brown suffix: p4est_2d_deg4
38130602db0SMatthew G. Knepley args: -petscspace_degree 4
382c4762a1bSJed Brown requires: p4est !single
383c4762a1bSJed Brown
384c4762a1bSJed Brown test:
385c4762a1bSJed Brown output_file: output/ex2_2d.out
386c4762a1bSJed Brown suffix: p4est_2d_deg8
38730602db0SMatthew G. Knepley args: -petscspace_degree 8
388c4762a1bSJed Brown requires: p4est !single
389c4762a1bSJed Brown
390c4762a1bSJed Brown test:
391c4762a1bSJed Brown output_file: output/ex2_steps2.out
392c4762a1bSJed Brown suffix: p4est_2d_deg2_steps2
39330602db0SMatthew G. Knepley args: -petscspace_degree 2 -coords -adapt_steps 2
394c4762a1bSJed Brown nsize: 3
395c4762a1bSJed Brown requires: p4est !single
396c4762a1bSJed Brown
397c4762a1bSJed Brown test:
398c4762a1bSJed Brown output_file: output/ex2_steps3.out
399c4762a1bSJed Brown suffix: p4est_2d_deg3_steps3
40030602db0SMatthew G. Knepley args: -petscspace_degree 3 -coords -adapt_steps 3 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1
401c4762a1bSJed Brown nsize: 3
402c4762a1bSJed Brown requires: p4est !single
403c4762a1bSJed Brown
404c4762a1bSJed Brown test:
405c4762a1bSJed Brown output_file: output/ex2_steps3.out
406c4762a1bSJed Brown suffix: p4est_2d_deg3_steps3_L2_periodic
40730602db0SMatthew 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
408c4762a1bSJed Brown nsize: 3
409c4762a1bSJed Brown requires: p4est !single
410c4762a1bSJed Brown
411c4762a1bSJed Brown test:
412c4762a1bSJed Brown output_file: output/ex2_steps3.out
413c4762a1bSJed Brown suffix: p4est_3d_deg2_steps3_L2_periodic
41430602db0SMatthew 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
415c4762a1bSJed Brown nsize: 3
416c4762a1bSJed Brown requires: p4est !single
417c4762a1bSJed Brown
418c4762a1bSJed Brown test:
419c4762a1bSJed Brown output_file: output/ex2_steps2.out
420c4762a1bSJed Brown suffix: p4est_3d_deg2_steps2
42130602db0SMatthew G. Knepley args: -dm_plex_dim 3 -petscspace_degree 2 -coords -adapt_steps 2
422c4762a1bSJed Brown nsize: 3
423c4762a1bSJed Brown requires: p4est !single
424c4762a1bSJed Brown
425c4762a1bSJed Brown test:
426c4762a1bSJed Brown output_file: output/ex2_steps3.out
427c4762a1bSJed Brown suffix: p4est_3d_deg3_steps3
42830602db0SMatthew G. Knepley args: -dm_plex_dim 3 -petscspace_degree 3 -coords -adapt_steps 3 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1
429c4762a1bSJed Brown nsize: 3
430c4762a1bSJed Brown requires: p4est !single
431c4762a1bSJed Brown
432c4762a1bSJed Brown test:
433c4762a1bSJed Brown output_file: output/ex2_3d.out
434c4762a1bSJed Brown suffix: p4est_3d
43530602db0SMatthew G. Knepley args: -dm_plex_dim 3 -petscspace_degree 1
436c4762a1bSJed Brown nsize: 3
437c4762a1bSJed Brown requires: p4est !single
438c4762a1bSJed Brown
439c4762a1bSJed Brown test:
440c4762a1bSJed Brown output_file: output/ex2_3d.out
441c4762a1bSJed Brown suffix: p4est_3d_deg3
44230602db0SMatthew G. Knepley args: -dm_plex_dim 3 -petscspace_degree 3
443c4762a1bSJed Brown nsize: 3
444c4762a1bSJed Brown requires: p4est !single
445c4762a1bSJed Brown
446c4762a1bSJed Brown test:
447c4762a1bSJed Brown output_file: output/ex2_2d.out
448c4762a1bSJed Brown suffix: p4est_2d_deg2_coords
44930602db0SMatthew G. Knepley args: -petscspace_degree 2 -coords
450c4762a1bSJed Brown nsize: 3
451c4762a1bSJed Brown requires: p4est !single
452c4762a1bSJed Brown
453c4762a1bSJed Brown test:
454c4762a1bSJed Brown output_file: output/ex2_3d.out
455c4762a1bSJed Brown suffix: p4est_3d_deg2_coords
45630602db0SMatthew G. Knepley args: -dm_plex_dim 3 -petscspace_degree 2 -coords
457c4762a1bSJed Brown nsize: 3
458c4762a1bSJed Brown requires: p4est !single
459c4762a1bSJed Brown
460c4762a1bSJed Brown test:
461c4762a1bSJed Brown suffix: p4est_3d_nans
46230602db0SMatthew G. Knepley args: -dm_plex_dim 3 -dm_forest_partition_overlap 1 -test_convert -petscspace_degree 1
463c4762a1bSJed Brown nsize: 2
464c4762a1bSJed Brown requires: p4est !single
465c4762a1bSJed Brown
466c4762a1bSJed Brown test:
467c4762a1bSJed Brown TODO: not broken, but the 3D case below is broken, so I do not trust this one
468c4762a1bSJed Brown output_file: output/ex2_steps2.out
469c4762a1bSJed Brown suffix: p4est_2d_tfb_distributed_nc
470e600fa54SMatthew 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
471c4762a1bSJed Brown nsize: 3
472c4762a1bSJed Brown requires: p4est !single
473c4762a1bSJed Brown
474c4762a1bSJed Brown test:
475c4762a1bSJed Brown TODO: broken
476c4762a1bSJed Brown output_file: output/ex2_steps2.out
477c4762a1bSJed Brown suffix: p4est_3d_tfb_distributed_nc
478e600fa54SMatthew 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
479c4762a1bSJed Brown nsize: 3
480c4762a1bSJed Brown requires: p4est !single
481c4762a1bSJed Brown
48230602db0SMatthew G. Knepley testset:
483012bc364SMatthew G. Knepley args: -petscspace_type tensor -dm_coord_space 0 -dm_plex_transform_type refine_tobox
48430602db0SMatthew G. Knepley
485c4762a1bSJed Brown test:
486c4762a1bSJed Brown TODO: broken
487c4762a1bSJed Brown output_file: output/ex2_3d.out
488c4762a1bSJed Brown suffix: p4est_3d_transfer_fails
489012bc364SMatthew 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
490c4762a1bSJed Brown requires: p4est !single
491c4762a1bSJed Brown
492c4762a1bSJed Brown test:
493c4762a1bSJed Brown TODO: broken
494c4762a1bSJed Brown output_file: output/ex2_steps2_notfb.out
495c4762a1bSJed Brown suffix: p4est_3d_transfer_fails_2
496012bc364SMatthew 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
497c4762a1bSJed Brown requires: p4est !single
498c4762a1bSJed Brown
499c4762a1bSJed Brown test:
500c4762a1bSJed Brown output_file: output/ex2_steps2.out
501c4762a1bSJed Brown suffix: p4est_3d_multi_transfer_s2t
502012bc364SMatthew 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
503012bc364SMatthew G. Knepley requires: p4est !single
504c4762a1bSJed Brown
505c4762a1bSJed Brown test:
506c4762a1bSJed Brown output_file: output/ex2_steps2.out
507c4762a1bSJed Brown suffix: p4est_3d_coords_transfer_s2t
508012bc364SMatthew 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
509012bc364SMatthew G. Knepley requires: p4est !single
51030602db0SMatthew G. Knepley
51130602db0SMatthew G. Knepley testset:
51230602db0SMatthew G. Knepley args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3
51330602db0SMatthew G. Knepley
51430602db0SMatthew G. Knepley test:
51530602db0SMatthew G. Knepley output_file: output/ex2_2d_fv.out
51630602db0SMatthew G. Knepley suffix: p4est_2d_fv
51730602db0SMatthew G. Knepley args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1
51830602db0SMatthew G. Knepley nsize: 3
51930602db0SMatthew G. Knepley requires: p4est !single
52030602db0SMatthew G. Knepley
52130602db0SMatthew G. Knepley test:
52230602db0SMatthew G. Knepley TODO: broken (codimension adjacency)
52330602db0SMatthew G. Knepley output_file: output/ex2_2d_fv.out
52430602db0SMatthew G. Knepley suffix: p4est_2d_fv_adjcodim
52530602db0SMatthew G. Knepley args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 -dm_forest_adjacency_codimension 1
52630602db0SMatthew G. Knepley nsize: 2
52730602db0SMatthew G. Knepley requires: p4est !single
52830602db0SMatthew G. Knepley
52930602db0SMatthew G. Knepley test:
53030602db0SMatthew G. Knepley TODO: broken (dimension adjacency)
53130602db0SMatthew G. Knepley output_file: output/ex2_2d_fv.out
53230602db0SMatthew G. Knepley suffix: p4est_2d_fv_adjdim
53330602db0SMatthew G. Knepley args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 -dm_forest_adjacency_dimension 1
53430602db0SMatthew G. Knepley nsize: 2
53530602db0SMatthew G. Knepley requires: p4est !single
53630602db0SMatthew G. Knepley
53730602db0SMatthew G. Knepley test:
53830602db0SMatthew G. Knepley output_file: output/ex2_2d_fv.out
53930602db0SMatthew G. Knepley suffix: p4est_2d_fv_zerocells
54030602db0SMatthew G. Knepley args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1
54130602db0SMatthew G. Knepley nsize: 10
54230602db0SMatthew G. Knepley requires: p4est !single
54330602db0SMatthew G. Knepley
54430602db0SMatthew G. Knepley test:
54530602db0SMatthew G. Knepley output_file: output/ex2_3d_fv.out
54630602db0SMatthew G. Knepley suffix: p4est_3d_fv
54730602db0SMatthew G. Knepley args: -dm_plex_dim 3 -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1
54830602db0SMatthew G. Knepley nsize: 3
549c4762a1bSJed Brown requires: p4est !single
550c4762a1bSJed Brown
551c4762a1bSJed Brown TEST*/
552