Actual source code: ex2.c
petsc-3.9.4 2018-09-11
1: static char help[] = "Create a mesh, refine and coarsen simultaneously, and transfer a field\n\n";
3: #include <petscds.h>
4: #include <petscdmplex.h>
5: #include <petscdmforest.h>
6: #include <petscoptions.h>
8: static PetscErrorCode AddIdentityLabel(DM dm)
9: {
10: PetscInt pStart,pEnd,p;
14: DMCreateLabel(dm, "identity");
15: DMPlexGetChart(dm, &pStart, &pEnd);
16: for (p = pStart; p < pEnd; p++) {DMSetLabelValue(dm, "identity", p, p);}
17: return(0);
18: }
20: static PetscErrorCode CreateAdaptivityLabel(DM forest,DMLabel *adaptLabel)
21: {
22: DMLabel identLabel;
23: PetscInt cStart, cEnd, c;
27: DMLabelCreate("adapt",adaptLabel);
28: DMLabelSetDefaultValue(*adaptLabel,DM_ADAPT_COARSEN);
29: DMGetLabel(forest,"identity",&identLabel);
30: DMForestGetCellChart(forest,&cStart,&cEnd);
31: for (c = cStart; c < cEnd; c++) {
32: PetscInt basePoint;
34: DMLabelGetValue(identLabel,c,&basePoint);
35: if (!basePoint) {DMLabelSetValue(*adaptLabel,c,DM_ADAPT_REFINE);}
36: }
37: return(0);
38: }
40: static PetscErrorCode LinearFunction(PetscInt dim,PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
41: {
43: u[0] = (x[0] * 2.0 + 1.) + (x[1] * 20.0 + 10.) + ((dim == 3) ? (x[2] * 200.0 + 100.) : 0.);
44: return(0);
45: }
47: static PetscErrorCode MultiaffineFunction(PetscInt dim,PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
48: {
50: u[0] = (x[0] * 1.0 + 2.0) * (x[1] * 3.0 - 4.0) * ((dim == 3) ? (x[2] * 5.0 + 6.0) : 1.);
51: return(0);
52: }
54: typedef struct _bc_func_ctx
55: {
56: PetscErrorCode (*func) (PetscInt,PetscReal,const PetscReal [], PetscInt, PetscScalar [], void *);
57: PetscInt dim;
58: PetscInt Nf;
59: void *ctx;
60: }
61: bc_func_ctx;
63: static PetscErrorCode bc_func_fv (PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
64: {
65: bc_func_ctx *bcCtx;
69: bcCtx = (bc_func_ctx *) ctx;
70: (bcCtx->func)(bcCtx->dim,time,c,bcCtx->Nf,xG,bcCtx->ctx);
71: return(0);
72: }
74: int main(int argc, char **argv)
75: {
76: MPI_Comm comm;
77: DM base, preForest, postForest;
78: PetscInt dim = 2;
79: PetscInt preCount, postCount;
80: Vec preVec, postVecTransfer, postVecExact;
81: PetscErrorCode (*funcs[1]) (PetscInt,PetscReal,const PetscReal [],PetscInt,PetscScalar [], void *) = {MultiaffineFunction};
82: void *ctxs[1] = {NULL};
83: const PetscInt cells[] = {3, 3, 3};
84: PetscReal diff, tol = PETSC_SMALL;
85: PetscBool linear = PETSC_FALSE;
86: PetscBool useFV = PETSC_FALSE;
87: PetscDS ds;
88: bc_func_ctx bcCtx;
89: DMLabel adaptLabel;
92: PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
93: comm = PETSC_COMM_WORLD;
94: PetscOptionsBegin(comm, "", "DMForestTransferVec() Test Options", "DMFOREST");
95: PetscOptionsInt("-dim", "The dimension (2 or 3)", "ex2.c", dim, &dim, NULL);
96: PetscOptionsBool("-linear","Transfer a simple linear function", "ex2.c", linear, &linear, NULL);
97: PetscOptionsBool("-use_fv","Use a finite volume approximation", "ex2.c", useFV, &useFV, NULL);
98: PetscOptionsEnd();
100: if (linear) {
101: funcs[0] = LinearFunction;
102: }
104: bcCtx.func = funcs[0];
105: bcCtx.dim = dim;
106: bcCtx.Nf = 1;
107: bcCtx.ctx = NULL;
109: /* the base mesh */
110: DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, cells, NULL, NULL, NULL, PETSC_TRUE, &base);
111: if (useFV) {
112: PetscFV fv;
113: PetscLimiter limiter;
114: DM baseFV;
116: DMPlexConstructGhostCells(base,NULL,NULL,&baseFV);
117: DMDestroy(&base);
118: base = baseFV;
119: PetscFVCreate(comm, &fv);
120: PetscFVSetSpatialDimension(fv,dim);
121: PetscFVSetType(fv,PETSCFVLEASTSQUARES);
122: PetscFVSetNumComponents(fv,1);
123: PetscLimiterCreate(comm,&limiter);
124: PetscLimiterSetType(limiter,PETSCLIMITERNONE);
125: PetscFVSetLimiter(fv,limiter);
126: PetscLimiterDestroy(&limiter);
127: PetscFVSetFromOptions(fv);
128: DMSetField(base,0,(PetscObject)fv);
129: PetscFVDestroy(&fv);
130: } else {
131: PetscFE fe;
132: PetscFECreateDefault(base,dim,1,PETSC_FALSE,NULL,PETSC_DEFAULT,&fe);
133: DMSetField(base,0,(PetscObject)fe);
134: PetscFEDestroy(&fe);
135: }
136: {
137: PetscDS prob;
138: PetscInt comps[] = {0};
139: PetscInt ids[] = {1, 2, 3, 4, 5, 6};
141: DMGetDS(base,&prob);
142: PetscDSAddBoundary(prob,PETSC_TRUE, "bc", "marker", 0, 1, comps, useFV ? (void(*)()) bc_func_fv : (void(*)()) funcs[0], 2 * dim, ids, useFV ? (void *) &bcCtx : NULL);
143: }
144: AddIdentityLabel(base);
145: DMViewFromOptions(base,NULL,"-dm_base_view");
147: /* the pre adaptivity forest */
148: DMCreate(comm,&preForest);
149: DMSetType(preForest,(dim == 2) ? DMP4EST : DMP8EST);
150: DMGetDS(base,&ds);
151: DMSetDS(preForest,ds);
152: DMForestSetBaseDM(preForest,base);
153: DMForestSetMinimumRefinement(preForest,1);
154: DMForestSetInitialRefinement(preForest,1);
155: DMSetFromOptions(preForest);
156: DMSetUp(preForest);
157: DMViewFromOptions(preForest,NULL,"-dm_pre_view");
159: /* the pre adaptivity field */
160: DMCreateGlobalVector(preForest,&preVec);
161: DMProjectFunction(preForest,0.,funcs,ctxs,INSERT_VALUES,preVec);
162: VecViewFromOptions(preVec,NULL,"-vec_pre_view");
164: PetscObjectGetReference((PetscObject)preForest,&preCount);
166: /* adapt */
167: CreateAdaptivityLabel(preForest,&adaptLabel);
168: DMForestTemplate(preForest,comm,&postForest);
169: DMForestSetMinimumRefinement(postForest,0);
170: DMForestSetInitialRefinement(postForest,0);
171: DMForestSetAdaptivityLabel(postForest,adaptLabel);
172: DMLabelDestroy(&adaptLabel);
173: DMSetUp(postForest);
174: DMViewFromOptions(postForest,NULL,"-dm_post_view");
176: /* transfer */
177: DMCreateGlobalVector(postForest,&postVecTransfer);
178: DMForestTransferVec(preForest,preVec,postForest,postVecTransfer,PETSC_TRUE,0.0);
179: VecViewFromOptions(postVecTransfer,NULL,"-vec_post_transfer_view");
181: /* the exact post adaptivity field */
182: DMCreateGlobalVector(postForest,&postVecExact);
183: DMProjectFunction(postForest,0.,funcs,ctxs,INSERT_VALUES,postVecExact);
184: VecViewFromOptions(postVecExact,NULL,"-vec_post_exact_view");
186: /* compare */
187: VecAXPY(postVecTransfer,-1.,postVecExact);
188: VecViewFromOptions(postVecTransfer,NULL,"-vec_diff_view");
189: VecNorm(postVecTransfer,NORM_2,&diff);
191: /* output */
192: if (diff < tol) {
193: PetscPrintf(comm,"DMForestTransferVec() passes.\n");
194: } else {
195: PetscPrintf(comm,"DMForestTransferVec() fails with error %g and tolerance %g\n",diff,tol);
196: }
198: /* disconnect preForest from postForest */
199: DMForestSetAdaptivityForest(postForest,NULL);
200: PetscObjectGetReference((PetscObject)preForest,&postCount);
201: if (postCount != preCount) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Adaptation not memory neutral: reference count increase from %d to %d\n",preCount,postCount);
203: /* cleanup */
204: VecDestroy(&postVecExact);
205: VecDestroy(&postVecTransfer);
206: DMDestroy(&postForest);
207: VecDestroy(&preVec);
208: DMDestroy(&preForest);
209: DMDestroy(&base);
210: PetscFinalize();
211: return ierr;
212: }