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 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 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 38d71ae5a4SJacob Faibussowitsch static PetscErrorCode LinearFunction(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *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 45d71ae5a4SJacob Faibussowitsch static PetscErrorCode MultiaffineFunction(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *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 52d71ae5a4SJacob Faibussowitsch static PetscErrorCode CoordsFunction(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *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 68d71ae5a4SJacob Faibussowitsch static PetscErrorCode bc_func_fv(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx) 69d71ae5a4SJacob Faibussowitsch { 70c4762a1bSJed Brown bc_func_ctx *bcCtx; 71c4762a1bSJed Brown 72c4762a1bSJed Brown PetscFunctionBegin; 73c4762a1bSJed Brown bcCtx = (bc_func_ctx *)ctx; 74*f4f49eeaSPierre Jolivet PetscCall(bcCtx->func(bcCtx->dim, time, c, bcCtx->Nf, xG, bcCtx->ctx)); 753ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 76c4762a1bSJed Brown } 77c4762a1bSJed Brown 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 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)); 2239566063dSJacob 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)); 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