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 11 #include <petscdmplex.h> 12 #include <petscsnes.h> 13 #include <petscds.h> 14 #include <petscconvest.h> 15 16 typedef enum {SOL_VLAP_QUADRATIC, SOL_ELAS_QUADRATIC, SOL_VLAP_TRIG, SOL_ELAS_TRIG, SOL_ELAS_AXIAL_DISP, SOL_ELAS_UNIFORM_STRAIN, NUM_SOLUTION_TYPES} SolutionType; 17 const char *solutionTypes[NUM_SOLUTION_TYPES+1] = {"vlap_quad", "elas_quad", "vlap_trig", "elas_trig", "elas_axial_disp", "elas_uniform_strain", "unknown"}; 18 19 typedef struct { 20 /* Domain and mesh definition */ 21 char dmType[256]; /* DM type for the solve */ 22 PetscBool shear; /* Shear the domain */ 23 /* Problem definition */ 24 SolutionType solType; /* Type of exact solution */ 25 /* Solver definition */ 26 PetscBool useNearNullspace; /* Use the rigid body modes as a near nullspace for AMG */ 27 } AppCtx; 28 29 static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 30 { 31 PetscInt d; 32 for (d = 0; d < dim; ++d) u[d] = 0.0; 33 return 0; 34 } 35 36 static PetscErrorCode quadratic_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 37 { 38 u[0] = x[0]*x[0]; 39 u[1] = x[1]*x[1] - 2.0*x[0]*x[1]; 40 return 0; 41 } 42 43 static PetscErrorCode quadratic_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 44 { 45 u[0] = x[0]*x[0]; 46 u[1] = x[1]*x[1] - 2.0*x[0]*x[1]; 47 u[2] = x[2]*x[2] - 2.0*x[1]*x[2]; 48 return 0; 49 } 50 51 /* 52 u = x^2 53 v = y^2 - 2xy 54 Delta <u,v> - f = <2, 2> - <2, 2> 55 56 u = x^2 57 v = y^2 - 2xy 58 w = z^2 - 2yz 59 Delta <u,v,w> - f = <2, 2, 2> - <2, 2, 2> 60 */ 61 static void f0_vlap_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 62 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 63 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 64 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 65 { 66 PetscInt d; 67 for (d = 0; d < dim; ++d) f0[d] += 2.0; 68 } 69 70 /* 71 u = x^2 72 v = y^2 - 2xy 73 \varepsilon = / 2x -y \ 74 \ -y 2y - 2x / 75 Tr(\varepsilon) = div u = 2y 76 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij} 77 = \lambda \partial_j (2y) + 2\mu < 2-1, 2 > 78 = \lambda < 0, 2 > + \mu < 2, 4 > 79 80 u = x^2 81 v = y^2 - 2xy 82 w = z^2 - 2yz 83 \varepsilon = / 2x -y 0 \ 84 | -y 2y - 2x -z | 85 \ 0 -z 2z - 2y/ 86 Tr(\varepsilon) = div u = 2z 87 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij} 88 = \lambda \partial_j (2z) + 2\mu < 2-1, 2-1, 2 > 89 = \lambda < 0, 0, 2 > + \mu < 2, 2, 4 > 90 */ 91 static void f0_elas_quadratic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 92 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 93 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 94 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 95 { 96 const PetscReal mu = 1.0; 97 const PetscReal lambda = 1.0; 98 PetscInt d; 99 100 for (d = 0; d < dim-1; ++d) f0[d] += 2.0*mu; 101 f0[dim-1] += 2.0*lambda + 4.0*mu; 102 } 103 104 static PetscErrorCode trig_2d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 105 { 106 u[0] = PetscSinReal(2.0*PETSC_PI*x[0]); 107 u[1] = PetscSinReal(2.0*PETSC_PI*x[1]) - 2.0*x[0]*x[1]; 108 return 0; 109 } 110 111 static PetscErrorCode trig_3d_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 112 { 113 u[0] = PetscSinReal(2.0*PETSC_PI*x[0]); 114 u[1] = PetscSinReal(2.0*PETSC_PI*x[1]) - 2.0*x[0]*x[1]; 115 u[2] = PetscSinReal(2.0*PETSC_PI*x[2]) - 2.0*x[1]*x[2]; 116 return 0; 117 } 118 119 /* 120 u = sin(2 pi x) 121 v = sin(2 pi y) - 2xy 122 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)> 123 124 u = sin(2 pi x) 125 v = sin(2 pi y) - 2xy 126 w = sin(2 pi z) - 2yz 127 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)> 128 */ 129 static void f0_vlap_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 130 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 131 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 132 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 133 { 134 PetscInt d; 135 for (d = 0; d < dim; ++d) f0[d] += -4.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[d]); 136 } 137 138 /* 139 u = sin(2 pi x) 140 v = sin(2 pi y) - 2xy 141 \varepsilon = / 2 pi cos(2 pi x) -y \ 142 \ -y 2 pi cos(2 pi y) - 2x / 143 Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y)) - 2 x 144 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij} 145 = \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) > 146 = \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) > 147 148 u = sin(2 pi x) 149 v = sin(2 pi y) - 2xy 150 w = sin(2 pi z) - 2yz 151 \varepsilon = / 2 pi cos(2 pi x) -y 0 \ 152 | -y 2 pi cos(2 pi y) - 2x -z | 153 \ 0 -z 2 pi cos(2 pi z) - 2y / 154 Tr(\varepsilon) = div u = 2 pi (cos(2 pi x) + cos(2 pi y) + cos(2 pi z)) - 2 x - 2 y 155 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij} 156 = \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) > 157 = \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) > 158 */ 159 static void f0_elas_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 160 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 161 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 162 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 163 { 164 const PetscReal mu = 1.0; 165 const PetscReal lambda = 1.0; 166 const PetscReal fact = 4.0*PetscSqr(PETSC_PI); 167 PetscInt d; 168 169 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); 170 } 171 172 static PetscErrorCode axial_disp_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 173 { 174 const PetscReal mu = 1.0; 175 const PetscReal lambda = 1.0; 176 const PetscReal N = 1.0; 177 PetscInt d; 178 179 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]; 180 u[1] = -0.25*lambda/mu/(lambda+mu)*N*x[1]; 181 for (d = 2; d < dim; ++d) u[d] = 0.0; 182 return 0; 183 } 184 185 /* 186 We will pull/push on the right side of a block of linearly elastic material. The uniform traction conditions on the 187 right side of the box will result in a uniform strain along x and y. The Neumann BC is given by 188 189 n_i \sigma_{ij} = t_i 190 191 u = (1/(2\mu) - 1) x 192 v = -y 193 f = 0 194 t = <4\mu/\lambda (\lambda + \mu), 0> 195 \varepsilon = / 1/(2\mu) - 1 0 \ 196 \ 0 -1 / 197 Tr(\varepsilon) = div u = 1/(2\mu) - 2 198 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij} 199 = \lambda \partial_j (1/(2\mu) - 2) + 2\mu < 0, 0 > 200 = \lambda < 0, 0 > + \mu < 0, 0 > = 0 201 NBC = <1,0> . <4\mu/\lambda (\lambda + \mu), 0> = 4\mu/\lambda (\lambda + \mu) 202 203 u = x - 1/2 204 v = 0 205 w = 0 206 \varepsilon = / x 0 0 \ 207 | 0 0 0 | 208 \ 0 0 0 / 209 Tr(\varepsilon) = div u = x 210 div \sigma = \partial_i \lambda \delta_{ij} \varepsilon_{kk} + \partial_i 2\mu\varepsilon_{ij} 211 = \lambda \partial_j x + 2\mu < 1, 0, 0 > 212 = \lambda < 1, 0, 0 > + \mu < 2, 0, 0 > 213 */ 214 static void f0_elas_axial_disp_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 215 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 216 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 217 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 218 { 219 const PetscReal N = -1.0; 220 221 f0[0] = N; 222 } 223 224 static PetscErrorCode uniform_strain_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 225 { 226 const PetscReal eps_xx = 0.1; 227 const PetscReal eps_xy = 0.3; 228 const PetscReal eps_yy = 0.25; 229 PetscInt d; 230 231 u[0] = eps_xx*x[0] + eps_xy*x[1]; 232 u[1] = eps_xy*x[0] + eps_yy*x[1]; 233 for (d = 2; d < dim; ++d) u[d] = 0.0; 234 return 0; 235 } 236 237 static void f1_vlap_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 238 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 239 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 240 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 241 { 242 const PetscInt Nc = dim; 243 PetscInt c, d; 244 245 for (c = 0; c < Nc; ++c) for (d = 0; d < dim; ++d) f1[c*dim+d] += u_x[c*dim+d]; 246 } 247 248 static void f1_elas_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 249 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 250 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 251 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 252 { 253 const PetscInt Nc = dim; 254 const PetscReal mu = 1.0; 255 const PetscReal lambda = 1.0; 256 PetscInt c, d; 257 258 for (c = 0; c < Nc; ++c) { 259 for (d = 0; d < dim; ++d) { 260 f1[c*dim+d] += mu*(u_x[c*dim+d] + u_x[d*dim+c]); 261 f1[c*dim+c] += lambda*u_x[d*dim+d]; 262 } 263 } 264 } 265 266 static void g3_vlap_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 267 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 268 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 269 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 270 { 271 const PetscInt Nc = dim; 272 PetscInt c, d; 273 274 for (c = 0; c < Nc; ++c) { 275 for (d = 0; d < dim; ++d) { 276 g3[((c*Nc + c)*dim + d)*dim + d] = 1.0; 277 } 278 } 279 } 280 281 /* 282 \partial_df \phi_fc g_{fc,gc,df,dg} \partial_dg \phi_gc 283 284 \partial_df \phi_fc \lambda \delta_{fc,df} \sum_gc \partial_dg \phi_gc \delta_{gc,dg} 285 = \partial_fc \phi_fc \sum_gc \partial_gc \phi_gc 286 */ 287 static void g3_elas_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 288 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 289 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 290 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 291 { 292 const PetscInt Nc = dim; 293 const PetscReal mu = 1.0; 294 const PetscReal lambda = 1.0; 295 PetscInt c, d; 296 297 for (c = 0; c < Nc; ++c) { 298 for (d = 0; d < dim; ++d) { 299 g3[((c*Nc + c)*dim + d)*dim + d] += mu; 300 g3[((c*Nc + d)*dim + d)*dim + c] += mu; 301 g3[((c*Nc + d)*dim + c)*dim + d] += lambda; 302 } 303 } 304 } 305 306 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 307 { 308 PetscInt sol; 309 PetscErrorCode ierr; 310 311 PetscFunctionBeginUser; 312 options->shear = PETSC_FALSE; 313 options->solType = SOL_VLAP_QUADRATIC; 314 options->useNearNullspace = PETSC_TRUE; 315 ierr = PetscStrncpy(options->dmType, DMPLEX, 256);CHKERRQ(ierr); 316 317 ierr = PetscOptionsBegin(comm, "", "Linear Elasticity Problem Options", "DMPLEX");CHKERRQ(ierr); 318 ierr = PetscOptionsBool("-shear", "Shear the domain", "ex17.c", options->shear, &options->shear, NULL);CHKERRQ(ierr); 319 sol = options->solType; 320 ierr = PetscOptionsEList("-sol_type", "Type of exact solution", "ex17.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL);CHKERRQ(ierr); 321 options->solType = (SolutionType) sol; 322 ierr = PetscOptionsBool("-near_nullspace", "Use the rigid body modes as an AMG near nullspace", "ex17.c", options->useNearNullspace, &options->useNearNullspace, NULL);CHKERRQ(ierr); 323 ierr = PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex17.c", DMList, options->dmType, options->dmType, 256, NULL);CHKERRQ(ierr); 324 ierr = PetscOptionsEnd(); 325 PetscFunctionReturn(0); 326 } 327 328 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 329 { 330 PetscErrorCode ierr; 331 332 PetscFunctionBeginUser; 333 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 334 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 335 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 336 if (user->shear) {ierr = DMPlexShearGeometry(*dm, DM_X, NULL);CHKERRQ(ierr);} 337 ierr = DMSetApplicationContext(*dm, user);CHKERRQ(ierr); 338 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 339 PetscFunctionReturn(0); 340 } 341 342 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 343 { 344 PetscErrorCode (*exact)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 345 PetscDS prob; 346 PetscWeakForm wf; 347 DMLabel label; 348 PetscInt id, bd; 349 PetscInt dim; 350 PetscErrorCode ierr; 351 352 PetscFunctionBeginUser; 353 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 354 ierr = PetscDSGetSpatialDimension(prob, &dim);CHKERRQ(ierr); 355 switch (user->solType) { 356 case SOL_VLAP_QUADRATIC: 357 ierr = PetscDSSetResidual(prob, 0, f0_vlap_quadratic_u, f1_vlap_u);CHKERRQ(ierr); 358 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_vlap_uu);CHKERRQ(ierr); 359 switch (dim) { 360 case 2: exact = quadratic_2d_u;break; 361 case 3: exact = quadratic_3d_u;break; 362 default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim); 363 } 364 break; 365 case SOL_ELAS_QUADRATIC: 366 ierr = PetscDSSetResidual(prob, 0, f0_elas_quadratic_u, f1_elas_u);CHKERRQ(ierr); 367 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);CHKERRQ(ierr); 368 switch (dim) { 369 case 2: exact = quadratic_2d_u;break; 370 case 3: exact = quadratic_3d_u;break; 371 default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim); 372 } 373 break; 374 case SOL_VLAP_TRIG: 375 ierr = PetscDSSetResidual(prob, 0, f0_vlap_trig_u, f1_vlap_u);CHKERRQ(ierr); 376 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_vlap_uu);CHKERRQ(ierr); 377 switch (dim) { 378 case 2: exact = trig_2d_u;break; 379 case 3: exact = trig_3d_u;break; 380 default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim); 381 } 382 break; 383 case SOL_ELAS_TRIG: 384 ierr = PetscDSSetResidual(prob, 0, f0_elas_trig_u, f1_elas_u);CHKERRQ(ierr); 385 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);CHKERRQ(ierr); 386 switch (dim) { 387 case 2: exact = trig_2d_u;break; 388 case 3: exact = trig_3d_u;break; 389 default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim); 390 } 391 break; 392 case SOL_ELAS_AXIAL_DISP: 393 ierr = PetscDSSetResidual(prob, 0, NULL, f1_elas_u);CHKERRQ(ierr); 394 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);CHKERRQ(ierr); 395 id = dim == 3 ? 5 : 2; 396 ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr); 397 ierr = DMAddBoundary(dm, DM_BC_NATURAL, "right", label, 1, &id, 0, 0, NULL, (void (*)(void)) NULL, NULL, user, &bd);CHKERRQ(ierr); 398 ierr = PetscDSGetBoundary(prob, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 399 ierr = PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_elas_axial_disp_bd_u, 0, NULL);CHKERRQ(ierr); 400 exact = axial_disp_u; 401 break; 402 case SOL_ELAS_UNIFORM_STRAIN: 403 ierr = PetscDSSetResidual(prob, 0, NULL, f1_elas_u);CHKERRQ(ierr); 404 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);CHKERRQ(ierr); 405 exact = uniform_strain_u; 406 break; 407 default: SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType); 408 } 409 ierr = PetscDSSetExactSolution(prob, 0, exact, user);CHKERRQ(ierr); 410 ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr); 411 if (user->solType == SOL_ELAS_AXIAL_DISP) { 412 PetscInt cmp; 413 414 id = dim == 3 ? 6 : 4; 415 cmp = 0; 416 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "left", label, 1, &id, 0, 1, &cmp, (void (*)(void)) zero, NULL, user, NULL);CHKERRQ(ierr); 417 cmp = dim == 3 ? 2 : 1; 418 id = dim == 3 ? 1 : 1; 419 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom", label, 1, &id, 0, 1, &cmp, (void (*)(void)) zero, NULL, user, NULL);CHKERRQ(ierr); 420 if (dim == 3) { 421 cmp = 1; 422 id = 3; 423 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "front", label, 1, &id, 0, 1, &cmp, (void (*)(void)) zero, NULL, user, NULL);CHKERRQ(ierr); 424 } 425 } else { 426 id = 1; 427 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) exact, NULL, user, NULL);CHKERRQ(ierr); 428 } 429 PetscFunctionReturn(0); 430 } 431 432 static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace) 433 { 434 PetscErrorCode ierr; 435 436 PetscFunctionBegin; 437 ierr = DMPlexCreateRigidBody(dm, origField, nullspace);CHKERRQ(ierr); 438 PetscFunctionReturn(0); 439 } 440 441 PetscErrorCode SetupFE(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx) 442 { 443 AppCtx *user = (AppCtx *) ctx; 444 DM cdm = dm; 445 PetscFE fe; 446 char prefix[PETSC_MAX_PATH_LEN]; 447 DMPolytopeType ct; 448 PetscBool simplex; 449 PetscInt dim, cStart; 450 PetscErrorCode ierr; 451 452 PetscFunctionBegin; 453 /* Create finite element */ 454 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 455 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr); 456 ierr = DMPlexGetCellType(dm, cStart, &ct);CHKERRQ(ierr); 457 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE; 458 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);CHKERRQ(ierr); 459 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, dim, simplex, name ? prefix : NULL, -1, &fe);CHKERRQ(ierr); 460 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 461 /* Set discretization and boundary conditions for each mesh */ 462 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 463 ierr = DMCreateDS(dm);CHKERRQ(ierr); 464 ierr = (*setup)(dm, user);CHKERRQ(ierr); 465 while (cdm) { 466 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr); 467 if (user->useNearNullspace) {ierr = DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace);CHKERRQ(ierr);} 468 /* TODO: Check whether the boundary of coarse meshes is marked */ 469 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 470 } 471 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 472 PetscFunctionReturn(0); 473 } 474 475 int main(int argc, char **argv) 476 { 477 DM dm; /* Problem specification */ 478 SNES snes; /* Nonlinear solver */ 479 Vec u; /* Solutions */ 480 AppCtx user; /* User-defined work context */ 481 PetscErrorCode ierr; 482 483 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 484 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 485 /* Primal system */ 486 ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr); 487 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 488 ierr = SNESSetDM(snes, dm);CHKERRQ(ierr); 489 ierr = SetupFE(dm, "displacement", SetupPrimalProblem, &user);CHKERRQ(ierr); 490 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 491 ierr = VecSet(u, 0.0);CHKERRQ(ierr); 492 ierr = PetscObjectSetName((PetscObject) u, "displacement");CHKERRQ(ierr); 493 ierr = DMPlexSetSNESLocalFEM(dm, &user, &user, &user);CHKERRQ(ierr); 494 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 495 ierr = DMSNESCheckFromOptions(snes, u);CHKERRQ(ierr); 496 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 497 ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr); 498 ierr = VecViewFromOptions(u, NULL, "-displacement_view");CHKERRQ(ierr); 499 /* Cleanup */ 500 ierr = VecDestroy(&u);CHKERRQ(ierr); 501 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 502 ierr = DMDestroy(&dm);CHKERRQ(ierr); 503 ierr = PetscFinalize(); 504 return ierr; 505 } 506 507 /*TEST 508 509 testset: 510 args: -dm_plex_box_faces 1,1,1 511 512 test: 513 suffix: 2d_p1_quad_vlap 514 requires: triangle 515 args: -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate 516 test: 517 suffix: 2d_p2_quad_vlap 518 requires: triangle 519 args: -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 520 test: 521 suffix: 2d_p3_quad_vlap 522 requires: triangle 523 args: -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 524 test: 525 suffix: 2d_q1_quad_vlap 526 args: -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate 527 test: 528 suffix: 2d_q2_quad_vlap 529 args: -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 530 test: 531 suffix: 2d_q3_quad_vlap 532 requires: !single 533 args: -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 534 test: 535 suffix: 2d_p1_quad_elas 536 requires: triangle 537 args: -sol_type elas_quad -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate 538 test: 539 suffix: 2d_p2_quad_elas 540 requires: triangle 541 args: -sol_type elas_quad -displacement_petscspace_degree 2 -dmsnes_check .0001 542 test: 543 suffix: 2d_p3_quad_elas 544 requires: triangle 545 args: -sol_type elas_quad -displacement_petscspace_degree 3 -dmsnes_check .0001 546 test: 547 suffix: 2d_q1_quad_elas 548 args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 549 test: 550 suffix: 2d_q1_quad_elas_shear 551 args: -sol_type elas_quad -dm_plex_simplex 0 -shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 552 test: 553 suffix: 2d_q2_quad_elas 554 args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dmsnes_check .0001 555 test: 556 suffix: 2d_q2_quad_elas_shear 557 args: -sol_type elas_quad -dm_plex_simplex 0 -shear -displacement_petscspace_degree 2 -dmsnes_check 558 test: 559 suffix: 2d_q3_quad_elas 560 args: -sol_type elas_quad -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dmsnes_check .0001 561 test: 562 suffix: 2d_q3_quad_elas_shear 563 requires: !single 564 args: -sol_type elas_quad -dm_plex_simplex 0 -shear -displacement_petscspace_degree 3 -dmsnes_check 565 566 test: 567 suffix: 3d_p1_quad_vlap 568 requires: ctetgen 569 args: -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 570 test: 571 suffix: 3d_p2_quad_vlap 572 requires: ctetgen 573 args: -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001 574 test: 575 suffix: 3d_p3_quad_vlap 576 requires: ctetgen 577 args: -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001 578 test: 579 suffix: 3d_q1_quad_vlap 580 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 581 test: 582 suffix: 3d_q2_quad_vlap 583 args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001 584 test: 585 suffix: 3d_q3_quad_vlap 586 args: -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001 587 test: 588 suffix: 3d_p1_quad_elas 589 requires: ctetgen 590 args: -sol_type elas_quad -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 591 test: 592 suffix: 3d_p2_quad_elas 593 requires: ctetgen 594 args: -sol_type elas_quad -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001 595 test: 596 suffix: 3d_p3_quad_elas 597 requires: ctetgen 598 args: -sol_type elas_quad -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001 599 test: 600 suffix: 3d_q1_quad_elas 601 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 602 test: 603 suffix: 3d_q2_quad_elas 604 args: -sol_type elas_quad -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001 605 test: 606 suffix: 3d_q3_quad_elas 607 requires: !single 608 args: -sol_type elas_quad -dm_plex_dim 3 -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001 609 610 test: 611 suffix: 2d_p1_trig_vlap 612 requires: triangle 613 args: -sol_type vlap_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 614 test: 615 suffix: 2d_p2_trig_vlap 616 requires: triangle 617 args: -sol_type vlap_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 618 test: 619 suffix: 2d_p3_trig_vlap 620 requires: triangle 621 args: -sol_type vlap_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 622 test: 623 suffix: 2d_q1_trig_vlap 624 args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 625 test: 626 suffix: 2d_q2_trig_vlap 627 args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 628 test: 629 suffix: 2d_q3_trig_vlap 630 args: -sol_type vlap_trig -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 631 test: 632 suffix: 2d_p1_trig_elas 633 requires: triangle 634 args: -sol_type elas_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 635 test: 636 suffix: 2d_p2_trig_elas 637 requires: triangle 638 args: -sol_type elas_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 639 test: 640 suffix: 2d_p3_trig_elas 641 requires: triangle 642 args: -sol_type elas_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 643 test: 644 suffix: 2d_q1_trig_elas 645 args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 646 test: 647 suffix: 2d_q1_trig_elas_shear 648 args: -sol_type elas_trig -dm_plex_simplex 0 -shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 649 test: 650 suffix: 2d_q2_trig_elas 651 args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 652 test: 653 suffix: 2d_q2_trig_elas_shear 654 args: -sol_type elas_trig -dm_plex_simplex 0 -shear -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 655 test: 656 suffix: 2d_q3_trig_elas 657 args: -sol_type elas_trig -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 658 test: 659 suffix: 2d_q3_trig_elas_shear 660 args: -sol_type elas_trig -dm_plex_simplex 0 -shear -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 661 662 test: 663 suffix: 3d_p1_trig_vlap 664 requires: ctetgen 665 args: -sol_type vlap_trig -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 666 test: 667 suffix: 3d_p2_trig_vlap 668 requires: ctetgen 669 args: -sol_type vlap_trig -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 670 test: 671 suffix: 3d_p3_trig_vlap 672 requires: ctetgen 673 args: -sol_type vlap_trig -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 674 test: 675 suffix: 3d_q1_trig_vlap 676 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 677 test: 678 suffix: 3d_q2_trig_vlap 679 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 680 test: 681 suffix: 3d_q3_trig_vlap 682 requires: !__float128 683 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 684 test: 685 suffix: 3d_p1_trig_elas 686 requires: ctetgen 687 args: -sol_type elas_trig -dm_plex_dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 688 test: 689 suffix: 3d_p2_trig_elas 690 requires: ctetgen 691 args: -sol_type elas_trig -dm_plex_dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 692 test: 693 suffix: 3d_p3_trig_elas 694 requires: ctetgen 695 args: -sol_type elas_trig -dm_plex_dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 696 test: 697 suffix: 3d_q1_trig_elas 698 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 699 test: 700 suffix: 3d_q2_trig_elas 701 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 702 test: 703 suffix: 3d_q3_trig_elas 704 requires: !__float128 705 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 706 707 test: 708 suffix: 2d_p1_axial_elas 709 requires: triangle 710 args: -sol_type elas_axial_disp -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 2 -dmsnes_check .0001 -pc_type lu 711 test: 712 suffix: 2d_p2_axial_elas 713 requires: triangle 714 args: -sol_type elas_axial_disp -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu 715 test: 716 suffix: 2d_p3_axial_elas 717 requires: triangle 718 args: -sol_type elas_axial_disp -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu 719 test: 720 suffix: 2d_q1_axial_elas 721 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 722 test: 723 suffix: 2d_q2_axial_elas 724 args: -sol_type elas_axial_disp -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu 725 test: 726 suffix: 2d_q3_axial_elas 727 args: -sol_type elas_axial_disp -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu 728 729 test: 730 suffix: 2d_p1_uniform_elas 731 requires: triangle 732 args: -sol_type elas_uniform_strain -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 733 test: 734 suffix: 2d_p2_uniform_elas 735 requires: triangle 736 args: -sol_type elas_uniform_strain -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 737 test: 738 suffix: 2d_p3_uniform_elas 739 requires: triangle 740 args: -sol_type elas_uniform_strain -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 741 test: 742 suffix: 2d_q1_uniform_elas 743 args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 744 test: 745 suffix: 2d_q2_uniform_elas 746 requires: !single 747 args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 748 test: 749 suffix: 2d_q3_uniform_elas 750 requires: !single 751 args: -sol_type elas_uniform_strain -dm_plex_simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 752 753 TEST*/ 754