[WIP] heart problem now with stokes

include/boundary_values.h

@@ -0,0 +1,58 @@
+#ifndef __boundary_values_h
+#define __boundary_values_h
+
+#include <deal.II/base/function.h>
+#include <deal.II/lac/vector.h>
+#include "heart_fe.h"
+
+using namespace dealii;
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues (int color, bool derivative=false) 
+    : 
+    Function<dim>(2*dim+1), 
+    color(color),
+    derivative(derivative)
+    {}
+  BoundaryValues (int color, 
+                  double timestep, 
+                  double dt, 
+                  bool side, 
+                  int degree, 
+                  bool derivative=false) 
+    : 
+    Function<dim>(2*dim+1),
+    color(color), 
+    timestep(timestep),
+    dt(dt),
+    heartstep(timestep/heartinterval),
+    derivative(derivative),
+    heart(side, degree)
+  {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const;
+  virtual void vector_value (const Point<dim> &p,
+                             Vector<double>   &value) const;
+private:
+  int color;
+  double timestep;
+  double dt;
+  double heartinterval = 0.005;
+  int heartstep;// = timestep / heartinterval;
+  bool derivative;
+  Heart<2,3> heart; 
+  void transform_to_polar_coord(const Point<3> &p, 
+                                double rot, 
+                                double &angle, 
+                                double &height) const;
+  void swap_coord(Point<3> &p) const;
+  void get_heartdelta (const Point<dim> &p, Vector<double> &value, int heartstep) const;
+  void get_values (const Point<dim> &p, Vector<double> &value) const;
+  void get_values_dt (const Point<dim> &p, Vector<double> &value) const;
+};
+
+#endif

include/elastic_problem.h

@@ -0,0 +1,60 @@
+#ifndef __elastic_problem_h
+#define __elastic_problem_h
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_q.h>
+
+using namespace dealii;
+
+
+template <int dim>
+class ElasticProblem
+{
+public:
+  ElasticProblem ();
+  ElasticProblem (Triangulation<dim> &tria, double timestep, double dt);
+  ~ElasticProblem ();
+  void run ();
+private:
+  void setup_system ();
+  void assemble_system ();
+  void solve ();
+  void refine_grid ();
+  Triangulation<dim>   triangulation;
+  DoFHandler<dim>      dof_handler;
+  FESystem<dim>        fe;
+  ConstraintMatrix     hanging_node_constraints;
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+  Vector<double>       solution;
+  Vector<double>       system_rhs;
+  double               timestep;
+  double               dt;
+};
+
+#endif

include/heart_fe.h

@@ -0,0 +1,33 @@
+#ifndef __heart_fe_h
+#define __heart_fe_h
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/dofs/dof_handler.h>
+
+using namespace dealii;
+
+template <int dim,int spacedim>
+class Heart
+{
+public:
+  Heart ();
+  Heart (bool, const int);
+  Point<spacedim> push_forward (const Point<dim>, const int) const;
+
+private:
+  void reinit_data ();
+  void setup_system ();
+  void run_side ();
+  void run_bottom ();
+
+  Triangulation<dim>		   		triangulation;
+  FESystem<dim> 			       	fe;
+  DoFHandler<dim>			      	dof_handler;
+  std::vector<Vector<double> >      solution;
+
+  bool side;
+};
+
+#endif

main/heart_stokes.cc

@@ -0,0 +1,194 @@
+#include "interfaces/stokes.h"
+#include "pidomus.h"
+
+#include "deal.II/base/numbers.h"
+
+#include "Teuchos_CommandLineProcessor.hpp"
+#include "Teuchos_GlobalMPISession.hpp"
+#include "Teuchos_oblackholestream.hpp"
+#include "Teuchos_StandardCatchMacros.hpp"
+#include "Teuchos_Version.hpp"
+
+#include "mpi.h"
+#include <iostream>
+#include <string>
+
+////////////////////////////////////////////////////////////////////////////////
+// Begin macros
+#define problem_Stokes(dim,spacedim,LAC) \
+  StokesInterface<dim,spacedim,LATrilinos> interface; \
+  piDoMUS<dim,spacedim,LAC> stokes ( \
+                                            "piDoMUS", \
+                                            interface); \
+  ParameterAcceptor::initialize(prm_file, pde_name+"_used.prm"); \
+  stokes.run ();
+// End macros
+////////////////////////////////////////////////////////////////////////////////
+
+void print_status(  std::string name,
+                    std::string prm_file,
+                    int dim,
+                    // int spacedim,
+                    int n_threads,
+                    const MPI_Comm &comm,
+                    bool check_prm);
+
+
+int main (int argc, char *argv[])
+{
+  using namespace dealii;
+  using namespace deal2lkit;
+
+  Teuchos::CommandLineProcessor My_CLP;
+  My_CLP.setDocString(
+    "                ______     ___  ____   _ __________ \n"
+    "                |  _  \\    |  \\/  | | | /  ______  \\ \n"
+    "  ______________| | | |____| .  . | |_| \\ `--.___| | \n"
+    " |__  __  ______| | | / _ \\| |\\/| | |_| |`--. \\____/ \n"
+    "   | | | |      | |/ | (_) | |  | | |_| /\\__/ / \n"
+    "   |_| |_|      |___/ \\___/\\_|  |_/\\___/\\____/ \n"
+  );
+
+  std::string pde_name="stokes";
+  //My_CLP.setOption("pde", &pde_name, "name of the PDE (stokes, NS for navier stokes, or ALE for ALE navier stokes)");
+
+  // int spacedim = 2;
+  // My_CLP.setOption("spacedim", &spacedim, "dimensione of the whole space");
+
+  int dim = 2;
+  My_CLP.setOption("dim", &dim, "dimension of the problem");
+
+  int n_threads = 0;
+  My_CLP.setOption("n_threads", &n_threads, "number of threads");
+
+  std::string prm_file=pde_name+".prm";
+  My_CLP.setOption("prm", &prm_file, "name of the parameter file");
+
+  bool check_prm = false;
+  My_CLP.setOption("check","no-check", &check_prm, "wait for a key press before starting the run");
+
+  // My_CLP.recogniseAllOptions(true);
+  My_CLP.throwExceptions(false);
+
+  Teuchos::CommandLineProcessor::EParseCommandLineReturn
+  parseReturn= My_CLP.parse( argc, argv );
+
+  if ( parseReturn == Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED )
+    {
+      return 0;
+    }
+  if ( parseReturn != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL   )
+    {
+      return 1; // Error!
+    }
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
+                                                      n_threads == 0 ? numbers::invalid_unsigned_int  : n_threads);
+
+  const MPI_Comm &comm  = MPI_COMM_WORLD;
+
+  Teuchos::oblackholestream blackhole;
+  std::ostream &out = ( Utilities::MPI::this_mpi_process(comm) == 0 ? std::cout : blackhole );
+
+
+  std::string string_pde_name="Stokes";
+
+
+  My_CLP.printHelpMessage(argv[0], out);
+
+  // deallog.depth_console (0);
+  try
+    {
+      print_status(   string_pde_name+" Equations",
+                      prm_file,
+                      dim,
+                      // spacedim,
+                      n_threads,
+                      comm,
+                      check_prm);
+
+      if (dim==2)
+        {
+          problem_Stokes(2,2,LATrilinos);
+        }
+      else
+        {
+          problem_Stokes(3,3,LATrilinos);
+        }
+
+
+      out << std::endl;
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}
+
+void print_status(  std::string name,
+                    std::string prm_file,
+                    int dim,
+                    // int spacedim,
+                    int n_threads,
+                    const MPI_Comm &comm,
+                    bool check_prm)
+{
+  int numprocs  = Utilities::MPI::n_mpi_processes(comm);
+  //  int myid      = Utilities::MPI::this_mpi_process(comm);
+
+  Teuchos::oblackholestream blackhole;
+  std::ostream &out = ( Utilities::MPI::this_mpi_process(comm) == 0 ? std::cout : blackhole );
+
+
+  out << std::endl
+      << "============================================================="
+      << std::endl
+      << "     Name:  " << name
+      << std::endl
+      << " Prm file:  " << prm_file
+      << std::endl
+      << "n threads:  " << n_threads
+      << std::endl
+      << "  process:  "  << getpid()
+      << std::endl
+      << " proc.tot:  "  << numprocs
+      << std::endl
+      // << " spacedim:  " << spacedim
+      // << std::endl
+      << "      dim:  " << dim
+      // << std::endl
+      // << "    codim:  " << spacedim-dim
+      << std::endl;
+  if (check_prm)
+    {
+      out<< "-------------------------------------------------------------"
+         << std::endl;
+      out << "Press [Enter] key to start...";
+      if (std::cin.get() == '\n') {};
+    }
+  out << "============================================================="
+      <<std::endl<<std::endl;
+
+}