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