From 11ef1804fb9e3d41a5a50f44d482bdc651e00c95 Mon Sep 17 00:00:00 2001
From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Mon, 28 Dec 2015 09:22:11 +0100
Subject: [PATCH] Added parameter file; code cleanup; fixed small bug in
 post-processing

---
 .../cook_membrane.cc                          | 24 ++++---
 .../parameters.prm                            | 68 +++++++++++++++++++
 2 files changed, 82 insertions(+), 10 deletions(-)
 create mode 100644 Quasi_static_Finite_strain_Compressible_Elasticity/parameters.prm

diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc b/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc
index 6dffb8b..fbfc824 100644
--- a/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc
+++ b/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc
@@ -40,7 +40,6 @@
 #include <deal.II/grid/grid_in.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/grid/grid_out.h> // TEMP
 
 #include <deal.II/fe/fe_dgp_monomial.h>
 #include <deal.II/fe/fe_q.h>
@@ -48,6 +47,7 @@
 #include <deal.II/fe/fe_tools.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping_q_eulerian.h>
+#include <deal.II/fe/mapping_q.h>
 
 #include <deal.II/lac/block_sparse_matrix.h>
 #include <deal.II/lac/block_vector.h>
@@ -1747,9 +1747,18 @@ Point<dim> grid_y_transform (const Point<dim> &pt_in)
       for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
         if (cell->vertex(v).distance(soln_pt) < 1e-6)
       {
-        // Create some object to help us extract the solution value at 
-        // the desired point
-        const Quadrature<dim> soln_qrule (soln_pt);
+        // Extract y-component of solution at the given point
+        // This point is coindicent with a vertex, so we can
+        // extract it directly as we're using FE_Q finite elements
+        // that have support at the vertices
+        vertical_tip_displacement = solution_n(cell->vertex_dof_index(v,u_dof+1));
+        
+        // Sanity check using alternate method to extract the solution
+        // at the given point. To do this, we must create an FEValues instance
+        // to help us extract the solution value at the desired point
+        const MappingQ<dim> mapping (parameters.poly_degree);
+        const Point<dim> qp_unit = mapping.transform_real_to_unit_cell(cell,soln_pt);
+        const Quadrature<dim> soln_qrule (qp_unit);
         AssertThrow(soln_qrule.size() == 1, ExcInternalError());
         FEValues<dim> fe_values_soln (fe, soln_qrule, update_values);
         fe_values_soln.reinit(cell);
@@ -1758,12 +1767,7 @@ Point<dim> grid_y_transform (const Point<dim> &pt_in)
         std::vector< Tensor<1,dim> > soln_values (soln_qrule.size());
         fe_values_soln[u_fe].get_function_values(solution_n,
                                                  soln_values);
-        vertical_tip_displacement = soln_values[0][u_dof+1];
-        
-        // Sanity Check
-        for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-          if (cell->vertex(v).distance(soln_pt) < 1e-6)
-            vertical_tip_displacement_check = solution_n(cell->vertex_dof_index(v,u_dof+1));
+        vertical_tip_displacement_check = soln_values[0][u_dof+1];
         
         break;
       }
diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/parameters.prm b/Quasi_static_Finite_strain_Compressible_Elasticity/parameters.prm
new file mode 100644
index 0000000..595f06a
--- /dev/null
+++ b/Quasi_static_Finite_strain_Compressible_Elasticity/parameters.prm
@@ -0,0 +1,68 @@
+# Listing of Parameters
+# ---------------------
+subsection Finite element system
+  # Displacement system polynomial order
+  set Polynomial degree = 1
+
+  # Gauss quadrature order
+  set Quadrature order  = 2
+end
+
+
+subsection Geometry
+  # Number of elements per long edge of the beam
+  set Elements per edge  = 32
+
+  # Global grid scaling factor
+  set Grid scale         = 1e-3
+end
+
+
+subsection Linear solver
+  # Linear solver iterations (multiples of the system matrix size)
+  set Max iteration multiplier  = 1
+
+  # Linear solver residual (scaled by residual norm)
+  set Residual                  = 1e-6
+
+  # Preconditioner type
+  set Preconditioner type        = ssor
+
+  # Preconditioner relaxation value
+  set Preconditioner relaxation  = 0.65
+
+  # Type of solver used to solve the linear system
+  set Solver type               = Direct
+end
+
+
+subsection Material properties
+  # Poisson's ratio
+  set Poisson's ratio = 0.3
+
+  # Shear modulus
+  set Shear modulus   = 0.4225e6
+end
+
+
+subsection Nonlinear solver
+  # Number of Newton-Raphson iterations allowed
+  set Max iterations Newton-Raphson = 10
+
+  # Displacement error tolerance
+  set Tolerance displacement        = 1.0e-6
+
+  # Force residual tolerance
+  set Tolerance force               = 1.0e-9
+end
+
+
+subsection Time
+  # End time
+  set End time       = 1
+
+  # Time step size
+  set Time step size = 0.1
+end
+
+
-- 
2.39.5