#880 Poromechanics: fix the parallel issues.

Data/Lua/demo_input_poromechanics.lua

@@ -111,6 +111,24 @@ NumberingSubset15 = {
 } -- NumberingSubset15
 
 
+NumberingSubset16 = {
+    
+    -- Name of the numbering subset (not really used; at the moment I just need one input parameter to ground
+    -- the possible values to choose elsewhere).
+    -- Expected format: "VALUE"
+    name = 'solid_velocity_on_fluid',
+    
+    -- Whether a vector defined on this numbering subset might be used to compute a movemesh. If true, a
+    -- FEltSpace featuring this numbering subset will compute additional quantities to enable fast computation.
+    -- This should be false for most numbering subsets, and when it's true the sole unknown involved should be a
+    -- displacement.
+    -- Expected format: 'true' or 'false' (without the quote)
+    do_move_mesh = false
+    
+} -- NumberingSubset16
+
+
+
 
 NumberingSubset20 = {
     
@@ -539,15 +557,15 @@ FiniteElementSpace10 = {
 	-- List of all unknowns defined in the finite element space. Unknowns here must be defined in this file as 
 	-- an 'Unknown' block; expected name/identifier is the name given there. 
 	-- Expected format: {"VALUE1", "VALUE2", ...}
-    unknown_list = { 'fluid_velocity', 'porosity', 'solid_displacement' },
+    unknown_list = { 'fluid_velocity', 'porosity', 'solid_displacement', 'solid_velocity' },
 
 	-- List of the shape function to use for each unknown;
 	-- Expected format: {"VALUE1", "VALUE2", ...}
-	shape_function_list = { 'P1b', 'P1', 'P1' },
+    shape_function_list = { 'P1b', 'P1', 'P1', 'P1' },
 
 	-- List of the numbering subset to use for each unknown;
 	-- Expected format: {VALUE1, VALUE2, ...}
-	numbering_subset_list = { 10, 13, 14 }
+	numbering_subset_list = { 10, 13, 14, 16 }
 
 } -- FiniteElementSpace10
 
@@ -1242,7 +1260,7 @@ TransientSource1 = {
 porosity = 0.1
 
 
--- Fluid
+-- Fluid/porosity
 InitVertexMatchingInterpolator10 = {
     
     -- Finite element space for which the dofs index will be associated to each vertex.
@@ -1256,7 +1274,7 @@ InitVertexMatchingInterpolator10 = {
 } -- InitVertexMatchingInterpolator10
 
 
--- Solid
+-- Solid/porosity
 InitVertexMatchingInterpolator20 = {
     
     -- Finite element space for which the dofs index will be associated to each vertex.
@@ -1269,6 +1287,34 @@ InitVertexMatchingInterpolator20 = {
     
 } -- InitVertexMatchingInterpolator20
 
+
+-- Fluid/velocity
+InitVertexMatchingInterpolator11 = {
+    
+    -- Finite element space for which the dofs index will be associated to each vertex.
+    -- Expected format: VALUE
+    finite_element_space = 10,
+    
+    -- Numbering subsetfor which the dofs index will be associated to each vertex.
+    -- Expected format: "VALUE"
+    numbering_subset = 16
+    
+} -- InitVertexMatchingInterpolator10
+
+
+-- Solid/velocity
+InitVertexMatchingInterpolator21 = {
+    
+    -- Finite element space for which the dofs index will be associated to each vertex.
+    -- Expected format: VALUE
+    finite_element_space = 20,
+    
+    -- Numbering subsetfor which the dofs index will be associated to each vertex.
+    -- Expected format: "VALUE"
+    numbering_subset = 21
+    
+} -- InitVertexMatchingInterpolator20
+
 interpolation_file = "${HOME}/Codes/HappyHeart/Data/Interpolation/Poromechanics/Identity_2.hhdata"
 
 

Sources/ModelInstances/UnderDevelopment/Poromechanics/InputParameterList.hpp

@@ -107,6 +107,7 @@ namespace HappyHeart
             porosity = 13,
             solid_displacement_in_fluid_mesh = 14,
             fluid_mass_velocity_pressure = 15,
+            solid_velocity_on_fluid = 16,
             
             solid_displacement = 20,
             solid_velocity = 21,
@@ -117,8 +118,11 @@ namespace HappyHeart
         
         enum class InitVertexMatchingInterpolator
         {
-            fluid = 10,
-            solid = 20
+            fluid_porosity = 10,
+            solid_porosity = 20,
+            
+            fluid_velocity = 11,
+            solid_velocity = 21
         };
             
              
@@ -135,7 +139,8 @@ namespace HappyHeart
             InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::porosity_on_solid)>,
             InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_on_solid)>,
             InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure)>,
-        
+            InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::solid_velocity_on_fluid)>,
+            
             InputParameter::Unknown<EnumUnderlyingType(UnknownIndex::fluid_velocity)>,
             InputParameter::Unknown<EnumUnderlyingType(UnknownIndex::fluid_pressure)>,
             InputParameter::Unknown<EnumUnderlyingType(UnknownIndex::fluid_mass)>,
@@ -179,8 +184,10 @@ namespace HappyHeart
             
             InputParameter::InterpolationFile,
             
-            InputParameter::InitVertexMatchingInterpolator<EnumUnderlyingType(InitVertexMatchingInterpolator::solid)>,
-            InputParameter::InitVertexMatchingInterpolator<EnumUnderlyingType(InitVertexMatchingInterpolator::fluid)>
+            InputParameter::InitVertexMatchingInterpolator<EnumUnderlyingType(InitVertexMatchingInterpolator::solid_porosity)>,
+            InputParameter::InitVertexMatchingInterpolator<EnumUnderlyingType(InitVertexMatchingInterpolator::fluid_porosity)>,
+            InputParameter::InitVertexMatchingInterpolator<EnumUnderlyingType(InitVertexMatchingInterpolator::solid_velocity)>,
+            InputParameter::InitVertexMatchingInterpolator<EnumUnderlyingType(InitVertexMatchingInterpolator::fluid_velocity)>
 
         >;
         

Sources/ModelInstances/UnderDevelopment/Poromechanics/Model.hxx

@@ -20,7 +20,6 @@ namespace HappyHeart
     {
         
         
-        
         template<class SolidVariationalFormulationPolicyT>
         void Model<SolidVariationalFormulationPolicyT>::Forward()
         {
@@ -62,7 +61,6 @@ namespace HappyHeart
         }
         
         
-        
     } // namespace PoromechanicsNS
     
     

Sources/ModelInstances/UnderDevelopment/Poromechanics/ModelInitialize.hxx

@@ -48,7 +48,8 @@ namespace HappyHeart
             decltype(auto) solid_mesh = solid_god_of_dof.GetGeometricMeshRegion();
             
             solid_ = std::make_unique<Solid>(input_parameter_data, solid_mesh);
-
+            
+      
             // \todo #820 Order below will have to be adapted probably: porosity that is used in almost all
             // variational formulations and operators is defined in the very first one.
             // So whatever the organization, we must make sure PorosityVariationalFormulation comes first.
@@ -93,8 +94,8 @@ namespace HappyHeart
             
             {
                 solid_on_fluid_mesh_ =
-                    std::make_unique<Private::SolidOnFluidMesh>(fluid_god_of_dof,
-                                                                GetSolidDisplacement(),
+                    std::make_unique<Private::SolidOnFluidMesh>(GetSolidDisplacement(),
+                                                                input_parameter_data,
                                                                 time_manager);
             }
 

Sources/ModelInstances/UnderDevelopment/Poromechanics/Private/PorosityParameter.cpp

@@ -94,15 +94,16 @@ namespace HappyHeart
                 }
                 
                 {
+                
                     using type = NonConformInterpolatorNS::FromVertexMatching;
                     
                     auto& init_vertex_matching_manager =
-                        ::HappyHeart::Private::DofProgramWiseIndexListPerVertexCoordIndexManager::GetInstance();
+                    ::HappyHeart::Private::DofProgramWiseIndexListPerVertexCoordIndexManager::GetInstance();
                     
                     from_solid_to_fluid_ =
-                        std::make_unique<type>(input_parameter_data,
-                                               init_vertex_matching_manager.GetDofProgramWiseIndexListPerVertexCoordIndex(EnumUnderlyingType(InitVertexMatchingInterpolator::solid)),
-                                               init_vertex_matching_manager.GetDofProgramWiseIndexListPerVertexCoordIndex(EnumUnderlyingType(InitVertexMatchingInterpolator::fluid)));
+                    std::make_unique<type>(input_parameter_data,
+                                           init_vertex_matching_manager.GetDofProgramWiseIndexListPerVertexCoordIndex(EnumUnderlyingType(InitVertexMatchingInterpolator::solid_porosity)),
+                                           init_vertex_matching_manager.GetDofProgramWiseIndexListPerVertexCoordIndex(EnumUnderlyingType(InitVertexMatchingInterpolator::fluid_porosity)));
                     
                 }
                 

Sources/ModelInstances/UnderDevelopment/Poromechanics/Private/PorosityParameter.hpp

@@ -115,6 +115,7 @@ namespace HappyHeart
                 //! Constant accessor to the solution of the PorosityVariationalFormulation.
                 const GlobalVector& GetSolutionPorosityVarf() const noexcept;
                 
+                
                 //! Non constant accessor to the porosity parameter on the solid mesh.
                 ParameterAtDof<ParameterNS::Type::scalar>::type& GetNonCstOnSolidMesh() noexcept;
             
@@ -130,12 +131,14 @@ namespace HappyHeart
                 //! Non constant accessor to the underlying global vector for porosity on fluid mesh.
                 GlobalVector& GetNonCstOnFluidMeshVector() noexcept;
                 
+                
                 /*!
                  * \brief Constant accessor to the underlying global vector for porosity on fluid mesh from previous
                  * time iteration.
                  */
                 const GlobalVector& GetOnFluidMeshVectorPreviousTimeIteration() const noexcept;
                 
+                
                 /*!
                  * \brief Non constant accessor to the underlying global vector for porosity on fluid mesh from
                  * previous time iteration.

Sources/ModelInstances/UnderDevelopment/Poromechanics/Private/PorosityParameter.hxx

@@ -101,8 +101,7 @@ namespace HappyHeart
             }
             
             
-            inline const NonConformInterpolatorNS::FromVertexMatching& PorosityParameter
-            ::GetFromSolidToFluid() const noexcept
+            inline const NonConformInterpolatorNS::FromVertexMatching& PorosityParameter::GetFromSolidToFluid() const noexcept
             {
                 assert(!(!from_solid_to_fluid_));
                 return *from_solid_to_fluid_;

Sources/ModelInstances/UnderDevelopment/Poromechanics/Private/SolidOnFluidMesh.cpp

@@ -13,6 +13,7 @@
 #include "Core/NumberingSubset.hpp"
 #include "Core/TimeManager/TimeManager.hpp"
 
+#include "FiniteElement/FiniteElementSpace/GodOfDofManager.hpp"
 #include "FiniteElement/FiniteElementSpace/GodOfDof.hpp"
 #include "FiniteElement/FiniteElementSpace/Private/DofProgramWiseIndexListPerVertexCoordIndexManager.hpp"
 
@@ -34,38 +35,45 @@ namespace HappyHeart
         {
             
           
-            SolidOnFluidMesh::SolidOnFluidMesh(const GodOfDof& god_of_dof,
-                                               const GlobalVector& solid_displacement_on_solid_mesh,
+            SolidOnFluidMesh::SolidOnFluidMesh(const GlobalVector& solid_displacement_on_solid_mesh,
+                                               const InputParameterList& input_parameter_data,
                                                const TimeManager& time_manager)
             : time_manager_(time_manager),
             solid_displacement_on_solid_mesh_(solid_displacement_on_solid_mesh)
             {
+                const auto& god_of_dof_manager = GodOfDofManager::GetInstance();
+                const auto& solid_god_of_dof = god_of_dof_manager.GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));
+                const auto& fluid_god_of_dof = god_of_dof_manager.GetGodOfDof(EnumUnderlyingType(MeshIndex::fluid));
+        
+                const auto& solid_vel_numbering_subset = solid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::solid_velocity));
+                
                 former_displacement_on_solid_mesh_ = std::make_unique<GlobalVector>(solid_displacement_on_solid_mesh);
                 former_displacement_on_solid_mesh_->ZeroEntries(__FILE__, __LINE__);
                 
-                half_sum_velocity_on_solid_mesh_ = std::make_unique<GlobalVector>(solid_displacement_on_solid_mesh);
+                half_sum_velocity_on_solid_mesh_ = std::make_unique<GlobalVector>(solid_vel_numbering_subset);
+                AllocateGlobalVector(solid_god_of_dof, *half_sum_velocity_on_solid_mesh_);
                 
                 {
                     const auto& numbering_subset =
-                        god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::solid_displacement_in_fluid_mesh));
+                        fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::solid_displacement_in_fluid_mesh));
                     
                     current_displacement_ = std::make_unique<GlobalVector>(numbering_subset);
                 }
                 
                 auto& current_displacement = GetNonCstCurrentDisplacement();
-                AllocateGlobalVector(god_of_dof, current_displacement);
+                AllocateGlobalVector(fluid_god_of_dof, current_displacement);
 
                 {
                     const auto& numbering_subset =
-                        god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity));
+                        fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity));
 
                     half_sum_velocity_ = std::make_unique<GlobalVector>(numbering_subset);
                 }
                 auto& half_sum_velocity = GetNonCstHalfSumVelocity();
-                AllocateGlobalVector(god_of_dof, half_sum_velocity);
+                AllocateGlobalVector(fluid_god_of_dof, half_sum_velocity);
                 
                 
-                InitInterpolationMatrix(god_of_dof);
+                InitInterpolationMatrix(input_parameter_data);
             }
             
             
@@ -104,44 +112,72 @@ namespace HappyHeart
                                              half_sum_vel,
                                              __FILE__, __LINE__);
                 }
-                
-                
             }
             
             
-            void SolidOnFluidMesh::InitInterpolationMatrix(const GodOfDof& god_of_dof)
+            void SolidOnFluidMesh::InitInterpolationMatrix( const InputParameterList& input_parameter_data)
             {
-                decltype(auto) felt_space = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::fluid));
+                decltype(auto) god_of_dof_manager = GodOfDofManager::GetInstance();
+                
+                // =============================================
+                // First determine the matrix that transforms a P1 solid velocity on solid mesh to a P1 velocity
+                // on fluid mesh. It might not be exactly identity in parallel!
+                // =============================================
+              
+                auto& init_vertex_matching_manager =
+                ::HappyHeart::Private::DofProgramWiseIndexListPerVertexCoordIndexManager::GetInstance();
+                
+                NonConformInterpolatorNS::FromVertexMatching
+                    from_solid_to_fluid(input_parameter_data,
+                                        init_vertex_matching_manager.GetDofProgramWiseIndexListPerVertexCoordIndex(EnumUnderlyingType(InitVertexMatchingInterpolator::solid_velocity)),
+                                        init_vertex_matching_manager.GetDofProgramWiseIndexListPerVertexCoordIndex(EnumUnderlyingType(InitVertexMatchingInterpolator::fluid_velocity)));
+
                 
-                decltype(auto) unknown_manager = UnknownManager::GetInstance();
+
                 
-                decltype(auto) fluid_vel_unknown_ptr =
-                    unknown_manager.GetUnknownPtr(EnumUnderlyingType(UnknownIndex::fluid_velocity));
+                // =============================================
+                // Then determine on the fluid mesh the matrix that transforms a P1 velocity into a P1b velocity
+                // =============================================
+                decltype(auto) fluid_god_of_dof =
+                    god_of_dof_manager.GetGodOfDof(EnumUnderlyingType(MeshIndex::fluid));
+                
+                decltype(auto) felt_space = fluid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::fluid));
+                
+                decltype(auto) unknown_manager = UnknownManager::GetInstance();
+                                
+                decltype(auto) vel_solid_numbering_subset =
+                    fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::solid_velocity_on_fluid));
                 
-                decltype(auto) solid_disp_on_fluid_unknown_ptr =
-                    unknown_manager.GetUnknownPtr(EnumUnderlyingType(UnknownIndex::solid_displacement));
+                decltype(auto) vel_fluid_numbering_subset =
+                    fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity));
                 
-                decltype(auto) solid_disp_numbering_subset =
-                    felt_space.GetNumberingSubset(*solid_disp_on_fluid_unknown_ptr);
+                decltype(auto) solid_velocity =
+                    unknown_manager.GetUnknownPtr(EnumUnderlyingType(UnknownIndex::solid_velocity));
                 
-                decltype(auto) fluid_vel_numbering_subset =
-                    felt_space.GetNumberingSubset(*fluid_vel_unknown_ptr);
+                decltype(auto) fluid_velocity =
+                    unknown_manager.GetUnknownPtr(EnumUnderlyingType(UnknownIndex::fluid_velocity));
                 
-                matrix_solid_to_fluid_ = std::make_unique<GlobalMatrix>(fluid_vel_numbering_subset,
-                                                                        solid_disp_numbering_subset);
+                matrix_solid_to_fluid_ = std::make_unique<GlobalMatrix>(vel_fluid_numbering_subset,
+                                                                        vel_solid_numbering_subset);
                 
-                auto pairing = { std::make_pair(solid_disp_on_fluid_unknown_ptr, fluid_vel_unknown_ptr) };
+                auto pairing = { std::make_pair(solid_velocity, fluid_velocity) };
                 
                 ConformInterpolatorNS::P1_to_P1b interpolator_P1_to_P1b(felt_space,
-                                                                        solid_disp_numbering_subset,
+                                                                        vel_solid_numbering_subset,
                                                                         felt_space,
-                                                                        fluid_vel_numbering_subset,
+                                                                        vel_fluid_numbering_subset,
                                                                         std::move(pairing));
                 
                 interpolator_P1_to_P1b.Init();
                 
-                GetNonCstMatrixSolidToFluid().CompleteCopy(interpolator_P1_to_P1b.GetInterpolationMatrix(),
-                                                           __FILE__, __LINE__);
+                // =============================================
+                // Sought matrix is the product of both matrices.
+                // =============================================
+
+                Wrappers::Petsc::MatMatMult(interpolator_P1_to_P1b.GetInterpolationMatrix(),
+                                            from_solid_to_fluid.GetInterpolationMatrix(),
+                                            GetNonCstMatrixSolidToFluid(),
+                                            __FILE__, __LINE__);
             }
             
             

Sources/ModelInstances/UnderDevelopment/Poromechanics/Private/SolidOnFluidMesh.hpp

@@ -16,6 +16,8 @@
 
 # include "Core/LinearAlgebra/GlobalVector.hpp"
 
+# include "Operators/NonConformInterpolator/FromVertexMatching.hpp"
+
 # include "ModelInstances/UnderDevelopment/Poromechanics/InputParameterList.hpp"
 
 
@@ -75,12 +77,11 @@ namespace HappyHeart
                 /*!
                  * \brief Constructor.
                  *
-                 * \param[in] god_of_dof God of dof related to the fluid mesh.
                  * \param[in] solid_displacement_on_solid_mesh Solid displacement as computed on the solid mesh.
-                 * \param[in] time_manager Time manager.
+                 * \copydetails doxygen_hide_input_parameter_data_arg
                  */
-                explicit SolidOnFluidMesh(const GodOfDof& god_of_dof,
-                                          const GlobalVector& solid_displacement_on_solid_mesh,
+                explicit SolidOnFluidMesh(const GlobalVector& solid_displacement_on_solid_mesh,
+                                          const InputParameterList& input_parameter_data,
                                           const TimeManager& time_manager);
                 
                 
@@ -129,9 +130,9 @@ namespace HappyHeart
                  * \brief Init the interpolation matrix that transforms solid displacement on solid mesh (P1) into 
                  * P1b displacement on fluid mesh.
                  *
-                 * \param[in] god_of_dof God of dof related to the fluid mesh.
+                 * \copydetails doxygen_hide_input_parameter_data_arg
                  */
-                void InitInterpolationMatrix(const GodOfDof& god_of_dof);
+                void InitInterpolationMatrix(const InputParameterList& input_parameter_data);
                 
                 
             private:
@@ -232,6 +233,8 @@ namespace HappyHeart
                 
                 //! Matrix that transforms a P1 displacement on solid mesh into a P1b displacement on fluid mesh.
                 GlobalMatrix::unique_ptr matrix_solid_to_fluid_ = nullptr;
+                
+                
             };
             
             

Sources/ModelInstances/UnderDevelopment/Poromechanics/Private/SolidOnFluidMesh.hxx

@@ -100,7 +100,7 @@ namespace HappyHeart
                 return const_cast<GlobalVector&>(GetHalfSumVelocityOnSolidMesh());
             }
 
-            
+ 
             
         } // namespace Private
         

Sources/Operators/NonConformInterpolator/FromVertexMatching.hpp

@@ -79,6 +79,9 @@ namespace HappyHeart
             //! Alias to unique pointer.
             using unique_ptr = std::unique_ptr<self>;
 
+            //! Alias to unique pointer to const object.
+            using const_unique_ptr = std::unique_ptr<const self>;
+
             
         public:
             

TAGS

@@ -4,6 +4,7 @@ means week 46 of year 2013. Should a second tag be given the same week, an index
 
 Next tag:
 
+- Bug #880: Poromechanics code was not working in parallel.
 - Bug #879: In Poromechanics, porosity was evaluated on the wrong god of dof. It is now properly evaluated on the solid one and only then interpolated to the fluid one.