From: Jean-Paul Pelteret Date: Sun, 1 Aug 2021 12:54:13 +0000 (+0200) Subject: Update Cook membrane problem X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F90%2Fhead;p=code-gallery.git Update Cook membrane problem - Further corrections for running in 3d - Correct evaluation point such that it corresponds to the reference text (Miehe 1994) - Set the default case to be the 3d one, which implements the plane strain condition in the reference text - Updates for newer versions of deal.II --- diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/CMakeLists.txt b/Quasi_static_Finite_strain_Compressible_Elasticity/CMakeLists.txt index 6f40451..3e33e48 100644 --- a/Quasi_static_Finite_strain_Compressible_Elasticity/CMakeLists.txt +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/CMakeLists.txt @@ -20,7 +20,7 @@ SET(CLEAN_UP_FILES CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) -FIND_PACKAGE(deal.II 9.1 QUIET +FIND_PACKAGE(deal.II 9.2 QUIET HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} ) IF(NOT ${deal.II_FOUND}) diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/README.md b/Quasi_static_Finite_strain_Compressible_Elasticity/README.md index 5a5e5cc..c786ef8 100644 --- a/Quasi_static_Finite_strain_Compressible_Elasticity/README.md +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/README.md @@ -31,10 +31,6 @@ formulation, the two-dimensional case corresponds to neither plane-strain nor plane-stress conditions. -## Requirements -* Version 8.5 or greater of `deal.II` - - ## Compiling and running Similar to the example programs, run ``` @@ -103,6 +99,8 @@ Below we briefly document the tip displacement as predicted for different discretisation levels and ansatz for the displacement field. A direct and, by visual inspection, favourable comparison of the following results can be made with those found in Miehe (1994). +Specifically, figure 5 of Miehe (1994) plots the vertical displacement at the +midway point on the traction surface for the compressible 3d case. Since the material is compressible, shear-locking is not exhibited by the beam for low-order elements. @@ -110,7 +108,6 @@ beam for low-order elements. Elements per edge | Q1 | Q2 :---------------: | :-------------: | :-------------: -1 | 24 | 81 2 | 54 | 225 4 | 150 | 729 8 | 486 | 2601 @@ -122,10 +119,9 @@ Elements per edge | Q1 | Q2 Elements per edge | Q1 | Q2 :---------------: | :-------------: | :-------------: -1 | 5.15 | 12.19 -2 | 8.72 | 13.83 -4 | 12.02 | 14.22 -8 | 13.61 | 14.30 -16 | 14.13 | 14.32 -32 | 14.28 | 14.33 -64 | 14.32 | 14.33 +2 | 8.638 | 14.30 +4 | 12.07 | 14.65 +8 | 13.86 | 14.71 +16 | 14.49 | 14.73 +32 | 14.67 | 14.74 +64 | 14.72 | 14.74 diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc b/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc index f95f839..2eb7d8d 100644 --- a/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc @@ -1101,7 +1101,7 @@ namespace Cook_Membrane cell->face(face)->set_boundary_id(1); // -X faces else if (std::abs(cell->face(face)->center()[0] - 48.0) < tol_boundary) cell->face(face)->set_boundary_id(11); // +X faces - else if (std::abs(std::abs(cell->face(face)->center()[2]) - 0.5) < tol_boundary) + else if (dim == 3 && std::abs(std::abs(cell->face(face)->center()[2]) - 0.5) < tol_boundary) cell->face(face)->set_boundary_id(2); // +Z and -Z faces } @@ -1133,8 +1133,7 @@ namespace Cook_Membrane dof_handler_ref.distribute_dofs(fe); DoFRenumbering::Cuthill_McKee(dof_handler_ref); DoFRenumbering::component_wise(dof_handler_ref, block_component); - DoFTools::count_dofs_per_block(dof_handler_ref, dofs_per_block, - block_component); + dofs_per_block = DoFTools::count_dofs_per_fe_block(dof_handler_ref, block_component); std::cout << "Triangulation:" << "\n\t Number of active cells: " << triangulation.n_active_cells() @@ -1369,11 +1368,11 @@ namespace Cook_Membrane std::cout << "_"; std::cout << std::endl; - Point soln_pt{}； - soln_pt[0] = 48.0*parameters.scale; - soln_pt[1] = 60.0*parameters.scale; - if (dim == 3) - soln_pt[2] = 0.5*parameters.scale; + // The measurement point, as stated in the reference paper, is at the midway + // point of the surface on which the traction is applied. + const Point soln_pt = (dim == 3 ? + Point(48.0*parameters.scale, 52.0*parameters.scale, 0.5*parameters.scale) : + Point(48.0*parameters.scale, 52.0*parameters.scale)); double vertical_tip_displacement = 0.0; double vertical_tip_displacement_check = 0.0; @@ -2060,7 +2059,6 @@ namespace Cook_Membrane // and we can simply skip the rebuilding step if we do not clear it. if (it_nr > 1) return; - constraints.clear(); const bool apply_dirichlet_bc = (it_nr == 0); // The boundary conditions for the indentation problem are as follows: On @@ -2080,32 +2078,45 @@ namespace Cook_Membrane // select. To this end we first set up such extractor objects and later // use it when generating the relevant component masks: - // Fixed left hand side of the beam + if (apply_dirichlet_bc) { - const int boundary_id = 1; + constraints.clear(); - if (apply_dirichlet_bc == true) + // Fixed left hand side of the beam + { + const int boundary_id = 1; VectorTools::interpolate_boundary_values(dof_handler_ref, - boundary_id, - ZeroFunction(n_components), - constraints, - fe.component_mask(u_fe)); - } + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(u_fe)); + } - // Zero Z-displacement through thickness direction - // This corresponds to a plane strain condition being imposed on the beam - if (dim == 3) + // Zero Z-displacement through thickness direction + // This corresponds to a plane strain condition being imposed on the beam + if (dim == 3) { const int boundary_id = 2; const FEValuesExtractors::Scalar z_displacement(2); - - if (apply_dirichlet_bc == true) - VectorTools::interpolate_boundary_values(dof_handler_ref, - boundary_id, - ZeroFunction(n_components), - constraints, - fe.component_mask(z_displacement)); + VectorTools::interpolate_boundary_values(dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(z_displacement)); + } + } + else + { + if (constraints.has_inhomogeneities()) + { + AffineConstraints homogeneous_constraints(constraints); + for (unsigned int dof = 0; dof != dof_handler_ref.n_dofs(); ++dof) + if (homogeneous_constraints.is_inhomogeneously_constrained(dof)) + homogeneous_constraints.set_inhomogeneity(dof, 0.0); + constraints.clear(); + constraints.copy_from(homogeneous_constraints); } + } constraints.close(); } @@ -2235,7 +2246,7 @@ int main (int argc, char *argv[]) using namespace dealii; using namespace Cook_Membrane; - const unsigned int dim = 2; + const unsigned int dim = 3; try { deallog.depth_console(0);