Actual source code: ex23.c
1: static char help[] = "Poisson Problem with a split domain.\n\
2: We solve the Poisson problem in two halves of a domain.\n\
3: In one half, we include an additional field.\n\n\n";
5: #include <petscdmplex.h>
6: #include <petscsnes.h>
7: #include <petscds.h>
8: #include <petscconvest.h>
10: typedef struct {
11: PetscInt dummy;
12: } AppCtx;
14: static PetscErrorCode quad_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
15: {
16: u[0] = x[0] * x[0] + x[1] * x[1];
17: return 0;
18: }
20: static PetscErrorCode quad_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
21: {
22: u[0] = 2.0;
23: return 0;
24: }
26: static void f0_quad_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
27: {
28: PetscInt d;
29: for (d = 0; d < dim; ++d) f0[0] += 2.0;
30: }
32: static void f0_quad_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
33: {
34: f0[0] = u[uOff[1]] - 2.0;
35: }
37: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
38: {
39: PetscInt d;
40: for (d = 0; d < dim; ++d) f1[d] = u_x[d];
41: }
43: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
44: {
45: PetscInt d;
46: for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
47: }
49: static void g0_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
50: {
51: g0[0] = 1.0;
52: }
54: static PetscErrorCode DivideDomain(DM dm, AppCtx *user)
55: {
56: DMLabel top, bottom;
57: PetscReal low[3], high[3], midy;
58: PetscInt cStart, cEnd, c;
61: DMCreateLabel(dm, "top");
62: DMCreateLabel(dm, "bottom");
63: DMGetLabel(dm, "top", &top);
64: DMGetLabel(dm, "bottom", &bottom);
65: DMGetCoordinatesLocalSetUp(dm);
66: DMGetBoundingBox(dm, low, high);
67: midy = 0.5 * (high[1] - low[1]);
68: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
69: for (c = cStart; c < cEnd; ++c) {
70: PetscReal centroid[3];
72: DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
73: if (centroid[1] > midy) DMLabelSetValue(top, c, 1);
74: else DMLabelSetValue(bottom, c, 1);
75: }
76: DMPlexLabelComplete(dm, top);
77: return 0;
78: }
80: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
81: {
83: DMCreate(comm, dm);
84: DMSetType(*dm, DMPLEX);
85: DMSetFromOptions(*dm);
86: DivideDomain(*dm, user);
87: DMSetApplicationContext(*dm, user);
88: DMViewFromOptions(*dm, NULL, "-dm_view");
89: return 0;
90: }
92: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
93: {
94: PetscDS ds;
95: PetscWeakForm wf;
96: DMLabel label;
97: const PetscInt id = 1;
100: DMGetRegionNumDS(dm, 0, &label, NULL, &ds);
101: PetscDSGetWeakForm(ds, &wf);
102: PetscWeakFormSetIndexResidual(wf, label, 1, 0, 0, 0, f0_quad_u, 0, f1_u);
103: PetscWeakFormSetIndexJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu);
104: PetscDSSetExactSolution(ds, 0, quad_u, user);
105: DMGetRegionNumDS(dm, 1, &label, NULL, &ds);
106: PetscDSGetWeakForm(ds, &wf);
107: PetscWeakFormSetIndexResidual(wf, label, 1, 0, 0, 0, f0_quad_u, 0, f1_u);
108: PetscWeakFormSetIndexJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu);
109: PetscWeakFormSetIndexResidual(wf, label, 1, 1, 0, 0, f0_quad_p, 0, NULL);
110: PetscWeakFormSetIndexJacobian(wf, label, 1, 1, 1, 0, 0, g0_pp, 0, NULL, 0, NULL, 0, NULL);
111: PetscDSSetExactSolution(ds, 0, quad_u, user);
112: PetscDSSetExactSolution(ds, 1, quad_p, user);
113: DMGetLabel(dm, "marker", &label);
114: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))quad_u, NULL, user, NULL);
115: return 0;
116: }
118: static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
119: {
120: DM cdm = dm;
121: DMLabel top;
122: PetscFE fe, feTop;
123: PetscQuadrature q;
124: PetscInt dim;
125: PetscBool simplex;
126: const char *nameTop = "pressure";
127: char prefix[PETSC_MAX_PATH_LEN];
130: DMGetDimension(dm, &dim);
131: DMPlexIsSimplex(dm, &simplex);
132: PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);
133: PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, name ? prefix : NULL, -1, &fe);
134: PetscObjectSetName((PetscObject)fe, name);
135: DMSetField(dm, 0, NULL, (PetscObject)fe);
136: PetscFEGetQuadrature(fe, &q);
137: PetscFEDestroy(&fe);
138: DMGetLabel(dm, "top", &top);
139: PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_top_", nameTop);
140: PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, name ? prefix : NULL, -1, &feTop);
141: PetscObjectSetName((PetscObject)feTop, nameTop);
142: PetscFESetQuadrature(feTop, q);
143: DMSetField(dm, 1, top, (PetscObject)feTop);
144: PetscFEDestroy(&feTop);
145: DMCreateDS(dm);
146: (*setup)(dm, user);
147: while (cdm) {
148: DMCopyDisc(dm, cdm);
149: DMGetCoarseDM(cdm, &cdm);
150: }
151: return 0;
152: }
154: int main(int argc, char **argv)
155: {
156: DM dm; /* Problem specification */
157: SNES snes; /* Nonlinear solver */
158: Vec u; /* Solutions */
159: AppCtx user; /* User-defined work context */
162: PetscInitialize(&argc, &argv, NULL, help);
163: /* Primal system */
164: SNESCreate(PETSC_COMM_WORLD, &snes);
165: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
166: SNESSetDM(snes, dm);
167: SetupDiscretization(dm, "potential", SetupPrimalProblem, &user);
168: DMCreateGlobalVector(dm, &u);
169: VecSet(u, 0.0);
170: PetscObjectSetName((PetscObject)u, "solution");
171: DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
172: SNESSetFromOptions(snes);
173: DMSNESCheckFromOptions(snes, u);
174: SNESSolve(snes, NULL, u);
175: SNESGetSolution(snes, &u);
176: VecViewFromOptions(u, NULL, "-sol_view");
177: /* Cleanup */
178: VecDestroy(&u);
179: SNESDestroy(&snes);
180: DMDestroy(&dm);
181: PetscFinalize();
182: return 0;
183: }
185: /*TEST
187: test:
188: suffix: 2d_p1_0
189: requires: triangle
190: args: -potential_petscspace_degree 2 -pressure_top_petscspace_degree 0 -dm_refine 0 -dmsnes_check
192: test:
193: suffix: 2d_p1_1
194: requires: triangle
195: args: -potential_petscspace_degree 1 -pressure_top_petscspace_degree 0 -dm_refine 0 -convest_num_refine 3 -snes_convergence_estimate
197: TEST*/