1 static char help[] = "Poisson Problem in 2d and 3d with simplicial finite elements in both primal and mixed form.\n\ 2 We solve the Poisson problem in a rectangular\n\ 3 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\ 4 This example solves mixed form equation to get the flux field to calculate flux norm. We then use that for adaptive mesh refinement. \n\n\n"; 5 6 /* 7 The primal (or original) Poisson problem, in the strong form, is given by, 8 9 \begin{align} 10 - \nabla \cdot ( \nabla u ) = f 11 \end{align} 12 where $u$ is potential. 13 14 The weak form of this equation is 15 16 \begin{align} 17 < \nabla v, \nabla u > - < v, \nabla u \cdot \hat{n} >_\Gamma - < v, f > = 0 18 \end{align} 19 20 The mixed Poisson problem, in the strong form, is given by, 21 22 \begin{align} 23 q - \nabla u &= 0 \\ 24 - \nabla \cdot q &= f 25 \end{align} 26 where $u$ is the potential and $q$ is the flux. 27 28 The weak form of this equation is 29 30 \begin{align} 31 < t, q > + < \nabla \cdot t, u > - < t \cdot \hat{n}, u >_\Gamma &= 0 \\ 32 <v, \nabla \cdot q> - < v, f > &= 0 33 \end{align} 34 35 We solve both primal and mixed problem and calculate the error in the flux norm, namely || e || = || q^m_h - \nabla u^p_h ||. Here superscript 'm' represents field from mixed form and 'p' represents field from the primal form. 36 37 The following boundary conditions are prescribed. 38 39 Primal problem: 40 \begin{align} 41 u = u_0 on \Gamma_D 42 \nabla u \cdot \hat{n} = g on \Gamma_N 43 \end{align} 44 45 Mixed problem: 46 \begin{align} 47 u = u_0 on \Gamma_D 48 q \cdot \hat{n} = g on \Gamma_N 49 \end{align} 50 __________\Gamma_D_____________ 51 | | 52 | | 53 | | 54 \Gamma_N \Gamma_N 55 | | 56 | | 57 | | 58 |_________\Gamma_D_____________| 59 60 To visualize the automated adaptation 61 62 -dm_adapt_pre_view draw -dm_adapt_view draw -draw_pause -1 -geometry 0,0,1024,1024 63 64 and to compare with a naice gradient estimator use 65 66 -adaptor_type gradient 67 68 To see a sequence of adaptations 69 70 -snes_adapt_sequence 8 -adaptor_monitor_error draw::draw_lg 71 -dm_adapt_pre_view draw -dm_adapt_iter_view draw -dm_adapt_view draw -draw_pause 1 -geometry 0,0,1024,1024 72 73 To get a better view of the by-hand process, use 74 75 -dm_view hdf5:${PWD}/mesh.h5 76 -primal_sol_vec_view hdf5:${PWD}/mesh.h5::append 77 -flux_error_vec_view hdf5:${PWD}/mesh.h5::append 78 -exact_error_vec_view hdf5:${PWD}/mesh.h5::append 79 -refine_vec_view hdf5:${PWD}/mesh.h5::append 80 -adapt_dm_view draw -draw_pause -1 81 82 This is also possible with the automated path 83 84 -dm_view hdf5:${PWD}/mesh.h5 85 -adapt_primal_sol_vec_view hdf5:${PWD}/mesh.h5::append 86 -adapt_error_vec_view hdf5:${PWD}/mesh.h5::append 87 -adapt_vec_view hdf5:${PWD}/mesh.h5::append 88 */ 89 90 #include <petsc/private/petscfeimpl.h> 91 #include <petscdmplex.h> 92 #include <petscsnes.h> 93 #include <petscdmadaptor.h> 94 #include <petscds.h> 95 #include <petscviewerhdf5.h> 96 #include <petscbag.h> 97 98 PETSC_EXTERN PetscErrorCode SetupMixed(DMAdaptor, DM); 99 100 typedef enum { 101 SOL_QUADRATIC, 102 SOL_TRIG, 103 SOL_SENSOR, 104 SOL_UNKNOWN, 105 NUM_SOL_TYPE 106 } SolType; 107 const char *SolTypeNames[NUM_SOL_TYPE + 4] = {"quadratic", "trig", "sensor", "unknown", "SolType", "SOL_", NULL}; 108 109 typedef struct { 110 PetscBag param; 111 SolType solType; 112 PetscBool byHand; 113 } AppCtx; 114 115 typedef struct { 116 PetscReal k; 117 } Parameter; 118 119 /* Exact solution: u = x^2 + y^2 */ 120 static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 121 { 122 PetscInt d; 123 124 u[0] = 0.0; 125 for (d = 0; d < dim; ++d) u[0] += x[d] * x[d]; 126 return PETSC_SUCCESS; 127 } 128 /* Exact solution: q = (2x, 2y) */ 129 static PetscErrorCode quadratic_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 130 { 131 PetscInt c; 132 for (c = 0; c < Nc; ++c) u[c] = 2.0 * x[c]; 133 return PETSC_SUCCESS; 134 } 135 136 /* Exact solution: u = sin( n \pi x ) * sin( n \pi y ) */ 137 static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 138 { 139 const PetscReal n = 5.5; 140 141 u[0] = 1.0; 142 for (PetscInt d = 0; d < dim; ++d) u[0] *= PetscSinReal(n * PETSC_PI * x[d]); 143 return PETSC_SUCCESS; 144 } 145 static PetscErrorCode trig_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 146 { 147 const PetscReal n = 5.5; 148 149 for (PetscInt c = 0; c < Nc; c++) u[c] = n * PETSC_PI * PetscCosReal(n * PETSC_PI * x[c]) * PetscSinReal(n * PETSC_PI * x[Nc - c - 1]); 150 return PETSC_SUCCESS; 151 } 152 153 /* 154 Classic hyperbolic sensor function for testing multi-scale anisotropic mesh adaptation: 155 156 f:[-1, 1]x[-1, 1] \to R, 157 f(x, y) = sin(50xy)/100 if |xy| > 2\pi/50 else sin(50xy) 158 159 (mapped to have domain [0,1] x [0,1] in this case). 160 */ 161 static PetscErrorCode sensor_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) 162 { 163 const PetscReal xref = 2. * x[0] - 1.; 164 const PetscReal yref = 2. * x[1] - 1.; 165 const PetscReal xy = xref * yref; 166 167 u[0] = PetscSinReal(50. * xy); 168 if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) u[0] *= 0.01; 169 170 return PETSC_SUCCESS; 171 } 172 173 /* Flux is (cos(50xy) * 50y/100, cos(50xy) * 50x/100) if |xy| > 2\pi/50 else (cos(50xy) * 50y, cos(50xy) * 50x) */ 174 static PetscErrorCode sensor_q(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) 175 { 176 const PetscReal xref = 2. * x[0] - 1.; 177 const PetscReal yref = 2. * x[1] - 1.; 178 const PetscReal xy = xref * yref; 179 180 u[0] = 50. * yref * PetscCosReal(50. * xy) * 2.0; 181 u[1] = 50. * xref * PetscCosReal(50. * xy) * 2.0; 182 if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) { 183 u[0] *= 0.01; 184 u[1] *= 0.01; 185 } 186 return PETSC_SUCCESS; 187 } 188 189 /* We set up residuals and Jacobians for the primal problem. */ 190 static void f0_quadratic_primal(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[]) 191 { 192 f0[0] = 4.0; 193 } 194 195 static void f0_trig_primal(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[]) 196 { 197 const PetscReal n = 5.5; 198 199 f0[0] = -2.0 * PetscSqr(n * PETSC_PI) * PetscSinReal(n * PETSC_PI * x[0]) * PetscSinReal(n * PETSC_PI * x[1]); 200 } 201 202 static void f0_sensor_primal(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[]) 203 { 204 const PetscReal xref = 2. * x[0] - 1.; 205 const PetscReal yref = 2. * x[1] - 1.; 206 const PetscReal xy = xref * yref; 207 208 f0[0] = -2500.0 * PetscSinReal(50. * xy) * (xref * xref + yref * yref) * 4.0; 209 if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) f0[0] *= 0.01; 210 } 211 212 static void f1_primal(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[]) 213 { 214 const PetscReal k = PetscRealPart(constants[0]); 215 216 for (PetscInt d = 0; d < dim; ++d) f1[d] = k * u_x[uOff_x[0] + d]; 217 } 218 219 static void f0_quadratic_bd_primal(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 220 { 221 const PetscReal k = PetscRealPart(constants[0]); 222 PetscScalar flux; 223 224 PetscCallAbort(PETSC_COMM_SELF, quadratic_q(dim, t, x, dim, &flux, NULL)); 225 for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * flux * n[d]; 226 } 227 228 static void f0_trig_bd_primal(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 229 { 230 const PetscReal k = PetscRealPart(constants[0]); 231 PetscScalar flux; 232 233 PetscCallAbort(PETSC_COMM_SELF, trig_q(dim, t, x, dim, &flux, NULL)); 234 for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * flux * n[d]; 235 } 236 237 static void f0_sensor_bd_primal(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 238 { 239 const PetscReal k = PetscRealPart(constants[0]); 240 PetscScalar flux; 241 242 PetscCallAbort(PETSC_COMM_SELF, sensor_q(dim, t, x, dim, &flux, NULL)); 243 for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * flux * n[d]; 244 } 245 246 static void g3_primal_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[]) 247 { 248 const PetscReal k = PetscRealPart(constants[0]); 249 250 for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = k; 251 } 252 253 /* Now we set up the residuals and Jacobians mixed problem. */ 254 static void f0_mixed_quadratic_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[]) 255 { 256 f0[0] = 4.0; 257 for (PetscInt d = 0; d < dim; ++d) f0[0] += -u_x[uOff_x[0] + d * dim + d]; 258 } 259 static void f0_mixed_trig_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[]) 260 { 261 PetscReal n = 5.5; 262 263 f0[0] = -2.0 * PetscSqr(n * PETSC_PI) * PetscSinReal(n * PETSC_PI * x[0]) * PetscSinReal(n * PETSC_PI * x[1]); 264 for (PetscInt d = 0; d < dim; ++d) f0[0] += -u_x[uOff_x[0] + d * dim + d]; 265 } 266 static void f0_mixed_sensor_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[]) 267 { 268 const PetscReal xref = 2. * x[0] - 1.; 269 const PetscReal yref = 2. * x[1] - 1.; 270 const PetscReal xy = xref * yref; 271 272 f0[0] = -2500.0 * PetscSinReal(50. * xy) * (xref * xref + yref * yref) * 4.0; 273 if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) f0[0] *= 0.01; 274 for (PetscInt d = 0; d < dim; ++d) f0[0] += -u_x[uOff_x[0] + d * dim + d]; 275 } 276 277 static void f0_mixed_q(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[]) 278 { 279 for (PetscInt d = 0; d < dim; d++) f0[d] = u[uOff[0] + d]; 280 } 281 282 static void f1_mixed_q(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[]) 283 { 284 const PetscReal k = PetscRealPart(constants[0]); 285 286 for (PetscInt d = 0; d < dim; d++) f1[d * dim + d] = k * u[uOff[1]]; 287 } 288 289 static void f0_quadratic_bd_mixed_q(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 290 { 291 const PetscReal k = PetscRealPart(constants[0]); 292 PetscScalar potential; 293 294 PetscCallAbort(PETSC_COMM_SELF, quadratic_u(dim, t, x, dim, &potential, NULL)); 295 for (PetscInt d = 0; d < dim; ++d) f0[d] = -k * potential * n[d]; 296 } 297 298 static void f0_trig_bd_mixed_q(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 299 { 300 const PetscReal k = PetscRealPart(constants[0]); 301 PetscScalar potential; 302 303 PetscCallAbort(PETSC_COMM_SELF, trig_u(dim, t, x, dim, &potential, NULL)); 304 for (PetscInt d = 0; d < dim; ++d) f0[d * dim + d] = -k * potential * n[d]; 305 } 306 307 static void f0_sensor_bd_mixed_q(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 308 { 309 const PetscReal k = PetscRealPart(constants[0]); 310 PetscScalar potential; 311 312 PetscCallAbort(PETSC_COMM_SELF, sensor_u(dim, t, x, dim, &potential, NULL)); 313 for (PetscInt d = 0; d < dim; ++d) f0[d * dim + d] = -k * potential * n[d]; 314 } 315 316 /* <v, \nabla\cdot q> */ 317 static void g1_mixed_uq(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 g1[]) 318 { 319 for (PetscInt d = 0; d < dim; d++) g1[d * dim + d] = -1.0; 320 } 321 322 /* < t, q> */ 323 static void g0_mixed_qq(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[]) 324 { 325 for (PetscInt d = 0; d < dim; d++) g0[d * dim + d] = 1.0; 326 } 327 328 /* <\nabla\cdot t, u> = <\nabla t, Iu> */ 329 static void g2_mixed_qu(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 g2[]) 330 { 331 const PetscReal k = PetscRealPart(constants[0]); 332 333 for (PetscInt d = 0; d < dim; d++) g2[d * dim + d] = k; 334 } 335 336 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *user) 337 { 338 PetscFunctionBeginUser; 339 PetscOptionsBegin(comm, "", "Flux norm error in primal poisson problem Options", "DMPLEX"); 340 user->byHand = PETSC_TRUE; 341 user->solType = SOL_QUADRATIC; 342 343 PetscCall(PetscOptionsGetBool(NULL, NULL, "-by_hand", &user->byHand, NULL)); 344 PetscCall(PetscOptionsEnum("-sol_type", "Type of exact solution", "ex27.c", SolTypeNames, (PetscEnum)user->solType, (PetscEnum *)&user->solType, NULL)); 345 PetscOptionsEnd(); 346 PetscFunctionReturn(PETSC_SUCCESS); 347 } 348 349 static PetscErrorCode SetupParameters(PetscBag bag, AppCtx *user) 350 { 351 Parameter *param; 352 353 PetscFunctionBeginUser; 354 PetscCall(PetscBagGetData(bag, (void **)¶m)); 355 PetscCall(PetscBagSetName(bag, "par", "Poisson parameters")); 356 PetscCall(PetscBagRegisterReal(bag, ¶m->k, 1.0, "k", "Thermal conductivity")); 357 PetscCall(PetscBagSetFromOptions(bag)); 358 PetscFunctionReturn(PETSC_SUCCESS); 359 } 360 361 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 362 { 363 PetscFunctionBeginUser; 364 PetscCall(DMCreate(comm, dm)); 365 PetscCall(DMSetType(*dm, DMPLEX)); 366 PetscCall(DMSetFromOptions(*dm)); 367 PetscCall(DMSetApplicationContext(*dm, &user)); 368 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 369 PetscFunctionReturn(PETSC_SUCCESS); 370 } 371 372 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 373 { 374 PetscDS ds; 375 DMLabel label; 376 PetscInt id, bd; 377 Parameter *param; 378 PetscWeakForm wf; 379 380 PetscFunctionBeginUser; 381 PetscCall(DMGetDS(dm, &ds)); 382 PetscCall(DMGetLabel(dm, "marker", &label)); 383 384 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_primal_uu)); 385 386 switch (user->solType) { 387 case SOL_QUADRATIC: 388 PetscCall(PetscDSSetResidual(ds, 0, f0_quadratic_primal, f1_primal)); 389 390 id = 1; 391 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (void (*)(void))quadratic_u, NULL, user, NULL)); 392 393 id = 2; 394 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 395 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 396 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_primal, 0, NULL)); 397 398 id = 3; 399 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (void (*)(void))quadratic_u, NULL, user, NULL)); 400 401 id = 4; 402 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 403 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 404 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_primal, 0, NULL)); 405 406 PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_u, user)); 407 break; 408 case SOL_TRIG: 409 PetscCall(PetscDSSetResidual(ds, 0, f0_trig_primal, f1_primal)); 410 411 id = 1; 412 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL)); 413 414 id = 2; 415 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 416 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 417 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_primal, 0, NULL)); 418 419 id = 3; 420 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL)); 421 422 id = 4; 423 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 424 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 425 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_primal, 0, NULL)); 426 427 PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user)); 428 break; 429 case SOL_SENSOR: 430 PetscCall(PetscDSSetResidual(ds, 0, f0_sensor_primal, f1_primal)); 431 432 id = 1; 433 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall primal potential", label, 1, &id, 0, 0, NULL, (void (*)(void))sensor_u, NULL, user, NULL)); 434 435 id = 2; 436 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 437 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 438 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_primal, 0, NULL)); 439 440 id = 3; 441 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall primal potential", label, 1, &id, 0, 0, NULL, (void (*)(void))sensor_u, NULL, user, NULL)); 442 443 id = 4; 444 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall flux", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 445 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 446 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_primal, 0, NULL)); 447 448 PetscCall(PetscDSSetExactSolution(ds, 0, sensor_u, user)); 449 break; 450 default: 451 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid exact solution type %s", SolTypeNames[PetscMin(user->solType, SOL_UNKNOWN)]); 452 } 453 454 /* Setup constants */ 455 { 456 PetscCall(PetscBagGetData(user->param, (void **)¶m)); 457 PetscScalar constants[1]; 458 459 constants[0] = param->k; 460 461 PetscCall(PetscDSSetConstants(ds, 1, constants)); 462 } 463 PetscFunctionReturn(PETSC_SUCCESS); 464 } 465 466 static PetscErrorCode SetupPrimalDiscretization(DM dm, AppCtx *user) 467 { 468 DM cdm = dm; 469 PetscFE fe[1]; 470 DMPolytopeType ct; 471 PetscInt dim, cStart; 472 MPI_Comm comm; 473 474 PetscFunctionBeginUser; 475 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 476 PetscCall(DMGetDimension(dm, &dim)); 477 478 /* Create finite element */ 479 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 480 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 481 PetscCall(PetscFECreateByCell(comm, dim, 1, ct, "primal_potential_", PETSC_DEFAULT, &fe[0])); 482 PetscCall(PetscObjectSetName((PetscObject)fe[0], "primal potential")); 483 484 /* Set discretization and boundary conditions for each mesh */ 485 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0])); 486 PetscCall(DMCreateDS(dm)); 487 while (cdm) { 488 PetscCall(DMCopyDisc(dm, cdm)); 489 PetscCall(DMGetCoarseDM(cdm, &cdm)); 490 } 491 492 PetscCall(PetscFEDestroy(&fe[0])); 493 PetscFunctionReturn(PETSC_SUCCESS); 494 } 495 496 static PetscErrorCode SetupMixedProblem(DM dm, AppCtx *user) 497 { 498 PetscDS ds; 499 DMLabel label; 500 PetscInt id, bd; 501 Parameter *param; 502 PetscWeakForm wf; 503 504 PetscFunctionBeginUser; 505 PetscCall(DMGetDS(dm, &ds)); 506 PetscCall(DMGetLabel(dm, "marker", &label)); 507 508 /* Residual terms */ 509 PetscCall(PetscDSSetResidual(ds, 0, f0_mixed_q, f1_mixed_q)); 510 511 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_mixed_qq, NULL, NULL, NULL)); 512 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_mixed_qu, NULL)); 513 514 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_mixed_uq, NULL, NULL)); 515 516 switch (user->solType) { 517 case SOL_QUADRATIC: 518 PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_quadratic_u, NULL)); 519 520 id = 1; 521 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral bottom wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 522 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 523 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_mixed_q, 0, NULL)); 524 525 id = 2; 526 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 0, 0, NULL, (void (*)(void))quadratic_q, NULL, user, NULL)); 527 528 id = 4; 529 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 0, 0, NULL, (void (*)(void))quadratic_q, NULL, user, NULL)); 530 531 id = 3; 532 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral top wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 533 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 534 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_quadratic_bd_mixed_q, 0, NULL)); 535 536 PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_q, user)); 537 PetscCall(PetscDSSetExactSolution(ds, 1, quadratic_u, user)); 538 break; 539 case SOL_TRIG: 540 PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_trig_u, NULL)); 541 542 id = 1; 543 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral bottom wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 544 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 545 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_mixed_q, 0, NULL)); 546 547 id = 2; 548 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 1, 0, NULL, (void (*)(void))trig_q, NULL, user, NULL)); 549 550 id = 4; 551 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 1, 0, NULL, (void (*)(void))trig_q, NULL, user, NULL)); 552 553 id = 3; 554 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral top wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 555 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 556 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_trig_bd_mixed_q, 0, NULL)); 557 558 PetscCall(PetscDSSetExactSolution(ds, 0, trig_q, user)); 559 PetscCall(PetscDSSetExactSolution(ds, 1, trig_u, user)); 560 break; 561 case SOL_SENSOR: 562 PetscCall(PetscDSSetResidual(ds, 1, f0_mixed_sensor_u, NULL)); 563 564 id = 1; 565 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral bottom wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 566 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 567 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_mixed_q, 0, NULL)); 568 569 id = 2; 570 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right wall flux", label, 1, &id, 1, 0, NULL, (void (*)(void))sensor_q, NULL, user, NULL)); 571 572 id = 4; 573 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left wall flux", label, 1, &id, 1, 0, NULL, (void (*)(void))sensor_q, NULL, user, NULL)); 574 575 id = 3; 576 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "Dirichlet Bd Integral top wall", label, 1, &id, 0, 0, NULL, NULL, NULL, user, &bd)); 577 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 578 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_sensor_bd_mixed_q, 0, NULL)); 579 580 PetscCall(PetscDSSetExactSolution(ds, 0, sensor_q, user)); 581 PetscCall(PetscDSSetExactSolution(ds, 1, sensor_u, user)); 582 break; 583 default: 584 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid exact solution type %s", SolTypeNames[PetscMin(user->solType, SOL_UNKNOWN)]); 585 } 586 587 /* Setup constants */ 588 { 589 PetscCall(PetscBagGetData(user->param, (void **)¶m)); 590 PetscScalar constants[1]; 591 592 constants[0] = param->k; 593 594 PetscCall(PetscDSSetConstants(ds, 1, constants)); 595 } 596 PetscFunctionReturn(PETSC_SUCCESS); 597 } 598 599 static PetscErrorCode SetupMixedDiscretization(DM dm, AppCtx *user) 600 { 601 DM cdm = dm; 602 PetscFE fe[2]; 603 DMPolytopeType ct; 604 PetscInt dim, cStart; 605 MPI_Comm comm; 606 607 PetscFunctionBeginUser; 608 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 609 PetscCall(DMGetDimension(dm, &dim)); 610 611 /* Create finite element */ 612 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 613 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 614 PetscCall(PetscFECreateByCell(comm, dim, dim, ct, "mixed_flux_", PETSC_DEFAULT, &fe[0])); 615 PetscCall(PetscObjectSetName((PetscObject)fe[0], "mixed flux")); 616 /* NOTE: Set the same quadrature order as the primal problem here or use the command line option. */ 617 618 PetscCall(PetscFECreateByCell(comm, dim, 1, ct, "mixed_potential_", PETSC_DEFAULT, &fe[1])); 619 PetscCall(PetscFECopyQuadrature(fe[0], fe[1])); 620 PetscCall(PetscObjectSetName((PetscObject)fe[1], "mixed potential")); 621 622 /* Set discretization and boundary conditions for each mesh */ 623 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0])); 624 PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe[1])); 625 PetscCall(DMCreateDS(dm)); 626 while (cdm) { 627 PetscCall(DMCopyDisc(dm, cdm)); 628 PetscCall(DMGetCoarseDM(cdm, &cdm)); 629 } 630 PetscCall(PetscFEDestroy(&fe[0])); 631 PetscCall(PetscFEDestroy(&fe[1])); 632 PetscFunctionReturn(PETSC_SUCCESS); 633 } 634 635 PetscErrorCode SetupMixed(DMAdaptor adaptor, DM mdm) 636 { 637 AppCtx *ctx; 638 639 PetscFunctionBeginUser; 640 PetscCall(DMGetApplicationContext(mdm, (void **)&ctx)); 641 PetscCall(SetupMixedDiscretization(mdm, ctx)); 642 PetscCall(SetupMixedProblem(mdm, ctx)); 643 PetscFunctionReturn(PETSC_SUCCESS); 644 } 645 646 int main(int argc, char **argv) 647 { 648 DM dm, mdm; /* problem specification */ 649 SNES snes, msnes; /* nonlinear solvers */ 650 Vec u, mu; /* solution vectors */ 651 Vec fluxError, exactError; /* Element wise error vector */ 652 PetscReal fluxNorm, exactNorm; /* Flux error norm */ 653 AppCtx user; /* user-defined work context */ 654 655 PetscFunctionBeginUser; 656 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 657 PetscCall(PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.param)); 658 PetscCall(SetupParameters(user.param, &user)); 659 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 660 661 // Set up and solve primal problem 662 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 663 PetscCall(DMSetApplicationContext(dm, &user)); 664 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 665 PetscCall(SNESSetDM(snes, dm)); 666 667 PetscCall(SetupPrimalDiscretization(dm, &user)); 668 PetscCall(SetupPrimalProblem(dm, &user)); 669 PetscCall(DMCreateGlobalVector(dm, &u)); 670 PetscCall(VecSet(u, 0.0)); 671 PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user)); 672 PetscCall(SNESSetFromOptions(snes)); 673 674 PetscCall(DMSNESCheckFromOptions(snes, u)); 675 PetscCall(SNESSolve(snes, NULL, u)); 676 677 // Needed if you allow SNES to refine 678 PetscCall(SNESGetSolution(snes, &u)); 679 PetscCall(VecGetDM(u, &dm)); 680 681 PetscCall(PetscObjectSetName((PetscObject)u, "Primal Solution ")); 682 PetscCall(VecViewFromOptions(u, NULL, "-primal_sol_vec_view")); 683 684 if (user.byHand) { 685 // Set up and solve mixed problem 686 PetscCall(DMClone(dm, &mdm)); 687 PetscCall(SNESCreate(PETSC_COMM_WORLD, &msnes)); 688 PetscCall(SNESSetOptionsPrefix(msnes, "mixed_")); 689 PetscCall(SNESSetDM(msnes, mdm)); 690 691 PetscCall(SetupMixedDiscretization(mdm, &user)); 692 PetscCall(SetupMixedProblem(mdm, &user)); 693 PetscCall(DMCreateGlobalVector(mdm, &mu)); 694 PetscCall(VecSet(mu, 0.0)); 695 PetscCall(DMPlexSetSNESLocalFEM(mdm, PETSC_FALSE, &user)); 696 PetscCall(SNESSetFromOptions(msnes)); 697 698 PetscCall(DMSNESCheckFromOptions(msnes, mu)); 699 PetscCall(SNESSolve(msnes, NULL, mu)); 700 PetscCall(PetscObjectSetName((PetscObject)mu, "Mixed Solution ")); 701 PetscCall(VecViewFromOptions(mu, NULL, "-mixed_sol_vec_view")); 702 703 // Create the error space of piecewise constants 704 DM dmErr; 705 PetscFE feErr; 706 DMPolytopeType ct; 707 PetscInt dim, cStart; 708 709 PetscCall(DMClone(dm, &dmErr)); 710 PetscCall(DMGetDimension(dmErr, &dim)); 711 PetscCall(DMPlexGetHeightStratum(dmErr, 0, &cStart, NULL)); 712 PetscCall(DMPlexGetCellType(dmErr, cStart, &ct)); 713 PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, 1, ct, 0, PETSC_DEFAULT, &feErr)); 714 PetscCall(PetscObjectSetName((PetscObject)feErr, "Error")); 715 PetscCall(DMSetField(dmErr, 0, NULL, (PetscObject)feErr)); 716 PetscCall(PetscFEDestroy(&feErr)); 717 PetscCall(DMCreateDS(dmErr)); 718 PetscCall(DMViewFromOptions(dmErr, NULL, "-dmerr_view")); 719 720 // Compute the flux norm 721 PetscCall(DMGetGlobalVector(dmErr, &fluxError)); 722 PetscCall(PetscObjectSetName((PetscObject)fluxError, "Flux Error")); 723 PetscCall(DMGetGlobalVector(dmErr, &exactError)); 724 PetscCall(PetscObjectSetName((PetscObject)exactError, "Analytical Error")); 725 PetscCall(DMPlexComputeL2FluxDiffVec(u, 0, mu, 0, fluxError)); 726 { 727 PetscDS ds; 728 PetscSimplePointFn *func[2] = {NULL, NULL}; 729 void *ctx[2] = {NULL, NULL}; 730 731 PetscCall(DMGetDS(mdm, &ds)); 732 PetscCall(PetscDSGetExactSolution(ds, 0, &func[0], &ctx[0])); 733 PetscCall(DMPlexComputeL2DiffVec(mdm, 0.0, func, ctx, mu, exactError)); 734 } 735 PetscCall(VecNorm(fluxError, NORM_2, &fluxNorm)); 736 PetscCall(VecNorm(exactError, NORM_2, &exactNorm)); 737 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Flux error norm = %g\t Exact flux error norm = %g\n", (double)fluxNorm, (double)exactNorm)); 738 PetscCall(VecViewFromOptions(fluxError, NULL, "-flux_error_vec_view")); 739 PetscCall(VecViewFromOptions(exactError, NULL, "-exact_error_vec_view")); 740 741 // Adaptive refinement based on calculated error 742 DM rdm; 743 VecTagger refineTag; 744 DMLabel adaptLabel; 745 IS refineIS; 746 Vec ref; 747 748 PetscCall(DMLabelCreate(PETSC_COMM_WORLD, "adapt", &adaptLabel)); 749 PetscCall(DMLabelSetDefaultValue(adaptLabel, DM_ADAPT_COARSEN)); 750 PetscCall(VecTaggerCreate(PETSC_COMM_WORLD, &refineTag)); 751 PetscCall(VecTaggerSetFromOptions(refineTag)); 752 PetscCall(VecTaggerSetUp(refineTag)); 753 PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view")); 754 755 PetscCall(VecTaggerComputeIS(refineTag, fluxError, &refineIS, NULL)); 756 PetscCall(VecTaggerDestroy(&refineTag)); 757 PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS)); 758 PetscCall(ISViewFromOptions(refineIS, NULL, "-refine_is_view")); 759 PetscCall(ISDestroy(&refineIS)); 760 761 PetscCall(DMPlexCreateLabelField(dm, adaptLabel, &ref)); 762 PetscCall(VecViewFromOptions(ref, NULL, "-refine_vec_view")); 763 PetscCall(VecDestroy(&ref)); 764 765 // Mark adaptation phase with prefix ref_ 766 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, "adapt_")); 767 PetscCall(DMAdaptLabel(dm, adaptLabel, &rdm)); 768 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, NULL)); 769 PetscCall(PetscObjectSetName((PetscObject)rdm, "Adaptively Refined DM")); 770 PetscCall(DMViewFromOptions(rdm, NULL, "-adapt_dm_view")); 771 PetscCall(DMDestroy(&rdm)); 772 PetscCall(DMLabelDestroy(&adaptLabel)); 773 774 // Destroy the error structures 775 PetscCall(DMRestoreGlobalVector(dmErr, &fluxError)); 776 PetscCall(DMRestoreGlobalVector(dmErr, &exactError)); 777 PetscCall(DMDestroy(&dmErr)); 778 779 // Destroy the mixed structures 780 PetscCall(VecDestroy(&mu)); 781 PetscCall(DMDestroy(&mdm)); 782 PetscCall(SNESDestroy(&msnes)); 783 } 784 785 // Destroy the primal structures 786 PetscCall(VecDestroy(&u)); 787 PetscCall(DMDestroy(&dm)); 788 PetscCall(SNESDestroy(&snes)); 789 PetscCall(PetscBagDestroy(&user.param)); 790 PetscCall(PetscFinalize()); 791 return 0; 792 } 793 794 /*TEST 795 796 # Tests using the explicit code above 797 testset: 798 suffix: 2d_p2_rt0p0_byhand 799 requires: triangle 800 args: -dm_plex_box_faces 1,1 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -dm_refine 3 \ 801 -primal_potential_petscspace_degree 2 \ 802 -mixed_potential_petscdualspace_lagrange_continuity 0 \ 803 -mixed_flux_petscspace_type ptrimmed \ 804 -mixed_flux_petscspace_components 2 \ 805 -mixed_flux_petscspace_ptrimmed_form_degree -1 \ 806 -mixed_flux_petscdualspace_order 1 \ 807 -mixed_flux_petscdualspace_form_degree -1 \ 808 -mixed_flux_petscdualspace_lagrange_trimmed true \ 809 -mixed_flux_petscfe_default_quadrature_order 2 \ 810 -vec_tagger_type cdf -vec_tagger_box 0.9,1.0 \ 811 -tag_view \ 812 -adapt_dm_adaptor cellrefiner -adapt_dm_plex_transform_type refine_sbr \ 813 -dmsnes_check 0.001 -mixed_dmsnes_check 0.001 -pc_type jacobi -mixed_pc_type jacobi 814 test: 815 suffix: quadratic 816 args: -sol_type quadratic 817 test: 818 suffix: trig 819 args: -sol_type trig 820 test: 821 suffix: sensor 822 args: -sol_type sensor 823 824 # Tests using the embedded adaptor in SNES 825 testset: 826 suffix: 2d_p2_rt0p0 827 requires: triangle defined(PETSC_HAVE_EXECUTABLE_EXPORT) 828 args: -dm_plex_box_faces 1,1 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -dm_refine 3 \ 829 -primal_potential_petscspace_degree 2 \ 830 -mixed_potential_petscdualspace_lagrange_continuity 0 \ 831 -mixed_flux_petscspace_type ptrimmed \ 832 -mixed_flux_petscspace_components 2 \ 833 -mixed_flux_petscspace_ptrimmed_form_degree -1 \ 834 -mixed_flux_petscdualspace_order 1 \ 835 -mixed_flux_petscdualspace_form_degree -1 \ 836 -mixed_flux_petscdualspace_lagrange_trimmed true \ 837 -mixed_flux_petscfe_default_quadrature_order 2 \ 838 -by_hand 0 \ 839 -refine_vec_tagger_type cdf -refine_vec_tagger_box 0.9,1.0 \ 840 -snes_adapt_view \ 841 -adapt_dm_adaptor cellrefiner -adapt_dm_plex_transform_type refine_sbr \ 842 -adaptor_criterion label -adaptor_type flux -adaptor_mixed_setup_function SetupMixed \ 843 -snes_adapt_sequence 1 -pc_type jacobi -mixed_pc_type jacobi 844 test: 845 suffix: quadratic 846 args: -sol_type quadratic -adaptor_monitor_error 847 test: 848 suffix: trig 849 args: -sol_type trig -adaptor_monitor_error 850 test: 851 suffix: sensor 852 args: -sol_type sensor 853 854 TEST*/ 855