]> https://gitweb.dealii.org/ - code-gallery.git/commitdiff
Documentation and code improvements 72/head
authorJean-Paul Pelteret <jppelteret@gmail.com>
Sat, 23 May 2020 19:18:24 +0000 (21:18 +0200)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Tue, 9 Jun 2020 20:45:49 +0000 (22:45 +0200)
15 files changed:
Nonlinear_PoroViscoelasticity/CMakeLists.txt
Nonlinear_PoroViscoelasticity/doc/images-results/Franceschini-consolidation.png [new file with mode: 0644]
Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-reac_p-poro-COMP.png [new file with mode: 0755]
Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-COMP.png [new file with mode: 0755]
Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-visco-COMP.png [new file with mode: 0755]
Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-COMP.png [new file with mode: 0755]
Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-visco-COMP.png [new file with mode: 0755]
Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-experiments-COMP.png [new file with mode: 0644]
Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-COMP.png [new file with mode: 0755]
Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-visco-COMP.png [new file with mode: 0755]
Nonlinear_PoroViscoelasticity/doc/images-results/consolidation-mesh.png [new file with mode: 0644]
Nonlinear_PoroViscoelasticity/doc/images-results/cyclic-loading-mesh.png [new file with mode: 0644]
Nonlinear_PoroViscoelasticity/doc/images-results/poro-visco-compression.gif [new file with mode: 0644]
Nonlinear_PoroViscoelasticity/nonlinear-poro-viscoelasticity.cc
Nonlinear_PoroViscoelasticity/readme.md

index a0a416dcf47a2367fd6017249df4a2070e1ff1e1..6c48e74263343b9d2a839c2390df52fd52a0b8e1 100644 (file)
@@ -26,16 +26,12 @@ ENDIF()
 
 
 # Are all dependencies fullfilled?
-IF(NOT DEAL_II_WITH_CXX11     OR
-   NOT DEAL_II_WITH_MPI       OR
-   NOT DEAL_II_WITH_METIS     OR
+IF(NOT DEAL_II_WITH_MPI       OR
    NOT DEAL_II_WITH_TRILINOS  OR
    NOT DEAL_II_TRILINOS_WITH_SACADO)
   MESSAGE(FATAL_ERROR "
 Error! The deal.II library found at ${DEAL_II_PATH} was not configured with
-    DEAL_II_WITH_CXX11 = ON
     DEAL_II_WITH_MPI = ON
-    DEAL_II_WITH_METIS = ON
     DEAL_II_WITH_TRILINOS = ON
     DEAL_II_TRILINOS_WITH_SACADO = ON
 One or all of these are OFF in your installation but are required for this code gallery example."
diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/Franceschini-consolidation.png b/Nonlinear_PoroViscoelasticity/doc/images-results/Franceschini-consolidation.png
new file mode 100644 (file)
index 0000000..18ce85e
Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/Franceschini-consolidation.png differ
diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-reac_p-poro-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-reac_p-poro-COMP.png
new file mode 100755 (executable)
index 0000000..0efb2ef
Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-reac_p-poro-COMP.png differ
diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-COMP.png
new file mode 100755 (executable)
index 0000000..61f472e
Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-COMP.png differ
diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-visco-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-visco-COMP.png
new file mode 100755 (executable)
index 0000000..4576b68
Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-visco-COMP.png differ
diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-COMP.png
new file mode 100755 (executable)
index 0000000..c446723
Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-COMP.png differ
diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-visco-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-visco-COMP.png
new file mode 100755 (executable)
index 0000000..aed8542
Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-visco-COMP.png differ
diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-experiments-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-experiments-COMP.png
new file mode 100644 (file)
index 0000000..75cceee
Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-experiments-COMP.png differ
diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-COMP.png
new file mode 100755 (executable)
index 0000000..1868863
Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-COMP.png differ
diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-visco-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-visco-COMP.png
new file mode 100755 (executable)
index 0000000..901553c
Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-visco-COMP.png differ
diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/consolidation-mesh.png b/Nonlinear_PoroViscoelasticity/doc/images-results/consolidation-mesh.png
new file mode 100644 (file)
index 0000000..6902ba6
Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/consolidation-mesh.png differ
diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/cyclic-loading-mesh.png b/Nonlinear_PoroViscoelasticity/doc/images-results/cyclic-loading-mesh.png
new file mode 100644 (file)
index 0000000..afae378
Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/cyclic-loading-mesh.png differ
diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/poro-visco-compression.gif b/Nonlinear_PoroViscoelasticity/doc/images-results/poro-visco-compression.gif
new file mode 100644 (file)
index 0000000..5b0c2cb
Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/poro-visco-compression.gif differ
index 74a431f9f329ee5c626387e2759f7a4f89e56969..c96f2b1aabfa9b5a14f1bf8a83b0dedecfff8906 100644 (file)
 #include <deal.II/base/work_stream.h>\r
 #include <deal.II/base/mpi.h>\r
 #include <deal.II/base/quadrature_point_data.h>\r
-#include <deal.II/base/std_cxx11/shared_ptr.h>\r
 \r
 #include <deal.II/differentiation/ad.h>\r
 \r
+#include <deal.II/distributed/shared_tria.h>\r
+\r
 #include <deal.II/dofs/dof_renumbering.h>\r
 #include <deal.II/dofs/dof_tools.h>\r
 #include <deal.II/dofs/dof_accessor.h>\r
@@ -47,7 +48,6 @@
 #include <deal.II/grid/grid_in.h>\r
 #include <deal.II/grid/grid_out.h>\r
 #include <deal.II/grid/manifold_lib.h>\r
-#include <deal.II/grid/tria.h>\r
 #include <deal.II/grid/tria_accessor.h>\r
 #include <deal.II/grid/tria_boundary_lib.h>\r
 #include <deal.II/grid/tria_iterator.h>\r
@@ -448,12 +448,12 @@ namespace NonLinearPoroViscoElasticity
           gravity_value = prm.get_double("gravity value");\r
 \r
           if ( (fluid_type == "Markert") && ((init_intrinsic_perm == 0.0) || (viscosity_FR == 0.0)) )\r
-              throw std::runtime_error ("Markert seepage velocity formulation requires the definition of "\r
-                                          "'initial intrinsic permeability' and 'fluid viscosity' greater than 0.0.");\r
+              AssertThrow(false, ExcMessage("Markert seepage velocity formulation requires the definition of "\r
+                                            "'initial intrinsic permeability' and 'fluid viscosity' greater than 0.0."));\r
 \r
           if ( (fluid_type == "Ehlers") && ((init_darcy_coef == 0.0) || (weight_FR == 0.0)) )\r
-              throw std::runtime_error ("Ehler seepage velocity formulation requires the definition of "\r
-                                          "'initial Darcy coefficient' and 'fluid weight' greater than 0.0.");\r
+              AssertThrow(false, ExcMessage("Ehler seepage velocity formulation requires the definition of "\r
+                                            "'initial Darcy coefficient' and 'fluid weight' greater than 0.0."));\r
 \r
           const std::string eigen_solver_type = prm.get("eigen solver");\r
           if (eigen_solver_type == "QL Implicit Shifts")\r
@@ -1009,8 +1009,8 @@ namespace NonLinearPoroViscoElasticity
                   }\r
                   iteration += 1;\r
                   if (iteration > 15 )\r
-                      throw std::runtime_error ("No convergence in local Newton iteration for the "\r
-                                                "viscoelastic exponential time integration algorithm.");\r
+                      AssertThrow(false, ExcMessage("No convergence in local Newton iteration for the "\r
+                                                    "viscoelastic exponential time integration algorithm."));\r
               }\r
 \r
               NumberType aux_J_e_1 = 1.0;\r
@@ -1222,8 +1222,9 @@ namespace NonLinearPoroViscoElasticity
                  permeability_term = get_darcy_flow_current(F) / weight_FR;\r
 \r
              else\r
-                 throw std::runtime_error ("Material_Darcy_Fluid --> Only Markert "\r
-                          "and Ehlers formulations have been implemented.");\r
+                 AssertThrow(false, ExcMessage(\r
+                   "Material_Darcy_Fluid --> Only Markert "\r
+                   "and Ehlers formulations have been implemented."));\r
 \r
              return ( -1.0 * permeability_term * det_F\r
                       * (grad_p_fluid - get_body_force_FR_current()) );\r
@@ -1250,7 +1251,9 @@ namespace NonLinearPoroViscoElasticity
                  seepage_velocity = get_seepage_velocity_current(F,grad_p_fluid);\r
              }\r
              else\r
-                 throw std::runtime_error ("Material_Darcy_Fluid --> Only Markert and Ehlers formulations have been implemented.");\r
+                 AssertThrow(false, ExcMessage(\r
+                   "Material_Darcy_Fluid --> Only Markert and Ehlers "\r
+                   "formulations have been implemented."));\r
 \r
              dissipation_term = ( invert(permeability_term) * seepage_velocity ) * seepage_velocity;\r
              dissipation_term *= 1.0/(det_F*det_F);\r
@@ -1423,8 +1426,8 @@ namespace NonLinearPoroViscoElasticity
                 return body_force;\r
             }\r
         private:\r
-            std_cxx11::shared_ptr< Material_Hyperelastic<dim, NumberType> > solid_material;\r
-            std_cxx11::shared_ptr< Material_Darcy_Fluid<dim, NumberType> > fluid_material;\r
+            std::shared_ptr< Material_Hyperelastic<dim, NumberType> > solid_material;\r
+            std::shared_ptr< Material_Darcy_Fluid<dim, NumberType> > fluid_material;\r
     };\r
 \r
 // @sect3{Nonlinear poro-viscoelastic solid}\r
@@ -1439,7 +1442,7 @@ namespace NonLinearPoroViscoElasticity
             void run();\r
 \r
           protected:\r
-            typedef Sacado::Fad::DFad<double> ADNumberType;\r
+            using ADNumberType = Sacado::Fad::DFad<double>;\r
 \r
             std::ofstream outfile;\r
             std::ofstream pointfile;\r
@@ -1525,11 +1528,11 @@ namespace NonLinearPoroViscoElasticity
             const unsigned int               this_mpi_process;\r
             mutable ConditionalOStream       pcout;\r
 \r
-            //A collection of the parameters used to describe the problem setup\r
+            // A collection of the parameters used to describe the problem setup\r
             const Parameters::AllParameters &parameters;\r
 \r
-            //Declare an instance of dealii Triangulation class (mesh)\r
-            Triangulation<dim>  triangulation;\r
+            // Declare an instance of dealii Triangulation class (mesh)\r
+            parallel::shared::Triangulation<dim>  triangulation;\r
 \r
             // Keep track of the current time and the time spent evaluating certain functions\r
             Time          time;\r
@@ -1674,7 +1677,7 @@ namespace NonLinearPoroViscoElasticity
         this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)),\r
         pcout(std::cout, this_mpi_process == 0),\r
         parameters(parameters),\r
-        triangulation(Triangulation<dim>::maximum_smoothing),\r
+        triangulation(mpi_communicator,Triangulation<dim>::maximum_smoothing),\r
         time(parameters.end_time, parameters.delta_t),\r
         timerconsole( mpi_communicator,\r
                       pcout,\r
@@ -1993,10 +1996,6 @@ namespace NonLinearPoroViscoElasticity
         timerconsole.enter_subsection("Setup system");\r
         timerfile.enter_subsection("Setup system");\r
 \r
-        // Partition triangulation\r
-        GridTools::partition_triangulation (n_mpi_processes, triangulation);\r
-\r
-\r
         //Determine number of components per block\r
         std::vector<unsigned int> block_component(n_components, u_block);\r
         block_component[p_fluid_component] = p_fluid_block;\r
@@ -2174,13 +2173,15 @@ namespace NonLinearPoroViscoElasticity
 \r
         //Setup the initial quadrature point data using the info stored in parameters\r
         FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>\r
-        cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+        cell (IteratorFilters::LocallyOwnedCell(),\r
               dof_handler_ref.begin_active()),\r
-        endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+        endc (IteratorFilters::LocallyOwnedCell(),\r
               dof_handler_ref.end());\r
         for (; cell!=endc; ++cell)\r
           {\r
-            Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());\r
+            Assert(cell->is_locally_owned(), ExcInternalError());\r
+            Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());\r
+\r
             const std::vector<std::shared_ptr<PointHistory<dim, ADNumberType> > >\r
                 lqph = quadrature_point_history.get_data(cell);\r
             Assert(lqph.size() == n_q_points, ExcInternalError());\r
@@ -2481,13 +2482,15 @@ namespace NonLinearPoroViscoElasticity
                                                    solution_total);\r
 \r
         FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>\r
-        cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+        cell (IteratorFilters::LocallyOwnedCell(),\r
               dof_handler_ref.begin_active()),\r
-        endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+        endc (IteratorFilters::LocallyOwnedCell(),\r
               dof_handler_ref.end());\r
         for (; cell != endc; ++cell)\r
           {\r
-            Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());\r
+            Assert(cell->is_locally_owned(), ExcInternalError());\r
+            Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());\r
+\r
             assemble_system_one_cell(cell, scratch_data, per_task_data);\r
             copy_local_to_global_system(per_task_data);\r
           }\r
@@ -2526,7 +2529,7 @@ namespace NonLinearPoroViscoElasticity
               ScratchData_ASM<ADNumberType>                        &scratch,\r
               PerTaskData_ASM                                      &data) const\r
     {\r
-        Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());\r
+        Assert(cell->is_locally_owned(), ExcInternalError());\r
 \r
         data.reset();\r
         scratch.reset();\r
@@ -2747,13 +2750,15 @@ namespace NonLinearPoroViscoElasticity
     void Solid<dim>::update_end_timestep()\r
     {\r
           FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>\r
-          cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+          cell (IteratorFilters::LocallyOwnedCell(),\r
                 dof_handler_ref.begin_active()),\r
-          endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+          endc (IteratorFilters::LocallyOwnedCell(),\r
                 dof_handler_ref.end());\r
           for (; cell!=endc; ++cell)\r
           {\r
-            Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());\r
+            Assert(cell->is_locally_owned(), ExcInternalError());\r
+            Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());\r
+\r
             const std::vector<std::shared_ptr<PointHistory<dim, ADNumberType> > >\r
                 lqph = quadrature_point_history.get_data(cell);\r
             Assert(lqph.size() == n_q_points, ExcInternalError());\r
@@ -2799,9 +2804,7 @@ namespace NonLinearPoroViscoElasticity
     class FilteredDataOut : public DataOut<dim, DH>\r
     {\r
         public:\r
-          FilteredDataOut (const unsigned int subdomain_id)\r
-            :\r
-            subdomain_id (subdomain_id)\r
+          FilteredDataOut ()\r
           {}\r
 \r
           virtual ~FilteredDataOut() {}\r
@@ -2812,7 +2815,7 @@ namespace NonLinearPoroViscoElasticity
             typename DataOut<dim, DH>::active_cell_iterator\r
             cell = this->dofs->begin_active();\r
             while ((cell != this->dofs->end()) &&\r
-                   (cell->subdomain_id() != subdomain_id))\r
+                   (!cell->is_locally_owned()))\r
               ++cell;\r
             return cell;\r
           }\r
@@ -2822,7 +2825,7 @@ namespace NonLinearPoroViscoElasticity
           {\r
             if (old_cell != this->dofs->end())\r
               {\r
-                const IteratorFilters::SubdomainEqualTo predicate(subdomain_id);\r
+                const IteratorFilters::LocallyOwnedCell predicate{};\r
                 return\r
                   ++(FilteredIterator<typename DataOut<dim, DH>::active_cell_iterator>\r
                      (predicate,old_cell));\r
@@ -2830,18 +2833,13 @@ namespace NonLinearPoroViscoElasticity
             else\r
               return old_cell;\r
           }\r
-\r
-        private:\r
-          const unsigned int subdomain_id;\r
     };\r
 \r
     template<int dim, class DH=DoFHandler<dim> >\r
      class FilteredDataOutFaces : public DataOutFaces<dim,DH>\r
      {\r
          public:\r
-           FilteredDataOutFaces (const unsigned int subdomain_id)\r
-             :\r
-             subdomain_id (subdomain_id)\r
+           FilteredDataOutFaces ()\r
            {}\r
 \r
            virtual ~FilteredDataOutFaces() {}\r
@@ -2851,7 +2849,7 @@ namespace NonLinearPoroViscoElasticity
            {\r
              typename DataOutFaces<dim,DH>::active_cell_iterator\r
              cell = this->dofs->begin_active();\r
-             while ((cell!=this->dofs->end()) && (cell->subdomain_id()!=subdomain_id))\r
+             while ((cell!=this->dofs->end()) && (!cell->is_locally_owned()))\r
                ++cell;\r
              return cell;\r
            }\r
@@ -2861,7 +2859,7 @@ namespace NonLinearPoroViscoElasticity
            {\r
              if (old_cell!=this->dofs->end())\r
              {\r
-                 const IteratorFilters::SubdomainEqualTo predicate(subdomain_id);\r
+                 const IteratorFilters::LocallyOwnedCell predicate{};\r
                  return\r
                    ++(FilteredIterator<typename DataOutFaces<dim,DH>::active_cell_iterator>\r
                       (predicate,old_cell));\r
@@ -2869,9 +2867,6 @@ namespace NonLinearPoroViscoElasticity
              else\r
                return old_cell;\r
            }\r
-\r
-         private:\r
-           const unsigned int subdomain_id;\r
      };\r
 \r
     //Class to compute gradient of the pressure\r
@@ -3028,20 +3023,22 @@ namespace NonLinearPoroViscoElasticity
 \r
         //Iterate through elements (cells) and Gauss Points\r
         FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>\r
-          cell(IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+          cell(IteratorFilters::LocallyOwnedCell(),\r
                dof_handler_ref.begin_active()),\r
-          endc(IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+          endc(IteratorFilters::LocallyOwnedCell(),\r
                dof_handler_ref.end()),\r
-          cell_v(IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+          cell_v(IteratorFilters::LocallyOwnedCell(),\r
                  vertex_handler_ref.begin_active()),\r
-          cell_v_vec(IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+          cell_v_vec(IteratorFilters::LocallyOwnedCell(),\r
                      vertex_vec_handler_ref.begin_active());\r
         //start cell loop\r
         for (; cell!=endc; ++cell, ++cell_v, ++cell_v_vec)\r
         {\r
-          if (cell->subdomain_id() != this_mpi_process) continue;\r
-                material_id(cell->active_cell_index())=\r
-                   static_cast<int>(cell->material_id());\r
+            Assert(cell->is_locally_owned(), ExcInternalError());\r
+            Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());\r
+\r
+            material_id(cell->active_cell_index())=\r
+              static_cast<int>(cell->material_id());\r
 \r
             fe_values_ref.reinit(cell);\r
 \r
@@ -3274,7 +3271,7 @@ namespace NonLinearPoroViscoElasticity
         }\r
 \r
         // Add the results to the solution to create the output file for Paraview\r
-        FilteredDataOut<dim> data_out(this_mpi_process);\r
+        FilteredDataOut<dim> data_out;\r
         std::vector<DataComponentInterpretation::DataComponentInterpretation>\r
           comp_type(dim,\r
                     DataComponentInterpretation::component_is_part_of_vector);\r
@@ -3532,15 +3529,15 @@ namespace NonLinearPoroViscoElasticity
 \r
         //Iterate through elements (cells) and Gauss Points\r
         FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>\r
-          cell(IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+          cell(IteratorFilters::LocallyOwnedCell(),\r
                dof_handler_ref.begin_active()),\r
-          endc(IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+          endc(IteratorFilters::LocallyOwnedCell(),\r
                dof_handler_ref.end());\r
         //start cell loop\r
         for (; cell!=endc; ++cell)\r
         {\r
-\r
-            if (cell->subdomain_id() != this_mpi_process) continue;\r
+            Assert(cell->is_locally_owned(), ExcInternalError());\r
+            Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());\r
 \r
             fe_values_ref.reinit(cell);\r
 \r
@@ -4092,7 +4089,7 @@ namespace NonLinearPoroViscoElasticity
               std::vector<double> displ_incr(dim, 0.0);\r
               (void)boundary_id;\r
               (void)direction;\r
-              throw std::runtime_error ("Displacement loading not implemented for Ehlers verification examples.");\r
+              AssertThrow(false, ExcMessage("Displacement loading not implemented for Ehlers verification examples."));\r
 \r
               return displ_incr;\r
           }\r
@@ -4320,7 +4317,7 @@ namespace NonLinearPoroViscoElasticity
               std::vector<double> displ_incr(dim, 0.0);\r
               (void)boundary_id;\r
               (void)direction;\r
-              throw std::runtime_error ("Displacement loading not implemented for Ehlers verification examples.");\r
+              AssertThrow(false, ExcMessage("Displacement loading not implemented for Ehlers verification examples."));\r
 \r
               return displ_incr;\r
           }\r
@@ -4480,7 +4477,7 @@ namespace NonLinearPoroViscoElasticity
               std::vector<double> displ_incr(dim, 0.0);\r
               (void)boundary_id;\r
               (void)direction;\r
-              throw std::runtime_error ("Displacement loading not implemented for Franceschini examples.");\r
+              AssertThrow(false, ExcMessage("Displacement loading not implemented for Franceschini examples."));\r
 \r
               return displ_incr;\r
           }\r
index 4981df818464e9de6c75ed967c3de7d566829a31..372c96d5d126ee27b431e707932cb32bf48f4cf4 100644 (file)
@@ -6,35 +6,34 @@ Overview
 
 We implemented a nonlinear poro-viscoelastic formulation with the aim of characterising brain tissue response to cyclic loading. Our model captures both experimentally observed fluid flow and conditioning aspects of brain tissue behavior in addition to its well-established nonlinear, preconditioning, hysteretic, and tension-compression asymmetric characteristics.
 
-The tissue is modelled as a biphasic material consisting in an immiscible aggregate of a nonlinear viscoelastic solid skeleton saturated with pore fluid. The governing equations are linearised using automatic differentiation and solved monolithically for the unknown solid displacements and fluid pore pressure values. 
+The tissue is modelled as a biphasic material consisting of an immiscible aggregate of a nonlinear viscoelastic solid skeleton saturated with pore fluid. The governing equations are linearised using automatic differentiation and solved monolithically for the unknown solid displacements and fluid pore pressure values. 
 
-A detailed description of the formulation, it's verification, and the results obtained can be found in:
-* E. Comellas, S. Budday, J.-P. Pelteret, G. A. Holzapfel and P. Steinamnn (2020), Modeling the porous and viscous responses of human brain tissue behavior, Computer Methods in Applied Mechanics and Engineering (accepted for publication).
+A detailed description of the formulation, its verification, and the results obtained can be found in:
+* E. Comellas, S. Budday, J.-P. Pelteret, G. A. Holzapfel and P. Steinmann (2020), Modeling the porous and viscous responses of human brain tissue behavior, Computer Methods in Applied Mechanics and Engineering, 113128. DOI: [10.1016/j.cma.2020.113128](https://doi.org/10.1016/j.cma.2020.113128);
 
 In this paper we show that nonlinear poroelasticity alone can reproduce consolidation experiments, yet it is insufficient to capture stress conditioning due to cyclic loading. We also discuss how the poroelastic response exhibits preconditioning and hysteresis in the fluid flow space, with porous and viscous effects being highly interrelated.
 
 
-Quick facts on the code
------------------------
+Quick facts about the code
+--------------------------
 
 * Biphasic material following the Theory of Porous Media
 * Nonlinear finite viscoelasticity built on Ogden hyperelasticity
 * Darcy-like fluid flow
-* Spatial discretisation with continuous Q2P1 Lagrangian finite elements
+* Spatial discretisation with continuous Q2-P1 Lagrangian finite elements
 * Temporal discretisation with a stable implicit one-step backward differentiation method
 * Newton-Raphson scheme to solve the nonlinear system of governing equations
-* Forward mode automatic differentiation with the number of derivative components chosen at run-time (Sacado library within Trilinos package) to linearise the equations
+* Forward mode automatic differentiation with the number of derivative components chosen at run-time (Sacado library within Trilinos package) to linearise the governing equations (and, implicitly, the constitutive laws)
 * Trilinos direct solver for the (non-symmetric) linear system of equations using a monolithic scheme
-* Parallelization through Threading Building Blocks and across nodes via MPI (Trilinos package)
-* Based on step-44 and code gallery contributions 'Quasi-Static Finite-Strain Compressible Elasticity' and 'Quasi-Static Finite-Strain Quasi-incompressible Visco-elasticity'
+* Parallelization through Threading Building Blocks and across nodes via MPI (using Trilinos linear algebra)
+* Based on step-44 and the code gallery contributions 'Quasi-Static Finite-Strain Compressible Elasticity' and 'Quasi-Static Finite-Strain Quasi-incompressible Visco-elasticity'
 * Only works in 3D
 
 
 Running the code
 ---------------- 
 ### Requirements
-* Version 9.1 or greater of deal.II
-* C++11, MPI, Metis and Trilinos with Sacado must be enabled
+* MPI and Trilinos (built with the Sacado library) must be enabled
 
 ### Compiling and running
 Similar to the example programs, run
@@ -70,10 +69,11 @@ python run-multi-calc.py
 The 'run-multi-calc.py' and 'runPoro.sh' files provided must both be in the main directory. This will automatically generate the required input files and run them in sequence.
 
 
+
 Reference for this work
 -----------------------
 
-If you use this program as a basis for your own work, please consider citing the paper referenced in the introduction. The initial version of this work was contributed to the deal.II project by E. Comellas and J.-P. Pelteret.
+If you use this program as a basis for your own work, please consider citing the paper referenced in the introduction. The initial version of this work was contributed to the deal.II project by E. Comellas and J-P. Pelteret.
 
 
 
@@ -87,3 +87,91 @@ Recommended literature
 * G.A. Holzapfel (2001), Nonlinear Solid Mechanics. A Continuum Approach for Engineering, John Wiley & Sons. ISBN: [978-0-471-82319-3](http://eu.wiley.com/WileyCDA/WileyTitle/productCd-0471823198.html);
 
 
+Results
+----------------------
+The results shown here are a selection of those presented and discussed in the paper referenced in the introduction.
+
+### Consolidation experiments
+We reproduce the uniaxial consolidation experiments of brain tissue from the seminal paper by Franceschini et al. (2006) using a reduction of our formulation to nonlinear poroelasticity without the viscous component. 
+
+The following geometry (320 cells and 9544 degrees of freedom), boundary conditions, and loading are used:
+
+![Consolidation geometry](./doc/images-results/consolidation-mesh.png)
+
+The material properties are:
+* First Lamé parameter: 334kPa
+* Ogden parameters:  μ<sub>1</sub>=1.044kPa,  μ<sub>2</sub>=1.183kPa, α<sub>1</sub>=4.309, and α<sub>2</sub>=7.736
+* Initial solid volume fraction: 0.80
+* Isotropic initial intrinsic permeability: 8.0e-11mm<sup>2</sup>
+* Deformation-dependency control parameter for the permeability: 40
+* Fluid viscosity: 0.89mPa·s
+
+We consider the effect of gravity in our formulation, with the simplifying assumption that the fluid and solid density are both 0.997mg/mm<sup>3</sup>.
+
+Using these parameters and simulation conditions, we obtain a reasonable fit to the experimental curve:
+
+![Consolidation results](./doc/images-results/Franceschini-consolidation.png)
+
+What is interesting about these results is that they show that nonlinear poroelasticity alone, i.e. with a purely hyperelastic solid component, can capture the consolidation behavior of brain tissue in response to an [oedometric test](https://en.wikipedia.org/wiki/Oedometer_test).
+
+However, as we will see in the next section, the viscous component is necessary to also capture stress conditioning due to cyclic loading.
+
+
+### Cyclic experiments under multiple loading modes
+Budday et al. (2017) performed an exhaustive set of cyclic loading experiments on brain tissue. They showed in this and subsequent publications that a finite viscoelastic Ogden model is able to predict essential features like nonlinearity, preconditioning, hysteresis, and tension-compression asymmetry. However, such a monophasic viscoelastic formulation can only implicitly capture the effects due to the extracellular fluid flow within the tissue. This was the motivation behind the formulation presented here: to provide a comprehensive  model for brain tissue response that will enable us to explore the porous and viscous effects in its mechanical response to a wide range of loading scenarios. 
+
+As an example, we show here the results for cyclic compressive loading. The following geometry (512 cells and 15468 degrees of freedom), boundary conditions, and loading are used:
+
+![Cyclic loading geometry](./doc/images-results/cyclic-loading-mesh.png)
+
+The material properties are:
+* First Lamé parameter: 24.5kPa
+* Ogden hyperelastic parameters: μ<sub>∞,1</sub>=-83.9Pa, and α<sub>∞,1</sub>=-11.92
+* Ogden viscous parameters: μ<sub>1</sub>=-2100Pa, and α<sub>1</sub>=-2.2
+* Deformation-independent viscosity of the solid: 14kPa·s
+* Initial solid volume fraction: 0.75
+* Isotropic initial intrinsic permeability: 1.0e-8mm<sup>2</sup>
+* Deformation-dependency control parameter for the permeability: 50
+* Fluid viscosity: 0.89mPa·s
+
+To simplify the problem, we neglect the effect of gravity for this example. Here is an animation of our results, visualised with Paraview:
+
+[![Cyclic loading video results](./doc/images-results/poro-visco-compression.gif)](https://youtu.be/1a4LbMXrnrs)
+
+To compare with the experimental results, our code computes the nominal stress on the loading surface and provides the values for each timestep in the "data-for-gnuplot.sol" output file. We can then plot compressive nominal stress versus the averaged vertical stretch:
+
+![Cyclic loading experimental results](./doc/images-results/brain_cube_stress-stretch-experiments-COMP.png)
+![Cyclic loading poroelastic results](./doc/images-results/brain_cube_stress-stretch-poro-COMP.png)
+![Cyclic loading poro-viscoelastic results](./doc/images-results/brain_cube_stress-stretch-poro-visco-COMP.png)
+
+We see that viscosity is required in the solid component to reproduce the preconditioning and hysteretic response seen in the experiments. 
+
+We were somewhat surprised by these results, because we were expecting to see preconditioning and hysteresis with poroelasticity alone. We discussed this in a conference talk:
+
+* E. Comellas, J.-P. Pelteret, S. Budday and P. Steinmann (2018). Unsolved issues in the numerical modelling of experimentally-observed porous effects in brain tissue behaviour, 6th European Conference on Computational Mechanics and 7th European Conference on Computational Fluid Dynamics (ECCM-ECFD 2018), Glasgow (UK), 11th–15th June 2018. DOI: [10.13140/RG.2.2.18553.29283](https://doi.org/10.13140/RG.2.2.18553.29283);
+
+So, we set out to explore why this is not the case. For that, we studied the porous and viscous dissipation, which are computed in the code and also provided in the "data-for-gnuplot.sol" output file. We determined that there is dissipation occurring in the poroelastic case, but we were barely seeing it in the stress-stretch plot. Even if we played around with the material parameters to increase the amount of dissipation, the slight hysteresis in our plot barely changed. *Where is all that energy going?* 
+
+We found the answer by looking closely into the thermodynamic basis of our constitutive equations. The viscous dissipation is a function of the viscous part of the stress tensor and the viscosity of the solid component. However, the porous dissipation depends on the porous material parameters and the seepage velocity, i.e. the relative velocity of the fluid with respect to the deforming solid component, which is proportional to the gradient of the pressure. Now, the plots we were studying are in the displacement-related space. We were looking for the porous dissipation in the wrong place! Our unknowns, the displacements and pressure, are slightly coupled. This is why we see a bit of porous dissipation in the stress-stretch plot for the poroelastic example. 
+
+Indeed, when we computed an equivalent to the nominal stress but in the pressure-related space, we found the "missing" dissipation. The "accumulated outgoing fluid volume" in the pressure-space is the equivalent to the nominal stress in the displacement-space. We plot it versus the fluid reaction pressure on the loading surface, which can be loosely seen as an equivalent to the averaged vertical stretch in the previous plots.
+
+![Cyclic loading poroelastic hysteresis](./doc/images-results/brain_cube_fluid-reac_p-poro-COMP.png)
+
+
+### Our takeaways
+In developing this formulation and running the examples provided in this code gallery, we have realised the importance of the loading and boundary conditions selected in the problems. Inhomogeneous displacements and pressure distributions require a fine enough discretization as well as a careful set-up of representative loading and boundary conditions. Without inhomogeneous stretch and pressure maps, the predicted viscous and porous behaviours of the material will not be realistic. 
+
+We have started to explore the relation between viscous and porous effects in response to different types of loading. Interestingly, the fluid accumulated outside the sample increases for each loading cycle in the poroelastic model, while it decreases for the poro-viscoelastic one. 
+
+![Cyclic loading poroelastic accum fluid](./doc/images-results/brain_cube_fluid-time-poro-COMP.png)
+![Cyclic loading poro-viscoelastic accum fluid](./doc/images-results/brain_cube_fluid-time-poro-visco-COMP.png)
+
+The evolution of pressure in the centre of the loading surface, where we observe a pressure concentration due to loading, has a noticeably different pattern for the two models. Note also the delay between the peak pressure values and peak displacement ones.
+
+![Cyclic loading poroelastic accum fluid](./doc/images-results/brain_cube_pres-displ-vs-time-poro-COMP.png)
+![Cyclic loading poro-viscoelastic accum fluid](./doc/images-results/brain_cube_pres-displ-vs-time-poro-visco-COMP.png)
+
+We are not sure yet how much of these differences is due to the incorporation of the viscous component, and how much can be explained by the choice of material parameters. (Note that the solid part of the poroelastic example has been adjusted with the Franceschini experiment, but the poro-viscoelastic one used the Budday experiments.) We will keep exploring with more numerical examples, but additional experiments are needed to adjust our material parameters in order to better distinguish the porous and viscous parts of brain tissue response to loading. 
+
+We hope the community finds this code useful, whether it is for brain mechanics or other applications, and look forward to learning and discussing about new insights gained about the interrelations of the viscous, porous, and elastic parts of the nonlinear poro-viscoelastic formulation.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.