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; 13*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateLabel(dm, "identity")); 14*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 15*5f80ce2aSJacob Faibussowitsch for (p = pStart; p < pEnd; p++) CHKERRQ(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; 25*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelCreate(PETSC_COMM_SELF,"adapt",adaptLabel)); 26*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelSetDefaultValue(*adaptLabel,DM_ADAPT_COARSEN)); 27*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(forest,"identity",&identLabel)); 28*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMForestGetCellChart(forest,&cStart,&cEnd)); 29c4762a1bSJed Brown for (c = cStart; c < cEnd; c++) { 30c4762a1bSJed Brown PetscInt basePoint; 31c4762a1bSJed Brown 32*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelGetValue(identLabel,c,&basePoint)); 33*5f80ce2aSJacob Faibussowitsch if (!basePoint) CHKERRQ(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; 76*5f80ce2aSJacob Faibussowitsch CHKERRQ((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; 89*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateLocalVector(dm, &vecLocal)); 90*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal)); 91*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal)); 92*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMConvert(dm ,DMPLEX, &dmplex)); 93*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetChart(dmplex, &pStart, &pEnd)); 94*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetDepthLabel(dmplex, &depthLabel)); 95*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLocalSection(dmplex, §ion)); 96*5f80ce2aSJacob Faibussowitsch CHKERRQ(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 103*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetValuesSection(vecLocal, section, p, &values)); 104*5f80ce2aSJacob Faibussowitsch CHKERRQ(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; 110*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "Bad point %D\n", p)); 111*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelGetValue(depthLabel, p, &s)); 112*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " Depth %D\n", s)); 113*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetTransitiveClosure(dmplex, p, PETSC_TRUE, &closureSize, &closure)); 114c4762a1bSJed Brown for (cl = 0; cl < closureSize; cl++) { 115c4762a1bSJed Brown PetscInt cp = closure[2 * cl]; 116*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetTreeParent(dmplex, cp, &parent, &childID)); 117c4762a1bSJed Brown if (parent != cp) { 118*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " Closure point %D (%D) child of %D (ID %D)\n", cl, cp, parent, childID)); 119c4762a1bSJed Brown } 120*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetTreeChildren(dmplex, cp, &numChildren, NULL)); 121c4762a1bSJed Brown if (numChildren) { 122*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " Closure point %D (%D) is parent\n", cl, cp)); 123c4762a1bSJed Brown } 124c4762a1bSJed Brown } 125*5f80ce2aSJacob Faibussowitsch CHKERRQ(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) { 129*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " Bad dof %D\n", c)); 130c4762a1bSJed Brown } 131c4762a1bSJed Brown } 132c4762a1bSJed Brown } 133*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(&dmplex)); 134*5f80ce2aSJacob Faibussowitsch CHKERRQ(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 PetscErrorCode ierr; 158c4762a1bSJed Brown 159c4762a1bSJed Brown ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 160c4762a1bSJed Brown comm = PETSC_COMM_WORLD; 161c4762a1bSJed Brown ierr = PetscOptionsBegin(comm, "", "DMForestTransferVec() Test Options", "DMFOREST");CHKERRQ(ierr); 162*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBool("-linear","Transfer a simple linear function", "ex2.c", linear, &linear, NULL)); 163*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBool("-coords","Transfer a simple coordinate function", "ex2.c", coords, &coords, NULL)); 164*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBool("-use_fv","Use a finite volume approximation", "ex2.c", useFV, &useFV, NULL)); 165*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBool("-test_convert","Test conversion to DMPLEX",NULL,conv,&conv,NULL)); 166*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBool("-transfer_from_base","Transfer a vector from base DM to DMForest", "ex2.c", transfer_from_base[0], &transfer_from_base[0], NULL)); 167c4762a1bSJed Brown transfer_from_base[1] = transfer_from_base[0]; 168*5f80ce2aSJacob Faibussowitsch CHKERRQ(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)); 169*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBool("-use_bcs","Use dirichlet boundary conditions", "ex2.c", use_bcs, &use_bcs, NULL)); 170*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBoundedInt("-adapt_steps","Number of adaptivity steps", "ex2.c", adaptSteps, &adaptSteps, NULL,0)); 171c4762a1bSJed Brown ierr = PetscOptionsEnd();CHKERRQ(ierr); 172c4762a1bSJed Brown 173c4762a1bSJed Brown tol = PetscMax(1.e-10,tol); /* XXX fix for quadruple precision -> why do I need to do this? */ 174c4762a1bSJed Brown 17530602db0SMatthew G. Knepley /* the base mesh */ 176*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreate(comm, &base)); 177*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetType(base, DMPLEX)); 178*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetFromOptions(base)); 17930602db0SMatthew G. Knepley 180*5f80ce2aSJacob Faibussowitsch CHKERRQ(AddIdentityLabel(base)); 181*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetDimension(base, &dim)); 182c4762a1bSJed Brown 183c4762a1bSJed Brown if (linear) { 184c4762a1bSJed Brown funcs[0] = LinearFunction; 185c4762a1bSJed Brown } 186c4762a1bSJed Brown if (coords) { 187c4762a1bSJed Brown funcs[0] = CoordsFunction; 188c4762a1bSJed Brown Nf = dim; 189c4762a1bSJed Brown } 190c4762a1bSJed Brown 191c4762a1bSJed Brown bcCtx.func = funcs[0]; 192c4762a1bSJed Brown bcCtx.dim = dim; 193c4762a1bSJed Brown bcCtx.Nf = Nf; 194c4762a1bSJed Brown bcCtx.ctx = NULL; 195c4762a1bSJed Brown 196c4762a1bSJed Brown if (useFV) { 197c4762a1bSJed Brown PetscFV fv; 198c4762a1bSJed Brown PetscLimiter limiter; 199c4762a1bSJed Brown DM baseFV; 200c4762a1bSJed Brown 201*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexConstructGhostCells(base,NULL,NULL,&baseFV)); 202*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(baseFV, NULL, "-fv_dm_view")); 203*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(&base)); 204c4762a1bSJed Brown base = baseFV; 205*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFVCreate(comm, &fv)); 206*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFVSetSpatialDimension(fv,dim)); 207*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFVSetType(fv,PETSCFVLEASTSQUARES)); 208*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFVSetNumComponents(fv,Nf)); 209*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscLimiterCreate(comm,&limiter)); 210*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscLimiterSetType(limiter,PETSCLIMITERNONE)); 211*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFVSetLimiter(fv,limiter)); 212*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscLimiterDestroy(&limiter)); 213*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFVSetFromOptions(fv)); 214*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetField(base,0,NULL,(PetscObject)fv)); 215*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFVDestroy(&fv)); 216c4762a1bSJed Brown } else { 217c4762a1bSJed Brown PetscFE fe; 218c4762a1bSJed Brown 219*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFECreateDefault(comm,dim,Nf,PETSC_FALSE,NULL,PETSC_DEFAULT,&fe)); 220*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetField(base,0,NULL,(PetscObject)fe)); 221*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFEDestroy(&fe)); 222c4762a1bSJed Brown } 223*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateDS(base)); 224c4762a1bSJed Brown 225c4762a1bSJed Brown if (use_bcs) { 226c4762a1bSJed Brown PetscInt ids[] = {1, 2, 3, 4, 5, 6}; 22745480ffeSMatthew G. Knepley DMLabel label; 228c4762a1bSJed Brown 229*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetLabel(base, "marker", &label)); 230*5f80ce2aSJacob Faibussowitsch CHKERRQ(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)); 231c4762a1bSJed Brown } 232*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(base,NULL,"-dm_base_view")); 233c4762a1bSJed Brown 234c4762a1bSJed Brown /* the pre adaptivity forest */ 235*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreate(comm,&preForest)); 236*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetType(preForest,(dim == 2) ? DMP4EST : DMP8EST)); 237*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCopyDisc(base,preForest)); 238*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMForestSetBaseDM(preForest,base)); 239*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMForestSetMinimumRefinement(preForest,0)); 240*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMForestSetInitialRefinement(preForest,1)); 241*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetFromOptions(preForest)); 242*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetUp(preForest)); 243*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(preForest,NULL,"-dm_pre_view")); 244c4762a1bSJed Brown 245c4762a1bSJed Brown /* the pre adaptivity field */ 246*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateGlobalVector(preForest,&preVec)); 247*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMProjectFunction(preForest,0.,funcs,ctxs,INSERT_VALUES,preVec)); 248*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(preVec,NULL,"-vec_pre_view")); 249c4762a1bSJed Brown 250c4762a1bSJed Brown /* communicate between base and pre adaptivity forest */ 251c4762a1bSJed Brown if (transfer_from_base[0]) { 252c4762a1bSJed Brown Vec baseVec, baseVecMapped; 253c4762a1bSJed Brown 254*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetGlobalVector(base,&baseVec)); 255*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMProjectFunction(base,0.,funcs,ctxs,INSERT_VALUES,baseVec)); 256*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetName((PetscObject)baseVec,"Function Base")); 257*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(baseVec,NULL,"-vec_base_view")); 258c4762a1bSJed Brown 259*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetGlobalVector(preForest,&baseVecMapped)); 260*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMForestTransferVecFromBase(preForest,baseVec,baseVecMapped)); 261*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(baseVecMapped,NULL,"-vec_map_base_view")); 262c4762a1bSJed Brown 263c4762a1bSJed Brown /* compare */ 264*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecAXPY(baseVecMapped,-1.,preVec)); 265*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(baseVecMapped,NULL,"-vec_map_diff_view")); 266*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecNorm(baseVecMapped,NORM_2,&diff)); 267c4762a1bSJed Brown 268c4762a1bSJed Brown /* output */ 269c4762a1bSJed Brown if (diff < tol) { 270*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(comm,"DMForestTransferVecFromBase() passes.\n")); 271c4762a1bSJed Brown } else { 272*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(comm,"DMForestTransferVecFromBase() fails with error %g and tolerance %g\n",(double)diff,(double)tol)); 273c4762a1bSJed Brown } 274c4762a1bSJed Brown 275*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMRestoreGlobalVector(base,&baseVec)); 276*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMRestoreGlobalVector(preForest,&baseVecMapped)); 277c4762a1bSJed Brown } 278c4762a1bSJed Brown 279c4762a1bSJed Brown for (step = 0; step < adaptSteps; ++step) { 280c4762a1bSJed Brown 281c4762a1bSJed Brown if (!transfer_from_base[1]) { 282*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectGetReference((PetscObject)preForest,&preCount)); 283c4762a1bSJed Brown } 284c4762a1bSJed Brown 285c4762a1bSJed Brown /* adapt */ 286*5f80ce2aSJacob Faibussowitsch CHKERRQ(CreateAdaptivityLabel(preForest,&adaptLabel)); 287*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMForestTemplate(preForest,comm,&postForest)); 288*5f80ce2aSJacob Faibussowitsch if (step) CHKERRQ(DMForestSetAdaptivityLabel(postForest,adaptLabel)); 289*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMLabelDestroy(&adaptLabel)); 290*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMSetUp(postForest)); 291*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(postForest,NULL,"-dm_post_view")); 292c4762a1bSJed Brown 293c4762a1bSJed Brown /* transfer */ 294*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateGlobalVector(postForest,&postVecTransfer)); 295*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMForestTransferVec(preForest,preVec,postForest,postVecTransfer,PETSC_TRUE,0.0)); 296*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(postVecTransfer,NULL,"-vec_post_transfer_view")); 297c4762a1bSJed Brown 298c4762a1bSJed Brown /* the exact post adaptivity field */ 299*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMCreateGlobalVector(postForest,&postVecExact)); 300*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMProjectFunction(postForest,0.,funcs,ctxs,INSERT_VALUES,postVecExact)); 301*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(postVecExact,NULL,"-vec_post_exact_view")); 302c4762a1bSJed Brown 303c4762a1bSJed Brown /* compare */ 304*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecAXPY(postVecExact,-1.,postVecTransfer)); 305*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(postVecExact,NULL,"-vec_diff_view")); 306*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecNorm(postVecExact,NORM_2,&diff)); 307c4762a1bSJed Brown 308c4762a1bSJed Brown /* output */ 309c4762a1bSJed Brown if (diff < tol) { 310*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(comm,"DMForestTransferVec() passes.\n")); 311c4762a1bSJed Brown } else { 312*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(comm,"DMForestTransferVec() fails with error %g and tolerance %g\n",(double)diff,(double)tol)); 313*5f80ce2aSJacob Faibussowitsch CHKERRQ(IdentifyBadPoints(postForest, postVecExact, tol)); 314c4762a1bSJed Brown } 315*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&postVecExact)); 316c4762a1bSJed Brown 317c4762a1bSJed Brown /* disconnect preForest from postForest if we don't test the transfer throughout the entire refinement process */ 318c4762a1bSJed Brown if (!transfer_from_base[1]) { 319*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMForestSetAdaptivityForest(postForest,NULL)); 320*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectGetReference((PetscObject)preForest,&postCount)); 3212c71b3e2SJacob Faibussowitsch PetscCheckFalse(postCount != preCount,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Adaptation not memory neutral: reference count increase from %d to %d",preCount,postCount); 322c4762a1bSJed Brown } 323c4762a1bSJed Brown 324c4762a1bSJed Brown if (conv) { 325c4762a1bSJed Brown DM dmConv; 326c4762a1bSJed Brown 327*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMConvert(postForest,DMPLEX,&dmConv)); 328*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMViewFromOptions(dmConv,NULL,"-dm_conv_view")); 329*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCheckCellShape(dmConv,PETSC_TRUE,PETSC_DETERMINE)); 330*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(&dmConv)); 331c4762a1bSJed Brown } 332c4762a1bSJed Brown 333*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&preVec)); 334*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(&preForest)); 335c4762a1bSJed Brown 336c4762a1bSJed Brown preVec = postVecTransfer; 337c4762a1bSJed Brown preForest = postForest; 338c4762a1bSJed Brown } 339c4762a1bSJed Brown 340c4762a1bSJed Brown if (transfer_from_base[1]) { 341c4762a1bSJed Brown Vec baseVec, baseVecMapped; 342c4762a1bSJed Brown 343c4762a1bSJed Brown /* communicate between base and last adapted forest */ 344*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetGlobalVector(base,&baseVec)); 345*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMProjectFunction(base,0.,funcs,ctxs,INSERT_VALUES,baseVec)); 346*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscObjectSetName((PetscObject)baseVec,"Function Base")); 347*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(baseVec,NULL,"-vec_base_view")); 348c4762a1bSJed Brown 349*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMGetGlobalVector(preForest,&baseVecMapped)); 350*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMForestTransferVecFromBase(preForest,baseVec,baseVecMapped)); 351*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(baseVecMapped,NULL,"-vec_map_base_view")); 352c4762a1bSJed Brown 353c4762a1bSJed Brown /* compare */ 354*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecAXPY(baseVecMapped,-1.,preVec)); 355*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecViewFromOptions(baseVecMapped,NULL,"-vec_map_diff_view")); 356*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecNorm(baseVecMapped,NORM_2,&diff)); 357c4762a1bSJed Brown 358c4762a1bSJed Brown /* output */ 359c4762a1bSJed Brown if (diff < tol) { 360*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(comm,"DMForestTransferVecFromBase() passes.\n")); 361c4762a1bSJed Brown } else { 362*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(comm,"DMForestTransferVecFromBase() fails with error %g and tolerance %g\n",(double)diff,(double)tol)); 363c4762a1bSJed Brown } 364c4762a1bSJed Brown 365*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMRestoreGlobalVector(base,&baseVec)); 366*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMRestoreGlobalVector(preForest,&baseVecMapped)); 367c4762a1bSJed Brown } 368c4762a1bSJed Brown 369c4762a1bSJed Brown /* cleanup */ 370*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&preVec)); 371*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(&preForest)); 372*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(&base)); 373c4762a1bSJed Brown ierr = PetscFinalize(); 374c4762a1bSJed Brown return ierr; 375c4762a1bSJed Brown } 376c4762a1bSJed Brown 377c4762a1bSJed Brown /*TEST 37830602db0SMatthew G. Knepley testset: 37930602db0SMatthew G. Knepley args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3 -petscspace_type tensor 380c4762a1bSJed Brown 381c4762a1bSJed Brown test: 382c4762a1bSJed Brown output_file: output/ex2_2d.out 383c4762a1bSJed Brown suffix: p4est_2d 38430602db0SMatthew G. Knepley args: -petscspace_degree 2 385c4762a1bSJed Brown nsize: 3 386c4762a1bSJed Brown requires: p4est !single 387c4762a1bSJed Brown 388c4762a1bSJed Brown test: 389c4762a1bSJed Brown output_file: output/ex2_2d.out 390c4762a1bSJed Brown suffix: p4est_2d_deg4 39130602db0SMatthew G. Knepley args: -petscspace_degree 4 392c4762a1bSJed Brown requires: p4est !single 393c4762a1bSJed Brown 394c4762a1bSJed Brown test: 395c4762a1bSJed Brown output_file: output/ex2_2d.out 396c4762a1bSJed Brown suffix: p4est_2d_deg8 39730602db0SMatthew G. Knepley args: -petscspace_degree 8 398c4762a1bSJed Brown requires: p4est !single 399c4762a1bSJed Brown 400c4762a1bSJed Brown test: 401c4762a1bSJed Brown output_file: output/ex2_steps2.out 402c4762a1bSJed Brown suffix: p4est_2d_deg2_steps2 40330602db0SMatthew G. Knepley args: -petscspace_degree 2 -coords -adapt_steps 2 404c4762a1bSJed Brown nsize: 3 405c4762a1bSJed Brown requires: p4est !single 406c4762a1bSJed Brown 407c4762a1bSJed Brown test: 408c4762a1bSJed Brown output_file: output/ex2_steps3.out 409c4762a1bSJed Brown suffix: p4est_2d_deg3_steps3 41030602db0SMatthew G. Knepley args: -petscspace_degree 3 -coords -adapt_steps 3 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1 411c4762a1bSJed Brown nsize: 3 412c4762a1bSJed Brown requires: p4est !single 413c4762a1bSJed Brown 414c4762a1bSJed Brown test: 415c4762a1bSJed Brown output_file: output/ex2_steps3.out 416c4762a1bSJed Brown suffix: p4est_2d_deg3_steps3_L2_periodic 41730602db0SMatthew 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 418c4762a1bSJed Brown nsize: 3 419c4762a1bSJed Brown requires: p4est !single 420c4762a1bSJed Brown 421c4762a1bSJed Brown test: 422c4762a1bSJed Brown output_file: output/ex2_steps3.out 423c4762a1bSJed Brown suffix: p4est_3d_deg2_steps3_L2_periodic 42430602db0SMatthew 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 425c4762a1bSJed Brown nsize: 3 426c4762a1bSJed Brown requires: p4est !single 427c4762a1bSJed Brown 428c4762a1bSJed Brown test: 429c4762a1bSJed Brown output_file: output/ex2_steps2.out 430c4762a1bSJed Brown suffix: p4est_3d_deg2_steps2 43130602db0SMatthew G. Knepley args: -dm_plex_dim 3 -petscspace_degree 2 -coords -adapt_steps 2 432c4762a1bSJed Brown nsize: 3 433c4762a1bSJed Brown requires: p4est !single 434c4762a1bSJed Brown 435c4762a1bSJed Brown test: 436c4762a1bSJed Brown output_file: output/ex2_steps3.out 437c4762a1bSJed Brown suffix: p4est_3d_deg3_steps3 43830602db0SMatthew G. Knepley args: -dm_plex_dim 3 -petscspace_degree 3 -coords -adapt_steps 3 -petscdualspace_lagrange_node_type equispaced -petscdualspace_lagrange_node_endpoints 1 439c4762a1bSJed Brown nsize: 3 440c4762a1bSJed Brown requires: p4est !single 441c4762a1bSJed Brown 442c4762a1bSJed Brown test: 443c4762a1bSJed Brown output_file: output/ex2_3d.out 444c4762a1bSJed Brown suffix: p4est_3d 44530602db0SMatthew G. Knepley args: -dm_plex_dim 3 -petscspace_degree 1 446c4762a1bSJed Brown nsize: 3 447c4762a1bSJed Brown requires: p4est !single 448c4762a1bSJed Brown 449c4762a1bSJed Brown test: 450c4762a1bSJed Brown output_file: output/ex2_3d.out 451c4762a1bSJed Brown suffix: p4est_3d_deg3 45230602db0SMatthew G. Knepley args: -dm_plex_dim 3 -petscspace_degree 3 453c4762a1bSJed Brown nsize: 3 454c4762a1bSJed Brown requires: p4est !single 455c4762a1bSJed Brown 456c4762a1bSJed Brown test: 457c4762a1bSJed Brown output_file: output/ex2_2d.out 458c4762a1bSJed Brown suffix: p4est_2d_deg2_coords 45930602db0SMatthew G. Knepley args: -petscspace_degree 2 -coords 460c4762a1bSJed Brown nsize: 3 461c4762a1bSJed Brown requires: p4est !single 462c4762a1bSJed Brown 463c4762a1bSJed Brown test: 464c4762a1bSJed Brown output_file: output/ex2_3d.out 465c4762a1bSJed Brown suffix: p4est_3d_deg2_coords 46630602db0SMatthew G. Knepley args: -dm_plex_dim 3 -petscspace_degree 2 -coords 467c4762a1bSJed Brown nsize: 3 468c4762a1bSJed Brown requires: p4est !single 469c4762a1bSJed Brown 470c4762a1bSJed Brown test: 471c4762a1bSJed Brown suffix: p4est_3d_nans 47230602db0SMatthew G. Knepley args: -dm_plex_dim 3 -dm_forest_partition_overlap 1 -test_convert -petscspace_degree 1 473c4762a1bSJed Brown nsize: 2 474c4762a1bSJed Brown requires: p4est !single 475c4762a1bSJed Brown 476c4762a1bSJed Brown test: 477c4762a1bSJed Brown TODO: not broken, but the 3D case below is broken, so I do not trust this one 478c4762a1bSJed Brown output_file: output/ex2_steps2.out 479c4762a1bSJed Brown suffix: p4est_2d_tfb_distributed_nc 480e600fa54SMatthew 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 481c4762a1bSJed Brown nsize: 3 482c4762a1bSJed Brown requires: p4est !single 483c4762a1bSJed Brown 484c4762a1bSJed Brown test: 485c4762a1bSJed Brown TODO: broken 486c4762a1bSJed Brown output_file: output/ex2_steps2.out 487c4762a1bSJed Brown suffix: p4est_3d_tfb_distributed_nc 488e600fa54SMatthew 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 489c4762a1bSJed Brown nsize: 3 490c4762a1bSJed Brown requires: p4est !single 491c4762a1bSJed Brown 49230602db0SMatthew G. Knepley testset: 493012bc364SMatthew G. Knepley args: -petscspace_type tensor -dm_coord_space 0 -dm_plex_transform_type refine_tobox 49430602db0SMatthew G. Knepley 495c4762a1bSJed Brown test: 496c4762a1bSJed Brown TODO: broken 497c4762a1bSJed Brown output_file: output/ex2_3d.out 498c4762a1bSJed Brown suffix: p4est_3d_transfer_fails 499012bc364SMatthew 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 500c4762a1bSJed Brown requires: p4est !single 501c4762a1bSJed Brown 502c4762a1bSJed Brown test: 503c4762a1bSJed Brown TODO: broken 504c4762a1bSJed Brown output_file: output/ex2_steps2_notfb.out 505c4762a1bSJed Brown suffix: p4est_3d_transfer_fails_2 506012bc364SMatthew 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 507c4762a1bSJed Brown requires: p4est !single 508c4762a1bSJed Brown 509c4762a1bSJed Brown test: 510c4762a1bSJed Brown output_file: output/ex2_steps2.out 511c4762a1bSJed Brown suffix: p4est_3d_multi_transfer_s2t 512012bc364SMatthew 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 513012bc364SMatthew G. Knepley requires: p4est !single 514c4762a1bSJed Brown 515c4762a1bSJed Brown test: 516c4762a1bSJed Brown output_file: output/ex2_steps2.out 517c4762a1bSJed Brown suffix: p4est_3d_coords_transfer_s2t 518012bc364SMatthew 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 519012bc364SMatthew G. Knepley requires: p4est !single 52030602db0SMatthew G. Knepley 52130602db0SMatthew G. Knepley testset: 52230602db0SMatthew G. Knepley args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3 52330602db0SMatthew G. Knepley 52430602db0SMatthew G. Knepley test: 52530602db0SMatthew G. Knepley output_file: output/ex2_2d_fv.out 52630602db0SMatthew G. Knepley suffix: p4est_2d_fv 52730602db0SMatthew G. Knepley args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 52830602db0SMatthew G. Knepley nsize: 3 52930602db0SMatthew G. Knepley requires: p4est !single 53030602db0SMatthew G. Knepley 53130602db0SMatthew G. Knepley test: 53230602db0SMatthew G. Knepley TODO: broken (codimension adjacency) 53330602db0SMatthew G. Knepley output_file: output/ex2_2d_fv.out 53430602db0SMatthew G. Knepley suffix: p4est_2d_fv_adjcodim 53530602db0SMatthew G. Knepley args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 -dm_forest_adjacency_codimension 1 53630602db0SMatthew G. Knepley nsize: 2 53730602db0SMatthew G. Knepley requires: p4est !single 53830602db0SMatthew G. Knepley 53930602db0SMatthew G. Knepley test: 54030602db0SMatthew G. Knepley TODO: broken (dimension adjacency) 54130602db0SMatthew G. Knepley output_file: output/ex2_2d_fv.out 54230602db0SMatthew G. Knepley suffix: p4est_2d_fv_adjdim 54330602db0SMatthew G. Knepley args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 -dm_forest_adjacency_dimension 1 54430602db0SMatthew G. Knepley nsize: 2 54530602db0SMatthew G. Knepley requires: p4est !single 54630602db0SMatthew G. Knepley 54730602db0SMatthew G. Knepley test: 54830602db0SMatthew G. Knepley output_file: output/ex2_2d_fv.out 54930602db0SMatthew G. Knepley suffix: p4est_2d_fv_zerocells 55030602db0SMatthew G. Knepley args: -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 55130602db0SMatthew G. Knepley nsize: 10 55230602db0SMatthew G. Knepley requires: p4est !single 55330602db0SMatthew G. Knepley 55430602db0SMatthew G. Knepley test: 55530602db0SMatthew G. Knepley output_file: output/ex2_3d_fv.out 55630602db0SMatthew G. Knepley suffix: p4est_3d_fv 55730602db0SMatthew G. Knepley args: -dm_plex_dim 3 -transfer_from_base 0 -use_fv -linear -dm_forest_partition_overlap 1 55830602db0SMatthew G. Knepley nsize: 3 559c4762a1bSJed Brown requires: p4est !single 560c4762a1bSJed Brown 561c4762a1bSJed Brown TEST*/ 562