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 /* Domain and mesh definition */ 12 PetscInt dim; /* The topological mesh dimension */ 13 PetscBool simplex; /* Simplicial mesh */ 14 } AppCtx; 15 16 static PetscErrorCode quad_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 17 { 18 u[0] = x[0]*x[0] + x[1]*x[1]; 19 return 0; 20 } 21 22 static PetscErrorCode quad_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 23 { 24 u[0] = 2.0; 25 return 0; 26 } 27 28 static void f0_quad_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 29 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 30 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 31 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 32 { 33 PetscInt d; 34 for (d = 0; d < dim; ++d) f0[0] += 2.0; 35 } 36 37 static void f0_quad_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, 38 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 39 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 40 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 41 { 42 f0[0] = u[uOff[1]] - 2.0; 43 } 44 45 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 46 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 47 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 48 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 49 { 50 PetscInt d; 51 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 52 } 53 54 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 55 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 56 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 57 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 58 { 59 PetscInt d; 60 for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0; 61 } 62 63 static void g0_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 64 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 65 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 66 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 67 { 68 g0[0] = 1.0; 69 } 70 71 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 72 { 73 PetscErrorCode ierr; 74 75 PetscFunctionBeginUser; 76 options->dim = 2; 77 options->simplex = PETSC_TRUE; 78 79 ierr = PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");CHKERRQ(ierr); 80 ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex23.c", options->dim, &options->dim, NULL);CHKERRQ(ierr); 81 ierr = PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex23.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr); 82 ierr = PetscOptionsEnd(); 83 PetscFunctionReturn(0); 84 } 85 86 static PetscErrorCode DivideDomain(DM dm, AppCtx *user) 87 { 88 DMLabel top, bottom; 89 PetscReal low[3], high[3], midy; 90 PetscInt cStart, cEnd, c; 91 PetscErrorCode ierr; 92 93 PetscFunctionBeginUser; 94 ierr = DMCreateLabel(dm, "top");CHKERRQ(ierr); 95 ierr = DMCreateLabel(dm, "bottom");CHKERRQ(ierr); 96 ierr = DMGetLabel(dm, "top", &top);CHKERRQ(ierr); 97 ierr = DMGetLabel(dm, "bottom", &bottom);CHKERRQ(ierr); 98 ierr = DMGetBoundingBox(dm, low, high);CHKERRQ(ierr); 99 midy = 0.5*(high[1] - low[1]); 100 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 101 for (c = cStart; c < cEnd; ++c) { 102 PetscReal centroid[3]; 103 104 ierr = DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);CHKERRQ(ierr); 105 if (centroid[1] > midy) {ierr = DMLabelSetValue(top, c, 1);CHKERRQ(ierr);} 106 else {ierr = DMLabelSetValue(bottom, c, 1);CHKERRQ(ierr);} 107 } 108 ierr = DMPlexLabelComplete(dm, top);CHKERRQ(ierr); 109 PetscFunctionReturn(0); 110 } 111 112 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 113 { 114 PetscErrorCode ierr; 115 116 PetscFunctionBeginUser; 117 /* Create box mesh */ 118 ierr = DMPlexCreateBoxMesh(comm, user->dim, user->simplex, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 119 /* Distribute mesh over processes */ 120 { 121 DM dmDist = NULL; 122 PetscPartitioner part; 123 124 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr); 125 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 126 ierr = DMPlexDistribute(*dm, 0, NULL, &dmDist);CHKERRQ(ierr); 127 if (dmDist) { 128 ierr = DMDestroy(dm);CHKERRQ(ierr); 129 *dm = dmDist; 130 } 131 } 132 /* TODO: This should be pulled into the library */ 133 { 134 char convType[256]; 135 PetscBool flg; 136 137 ierr = PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");CHKERRQ(ierr); 138 ierr = PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex12",DMList,DMPLEX,convType,256,&flg);CHKERRQ(ierr); 139 ierr = PetscOptionsEnd(); 140 if (flg) { 141 DM dmConv; 142 143 ierr = DMConvert(*dm,convType,&dmConv);CHKERRQ(ierr); 144 if (dmConv) { 145 ierr = DMDestroy(dm);CHKERRQ(ierr); 146 *dm = dmConv; 147 } 148 } 149 } 150 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 151 152 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr); 153 ierr = DMSetApplicationContext(*dm, user);CHKERRQ(ierr); 154 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 155 ierr = DivideDomain(*dm, user);CHKERRQ(ierr); 156 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 157 PetscFunctionReturn(0); 158 } 159 160 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 161 { 162 PetscDS ds; 163 PetscWeakForm wf; 164 DMLabel label; 165 const PetscInt id = 1; 166 PetscErrorCode ierr; 167 168 PetscFunctionBeginUser; 169 ierr = DMGetRegionNumDS(dm, 0, &label, NULL, &ds);CHKERRQ(ierr); 170 ierr = PetscDSGetWeakForm(ds, &wf);CHKERRQ(ierr); 171 ierr = PetscWeakFormSetIndexResidual(wf, label, 1, 0, 0, f0_quad_u, 0, f1_u);CHKERRQ(ierr); 172 ierr = PetscWeakFormSetIndexJacobian(wf, label, 1, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu);CHKERRQ(ierr); 173 ierr = PetscDSSetExactSolution(ds, 0, quad_u, user);CHKERRQ(ierr); 174 ierr = DMGetRegionNumDS(dm, 1, &label, NULL, &ds);CHKERRQ(ierr); 175 ierr = PetscDSGetWeakForm(ds, &wf);CHKERRQ(ierr); 176 ierr = PetscWeakFormSetIndexResidual(wf, label, 1, 0, 0, f0_quad_u, 0, f1_u);CHKERRQ(ierr); 177 ierr = PetscWeakFormSetIndexJacobian(wf, label, 1, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu);CHKERRQ(ierr); 178 ierr = PetscWeakFormSetIndexResidual(wf, label, 1, 1, 0, f0_quad_p, 0, NULL);CHKERRQ(ierr); 179 ierr = PetscWeakFormSetIndexJacobian(wf, label, 1, 1, 1, 0, g0_pp, 0, NULL, 0, NULL, 0, NULL);CHKERRQ(ierr); 180 ierr = PetscDSSetExactSolution(ds, 0, quad_u, user);CHKERRQ(ierr); 181 ierr = PetscDSSetExactSolution(ds, 1, quad_p, user);CHKERRQ(ierr); 182 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) quad_u, NULL, 1, &id, user);CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185 186 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) 187 { 188 DM cdm = dm; 189 DMLabel top; 190 PetscFE fe, feTop; 191 PetscQuadrature q; 192 const char *nameTop = "pressure"; 193 char prefix[PETSC_MAX_PATH_LEN]; 194 PetscErrorCode ierr; 195 196 PetscFunctionBeginUser; 197 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);CHKERRQ(ierr); 198 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), user->dim, 1, user->simplex, name ? prefix : NULL, -1, &fe);CHKERRQ(ierr); 199 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 200 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 201 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 202 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 203 ierr = DMGetLabel(dm, "top", &top);CHKERRQ(ierr); 204 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_top_", nameTop);CHKERRQ(ierr); 205 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), user->dim, 1, user->simplex, name ? prefix : NULL, -1, &feTop);CHKERRQ(ierr); 206 ierr = PetscObjectSetName((PetscObject) feTop, nameTop);CHKERRQ(ierr); 207 ierr = PetscFESetQuadrature(feTop, q);CHKERRQ(ierr); 208 ierr = DMSetField(dm, 1, top, (PetscObject) feTop);CHKERRQ(ierr); 209 ierr = PetscFEDestroy(&feTop);CHKERRQ(ierr); 210 ierr = DMCreateDS(dm);CHKERRQ(ierr); 211 ierr = (*setup)(dm, user);CHKERRQ(ierr); 212 while (cdm) { 213 ierr = DMCopyDisc(dm,cdm);CHKERRQ(ierr); 214 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 215 } 216 PetscFunctionReturn(0); 217 } 218 219 int main(int argc, char **argv) 220 { 221 DM dm; /* Problem specification */ 222 SNES snes; /* Nonlinear solver */ 223 Vec u; /* Solutions */ 224 AppCtx user; /* User-defined work context */ 225 PetscErrorCode ierr; 226 227 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 228 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 229 /* Primal system */ 230 ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr); 231 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 232 ierr = SNESSetDM(snes, dm);CHKERRQ(ierr); 233 ierr = SetupDiscretization(dm, "potential", SetupPrimalProblem, &user);CHKERRQ(ierr); 234 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 235 ierr = VecSet(u, 0.0);CHKERRQ(ierr); 236 ierr = PetscObjectSetName((PetscObject) u, "solution");CHKERRQ(ierr); 237 ierr = DMPlexSetSNESLocalFEM(dm, &user, &user, &user);CHKERRQ(ierr); 238 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 239 ierr = DMSNESCheckFromOptions(snes, u);CHKERRQ(ierr); 240 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 241 ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr); 242 ierr = VecViewFromOptions(u, NULL, "-sol_view");CHKERRQ(ierr); 243 /* Cleanup */ 244 ierr = VecDestroy(&u);CHKERRQ(ierr); 245 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 246 ierr = DMDestroy(&dm);CHKERRQ(ierr); 247 ierr = PetscFinalize(); 248 return ierr; 249 } 250 251 /*TEST 252 253 test: 254 suffix: 2d_p1_0 255 requires: triangle 256 args: -potential_petscspace_degree 2 -pressure_top_petscspace_degree 0 -dm_refine 0 -dmsnes_check 257 258 test: 259 suffix: 2d_p1_1 260 requires: triangle 261 args: -potential_petscspace_degree 1 -pressure_top_petscspace_degree 0 -dm_refine 0 -convest_num_refine 3 -snes_convergence_estimate 262 263 TEST*/ 264