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"; 4 5 #include <petscdmplex.h> 6 #include <petscsnes.h> 7 #include <petscds.h> 8 #include <petscconvest.h> 9 10 typedef struct { 11 PetscInt dummy; 12 } AppCtx; 13 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 } 19 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 } 25 26 static void f0_quad_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 27 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 28 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 29 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 30 { 31 PetscInt d; 32 for (d = 0; d < dim; ++d) f0[0] += 2.0; 33 } 34 35 static void f0_quad_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, 36 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 37 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 38 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 39 { 40 f0[0] = u[uOff[1]] - 2.0; 41 } 42 43 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 44 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 45 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 46 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 47 { 48 PetscInt d; 49 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 50 } 51 52 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 53 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 54 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 55 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 56 { 57 PetscInt d; 58 for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0; 59 } 60 61 static void g0_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 62 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 63 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 64 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 65 { 66 g0[0] = 1.0; 67 } 68 69 static PetscErrorCode DivideDomain(DM dm, AppCtx *user) 70 { 71 DMLabel top, bottom; 72 PetscReal low[3], high[3], midy; 73 PetscInt cStart, cEnd, c; 74 PetscErrorCode ierr; 75 76 PetscFunctionBeginUser; 77 ierr = DMCreateLabel(dm, "top");CHKERRQ(ierr); 78 ierr = DMCreateLabel(dm, "bottom");CHKERRQ(ierr); 79 ierr = DMGetLabel(dm, "top", &top);CHKERRQ(ierr); 80 ierr = DMGetLabel(dm, "bottom", &bottom);CHKERRQ(ierr); 81 ierr = DMGetBoundingBox(dm, low, high);CHKERRQ(ierr); 82 midy = 0.5*(high[1] - low[1]); 83 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 84 for (c = cStart; c < cEnd; ++c) { 85 PetscReal centroid[3]; 86 87 ierr = DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);CHKERRQ(ierr); 88 if (centroid[1] > midy) {ierr = DMLabelSetValue(top, c, 1);CHKERRQ(ierr);} 89 else {ierr = DMLabelSetValue(bottom, c, 1);CHKERRQ(ierr);} 90 } 91 ierr = DMPlexLabelComplete(dm, top);CHKERRQ(ierr); 92 PetscFunctionReturn(0); 93 } 94 95 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 96 { 97 PetscErrorCode ierr; 98 99 PetscFunctionBeginUser; 100 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 101 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 102 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 103 ierr = DivideDomain(*dm, user);CHKERRQ(ierr); 104 ierr = DMSetApplicationContext(*dm, user);CHKERRQ(ierr); 105 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 106 PetscFunctionReturn(0); 107 } 108 109 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 110 { 111 PetscDS ds; 112 PetscWeakForm wf; 113 DMLabel label; 114 const PetscInt id = 1; 115 PetscErrorCode ierr; 116 117 PetscFunctionBeginUser; 118 ierr = DMGetRegionNumDS(dm, 0, &label, NULL, &ds);CHKERRQ(ierr); 119 ierr = PetscDSGetWeakForm(ds, &wf);CHKERRQ(ierr); 120 ierr = PetscWeakFormSetIndexResidual(wf, label, 1, 0, 0, 0, f0_quad_u, 0, f1_u);CHKERRQ(ierr); 121 ierr = PetscWeakFormSetIndexJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu);CHKERRQ(ierr); 122 ierr = PetscDSSetExactSolution(ds, 0, quad_u, user);CHKERRQ(ierr); 123 ierr = DMGetRegionNumDS(dm, 1, &label, NULL, &ds);CHKERRQ(ierr); 124 ierr = PetscDSGetWeakForm(ds, &wf);CHKERRQ(ierr); 125 ierr = PetscWeakFormSetIndexResidual(wf, label, 1, 0, 0, 0, f0_quad_u, 0, f1_u);CHKERRQ(ierr); 126 ierr = PetscWeakFormSetIndexJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu);CHKERRQ(ierr); 127 ierr = PetscWeakFormSetIndexResidual(wf, label, 1, 1, 0, 0, f0_quad_p, 0, NULL);CHKERRQ(ierr); 128 ierr = PetscWeakFormSetIndexJacobian(wf, label, 1, 1, 1, 0, 0, g0_pp, 0, NULL, 0, NULL, 0, NULL);CHKERRQ(ierr); 129 ierr = PetscDSSetExactSolution(ds, 0, quad_u, user);CHKERRQ(ierr); 130 ierr = PetscDSSetExactSolution(ds, 1, quad_p, user);CHKERRQ(ierr); 131 ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr); 132 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) quad_u, NULL, user, NULL);CHKERRQ(ierr); 133 PetscFunctionReturn(0); 134 } 135 136 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) 137 { 138 DM cdm = dm; 139 DMLabel top; 140 PetscFE fe, feTop; 141 PetscQuadrature q; 142 PetscInt dim; 143 PetscBool simplex; 144 const char *nameTop = "pressure"; 145 char prefix[PETSC_MAX_PATH_LEN]; 146 PetscErrorCode ierr; 147 148 PetscFunctionBeginUser; 149 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 150 ierr = DMPlexIsSimplex(dm, &simplex);CHKERRQ(ierr); 151 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);CHKERRQ(ierr); 152 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, simplex, name ? prefix : NULL, -1, &fe);CHKERRQ(ierr); 153 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 154 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 155 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 156 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 157 ierr = DMGetLabel(dm, "top", &top);CHKERRQ(ierr); 158 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_top_", nameTop);CHKERRQ(ierr); 159 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, simplex, name ? prefix : NULL, -1, &feTop);CHKERRQ(ierr); 160 ierr = PetscObjectSetName((PetscObject) feTop, nameTop);CHKERRQ(ierr); 161 ierr = PetscFESetQuadrature(feTop, q);CHKERRQ(ierr); 162 ierr = DMSetField(dm, 1, top, (PetscObject) feTop);CHKERRQ(ierr); 163 ierr = PetscFEDestroy(&feTop);CHKERRQ(ierr); 164 ierr = DMCreateDS(dm);CHKERRQ(ierr); 165 ierr = (*setup)(dm, user);CHKERRQ(ierr); 166 while (cdm) { 167 ierr = DMCopyDisc(dm,cdm);CHKERRQ(ierr); 168 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 169 } 170 PetscFunctionReturn(0); 171 } 172 173 int main(int argc, char **argv) 174 { 175 DM dm; /* Problem specification */ 176 SNES snes; /* Nonlinear solver */ 177 Vec u; /* Solutions */ 178 AppCtx user; /* User-defined work context */ 179 PetscErrorCode ierr; 180 181 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 182 /* Primal system */ 183 ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr); 184 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 185 ierr = SNESSetDM(snes, dm);CHKERRQ(ierr); 186 ierr = SetupDiscretization(dm, "potential", SetupPrimalProblem, &user);CHKERRQ(ierr); 187 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 188 ierr = VecSet(u, 0.0);CHKERRQ(ierr); 189 ierr = PetscObjectSetName((PetscObject) u, "solution");CHKERRQ(ierr); 190 ierr = DMPlexSetSNESLocalFEM(dm, &user, &user, &user);CHKERRQ(ierr); 191 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 192 ierr = DMSNESCheckFromOptions(snes, u);CHKERRQ(ierr); 193 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 194 ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr); 195 ierr = VecViewFromOptions(u, NULL, "-sol_view");CHKERRQ(ierr); 196 /* Cleanup */ 197 ierr = VecDestroy(&u);CHKERRQ(ierr); 198 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 199 ierr = DMDestroy(&dm);CHKERRQ(ierr); 200 ierr = PetscFinalize(); 201 return ierr; 202 } 203 204 /*TEST 205 206 test: 207 suffix: 2d_p1_0 208 requires: triangle 209 args: -potential_petscspace_degree 2 -pressure_top_petscspace_degree 0 -dm_refine 0 -dmsnes_check 210 211 test: 212 suffix: 2d_p1_1 213 requires: triangle 214 args: -potential_petscspace_degree 1 -pressure_top_petscspace_degree 0 -dm_refine 0 -convest_num_refine 3 -snes_convergence_estimate 215 216 TEST*/ 217