Merge pull request #52 from mathLab/andrea_new

Andrea new

Author: luca-heltai

CMakeLists.txt

@@ -38,7 +38,7 @@ ENDIF()
 
 CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
 
-FIND_PACKAGE(deal.II 8.3 REQUIRED
+FIND_PACKAGE(deal.II 9.5 REQUIRED
   HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
 )
 
@@ -62,10 +62,6 @@ set(CMAKE_MACOSX_RPATH 1)
 endif()
 
 
-FIND_PACKAGE(deal2lkit 1.0 REQUIRED
-  HINTS ${D2K_DIR} $ENV{D2K_DIR} $ENV{DEAL2LKIT_DIR}
-  )
-D2K_INITIALIZE_CACHED_VARIABLES()
 
 # We one library and one target for each type of deal.II library
 # we found. If you compiled deal.II with both Release and Debug
@@ -80,10 +76,8 @@ FOREACH(_build_type ${_d2_build_types})
     SET(_p "${${_build_type}_postfix}")
     # Only build this type, if deal.II was compiled with it.
     IF(CMAKE_BUILD_TYPE MATCHES "${_build_type}" AND
-       DEAL_II_BUILD_TYPE MATCHES "${_build_type}"  AND
-       D2K_BUILD_TYPE MATCHES "${_build_type}")
+       DEAL_II_BUILD_TYPE MATCHES "${_build_type}")
 	MESSAGE("-- Found ${_build_type} version of deal.II.")
-	MESSAGE("-- Found ${_build_type} version of deal2lkit.")
 
 	SET(_lib "${TARGET}-lib${_p}")
 	MESSAGE("-- Configuring library ${_lib}")
@@ -98,24 +92,16 @@ FOREACH(_build_type ${_d2_build_types})
 		ADD_EXECUTABLE(${_exe} ${_main})
 		target_compile_definitions(${_exe} PUBLIC DEAL_II_DIMENSION=${_dim})
 		TARGET_LINK_LIBRARIES(${_exe} ${_lib})
-		D2K_SETUP_TARGET(${_exe} ${_BUILD_TYPE})
+		DEAL_II_SETUP_TARGET(${_exe} ${_BUILD_TYPE})
     	ENDFOREACH()
 
-	# Compile an executable only if required by the user
-	#IF(NOT "${_main}" STREQUAL "")
-	#    SET(_exe "${TARGET}${${_build_type}_postfix}")
-	#    MESSAGE("-- Configuring executable ${_exe}")
-	#    ADD_EXECUTABLE(${_exe} ${_main})
-	#    TARGET_LINK_LIBRARIES(${_exe} ${_lib})
-	#    D2K_SETUP_TARGET(${_exe} ${_BUILD_TYPE})
-	#ENDIF()
 
-	D2K_SETUP_TARGET(${_lib} ${_BUILD_TYPE})
+
+	DEAL_II_SETUP_TARGET(${_lib} ${_BUILD_TYPE})
   SET(TEST_LIBRARIES_${_BUILD_TYPE} ${_lib})
   ENDIF()
 ENDFOREACH()
 
-INCLUDE_DIRECTORIES(${D2K_DIR}/include)
 INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}/include)
 
 if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/tests AND BEM_ENABLE_TESTING)
@@ -138,6 +124,8 @@ if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in)
 endif()
 
 
+message(STATUS "DEAL_II_INCLUDE_DIRECTORIES=${deal.II_INCLUDE_DIRECTORIES}")
+
 ADD_CUSTOM_TARGET(indent
   WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
   COMMAND ./scripts/indent

Doxyfile.in

@@ -32,7 +32,7 @@ DOXYFILE_ENCODING      = UTF-8
 # title of most generated pages and in a few other places.
 # The default value is: My Project.
 
-PROJECT_NAME           = "pi-BEM: a parallel BEM solver for deal.II"
+PROJECT_NAME           = "Heat Shell: a heat transmission model for shells in deal.II"
 
 # The PROJECT_NUMBER tag can be used to enter a project or revision number. This
 # could be handy for archiving the generated documentation or if some version

grids/coarse_straight_cube.inp

@@ -0,0 +1,32 @@
+24 6 0 0 0
+1	0 0 0
+2	0 1 0
+3	1 1 0
+4	1 0 0 
+5	0 0 1
+6	0 1 1
+7	1 1 1
+8	1 0 1
+9	0 0 0
+10	1 0 0
+11	1 0 1
+12	0 0 1
+13	0 1 0
+14	1 1 0
+15	1 1 1
+16	0 1 1
+17	0 0 0
+18	0 0 1
+19	0 1 1
+20	0 1 0
+21	1 0 0
+22	1 0 1
+23	1 1 1
+24	1 1 0
+
+1 1 quad 1 2 3 4
+2 0 quad 5 8 7 6
+3 1 quad 9 10 11 12 
+4 1 quad 13 16 15 14
+5 1 quad 17 18 19 20
+6 1 quad 21 24 23 22

include/ass_leg_function.h

@@ -20,7 +20,6 @@
 #include <deal.II/lac/precondition.h>
 #include <deal.II/lac/solver_control.h>
 #include <deal.II/lac/solver_gmres.h>
-#include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/vector.h>
 // #include <deal.II/grid/tria_boundary_lib.h>
 

include/bem_fma.h

@@ -24,7 +24,7 @@
 #ifndef bem_fma_h
 #define bem_fma_h
 
-#include <deal.II/base/std_cxx11/bind.h>
+#include <deal.II/base/parameter_acceptor.h>
 #include <deal.II/base/types.h>
 #include <deal.II/base/work_stream.h>
 
@@ -33,14 +33,13 @@
 #include <deal.II/lac/trilinos_sparse_matrix.h>
 #include <deal.II/lac/trilinos_vector.h>
 
-#include <deal2lkit/parameter_acceptor.h>
-#include <deal2lkit/utilities.h>
 #include <mpi.h>
 
 #include <cmath>
 #include <fstream>
 #include <iostream>
 #include <map>
+#include <memory>
 #include <set>
 #include <string>
 
@@ -57,7 +56,6 @@ namespace Operator
 }
 
 using namespace dealii;
-using namespace deal2lkit;
 
 /**
 * A class for the handling of the Fast Multiple Method coupled with the Bundary
@@ -70,7 +68,7 @@ range interactions
 
 */
 template <int dim> //, Type V>
-class BEMFMA : public deal2lkit::ParameterAcceptor
+class BEMFMA : public ParameterAcceptor
 {
 public:
   /// Function to be used in the tests, it can access everything inside bemfma
@@ -100,10 +98,10 @@ class BEMFMA : public deal2lkit::ParameterAcceptor
   /// Mapping in the BEMFMA class. It also sets up some useful vector for,
   /// mixed boundary conditions and double nodes handling.
   void
-  init_fma(const DoFHandler<dim - 1, dim> &                      input_dh,
+  init_fma(const DoFHandler<dim - 1, dim>                       &input_dh,
            const std::vector<std::set<types::global_dof_index>> &db_in,
-           const TrilinosWrappers::MPI::Vector &                 input_sn,
-           const Mapping<dim - 1, dim> &                         input_mapping =
+           const TrilinosWrappers::MPI::Vector                  &input_sn,
+           const Mapping<dim - 1, dim>                          &input_mapping =
              StaticMappingQ1<dim - 1, dim>::mapping,
            unsigned int quad_order      = 4,
            unsigned int sing_quad_order = 5);
@@ -175,8 +173,8 @@ class BEMFMA : public deal2lkit::ParameterAcceptor
   multipole_matr_vect_products(
     const TrilinosWrappers::MPI::Vector &phi_values,
     const TrilinosWrappers::MPI::Vector &dphi_dn_values,
-    TrilinosWrappers::MPI::Vector &      matrVectProdN,
-    TrilinosWrappers::MPI::Vector &      matrVectProdD) const;
+    TrilinosWrappers::MPI::Vector       &matrVectProdN,
+    TrilinosWrappers::MPI::Vector       &matrVectProdD) const;
 
 
   // void compute_m2l_flags();
@@ -211,7 +209,7 @@ class BEMFMA : public deal2lkit::ParameterAcceptor
   /// architectures.
   TrilinosWrappers::PreconditionILU &
   FMA_preconditioner(const TrilinosWrappers::MPI::Vector &alpha,
-                     AffineConstraints<double> &          c);
+                     AffineConstraints<double>           &c);
 
 protected:
   /// Three pointers to the problem parameters to be set equal to
@@ -449,8 +447,8 @@ class BEMFMA : public deal2lkit::ParameterAcceptor
 
 
   /// TODO parsed quadrature?
-  shared_ptr<Quadrature<dim - 1>>    quadrature;
-  SmartPointer<const Vector<double>> dirichlet_nodes;
+  std::shared_ptr<Quadrature<dim - 1>> quadrature;
+  SmartPointer<const Vector<double>>   dirichlet_nodes;
   /// This should be erased by the usage of the constraint matrix.
   const std::vector<std::set<types::global_dof_index>> *double_nodes_set;
 

include/bem_problem.h

@@ -23,6 +23,7 @@
 
 #include <deal.II/base/conditional_ostream.h>
 #include <deal.II/base/convergence_table.h>
+#include <deal.II/base/parameter_acceptor.h>
 #include <deal.II/base/parsed_function.h>
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/quadrature_selector.h>
@@ -59,6 +60,7 @@
 #include <deal.II/dofs/dof_renumbering.h>
 #include <deal.II/dofs/dof_tools.h>
 
+#include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_system.h>
 #include <deal.II/fe/fe_values.h>
@@ -69,6 +71,7 @@
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_in.h>
 #include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
 #include <deal.II/grid/manifold_lib.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_accessor.h>
@@ -82,10 +85,7 @@
 
 // And here are a few C++ standard header
 // files that we will need:
-#include <deal2lkit/parameter_acceptor.h>
-#include <deal2lkit/parsed_finite_element.h>
-#include <deal2lkit/parsed_grid_refinement.h>
-#include <deal2lkit/utilities.h>
+
 #include <mpi.h>
 
 #include <cmath>
@@ -104,7 +104,6 @@
 #include "../include/octree_block.h"
 
 using namespace dealii;
-using namespace deal2lkit;
 
 // using namespace TrilinosWrappers;
 // using namespace TrilinosWrappers::MPI;
@@ -118,7 +117,7 @@ using namespace deal2lkit;
  *   - it eventually interacts with the FMM accelerator.
  */
 template <int dim>
-class BEMProblem : public deal2lkit::ParameterAcceptor
+class BEMProblem : public ParameterAcceptor
 {
 public:
   typedef typename DoFHandler<dim - 1, dim>::active_cell_iterator cell_it;
@@ -127,8 +126,8 @@ class BEMProblem : public deal2lkit::ParameterAcceptor
              const MPI_Comm            comm = MPI_COMM_WORLD);
 
   void
-  solve(TrilinosWrappers::MPI::Vector &      phi,
-        TrilinosWrappers::MPI::Vector &      dphi_dn,
+  solve(TrilinosWrappers::MPI::Vector       &phi,
+        TrilinosWrappers::MPI::Vector       &dphi_dn,
         const TrilinosWrappers::MPI::Vector &tmp_rhs);
 
   /// This function takes care of the proper initialization of all the elements
@@ -152,8 +151,8 @@ class BEMProblem : public deal2lkit::ParameterAcceptor
   /// have kept this function serial. We stress that it needs to be called only
   /// once.
   void
-  compute_constraints(IndexSet &                           c_cpu_set,
-                      AffineConstraints<double> &          constraints,
+  compute_constraints(IndexSet                            &c_cpu_set,
+                      AffineConstraints<double>           &constraints,
                       const TrilinosWrappers::MPI::Vector &tmp_rhs);
 
   //  private:
@@ -170,6 +169,10 @@ class BEMProblem : public deal2lkit::ParameterAcceptor
   virtual void
   parse_parameters(ParameterHandler &prm);
 
+  /// This function computes the free coefficients appearing  in the
+  /// hypersingular BIE.
+  void
+  compute_hypersingular_free_coeffs();
 
   /// This function computes the fraction of solid angles seen by our domain. We
   /// use the Double Layer Operator (through the Neumann matrix) to determine
@@ -195,15 +198,15 @@ class BEMProblem : public deal2lkit::ParameterAcceptor
   /// vector src. The result is stored
   /// in the vector dst.
   void
-  vmult(TrilinosWrappers::MPI::Vector &      dst,
+  vmult(TrilinosWrappers::MPI::Vector       &dst,
         const TrilinosWrappers::MPI::Vector &src) const;
 
   /// The second method computes the
   /// right hand side vector of the
   /// system.
 
   void
-  compute_rhs(TrilinosWrappers::MPI::Vector &      dst,
+  compute_rhs(TrilinosWrappers::MPI::Vector       &dst,
               const TrilinosWrappers::MPI::Vector &src) const;
 
   /// The third method computes the
@@ -220,8 +223,8 @@ class BEMProblem : public deal2lkit::ParameterAcceptor
   /// Depending on the resolution stategy we go whether for the direct or fma
   /// strategy.
   void
-  solve_system(TrilinosWrappers::MPI::Vector &      phi,
-               TrilinosWrappers::MPI::Vector &      dphi_dn,
+  solve_system(TrilinosWrappers::MPI::Vector       &phi,
+               TrilinosWrappers::MPI::Vector       &dphi_dn,
                const TrilinosWrappers::MPI::Vector &tmp_rhs);
 
 
@@ -247,6 +250,14 @@ class BEMProblem : public deal2lkit::ParameterAcceptor
   compute_gradients(const TrilinosWrappers::MPI::Vector &phi,
                     const TrilinosWrappers::MPI::Vector &dphi_dn);
 
+
+  /// We compute the potential gradients also in an alternative way.
+  /// Here we make use of the hypersingular integrals computed with the
+  /// SingularKernelIntegral class
+  void
+  compute_gradients_hypersingular(const TrilinosWrappers::MPI::Vector &phi,
+                                  const TrilinosWrappers::MPI::Vector &dphi_dn);
+
   /// We have parallelised the computation of the L2 projection of the normal
   /// vector. We need a solution vector that has also ghost cells. for this
   /// reason we made use of a ghosted IndexSet that we have computed in the
@@ -286,9 +297,8 @@ class BEMProblem : public deal2lkit::ParameterAcceptor
   ConditionalOStream        pcout;
   ComputationalDomain<dim> &comp_dom;
 
-
-  ParsedFiniteElement<dim - 1, dim>            parsed_fe;
-  ParsedFiniteElement<dim - 1, dim>            parsed_gradient_fe;
+  std::string                                  scalar_fe_type, vector_fe_type;
+  unsigned int                                 scalar_fe_order, vector_fe_order;
   std::unique_ptr<FiniteElement<dim - 1, dim>> fe;
   std::unique_ptr<FiniteElement<dim - 1, dim>> gradient_fe;
   DoFHandler<dim - 1, dim>                     dh;
@@ -297,16 +307,16 @@ class BEMProblem : public deal2lkit::ParameterAcceptor
   // FE_Q<dim-1,dim>                   fe;
   // FESystem<dim-1,dim>      gradient_fe;
 
-  ParsedGridRefinement pgr;
+  double refinement_threshold, coarsening_threshold;
 
   /// An Eulerian Mapping is created to deal
   /// with the free surface and boat mesh
   /// deformation
 
-  Vector<double>                    map_vector;
-  shared_ptr<Mapping<dim - 1, dim>> mapping;
-  unsigned int                      mapping_degree;
-  Vector<double>                    map_points;
+  Vector<double>                         map_vector;
+  std::shared_ptr<Mapping<dim - 1, dim>> mapping;
+  unsigned int                           mapping_degree;
+  Vector<double>                         map_points;
 
 
   /// these are the std::vectors of std::sets
@@ -338,8 +348,18 @@ class BEMProblem : public deal2lkit::ParameterAcceptor
 
   TrilinosWrappers::MPI::Vector system_rhs;
 
+  /// solution and alpha vectors
   TrilinosWrappers::MPI::Vector sol;
   TrilinosWrappers::MPI::Vector alpha;
+  /// an alternatively computed alpha vector (obtained with geometric
+  /// computations)
+  TrilinosWrappers::MPI::Vector hyp_alpha;
+  /// a set of distributed vectors which contain all the entries of the
+  /// C_ij tensor appearing in the hypersingular BIE
+  std::vector<TrilinosWrappers::MPI::Vector> C_ij;
+  /// a set of distributed vectors which contain all the entries of the
+  /// b_i vector appearing in the hypersingular BIE
+  std::vector<TrilinosWrappers::MPI::Vector> b_i;
 
   mutable TrilinosWrappers::MPI::Vector serv_phi;
   mutable TrilinosWrappers::MPI::Vector serv_dphi_dn;