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 PetscInt id; 383 PetscInt dim; 384 PetscErrorCode ierr; 385 386 PetscFunctionBeginUser; 387 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 388 ierr = PetscDSGetSpatialDimension(prob, &dim);CHKERRQ(ierr); 389 switch (user->solType) { 390 case SOL_VLAP_QUADRATIC: 391 ierr = PetscDSSetResidual(prob, 0, f0_vlap_quadratic_u, f1_vlap_u);CHKERRQ(ierr); 392 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_vlap_uu);CHKERRQ(ierr); 393 switch (dim) { 394 case 2: exact = quadratic_2d_u;break; 395 case 3: exact = quadratic_3d_u;break; 396 default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim); 397 } 398 break; 399 case SOL_ELAS_QUADRATIC: 400 ierr = PetscDSSetResidual(prob, 0, f0_elas_quadratic_u, f1_elas_u);CHKERRQ(ierr); 401 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);CHKERRQ(ierr); 402 switch (dim) { 403 case 2: exact = quadratic_2d_u;break; 404 case 3: exact = quadratic_3d_u;break; 405 default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim); 406 } 407 break; 408 case SOL_VLAP_TRIG: 409 ierr = PetscDSSetResidual(prob, 0, f0_vlap_trig_u, f1_vlap_u);CHKERRQ(ierr); 410 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_vlap_uu);CHKERRQ(ierr); 411 switch (dim) { 412 case 2: exact = trig_2d_u;break; 413 case 3: exact = trig_3d_u;break; 414 default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim); 415 } 416 break; 417 case SOL_ELAS_TRIG: 418 ierr = PetscDSSetResidual(prob, 0, f0_elas_trig_u, f1_elas_u);CHKERRQ(ierr); 419 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);CHKERRQ(ierr); 420 switch (dim) { 421 case 2: exact = trig_2d_u;break; 422 case 3: exact = trig_3d_u;break; 423 default: SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid dimension: %D", dim); 424 } 425 break; 426 case SOL_ELAS_AXIAL_DISP: 427 ierr = PetscDSSetResidual(prob, 0, NULL, f1_elas_u);CHKERRQ(ierr); 428 ierr = PetscDSSetBdResidual(prob, 0, f0_elas_axial_disp_bd_u, NULL);CHKERRQ(ierr); 429 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);CHKERRQ(ierr); 430 exact = axial_disp_u; 431 break; 432 case SOL_ELAS_UNIFORM_STRAIN: 433 ierr = PetscDSSetResidual(prob, 0, NULL, f1_elas_u);CHKERRQ(ierr); 434 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_elas_uu);CHKERRQ(ierr); 435 exact = uniform_strain_u; 436 break; 437 default: SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(user->solType, NUM_SOLUTION_TYPES)], user->solType); 438 } 439 ierr = PetscDSSetExactSolution(prob, 0, exact, user);CHKERRQ(ierr); 440 if (user->solType == SOL_ELAS_AXIAL_DISP) { 441 PetscInt cmp; 442 443 id = dim == 3 ? 5 : 2; 444 ierr = DMAddBoundary(dm, DM_BC_NATURAL, "right", "marker", 0, 0, NULL, (void (*)(void)) zero, NULL, 1, &id, user);CHKERRQ(ierr); 445 id = dim == 3 ? 6 : 4; 446 cmp = 0; 447 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "left", "marker", 0, 1, &cmp, (void (*)(void)) zero, NULL, 1, &id, user);CHKERRQ(ierr); 448 cmp = dim == 3 ? 2 : 1; 449 id = dim == 3 ? 1 : 1; 450 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom", "marker", 0, 1, &cmp, (void (*)(void)) zero, NULL, 1, &id, user);CHKERRQ(ierr); 451 if (dim == 3) { 452 cmp = 1; 453 id = 3; 454 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "front", "marker", 0, 1, &cmp, (void (*)(void)) zero, NULL, 1, &id, user);CHKERRQ(ierr); 455 } 456 } else { 457 id = 1; 458 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) exact, NULL, 1, &id, user);CHKERRQ(ierr); 459 } 460 PetscFunctionReturn(0); 461 } 462 463 static PetscErrorCode CreateElasticityNullSpace(DM dm, PetscInt origField, PetscInt field, MatNullSpace *nullspace) 464 { 465 PetscErrorCode ierr; 466 467 PetscFunctionBegin; 468 ierr = DMPlexCreateRigidBody(dm, origField, nullspace);CHKERRQ(ierr); 469 PetscFunctionReturn(0); 470 } 471 472 PetscErrorCode SetupFE(DM dm, PetscInt Nc, PetscBool simplex, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), void *ctx) 473 { 474 AppCtx *user = (AppCtx *) ctx; 475 DM cdm = dm; 476 PetscFE fe; 477 char prefix[PETSC_MAX_PATH_LEN]; 478 PetscInt dim; 479 PetscErrorCode ierr; 480 481 PetscFunctionBegin; 482 /* Create finite element */ 483 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 484 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);CHKERRQ(ierr); 485 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, Nc, simplex, name ? prefix : NULL, -1, &fe);CHKERRQ(ierr); 486 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 487 /* Set discretization and boundary conditions for each mesh */ 488 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 489 ierr = DMCreateDS(dm);CHKERRQ(ierr); 490 ierr = (*setup)(dm, user);CHKERRQ(ierr); 491 while (cdm) { 492 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr); 493 if (user->useNearNullspace) {ierr = DMSetNearNullSpaceConstructor(cdm, 0, CreateElasticityNullSpace);CHKERRQ(ierr);} 494 /* TODO: Check whether the boundary of coarse meshes is marked */ 495 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 496 } 497 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 498 PetscFunctionReturn(0); 499 } 500 501 int main(int argc, char **argv) 502 { 503 DM dm; /* Problem specification */ 504 SNES snes; /* Nonlinear solver */ 505 Vec u; /* Solutions */ 506 AppCtx user; /* User-defined work context */ 507 PetscErrorCode ierr; 508 509 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 510 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 511 /* Primal system */ 512 ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr); 513 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 514 ierr = SNESSetDM(snes, dm);CHKERRQ(ierr); 515 ierr = SetupFE(dm, user.dim, user.simplex, "displacement", SetupPrimalProblem, &user);CHKERRQ(ierr); 516 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 517 ierr = VecSet(u, 0.0);CHKERRQ(ierr); 518 ierr = PetscObjectSetName((PetscObject) u, "displacement");CHKERRQ(ierr); 519 ierr = DMPlexSetSNESLocalFEM(dm, &user, &user, &user);CHKERRQ(ierr); 520 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 521 ierr = DMSNESCheckFromOptions(snes, u);CHKERRQ(ierr); 522 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 523 ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr); 524 ierr = VecViewFromOptions(u, NULL, "-displacement_view");CHKERRQ(ierr); 525 /* Cleanup */ 526 ierr = VecDestroy(&u);CHKERRQ(ierr); 527 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 528 ierr = DMDestroy(&dm);CHKERRQ(ierr); 529 ierr = PetscFinalize(); 530 return ierr; 531 } 532 533 /*TEST 534 535 test: 536 suffix: 2d_p1_quad_vlap 537 requires: triangle 538 args: -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate 539 test: 540 suffix: 2d_p2_quad_vlap 541 requires: triangle 542 args: -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 543 test: 544 suffix: 2d_p3_quad_vlap 545 requires: triangle 546 args: -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 547 test: 548 suffix: 2d_q1_quad_vlap 549 args: -simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate 550 test: 551 suffix: 2d_q2_quad_vlap 552 args: -simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 553 test: 554 suffix: 2d_q3_quad_vlap 555 requires: !single 556 args: -simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 557 test: 558 suffix: 2d_p1_quad_elas 559 requires: triangle 560 args: -sol_type elas_quad -displacement_petscspace_degree 1 -dm_refine 2 -convest_num_refine 3 -snes_convergence_estimate 561 test: 562 suffix: 2d_p2_quad_elas 563 requires: triangle 564 args: -sol_type elas_quad -displacement_petscspace_degree 2 -dmsnes_check .0001 565 test: 566 suffix: 2d_p3_quad_elas 567 requires: triangle 568 args: -sol_type elas_quad -displacement_petscspace_degree 3 -dmsnes_check .0001 569 test: 570 suffix: 2d_q1_quad_elas 571 args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 572 test: 573 suffix: 2d_q1_quad_elas_shear 574 args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 575 test: 576 suffix: 2d_q2_quad_elas 577 args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 2 -dmsnes_check .0001 578 test: 579 suffix: 2d_q2_quad_elas_shear 580 args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 2 -dmsnes_check 581 test: 582 suffix: 2d_q3_quad_elas 583 args: -sol_type elas_quad -simplex 0 -displacement_petscspace_degree 3 -dmsnes_check .0001 584 test: 585 suffix: 2d_q3_quad_elas_shear 586 requires: !single 587 args: -sol_type elas_quad -simplex 0 -shear -displacement_petscspace_degree 3 -dmsnes_check 588 589 test: 590 suffix: 3d_p1_quad_vlap 591 requires: ctetgen 592 args: -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 593 test: 594 suffix: 3d_p2_quad_vlap 595 requires: ctetgen 596 args: -dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001 597 test: 598 suffix: 3d_p3_quad_vlap 599 requires: ctetgen 600 args: -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001 601 test: 602 suffix: 3d_q1_quad_vlap 603 args: -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 604 test: 605 suffix: 3d_q2_quad_vlap 606 args: -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001 607 test: 608 suffix: 3d_q3_quad_vlap 609 args: -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001 610 test: 611 suffix: 3d_p1_quad_elas 612 requires: ctetgen 613 args: -sol_type elas_quad -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 614 test: 615 suffix: 3d_p2_quad_elas 616 requires: ctetgen 617 args: -sol_type elas_quad -dim 3 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001 618 test: 619 suffix: 3d_p3_quad_elas 620 requires: ctetgen 621 args: -sol_type elas_quad -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001 622 test: 623 suffix: 3d_q1_quad_elas 624 args: -sol_type elas_quad -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 625 test: 626 suffix: 3d_q2_quad_elas 627 args: -sol_type elas_quad -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -dmsnes_check .0001 628 test: 629 suffix: 3d_q3_quad_elas 630 requires: !single 631 args: -sol_type elas_quad -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -dmsnes_check .0001 632 633 test: 634 suffix: 2d_p1_trig_vlap 635 requires: triangle 636 args: -sol_type vlap_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 637 test: 638 suffix: 2d_p2_trig_vlap 639 requires: triangle 640 args: -sol_type vlap_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 641 test: 642 suffix: 2d_p3_trig_vlap 643 requires: triangle 644 args: -sol_type vlap_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 645 test: 646 suffix: 2d_q1_trig_vlap 647 args: -sol_type vlap_trig -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 648 test: 649 suffix: 2d_q2_trig_vlap 650 args: -sol_type vlap_trig -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 651 test: 652 suffix: 2d_q3_trig_vlap 653 args: -sol_type vlap_trig -simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 654 test: 655 suffix: 2d_p1_trig_elas 656 requires: triangle 657 args: -sol_type elas_trig -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 658 test: 659 suffix: 2d_p2_trig_elas 660 requires: triangle 661 args: -sol_type elas_trig -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 662 test: 663 suffix: 2d_p3_trig_elas 664 requires: triangle 665 args: -sol_type elas_trig -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 666 test: 667 suffix: 2d_q1_trig_elas 668 args: -sol_type elas_trig -simplex 0 -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 669 test: 670 suffix: 2d_q1_trig_elas_shear 671 args: -sol_type elas_trig -simplex 0 -shear -displacement_petscspace_degree 1 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 672 test: 673 suffix: 2d_q2_trig_elas 674 args: -sol_type elas_trig -simplex 0 -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 675 test: 676 suffix: 2d_q2_trig_elas_shear 677 args: -sol_type elas_trig -simplex 0 -shear -displacement_petscspace_degree 2 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 678 test: 679 suffix: 2d_q3_trig_elas 680 args: -sol_type elas_trig -simplex 0 -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 681 test: 682 suffix: 2d_q3_trig_elas_shear 683 args: -sol_type elas_trig -simplex 0 -shear -displacement_petscspace_degree 3 -dm_refine 1 -convest_num_refine 3 -snes_convergence_estimate 684 685 test: 686 suffix: 3d_p1_trig_vlap 687 requires: ctetgen 688 args: -sol_type vlap_trig -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 689 test: 690 suffix: 3d_p2_trig_vlap 691 requires: ctetgen 692 args: -sol_type vlap_trig -dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 693 test: 694 suffix: 3d_p3_trig_vlap 695 requires: ctetgen 696 args: -sol_type vlap_trig -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 697 test: 698 suffix: 3d_q1_trig_vlap 699 args: -sol_type vlap_trig -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 700 test: 701 suffix: 3d_q2_trig_vlap 702 args: -sol_type vlap_trig -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 703 test: 704 suffix: 3d_q3_trig_vlap 705 requires: !__float128 706 args: -sol_type vlap_trig -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 707 test: 708 suffix: 3d_p1_trig_elas 709 requires: ctetgen 710 args: -sol_type elas_trig -dim 3 -dm_refine 1 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 711 test: 712 suffix: 3d_p2_trig_elas 713 requires: ctetgen 714 args: -sol_type elas_trig -dim 3 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 715 test: 716 suffix: 3d_p3_trig_elas 717 requires: ctetgen 718 args: -sol_type elas_trig -dim 3 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 719 test: 720 suffix: 3d_q1_trig_elas 721 args: -sol_type elas_trig -dim 3 -cells 2,2,2 -simplex 0 -displacement_petscspace_degree 1 -convest_num_refine 2 -snes_convergence_estimate 722 test: 723 suffix: 3d_q2_trig_elas 724 args: -sol_type elas_trig -dim 3 -simplex 0 -displacement_petscspace_degree 2 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 725 test: 726 suffix: 3d_q3_trig_elas 727 requires: !__float128 728 args: -sol_type elas_trig -dim 3 -simplex 0 -displacement_petscspace_degree 3 -dm_refine 0 -convest_num_refine 1 -snes_convergence_estimate 729 730 test: 731 suffix: 2d_p1_axial_elas 732 requires: triangle 733 args: -sol_type elas_axial_disp -displacement_petscspace_degree 1 -dm_plex_separate_marker -dm_refine 2 -dmsnes_check .0001 -pc_type lu 734 test: 735 suffix: 2d_p2_axial_elas 736 requires: triangle 737 args: -sol_type elas_axial_disp -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu 738 test: 739 suffix: 2d_p3_axial_elas 740 requires: triangle 741 args: -sol_type elas_axial_disp -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu 742 test: 743 suffix: 2d_q1_axial_elas 744 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 745 test: 746 suffix: 2d_q2_axial_elas 747 args: -sol_type elas_axial_disp -simplex 0 -displacement_petscspace_degree 2 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu 748 test: 749 suffix: 2d_q3_axial_elas 750 args: -sol_type elas_axial_disp -simplex 0 -displacement_petscspace_degree 3 -dm_plex_separate_marker -dmsnes_check .0001 -pc_type lu 751 752 test: 753 suffix: 2d_p1_uniform_elas 754 requires: triangle 755 args: -sol_type elas_uniform_strain -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 756 test: 757 suffix: 2d_p2_uniform_elas 758 requires: triangle 759 args: -sol_type elas_uniform_strain -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 760 test: 761 suffix: 2d_p3_uniform_elas 762 requires: triangle 763 args: -sol_type elas_uniform_strain -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 764 test: 765 suffix: 2d_q1_uniform_elas 766 args: -sol_type elas_uniform_strain -simplex 0 -displacement_petscspace_degree 1 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 767 test: 768 suffix: 2d_q2_uniform_elas 769 requires: !single 770 args: -sol_type elas_uniform_strain -simplex 0 -displacement_petscspace_degree 2 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 771 test: 772 suffix: 2d_q3_uniform_elas 773 requires: !single 774 args: -sol_type elas_uniform_strain -simplex 0 -displacement_petscspace_degree 3 -dm_refine 2 -dmsnes_check .0001 -pc_type lu 775 776 TEST*/ 777