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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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;

60:   PetscFunctionBeginUser;
61:   PetscCall(DMCreateLabel(dm, "top"));
62:   PetscCall(DMCreateLabel(dm, "bottom"));
63:   PetscCall(DMGetLabel(dm, "top", &top));
64:   PetscCall(DMGetLabel(dm, "bottom", &bottom));
65:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
66:   PetscCall(DMGetBoundingBox(dm, low, high));
67:   midy = 0.5 * (high[1] - low[1]);
68:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
69:   for (c = cStart; c < cEnd; ++c) {
70:     PetscReal centroid[3];

72:     PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL));
73:     if (centroid[1] > midy) PetscCall(DMLabelSetValue(top, c, 1));
74:     else PetscCall(DMLabelSetValue(bottom, c, 1));
75:   }
76:   PetscCall(DMPlexLabelComplete(dm, top));
77:   PetscFunctionReturn(PETSC_SUCCESS);
78: }

80: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
81: {
82:   PetscFunctionBeginUser;
83:   PetscCall(DMCreate(comm, dm));
84:   PetscCall(DMSetType(*dm, DMPLEX));
85:   PetscCall(DMSetFromOptions(*dm));
86:   PetscCall(DivideDomain(*dm, user));
87:   PetscCall(DMSetApplicationContext(*dm, user));
88:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
89:   PetscFunctionReturn(PETSC_SUCCESS);
90: }

92: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
93: {
94:   PetscDS        ds;
95:   PetscWeakForm  wf;
96:   DMLabel        label;
97:   const PetscInt id = 1;

99:   PetscFunctionBeginUser;
100:   PetscCall(DMGetRegionNumDS(dm, 0, &label, NULL, &ds, NULL));
101:   PetscCall(PetscDSGetWeakForm(ds, &wf));
102:   PetscCall(PetscWeakFormSetIndexResidual(wf, label, 1, 0, 0, 0, f0_quad_u, 0, f1_u));
103:   PetscCall(PetscWeakFormSetIndexJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu));
105:   PetscCall(DMGetRegionNumDS(dm, 1, &label, NULL, &ds, NULL));
106:   PetscCall(PetscDSGetWeakForm(ds, &wf));
107:   PetscCall(PetscWeakFormSetIndexResidual(wf, label, 1, 0, 0, 0, f0_quad_u, 0, f1_u));
108:   PetscCall(PetscWeakFormSetIndexJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu));
109:   PetscCall(PetscWeakFormSetIndexResidual(wf, label, 1, 1, 0, 0, f0_quad_p, 0, NULL));
110:   PetscCall(PetscWeakFormSetIndexJacobian(wf, label, 1, 1, 1, 0, 0, g0_pp, 0, NULL, 0, NULL, 0, NULL));
113:   PetscCall(DMGetLabel(dm, "marker", &label));
114:   PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))quad_u, NULL, user, NULL));
115:   PetscFunctionReturn(PETSC_SUCCESS);
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;
124:   PetscInt        dim;
125:   PetscBool       simplex;
126:   const char     *nameTop = "pressure";
127:   char            prefix[PETSC_MAX_PATH_LEN];

129:   PetscFunctionBeginUser;
130:   PetscCall(DMGetDimension(dm, &dim));
131:   PetscCall(DMPlexIsSimplex(dm, &simplex));
132:   PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name));
133:   PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, name ? prefix : NULL, -1, &fe));
134:   PetscCall(PetscObjectSetName((PetscObject)fe, name));
135:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
137:   PetscCall(PetscFEDestroy(&fe));
138:   PetscCall(DMGetLabel(dm, "top", &top));
139:   PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_top_", nameTop));
140:   PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, simplex, name ? prefix : NULL, -1, &feTop));
141:   PetscCall(PetscObjectSetName((PetscObject)feTop, nameTop));
143:   PetscCall(DMSetField(dm, 1, top, (PetscObject)feTop));
144:   PetscCall(PetscFEDestroy(&feTop));
145:   PetscCall(DMCreateDS(dm));
146:   PetscCall((*setup)(dm, user));
147:   while (cdm) {
148:     PetscCall(DMCopyDisc(dm, cdm));
149:     PetscCall(DMGetCoarseDM(cdm, &cdm));
150:   }
151:   PetscFunctionReturn(PETSC_SUCCESS);
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 */

161:   PetscFunctionBeginUser;
162:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
163:   /* Primal system */
164:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
165:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
166:   PetscCall(SNESSetDM(snes, dm));
167:   PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user));
168:   PetscCall(DMCreateGlobalVector(dm, &u));
169:   PetscCall(VecSet(u, 0.0));
170:   PetscCall(PetscObjectSetName((PetscObject)u, "solution"));
171:   PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
172:   PetscCall(SNESSetFromOptions(snes));
173:   PetscCall(DMSNESCheckFromOptions(snes, u));
174:   PetscCall(SNESSolve(snes, NULL, u));
175:   PetscCall(SNESGetSolution(snes, &u));
176:   PetscCall(VecViewFromOptions(u, NULL, "-sol_view"));
177:   /* Cleanup */
178:   PetscCall(VecDestroy(&u));
179:   PetscCall(SNESDestroy(&snes));
180:   PetscCall(DMDestroy(&dm));
181:   PetscCall(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*/
```