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 6f40451d903fb8f9f9f623317462b8a12477cd73..3e33e485c48582d7933cf6191f62e8f237f1e83d 100644 (file)
--- 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 5a5e5cc1b6ec382ebd457d7ce818335c3183433f..c786ef879fb119664ea307200f13d57170922ba3 100644 (file)
--- 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 f95f8390235c3942ce9f9ca3eecebeed8147b681..2eb7d8d9c1565357af8e7d712cefc23c57799356 100644 (file)
--- 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<dim> 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<dim> soln_pt = (dim == 3 ? 
+                                Point<dim>(48.0*parameters.scale, 52.0*parameters.scale, 0.5*parameters.scale) : 
+                                Point<dim>(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<dim>(n_components),
-                                                 constraints,
-                                                 fe.component_mask(u_fe));
-    }
+                                                boundary_id,
+                                                ZeroFunction<dim>(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<dim>(n_components),
-                                                   constraints,
-                                                   fe.component_mask(z_displacement));
+        VectorTools::interpolate_boundary_values(dof_handler_ref,
+                                                boundary_id,
+                                                ZeroFunction<dim>(n_components),
+                                                constraints,
+                                                fe.component_mask(z_displacement));
+      }
+    }
+    else
+    {
+      if (constraints.has_inhomogeneities())
+      {
+        AffineConstraints<double> 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);