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