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