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 75 PetscFunctionBeginUser; 76 CHKERRQ(DMCreateLabel(dm, "top")); 77 CHKERRQ(DMCreateLabel(dm, "bottom")); 78 CHKERRQ(DMGetLabel(dm, "top", &top)); 79 CHKERRQ(DMGetLabel(dm, "bottom", &bottom)); 80 CHKERRQ(DMGetBoundingBox(dm, low, high)); 81 midy = 0.5*(high[1] - low[1]); 82 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 83 for (c = cStart; c < cEnd; ++c) { 84 PetscReal centroid[3]; 85 86 CHKERRQ(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL)); 87 if (centroid[1] > midy) CHKERRQ(DMLabelSetValue(top, c, 1)); 88 else CHKERRQ(DMLabelSetValue(bottom, c, 1)); 89 } 90 CHKERRQ(DMPlexLabelComplete(dm, top)); 91 PetscFunctionReturn(0); 92 } 93 94 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 95 { 96 PetscFunctionBeginUser; 97 CHKERRQ(DMCreate(comm, dm)); 98 CHKERRQ(DMSetType(*dm, DMPLEX)); 99 CHKERRQ(DMSetFromOptions(*dm)); 100 CHKERRQ(DivideDomain(*dm, user)); 101 CHKERRQ(DMSetApplicationContext(*dm, user)); 102 CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view")); 103 PetscFunctionReturn(0); 104 } 105 106 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 107 { 108 PetscDS ds; 109 PetscWeakForm wf; 110 DMLabel label; 111 const PetscInt id = 1; 112 113 PetscFunctionBeginUser; 114 CHKERRQ(DMGetRegionNumDS(dm, 0, &label, NULL, &ds)); 115 CHKERRQ(PetscDSGetWeakForm(ds, &wf)); 116 CHKERRQ(PetscWeakFormSetIndexResidual(wf, label, 1, 0, 0, 0, f0_quad_u, 0, f1_u)); 117 CHKERRQ(PetscWeakFormSetIndexJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu)); 118 CHKERRQ(PetscDSSetExactSolution(ds, 0, quad_u, user)); 119 CHKERRQ(DMGetRegionNumDS(dm, 1, &label, NULL, &ds)); 120 CHKERRQ(PetscDSGetWeakForm(ds, &wf)); 121 CHKERRQ(PetscWeakFormSetIndexResidual(wf, label, 1, 0, 0, 0, f0_quad_u, 0, f1_u)); 122 CHKERRQ(PetscWeakFormSetIndexJacobian(wf, label, 1, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu)); 123 CHKERRQ(PetscWeakFormSetIndexResidual(wf, label, 1, 1, 0, 0, f0_quad_p, 0, NULL)); 124 CHKERRQ(PetscWeakFormSetIndexJacobian(wf, label, 1, 1, 1, 0, 0, g0_pp, 0, NULL, 0, NULL, 0, NULL)); 125 CHKERRQ(PetscDSSetExactSolution(ds, 0, quad_u, user)); 126 CHKERRQ(PetscDSSetExactSolution(ds, 1, quad_p, user)); 127 CHKERRQ(DMGetLabel(dm, "marker", &label)); 128 CHKERRQ(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) quad_u, NULL, user, NULL)); 129 PetscFunctionReturn(0); 130 } 131 132 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) 133 { 134 DM cdm = dm; 135 DMLabel top; 136 PetscFE fe, feTop; 137 PetscQuadrature q; 138 PetscInt dim; 139 PetscBool simplex; 140 const char *nameTop = "pressure"; 141 char prefix[PETSC_MAX_PATH_LEN]; 142 143 PetscFunctionBeginUser; 144 CHKERRQ(DMGetDimension(dm, &dim)); 145 CHKERRQ(DMPlexIsSimplex(dm, &simplex)); 146 CHKERRQ(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name)); 147 CHKERRQ(PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, simplex, name ? prefix : NULL, -1, &fe)); 148 CHKERRQ(PetscObjectSetName((PetscObject) fe, name)); 149 CHKERRQ(DMSetField(dm, 0, NULL, (PetscObject) fe)); 150 CHKERRQ(PetscFEGetQuadrature(fe, &q)); 151 CHKERRQ(PetscFEDestroy(&fe)); 152 CHKERRQ(DMGetLabel(dm, "top", &top)); 153 CHKERRQ(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_top_", nameTop)); 154 CHKERRQ(PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, simplex, name ? prefix : NULL, -1, &feTop)); 155 CHKERRQ(PetscObjectSetName((PetscObject) feTop, nameTop)); 156 CHKERRQ(PetscFESetQuadrature(feTop, q)); 157 CHKERRQ(DMSetField(dm, 1, top, (PetscObject) feTop)); 158 CHKERRQ(PetscFEDestroy(&feTop)); 159 CHKERRQ(DMCreateDS(dm)); 160 CHKERRQ((*setup)(dm, user)); 161 while (cdm) { 162 CHKERRQ(DMCopyDisc(dm,cdm)); 163 CHKERRQ(DMGetCoarseDM(cdm, &cdm)); 164 } 165 PetscFunctionReturn(0); 166 } 167 168 int main(int argc, char **argv) 169 { 170 DM dm; /* Problem specification */ 171 SNES snes; /* Nonlinear solver */ 172 Vec u; /* Solutions */ 173 AppCtx user; /* User-defined work context */ 174 PetscErrorCode ierr; 175 176 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 177 /* Primal system */ 178 CHKERRQ(SNESCreate(PETSC_COMM_WORLD, &snes)); 179 CHKERRQ(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 180 CHKERRQ(SNESSetDM(snes, dm)); 181 CHKERRQ(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user)); 182 CHKERRQ(DMCreateGlobalVector(dm, &u)); 183 CHKERRQ(VecSet(u, 0.0)); 184 CHKERRQ(PetscObjectSetName((PetscObject) u, "solution")); 185 CHKERRQ(DMPlexSetSNESLocalFEM(dm, &user, &user, &user)); 186 CHKERRQ(SNESSetFromOptions(snes)); 187 CHKERRQ(DMSNESCheckFromOptions(snes, u)); 188 CHKERRQ(SNESSolve(snes, NULL, u)); 189 CHKERRQ(SNESGetSolution(snes, &u)); 190 CHKERRQ(VecViewFromOptions(u, NULL, "-sol_view")); 191 /* Cleanup */ 192 CHKERRQ(VecDestroy(&u)); 193 CHKERRQ(SNESDestroy(&snes)); 194 CHKERRQ(DMDestroy(&dm)); 195 ierr = PetscFinalize(); 196 return ierr; 197 } 198 199 /*TEST 200 201 test: 202 suffix: 2d_p1_0 203 requires: triangle 204 args: -potential_petscspace_degree 2 -pressure_top_petscspace_degree 0 -dm_refine 0 -dmsnes_check 205 206 test: 207 suffix: 2d_p1_1 208 requires: triangle 209 args: -potential_petscspace_degree 1 -pressure_top_petscspace_degree 0 -dm_refine 0 -convest_num_refine 3 -snes_convergence_estimate 210 211 TEST*/ 212