1 static char help[] = "Linear elasticity in 2d and 3d with finite elements.\n\ 2 We solve the elasticity problem in a rectangular\n\ 3 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\ 4 This example supports automatic convergence estimation\n\ 5 and eventually adaptivity.\n\n\n"; 6 7 /* 8 https://en.wikipedia.org/wiki/Linear_elasticity 9 10 Converting elastic constants: 11 lambda = E nu / ((1 + nu) (1 - 2 nu)) 12 mu = E / (2 (1 + nu)) 13 */ 14 15 #include <petscdmplex.h> 16 #include <petscsnes.h> 17 #include <petscds.h> 18 #include <petscbag.h> 19 #include <petscconvest.h> 20 21 typedef enum { 22 SOL_VLAP_QUADRATIC, 23 SOL_ELAS_QUADRATIC, 24 SOL_VLAP_TRIG, 25 SOL_ELAS_TRIG, 26 SOL_ELAS_AXIAL_DISP, 27 SOL_ELAS_UNIFORM_STRAIN, 28 SOL_ELAS_GE, 29 SOL_MASS_QUADRATIC, 30 NUM_SOLUTION_TYPES 31 } SolutionType; 32 const char *solutionTypes[NUM_SOLUTION_TYPES + 1] = {"vlap_quad", "elas_quad", "vlap_trig", "elas_trig", "elas_axial_disp", "elas_uniform_strain", "elas_ge", "mass_quad", "unknown"}; 33 34 typedef enum { 35 DEFORM_NONE, 36 DEFORM_SHEAR, 37 DEFORM_STEP, 38 NUM_DEFORM_TYPES 39 } DeformType; 40 const char *deformTypes[NUM_DEFORM_TYPES + 1] = {"none", "shear", "step", "unknown"}; 41 42 typedef struct { 43 PetscScalar mu; /* shear modulus */ 44 PetscScalar lambda; /* Lame's first parameter */ 45 PetscScalar N; /* Tension force on right wall */ 46 } Parameter; 47 48 typedef struct { 49 /* Domain and mesh definition */ 50 char dmType[256]; /* DM type for the solve */ 51 DeformType deform; /* Domain deformation type */ 52 /* Problem definition */ 53 SolutionType solType; /* Type of exact solution */ 54 PetscBag bag; /* Problem parameters */ 55 /* Solver definition */ 56 PetscBool useNearNullspace; /* Use the rigid body modes as a near nullspace for AMG */ 57 } AppCtx; 58 59 static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 60 { 61 PetscInt d; 62 for (d = 0; d < dim; ++d) u[d] = 0.0; 63 return PETSC_SUCCESS; 64 } 65 66 static PetscErrorCode ge_shift(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 67 { 68 PetscInt d; 69 u[0] = 0.1; 70 for (d = 1; d < dim; ++d) u[d] = 0.0; 71 return PETSC_SUCCESS; 72 } 73 74 static PetscErrorCode quadratic_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 75 { 76 u[0] = x[0] * x[0]; 77 u[1] = x[1] * x[1] - 2.0 * x[0] * x[1]; 78 return PETSC_SUCCESS; 79 } 80 81 static PetscErrorCode quadratic_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 82 { 83 u[0] = x[0] * x[0]; 84 u[1] = x[1] * x[1] - 2.0 * x[0] * x[1]; 85 u[2] = x[2] * x[2] - 2.0 * x[1] * x[2]; 86 return PETSC_SUCCESS; 87 } 88 89 /* 90 u = x^2 91 v = y^2 - 2xy 92 Delta <u,v> - f = <2, 2> - <2, 2> 93 94 u = x^2 95 v = y^2 - 2xy 96 w = z^2 - 2yz 97 Delta <u,v,w> - f = <2, 2, 2> - <2, 2, 2> 98 */ 99 static void f0_vlap_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[]) 100 { 101 PetscInt d; 102 for (d = 0; d < dim; ++d) f0[d] += 2.0; 103 } 104 105 /* 106 u = x^2 107 v = y^2 - 2xy 108 \varepsilon = / 2x -y \ 109 \ -y 2y - 2x / 110 Tr(\varepsilon) = div u = 2y 111 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij} 112 = \lambda \partial_j (2y) + 2\mu < 2-1, 2 > 113 = \lambda < 0, 2 > + \mu < 2, 4 > 114 115 u = x^2 116 v = y^2 - 2xy 117 w = z^2 - 2yz 118 \varepsilon = / 2x -y 0 \ 119 | -y 2y - 2x -z | 120 \ 0 -z 2z - 2y/ 121 Tr(\varepsilon) = div u = 2z 122 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij} 123 = \lambda \partial_j (2z) + 2\mu < 2-1, 2-1, 2 > 124 = \lambda < 0, 0, 2 > + \mu < 2, 2, 4 > 125 */ 126 static void f0_elas_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[]) 127 { 128 const PetscReal mu = PetscRealPart(constants[0]); 129 const PetscReal lambda = PetscRealPart(constants[1]); 130 131 for (PetscInt d = 0; d < dim - 1; ++d) f0[d] += 2.0 * mu; 132 f0[dim - 1] += 2.0 * lambda + 4.0 * mu; 133 } 134 135 static void f0_mass_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[]) 136 { 137 if (dim == 2) { 138 f0[0] -= x[0] * x[0]; 139 f0[1] -= x[1] * x[1] - 2.0 * x[0] * x[1]; 140 } else { 141 f0[0] -= x[0] * x[0]; 142 f0[1] -= x[1] * x[1] - 2.0 * x[0] * x[1]; 143 f0[2] -= x[2] * x[2] - 2.0 * x[1] * x[2]; 144 } 145 } 146 147 static PetscErrorCode trig_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 148 { 149 u[0] = PetscSinReal(2.0 * PETSC_PI * x[0]); 150 u[1] = PetscSinReal(2.0 * PETSC_PI * x[1]) - 2.0 * x[0] * x[1]; 151 return PETSC_SUCCESS; 152 } 153 154 static PetscErrorCode trig_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 155 { 156 u[0] = PetscSinReal(2.0 * PETSC_PI * x[0]); 157 u[1] = PetscSinReal(2.0 * PETSC_PI * x[1]) - 2.0 * x[0] * x[1]; 158 u[2] = PetscSinReal(2.0 * PETSC_PI * x[2]) - 2.0 * x[1] * x[2]; 159 return PETSC_SUCCESS; 160 } 161 162 /* 163 u = sin(2 pi x) 164 v = sin(2 pi y) - 2xy 165 Delta <u,v> - f = <-4 pi^2 u, -4 pi^2 v> - <-4 pi^2 sin(2 pi x), -4 pi^2 sin(2 pi y)> 166 167 u = sin(2 pi x) 168 v = sin(2 pi y) - 2xy 169 w = sin(2 pi z) - 2yz 170 Delta <u,v,2> - f = <-4 pi^2 u, -4 pi^2 v, -4 pi^2 w> - <-4 pi^2 sin(2 pi x), -4 pi^2 sin(2 pi y), -4 pi^2 sin(2 pi z)> 171 */ 172 static void f0_vlap_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[]) 173 { 174 PetscInt d; 175 for (d = 0; d < dim; ++d) f0[d] += -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]); 176 } 177 178 /* 179 u = sin(2 pi x) 180 v = sin(2 pi y) - 2xy 181 \varepsilon = / 2 pi cos(2 pi x) -y \ 182 \ -y 2 pi cos(2 pi y) - 2x / 183 Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x 184 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij} 185 = \lambda \partial_j 2 pi (cos(2 pi x) + cos(2 pi y)) + 2\mu < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) > 186 = \lambda < -4 pi^2 sin(2 pi x) - 2, -4 pi^2 sin(2 pi y) > + \mu < -8 pi^2 sin(2 pi x) - 2, -8 pi^2 sin(2 pi y) > 187 188 u = sin(2 pi x) 189 v = sin(2 pi y) - 2xy 190 w = sin(2 pi z) - 2yz 191 \varepsilon = / 2 pi cos(2 pi x) -y 0 \ 192 | -y 2 pi cos(2 pi y) - 2x -z | 193 \ 0 -z 2 pi cos(2 pi z) - 2y / 194 Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2 y 195 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij} 196 = \lambda \partial_j (2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2 y) + 2\mu < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) - 1, -4 pi^2 sin(2 pi z) > 197 = \lambda < -4 pi^2 sin(2 pi x) - 2, -4 pi^2 sin(2 pi y) - 2, -4 pi^2 sin(2 pi z) > + 2\mu < -4 pi^2 sin(2 pi x) - 1, -4 pi^2 sin(2 pi y) - 1, -4 pi^2 sin(2 pi z) > 198 */ 199 static void f0_elas_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[]) 200 { 201 const PetscReal mu = PetscRealPart(constants[0]); 202 const PetscReal lambda = PetscRealPart(constants[1]); 203 const PetscReal fact = 4.0 * PetscSqr(PETSC_PI); 204 205 for (PetscInt d = 0; d < dim; ++d) f0[d] += -(2.0 * mu + lambda) * fact * PetscSinReal(2.0 * PETSC_PI * x[d]) - (d < dim - 1 ? 2.0 * (mu + lambda) : 0.0); 206 } 207 208 static PetscErrorCode axial_disp_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 209 { 210 AppCtx *user = (AppCtx *)ctx; 211 Parameter *param; 212 213 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 214 { 215 const PetscReal mu = PetscRealPart(param->mu); 216 const PetscReal lambda = PetscRealPart(param->lambda); 217 const PetscReal N = PetscRealPart(param->N); 218 219 u[0] = (3. * lambda * lambda + 8. * lambda * mu + 4 * mu * mu) / (4 * mu * (3 * lambda * lambda + 5. * lambda * mu + 2 * mu * mu)) * N * x[0]; 220 u[1] = 0.25 * lambda / mu / (lambda + mu) * N * x[1]; 221 for (PetscInt d = 2; d < dim; ++d) u[d] = 0.0; 222 } 223 return PETSC_SUCCESS; 224 } 225 226 /* 227 We will pull/push on the right side of a block of linearly elastic material. The uniform traction conditions on the 228 right side of the box will result in a uniform strain along x and y. The Neumann BC is given by 229 230 n_i \sigma_{ij} = t_i 231 232 u = (1/(2\mu) - 1) x 233 v = -y 234 f = 0 235 t = <4\mu/\lambda (\lambda + \mu), 0> 236 \varepsilon = / 1/(2\mu) - 1 0 \ 237 \ 0 -1 / 238 Tr(\varepsilon) = div u = 1/(2\mu) - 2 239 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij} 240 = \lambda \partial_j (1/(2\mu) - 2) + 2\mu < 0, 0 > 241 = \lambda < 0, 0 > + \mu < 0, 0 > = 0 242 NBC = <1,0> . <4\mu/\lambda (\lambda + \mu), 0> = 4\mu/\lambda (\lambda + \mu) 243 244 u = x - 1/2 245 v = 0 246 w = 0 247 \varepsilon = / x 0 0 \ 248 | 0 0 0 | 249 \ 0 0 0 / 250 Tr(\varepsilon) = div u = x 251 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij} 252 = \lambda \partial_j x + 2\mu < 1, 0, 0 > 253 = \lambda < 1, 0, 0 > + \mu < 2, 0, 0 > 254 */ 255 static void f0_elas_axial_disp_bd_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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 256 { 257 const PetscReal N = PetscRealPart(constants[2]); 258 259 f0[0] = N; 260 } 261 262 static PetscErrorCode uniform_strain_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 263 { 264 const PetscReal eps_xx = 0.1; 265 const PetscReal eps_xy = 0.3; 266 const PetscReal eps_yy = 0.25; 267 PetscInt d; 268 269 u[0] = eps_xx * x[0] + eps_xy * x[1]; 270 u[1] = eps_xy * x[0] + eps_yy * x[1]; 271 for (d = 2; d < dim; ++d) u[d] = 0.0; 272 return PETSC_SUCCESS; 273 } 274 275 static void f0_mass_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[]) 276 { 277 const PetscInt Nc = dim; 278 PetscInt c; 279 280 for (c = 0; c < Nc; ++c) f0[c] = u[c]; 281 } 282 283 static void f1_vlap_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 f1[]) 284 { 285 const PetscInt Nc = dim; 286 PetscInt c, d; 287 288 for (c = 0; c < Nc; ++c) 289 for (d = 0; d < dim; ++d) f1[c * dim + d] += u_x[c * dim + d]; 290 } 291 292 static void f1_elas_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 f1[]) 293 { 294 const PetscReal mu = PetscRealPart(constants[0]); 295 const PetscReal lambda = PetscRealPart(constants[1]); 296 const PetscInt Nc = dim; 297 298 for (PetscInt c = 0; c < Nc; ++c) { 299 for (PetscInt d = 0; d < dim; ++d) { 300 f1[c * dim + d] += mu * (u_x[c * dim + d] + u_x[d * dim + c]); 301 f1[c * dim + c] += lambda * u_x[d * dim + d]; 302 } 303 } 304 } 305 306 static void g0_mass_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 g0[]) 307 { 308 const PetscInt Nc = dim; 309 PetscInt c; 310 311 for (c = 0; c < Nc; ++c) g0[c * Nc + c] = 1.0; 312 } 313 314 static void g3_vlap_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[]) 315 { 316 const PetscInt Nc = dim; 317 PetscInt c, d; 318 319 for (c = 0; c < Nc; ++c) { 320 for (d = 0; d < dim; ++d) g3[((c * Nc + c) * dim + d) * dim + d] = 1.0; 321 } 322 } 323 324 /* 325 \partial_df \phi_fc g_{fc,gc,df,dg} \partial_dg \phi_gc 326 327 \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg} 328 = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc 329 */ 330 static void g3_elas_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[]) 331 { 332 const PetscReal mu = PetscRealPart(constants[0]); 333 const PetscReal lambda = PetscRealPart(constants[1]); 334 const PetscInt Nc = dim; 335 336 for (PetscInt c = 0; c < Nc; ++c) { 337 for (PetscInt d = 0; d < dim; ++d) { 338 g3[((c * Nc + c) * dim + d) * dim + d] += mu; 339 g3[((c * Nc + d) * dim + d) * dim + c] += mu; 340 g3[((c * Nc + d) * dim + c) * dim + d] += lambda; 341 } 342 } 343 } 344 345 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 346 { 347 PetscInt sol = 0, def = 0; 348 349 PetscFunctionBeginUser; 350 options->deform = DEFORM_NONE; 351 options->solType = SOL_VLAP_QUADRATIC; 352 options->useNearNullspace = PETSC_TRUE; 353 PetscCall(PetscStrncpy(options->dmType, DMPLEX, 256)); 354 355 PetscOptionsBegin(comm, "", "Linear Elasticity Problem Options", "DMPLEX"); 356 PetscCall(PetscOptionsEList("-deform_type", "Type of domain deformation", "ex17.c", deformTypes, NUM_DEFORM_TYPES, deformTypes[options->deform], &def, NULL)); 357 options->deform = (DeformType)def; 358 PetscCall(PetscOptionsEList("-sol_type", "Type of exact solution", "ex17.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL)); 359 options->solType = (SolutionType)sol; 360 PetscCall(PetscOptionsBool("-near_nullspace", "Use the rigid body modes as an AMG near nullspace", "ex17.c", options->useNearNullspace, &options->useNearNullspace, NULL)); 361 PetscCall(PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex17.c", DMList, options->dmType, options->dmType, 256, NULL)); 362 PetscOptionsEnd(); 363 PetscFunctionReturn(PETSC_SUCCESS); 364 } 365 366 static PetscErrorCode SetupParameters(MPI_Comm comm, AppCtx *ctx) 367 { 368 PetscBag bag; 369 Parameter *p; 370 371 PetscFunctionBeginUser; 372 /* setup PETSc parameter bag */ 373 PetscCall(PetscBagGetData(ctx->bag, (void **)&p)); 374 PetscCall(PetscBagSetName(ctx->bag, "par", "Elastic Parameters")); 375 bag = ctx->bag; 376 PetscCall(PetscBagRegisterScalar(bag, &p->mu, 1.0, "mu", "Shear Modulus, Pa")); 377 PetscCall(PetscBagRegisterScalar(bag, &p->lambda, 1.0, "lambda", "Lame's first parameter, Pa")); 378 PetscCall(PetscBagRegisterScalar(bag, &p->N, -1.0, "N", "Tension on right wall, Pa")); 379 PetscCall(PetscBagSetFromOptions(bag)); 380 { 381 PetscViewer viewer; 382 PetscViewerFormat format; 383 PetscBool flg; 384 385 PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-param_view", &viewer, &format, &flg)); 386 if (flg) { 387 PetscCall(PetscViewerPushFormat(viewer, format)); 388 PetscCall(PetscBagView(bag, viewer)); 389 PetscCall(PetscViewerFlush(viewer)); 390 PetscCall(PetscViewerPopFormat(viewer)); 391 PetscCall(PetscViewerDestroy(&viewer)); 392 } 393 } 394 PetscFunctionReturn(PETSC_SUCCESS); 395 } 396 397 static PetscErrorCode DMPlexDistortGeometry(DM dm) 398 { 399 DM cdm; 400 DMLabel label; 401 Vec coordinates; 402 PetscScalar *coords; 403 PetscReal mid = 0.5; 404 PetscInt cdim, d, vStart, vEnd, v; 405 406 PetscFunctionBeginUser; 407 PetscCall(DMGetCoordinateDM(dm, &cdm)); 408 PetscCall(DMGetCoordinateDim(dm, &cdim)); 409 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 410 PetscCall(DMGetLabel(dm, "marker", &label)); 411 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 412 PetscCall(VecGetArrayWrite(coordinates, &coords)); 413 for (v = vStart; v < vEnd; ++v) { 414 PetscScalar *pcoords, shift; 415 PetscInt val; 416 417 PetscCall(DMLabelGetValue(label, v, &val)); 418 if (val >= 0) continue; 419 PetscCall(DMPlexPointLocalRef(cdm, v, coords, &pcoords)); 420 shift = 0.2 * PetscAbsScalar(pcoords[0] - mid); 421 shift = PetscRealPart(pcoords[0]) > mid ? shift : -shift; 422 for (d = 1; d < cdim; ++d) pcoords[d] += shift; 423 } 424 PetscCall(VecRestoreArrayWrite(coordinates, &coords)); 425 PetscFunctionReturn(PETSC_SUCCESS); 426 } 427 428 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 429 { 430 PetscFunctionBeginUser; 431 PetscCall(DMCreate(comm, dm)); 432 PetscCall(DMSetType(*dm, DMPLEX)); 433 PetscCall(DMSetFromOptions(*dm)); 434 switch (user->deform) { 435 case DEFORM_NONE: 436 break; 437 case DEFORM_SHEAR: 438 PetscCall(DMPlexShearGeometry(*dm, DM_X, NULL)); 439 break; 440 case DEFORM_STEP: 441 PetscCall(DMPlexDistortGeometry(*dm)); 442 break; 443 default: 444 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid deformation type: %s (%d)", deformTypes[PetscMin(user->deform, NUM_DEFORM_TYPES)], user->deform); 445 } 446 PetscCall(DMSetApplicationContext(*dm, user)); 447 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 448 PetscFunctionReturn(PETSC_SUCCESS); 449 } 450 451 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 452 { 453 PetscErrorCode (*exact)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 454 Parameter *param; 455 PetscDS ds; 456 PetscWeakForm wf; 457 DMLabel label; 458 PetscInt id, bd; 459 PetscInt dim; 460 461 PetscFunctionBeginUser; 462 PetscCall(DMGetDS(dm, &ds)); 463 PetscCall(PetscDSGetWeakForm(ds, &wf)); 464 PetscCall(PetscDSGetSpatialDimension(ds, &dim)); 465 PetscCall(PetscBagGetData(user->bag, (void **)¶m)); 466 switch (user->solType) { 467 case SOL_MASS_QUADRATIC: 468 PetscCall(PetscDSSetResidual(ds, 0, f0_mass_u, NULL)); 469 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_mass_uu, NULL, NULL, NULL)); 470 PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 1, f0_mass_quadratic_u, 0, NULL)); 471 switch (dim) { 472 case 2: 473 exact = quadratic_2d_u; 474 break; 475 case 3: 476 exact = quadratic_3d_u; 477 break; 478 default: 479 SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %" PetscInt_FMT, dim); 480 } 481 break; 482 case SOL_VLAP_QUADRATIC: 483 PetscCall(PetscDSSetResidual(ds, 0, f0_vlap_quadratic_u, f1_vlap_u)); 484 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_vlap_uu)); 485 switch (dim) { 486 case 2: 487 exact = quadratic_2d_u; 488 break; 489 case 3: 490 exact = quadratic_3d_u; 491 break; 492 default: 493 SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %" PetscInt_FMT, dim); 494 } 495 break; 496 case SOL_ELAS_QUADRATIC: 497 PetscCall(PetscDSSetResidual(ds, 0, f0_elas_quadratic_u, f1_elas_u)); 498 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu)); 499 switch (dim) { 500 case 2: 501 exact = quadratic_2d_u; 502 break; 503 case 3: 504 exact = quadratic_3d_u; 505 break; 506 default: 507 SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %" PetscInt_FMT, dim); 508 } 509 break; 510 case SOL_VLAP_TRIG: 511 PetscCall(PetscDSSetResidual(ds, 0, f0_vlap_trig_u, f1_vlap_u)); 512 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_vlap_uu)); 513 switch (dim) { 514 case 2: 515 exact = trig_2d_u; 516 break; 517 case 3: 518 exact = trig_3d_u; 519 break; 520 default: 521 SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %" PetscInt_FMT, dim); 522 } 523 break; 524 case SOL_ELAS_TRIG: 525 PetscCall(PetscDSSetResidual(ds, 0, f0_elas_trig_u, f1_elas_u)); 526 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu)); 527 switch (dim) { 528 case 2: 529 exact = trig_2d_u; 530 break; 531 case 3: 532 exact = trig_3d_u; 533 break; 534 default: 535 SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid dimension: %" PetscInt_FMT, dim); 536 } 537 break; 538 case SOL_ELAS_AXIAL_DISP: 539 PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_elas_u)); 540 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu)); 541 id = dim == 3 ? 5 : 2; 542 PetscCall(DMGetLabel(dm, "marker", &label)); 543 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right", label, 1, &id, 0, 0, NULL, (void (*)(void))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_elas_axial_disp_bd_u, 0, NULL)); 546 exact = axial_disp_u; 547 break; 548 case SOL_ELAS_UNIFORM_STRAIN: 549 PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_elas_u)); 550 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu)); 551 exact = uniform_strain_u; 552 break; 553 case SOL_ELAS_GE: 554 PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_elas_u)); 555 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_elas_uu)); 556 exact = zero; /* No exact solution available */ 557 break; 558 default: 559 SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType); 560 } 561 PetscCall(PetscDSSetExactSolution(ds, 0, exact, user)); 562 PetscCall(DMGetLabel(dm, "marker", &label)); 563 if (user->solType == SOL_ELAS_AXIAL_DISP) { 564 PetscInt cmp; 565 566 id = dim == 3 ? 6 : 4; 567 cmp = 0; 568 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left", label, 1, &id, 0, 1, &cmp, (void (*)(void))zero, NULL, user, NULL)); 569 cmp = dim == 3 ? 2 : 1; 570 id = dim == 3 ? 1 : 1; 571 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom", label, 1, &id, 0, 1, &cmp, (void (*)(void))zero, NULL, user, NULL)); 572 if (dim == 3) { 573 cmp = 1; 574 id = 3; 575 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "front", label, 1, &id, 0, 1, &cmp, (void (*)(void))zero, NULL, user, NULL)); 576 } 577 } else if (user->solType == SOL_ELAS_GE) { 578 PetscInt cmp; 579 580 id = dim == 3 ? 6 : 4; 581 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "left", label, 1, &id, 0, 0, NULL, (void (*)(void))zero, NULL, user, NULL)); 582 id = dim == 3 ? 5 : 2; 583 cmp = 0; 584 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "right", label, 1, &id, 0, 1, &cmp, (void (*)(void))ge_shift, NULL, user, NULL)); 585 } else { 586 id = 1; 587 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))exact, NULL, user, NULL)); 588 } 589 /* Setup constants */ 590 { 591 PetscScalar constants[3]; 592 593 constants[0] = param->mu; /* shear modulus, Pa */ 594 constants[1] = param->lambda; /* Lame's first parameter, Pa */ 595 constants[2] = param->N; /* Tension on right wall, Pa */ 596 PetscCall(PetscDSSetConstants(ds, 3, constants)); 597 } 598 PetscFunctionReturn(PETSC_SUCCESS); 599 } 600 601 static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace) 602 { 603 PetscFunctionBegin; 604 PetscCall(DMPlexCreateRigidBody(dm, origField, nullspace)); 605 PetscFunctionReturn(PETSC_SUCCESS); 606 } 607 608 PetscErrorCode SetupFE(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx) 609 { 610 AppCtx *user = (AppCtx *)ctx; 611 DM cdm = dm; 612 PetscFE fe; 613 char prefix[PETSC_MAX_PATH_LEN]; 614 DMPolytopeType ct; 615 PetscInt dim, cStart; 616 617 PetscFunctionBegin; 618 /* Create finite element */ 619 PetscCall(DMGetDimension(dm, &dim)); 620 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 621 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 622 PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name)); 623 PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, dim, ct, name ? prefix : NULL, -1, &fe)); 624 PetscCall(PetscObjectSetName((PetscObject)fe, name)); 625 /* Set discretization and boundary conditions for each mesh */ 626 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 627 PetscCall(DMCreateDS(dm)); 628 PetscCall((*setup)(dm, user)); 629 while (cdm) { 630 PetscCall(DMCopyDisc(dm, cdm)); 631 if (user->useNearNullspace) PetscCall(DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace)); 632 PetscCall(DMGetCoarseDM(cdm, &cdm)); 633 } 634 PetscCall(PetscFEDestroy(&fe)); 635 PetscFunctionReturn(PETSC_SUCCESS); 636 } 637 638 int main(int argc, char **argv) 639 { 640 DM dm; /* Problem specification */ 641 SNES snes; /* Nonlinear solver */ 642 Vec u; /* Solutions */ 643 AppCtx user; /* User-defined work context */ 644 645 PetscFunctionBeginUser; 646 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 647 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 648 PetscCall(PetscBagCreate(PETSC_COMM_SELF, sizeof(Parameter), &user.bag)); 649 PetscCall(SetupParameters(PETSC_COMM_WORLD, &user)); 650 /* Primal system */ 651 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 652 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 653 PetscCall(SNESSetDM(snes, dm)); 654 PetscCall(SetupFE(dm, "displacement", SetupPrimalProblem, &user)); 655 PetscCall(DMCreateGlobalVector(dm, &u)); 656 PetscCall(VecSet(u, 0.0)); 657 PetscCall(PetscObjectSetName((PetscObject)u, "displacement")); 658 PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user)); 659 PetscCall(SNESSetFromOptions(snes)); 660 PetscCall(DMSNESCheckFromOptions(snes, u)); 661 PetscCall(SNESSolve(snes, NULL, u)); 662 PetscCall(SNESGetSolution(snes, &u)); 663 PetscCall(VecViewFromOptions(u, NULL, "-displacement_view")); 664 /* Cleanup */ 665 PetscCall(VecDestroy(&u)); 666 PetscCall(SNESDestroy(&snes)); 667 PetscCall(DMDestroy(&dm)); 668 PetscCall(PetscBagDestroy(&user.bag)); 669 PetscCall(PetscFinalize()); 670 return 0; 671 } 672 673 /*TEST 674 675 testset: 676 args: -dm_plex_box_faces 1,1,1 677 678 test: 679 suffix: 2d_p1_quad_vlap 680 requires: triangle 681 args: -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate 682 test: 683 suffix: 2d_p2_quad_vlap 684 requires: triangle 685 args: -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 686 test: 687 suffix: 2d_p3_quad_vlap 688 requires: triangle 689 args: -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 690 test: 691 suffix: 2d_q1_quad_vlap 692 args: -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate 693 test: 694 suffix: 2d_q2_quad_vlap 695 args: -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 696 test: 697 suffix: 2d_q3_quad_vlap 698 requires: !single 699 args: -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 700 test: 701 suffix: 2d_p1_quad_elas 702 requires: triangle 703 args: -sol_type elas_quad -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate 704 test: 705 suffix: 2d_p2_quad_elas 706 requires: triangle 707 args: -sol_type elas_quad -displacement_petscspace_degree 2 -dmsnes_check .0001 708 test: 709 suffix: 2d_p3_quad_elas 710 requires: triangle 711 args: -sol_type elas_quad -displacement_petscspace_degree 3 -dmsnes_check .0001 712 test: 713 suffix: 2d_q1_quad_elas 714 args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 715 test: 716 suffix: 2d_q1_quad_elas_shear 717 args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 718 test: 719 suffix: 2d_q2_quad_elas 720 args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dmsnes_check .0001 721 test: 722 suffix: 2d_q2_quad_elas_shear 723 args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 2 -dmsnes_check 724 test: 725 suffix: 2d_q3_quad_elas 726 args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dmsnes_check .0001 727 test: 728 suffix: 2d_q3_quad_elas_shear 729 requires: !single 730 args: -sol_type elas_quad -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 3 -dmsnes_check 731 732 test: 733 suffix: 3d_p1_quad_vlap 734 requires: ctetgen 735 args: -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 736 test: 737 suffix: 3d_p2_quad_vlap 738 requires: ctetgen 739 args: -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001 740 test: 741 suffix: 3d_p3_quad_vlap 742 requires: ctetgen 743 args: -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001 744 test: 745 suffix: 3d_q1_quad_vlap 746 args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 747 test: 748 suffix: 3d_q2_quad_vlap 749 args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001 750 test: 751 suffix: 3d_q3_quad_vlap 752 args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001 753 test: 754 suffix: 3d_p1_quad_elas 755 requires: ctetgen 756 args: -sol_type elas_quad -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 757 test: 758 suffix: 3d_p2_quad_elas 759 requires: ctetgen 760 args: -sol_type elas_quad -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001 761 test: 762 suffix: 3d_p3_quad_elas 763 requires: ctetgen 764 args: -sol_type elas_quad -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001 765 test: 766 suffix: 3d_q1_quad_elas 767 args: -sol_type elas_quad -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 768 test: 769 suffix: 3d_q2_quad_elas 770 args: -sol_type elas_quad -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001 771 test: 772 suffix: 3d_q3_quad_elas 773 requires: !single 774 args: -sol_type elas_quad -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001 775 776 test: 777 suffix: 2d_p1_trig_vlap 778 requires: triangle 779 args: -sol_type vlap_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 780 test: 781 suffix: 2d_p2_trig_vlap 782 requires: triangle 783 args: -sol_type vlap_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 784 test: 785 suffix: 2d_p3_trig_vlap 786 requires: triangle 787 args: -sol_type vlap_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 788 test: 789 suffix: 2d_q1_trig_vlap 790 args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 791 test: 792 suffix: 2d_q2_trig_vlap 793 args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 794 test: 795 suffix: 2d_q3_trig_vlap 796 args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 797 test: 798 suffix: 2d_p1_trig_elas 799 requires: triangle 800 args: -sol_type elas_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 801 test: 802 suffix: 2d_p2_trig_elas 803 requires: triangle 804 args: -sol_type elas_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 805 test: 806 suffix: 2d_p3_trig_elas 807 requires: triangle 808 args: -sol_type elas_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 809 test: 810 suffix: 2d_q1_trig_elas 811 args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 812 test: 813 suffix: 2d_q1_trig_elas_shear 814 args: -sol_type elas_trig -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 815 test: 816 suffix: 2d_q2_trig_elas 817 args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 818 test: 819 suffix: 2d_q2_trig_elas_shear 820 args: -sol_type elas_trig -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 821 test: 822 suffix: 2d_q3_trig_elas 823 args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 824 test: 825 suffix: 2d_q3_trig_elas_shear 826 args: -sol_type elas_trig -dm_plex_simplex 0 -deform_type shear -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 827 828 test: 829 suffix: 3d_p1_trig_vlap 830 requires: ctetgen 831 args: -sol_type vlap_trig -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 832 test: 833 suffix: 3d_p2_trig_vlap 834 requires: ctetgen 835 args: -sol_type vlap_trig -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 836 test: 837 suffix: 3d_p3_trig_vlap 838 requires: ctetgen 839 args: -sol_type vlap_trig -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 840 test: 841 suffix: 3d_q1_trig_vlap 842 args: -sol_type vlap_trig -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 843 test: 844 suffix: 3d_q2_trig_vlap 845 args: -sol_type vlap_trig -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 846 test: 847 suffix: 3d_q3_trig_vlap 848 requires: !__float128 849 args: -sol_type vlap_trig -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 850 test: 851 suffix: 3d_p1_trig_elas 852 requires: ctetgen 853 args: -sol_type elas_trig -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 854 test: 855 suffix: 3d_p2_trig_elas 856 requires: ctetgen 857 args: -sol_type elas_trig -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 858 test: 859 suffix: 3d_p3_trig_elas 860 requires: ctetgen 861 args: -sol_type elas_trig -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 862 test: 863 suffix: 3d_q1_trig_elas 864 args: -sol_type elas_trig -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 865 test: 866 suffix: 3d_q2_trig_elas 867 args: -sol_type elas_trig -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 868 test: 869 suffix: 3d_q3_trig_elas 870 requires: !__float128 871 args: -sol_type elas_trig -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 872 873 test: 874 suffix: 2d_p1_axial_elas 875 requires: triangle 876 args: -sol_type elas_axial_disp -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 2 -dmsnes_check .0001 -pc_type lu 877 test: 878 suffix: 2d_p2_axial_elas 879 requires: triangle 880 args: -sol_type elas_axial_disp -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu 881 test: 882 suffix: 2d_p3_axial_elas 883 requires: triangle 884 args: -sol_type elas_axial_disp -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu 885 test: 886 suffix: 2d_q1_axial_elas 887 args: -sol_type elas_axial_disp -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 1 -dmsnes_check .0001 -pc_type lu 888 test: 889 suffix: 2d_q2_axial_elas 890 args: -sol_type elas_axial_disp -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu 891 test: 892 suffix: 2d_q3_axial_elas 893 args: -sol_type elas_axial_disp -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu 894 895 test: 896 suffix: 2d_p1_uniform_elas 897 requires: triangle 898 args: -sol_type elas_uniform_strain -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 899 test: 900 suffix: 2d_p2_uniform_elas 901 requires: triangle 902 args: -sol_type elas_uniform_strain -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 903 test: 904 suffix: 2d_p3_uniform_elas 905 requires: triangle 906 args: -sol_type elas_uniform_strain -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 907 test: 908 suffix: 2d_q1_uniform_elas 909 args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 910 test: 911 suffix: 2d_q2_uniform_elas 912 requires: !single 913 args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 914 test: 915 suffix: 2d_q3_uniform_elas 916 requires: !single 917 args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 918 test: 919 suffix: 2d_p1_uniform_elas_step 920 requires: triangle 921 args: -sol_type elas_uniform_strain -deform_type step -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 922 923 testset: 924 args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -deform_type step -displacement_petscspace_degree 1 -dmsnes_check .0001 -pc_type lu 925 926 test: 927 suffix: 2d_q1_uniform_elas_step 928 args: -sol_type elas_uniform_strain -dm_refine 2 929 test: 930 suffix: 2d_q1_quad_vlap_step 931 args: 932 test: 933 suffix: 2d_q2_quad_vlap_step 934 args: -displacement_petscspace_degree 2 935 test: 936 suffix: 2d_q1_quad_mass_step 937 args: -sol_type mass_quad 938 939 testset: 940 filter: grep -v "variant HERMITIAN" 941 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_lower -5,-5,-0.25 -dm_plex_box_upper 5,5,0.25 \ 942 -dm_plex_box_faces 5,5,2 -dm_plex_separate_marker -dm_refine 0 -petscpartitioner_type simple \ 943 -sol_type elas_ge 944 945 test: 946 suffix: ge_q1_0 947 args: -displacement_petscspace_degree 1 \ 948 -snes_max_it 2 -snes_rtol 1.e-10 \ 949 -ksp_type cg -ksp_rtol 1.e-10 -ksp_max_it 100 -ksp_norm_type unpreconditioned \ 950 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 \ 951 -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true \ 952 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 \ 953 -mg_levels_ksp_max_it 2 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.1 -mg_levels_pc_type jacobi \ 954 -matptap_via scalable 955 test: 956 suffix: ge_q1_gmg 957 args: -displacement_petscspace_degree 1 \ 958 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 959 -snes_max_it 2 -snes_rtol 1.e-10 \ 960 -ksp_type cg -ksp_rtol 1.e-10 -ksp_max_it 100 -ksp_norm_type unpreconditioned \ 961 -pc_type mg -pc_mg_type full \ 962 -mg_levels_ksp_max_it 4 -mg_levels_esteig_ksp_type cg \ 963 -mg_levels_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \ 964 -mg_levels_pc_type jacobi 965 test: 966 nsize: 5 967 suffix: ge_q1_gdsw 968 args: -snes_max_it 1 -ksp_type cg -ksp_norm_type natural -displacement_petscspace_degree 1 -snes_monitor_short -ksp_monitor_short -pc_type mg -pc_mg_adapt_interp_coarse_space gdsw -pc_mg_levels 2 -pc_mg_galerkin -mg_levels_pc_type bjacobi -mg_levels_esteig_ksp_type cg -mg_levels_sub_pc_type icc -mg_coarse_redundant_pc_type cholesky -ksp_view 969 970 TEST*/ 971