1 static const char help[] = "\ 2 Minimum surface area problem in 2D using DMPLEX.\n\ 3 It solves the unconstrained minimization problem \n\ 4 \n\ 5 argmin \\int_\\Omega (1 + ||u||^2), u = g on \\partial\\Omega.\n\ 6 \n\ 7 This example shows how to setup an optimization problem using DMPLEX FEM routines.\n\ 8 It supports nonlinear domain decomposition and multilevel solvers.\n\ 9 \n\n"; 10 11 #include <petscdmplex.h> 12 #include <petscsnes.h> 13 #include <petscds.h> 14 15 /* The pointwise evaluation routine for the objective function */ 16 static void objective_2d(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 obj[]) 17 { 18 const PetscScalar ux2 = PetscSqr(u_x[0]); 19 const PetscScalar uy2 = PetscSqr(u_x[1]); 20 21 obj[0] = PetscSqrtReal(PetscAbsScalar(1.0 + ux2 + uy2)); 22 } 23 24 /* The pointwise evaluation routine for the gradient wrt the gradients of the FEM basis */ 25 static void gradient_1_2d(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 f[]) 26 { 27 const PetscScalar ux2 = PetscSqr(u_x[0]); 28 const PetscScalar uy2 = PetscSqr(u_x[1]); 29 const PetscScalar v = 1. / PetscSqrtReal(PetscAbsScalar(1.0 + ux2 + uy2)); 30 31 f[0] = v * u_x[0]; 32 f[1] = v * u_x[1]; 33 } 34 35 /* The pointwise evaluation routine for the hessian wrt the gradients of the FEM basis */ 36 static void hessian_11_2d(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 jac[]) 37 { 38 const PetscScalar ux2 = PetscSqr(u_x[0]); 39 const PetscScalar uy2 = PetscSqr(u_x[1]); 40 const PetscScalar uxy = u_x[0] * u_x[1]; 41 const PetscScalar v1 = 1. / PetscSqrtReal(PetscAbsScalar(1.0 + ux2 + uy2)); 42 const PetscScalar v2 = v1 / (1.0 + ux2 + uy2); 43 44 jac[0] = v1 - v2 * ux2; 45 jac[1] = -v2 * uxy; 46 jac[2] = -v2 * uxy; 47 jac[3] = v1 - v2 * uy2; 48 } 49 50 /* The boundary conditions we impose */ 51 static PetscErrorCode sins_2d(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx) 52 { 53 PetscFunctionBeginUser; 54 const PetscReal pi2 = PETSC_PI * 2; 55 const PetscReal x = xx[0]; 56 const PetscReal y = xx[1]; 57 58 *u = (y - 0.5) * PetscSinReal(pi2 * x) + (x - 0.5) * PetscSinReal(pi2 * y); 59 PetscFunctionReturn(PETSC_SUCCESS); 60 } 61 62 static PetscErrorCode CreateBCLabel(DM dm, const char name[]) 63 { 64 DM plex; 65 DMLabel label; 66 67 PetscFunctionBeginUser; 68 PetscCall(DMCreateLabel(dm, name)); 69 PetscCall(DMGetLabel(dm, name, &label)); 70 PetscCall(DMConvert(dm, DMPLEX, &plex)); 71 PetscCall(DMPlexMarkBoundaryFaces(plex, 1, label)); 72 PetscCall(DMDestroy(&plex)); 73 PetscFunctionReturn(PETSC_SUCCESS); 74 } 75 76 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm) 77 { 78 PetscFunctionBeginUser; 79 PetscCall(DMCreate(comm, dm)); 80 PetscCall(DMSetType(*dm, DMPLEX)); 81 PetscCall(DMSetFromOptions(*dm)); 82 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 83 PetscFunctionReturn(PETSC_SUCCESS); 84 } 85 86 static PetscErrorCode SetupProblem(DM dm) 87 { 88 PetscDS ds; 89 DMLabel label; 90 PetscInt dim; 91 const PetscInt id = 1; 92 93 PetscFunctionBeginUser; 94 PetscCall(DMGetDS(dm, &ds)); 95 PetscCall(DMGetDimension(dm, &dim)); 96 PetscCheck(dim == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only for 2-D"); 97 PetscCall(PetscDSSetObjective(ds, 0, objective_2d)); 98 PetscCall(PetscDSSetResidual(ds, 0, NULL, gradient_1_2d)); 99 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, hessian_11_2d)); 100 PetscCall(DMGetLabel(dm, "marker", &label)); 101 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "data", label, 1, &id, 0, 0, NULL, (void (*)(void))sins_2d, NULL, NULL, NULL)); 102 PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_TRUE, NULL)); 103 PetscFunctionReturn(PETSC_SUCCESS); 104 } 105 106 static PetscErrorCode SetupDiscretization(DM dm) 107 { 108 DM plex, cdm = dm; 109 PetscFE fe; 110 PetscBool simplex; 111 PetscInt dim; 112 113 PetscFunctionBeginUser; 114 PetscCall(DMGetDimension(dm, &dim)); 115 PetscCall(DMConvert(dm, DMPLEX, &plex)); 116 PetscCall(DMPlexIsSimplex(plex, &simplex)); 117 PetscCall(DMDestroy(&plex)); 118 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, -1, &fe)); 119 PetscCall(PetscObjectSetName((PetscObject)fe, "potential")); 120 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 121 PetscCall(DMCreateDS(dm)); 122 PetscCall(SetupProblem(dm)); 123 while (cdm) { 124 PetscBool hasLabel; 125 126 PetscCall(DMHasLabel(cdm, "marker", &hasLabel)); 127 if (!hasLabel) PetscCall(CreateBCLabel(cdm, "marker")); 128 PetscCall(DMCopyDisc(dm, cdm)); 129 PetscCall(DMGetCoarseDM(cdm, &cdm)); 130 } 131 PetscCall(PetscFEDestroy(&fe)); 132 PetscFunctionReturn(PETSC_SUCCESS); 133 } 134 135 int main(int argc, char **argv) 136 { 137 DM dm; /* Problem specification */ 138 SNES snes; /* nonlinear solver */ 139 Vec u; /* solution vector */ 140 141 PetscFunctionBeginUser; 142 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 143 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 144 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm)); 145 PetscCall(SNESSetDM(snes, dm)); 146 147 PetscCall(SetupDiscretization(dm)); 148 149 PetscCall(SNESSetFromOptions(snes)); 150 151 PetscCall(DMCreateGlobalVector(dm, &u)); 152 PetscCall(SNESSolve(snes, NULL, u)); 153 154 PetscCall(VecDestroy(&u)); 155 PetscCall(SNESDestroy(&snes)); 156 PetscCall(DMDestroy(&dm)); 157 PetscCall(PetscFinalize()); 158 return 0; 159 } 160 161 /*TEST 162 163 test: 164 requires: !single 165 suffix: qn_nasm 166 filter: sed -e "s/CONVERGED_FNORM_ABS/CONVERGED_FNORM_RELATIVE/g" 167 nsize: 4 168 args: -petscspace_degree 1 -dm_refine 2 -snes_type qn -snes_npc_side {{left right}separate output} -npc_snes_type nasm -npc_snes_nasm_type restrict -npc_sub_snes_linesearch_order 1 -npc_sub_snes_linesearch_type bt -dm_plex_dd_overlap 1 -snes_linesearch_type bt -snes_linesearch_order 1 -npc_sub_pc_type lu -npc_sub_ksp_type preonly -snes_converged_reason -snes_monitor_short -petscpartitioner_type simple -npc_sub_snes_max_it 1 -dm_plex_simplex 0 -snes_rtol 1.e-6 169 170 TEST*/ 171