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