#782 Problem with the fibers in parallel.

Data/Lua/demo_input_bidomain_2D_MS.lua

@@ -21,7 +21,7 @@ transient = {
 	-- Maximum time, if set to zero run a static case.
 	-- Expected format: VALUE
 	-- Constraint: v >= 0.
-	timeMax = 0.1
+	timeMax = 10
 
 } -- transient
 

Data/Lua/demo_input_cardiac_mechanics.lua

@@ -715,7 +715,7 @@ Fiber_vector_1 = {
 Result = {
     -- Directory in which all the results will be written.
     -- Expected format: "VALUE"
-    output_directory = "/Volumes/Data/${USER}/HappyHeart/Results/CardiacMechanicsVentriHyperViscoActive",
+    output_directory = "/Volumes/Data/${USER}/HappyHeart/Results/CardiacMechanicsVentriHyperViscoActive_MPI",
     
     -- Enables to skip some printing in the console. Can be used to WriteSolution every n time.
     -- Expected format: VALUE

HappyHeart.xcodeproj/xcshareddata/xcschemes/Bidomain_MS.xcscheme

@@ -42,7 +42,7 @@
       </AdditionalOptions>
    </TestAction>
    <LaunchAction
-      buildConfiguration = "Debug"
+      buildConfiguration = "Release"
       selectedDebuggerIdentifier = "Xcode.DebuggerFoundation.Debugger.LLDB"
       selectedLauncherIdentifier = "Xcode.DebuggerFoundation.Launcher.LLDB"
       launchStyle = "0"

HappyHeart.xcodeproj/xcshareddata/xcschemes/CardiacMechanics-mpi.xcscheme

@@ -42,9 +42,9 @@
       </AdditionalOptions>
    </TestAction>
    <LaunchAction
-      buildConfiguration = "Debug"
-      selectedDebuggerIdentifier = "Xcode.DebuggerFoundation.Debugger.LLDB"
-      selectedLauncherIdentifier = "Xcode.DebuggerFoundation.Launcher.LLDB"
+      buildConfiguration = "Release"
+      selectedDebuggerIdentifier = ""
+      selectedLauncherIdentifier = "Xcode.IDEFoundation.Launcher.PosixSpawn"
       launchStyle = "0"
       useCustomWorkingDirectory = "NO"
       ignoresPersistentStateOnLaunch = "NO"

Sources/ModelInstances/UnderDevelopment/CardiacMechanics/Model.cpp

@@ -82,8 +82,9 @@ namespace HappyHeart
             }
                 
             formulation.WriteSolution(GetTimeManager(), displacement_numbering_subset);
+
             formulation.PrepareDynamicRuns();
-            formulation.ComputeEnergy();
+                        formulation.ComputeEnergy();
         }
         
         

Sources/ModelInstances/UnderDevelopment/CardiacMechanics/VariationalFormulation.cpp

@@ -24,7 +24,7 @@ namespace HappyHeart
     {
         
         
-        //int increment_file = 0;
+        int increment_file = 0;
         
         VariationalFormulation::VariationalFormulation(const Wrappers::Mpi& mpi,
                                                        const NumberingSubset& displacement_numbering_subset,
@@ -100,6 +100,7 @@ namespace HappyHeart
             }
             
             DefineStaticOperators(input_parameter_data);
+            ComputeApexElevation();
             AssembleStaticOperators();
         }
 
@@ -179,6 +180,9 @@ namespace HappyHeart
             
             MatchDofInNumberingSubset* raw_match_dof_in_numbering_subset_operator = match_dof_in_numbering_subset_operator_.get();
             
+            find_coords_of_global_vector_operator_ = std::make_unique<FindCoordsOfGlobalVector>(felt_space_highest_dimension,
+                                                                                                GetSystemRhs(electrical_activation_numbering_subset));
+            
             stiffness_operator_ =
                 std::make_unique<StiffnessOperatorType>(felt_space_highest_dimension,
                                                         displacement,
@@ -274,7 +278,7 @@ namespace HappyHeart
             // Get The list of dofs of the BC.
 //            const auto& displacement_numbering_subset = GetDisplacementNumberingSubset();
 //            const auto& dirichlet_bc_list = GetEssentialBoundaryConditionList();
-            
+//            
 //            std::vector<unsigned int> dofs_bc;
 //            
 //            for (const auto& dirichlet_bc_ptr : dirichlet_bc_list)
@@ -328,7 +332,7 @@ namespace HappyHeart
             AssembleDynamicOperators();
             
             UpdateForNextTimeStep();
-            
+
             InitializeElectricalActivation();
         }
         
@@ -350,11 +354,35 @@ namespace HappyHeart
                                                          mesh_dimension,
                                                          default_quadrature_rule_set,
                                                          DoComputeProcessorWiseLocal2Global::no);
+        }
+        
+        
+        void VariationalFormulation::ComputeApexElevation()
+        {
+            const auto& find_coords_of_global_vector_operator = GetFindCoordsOfGlobalVectorOperator();
             
-            const auto& electrical_activation_numbering_subset = GetElectricalActivationNumberingSubset();
+            const auto coords_list_ptr = find_coords_of_global_vector_operator.GetCoordsList();
+    
+            const auto it = std::max_element(coords_list_ptr.cbegin(),
+                                             coords_list_ptr.cend(),
+                                             [](const auto& lhs, const auto& rhs)
+                                             {
+                                                 assert(!(!lhs));
+                                                 assert(!(!rhs));
+                                                 
+                                                 return std::fabs(lhs->z()) < std::fabs(rhs->z());
+                                             });
             
-            find_coords_of_global_vector_operator_ = std::make_unique<FindCoordsOfGlobalVector>(felt_space_highest_dimension,
-                                                                                                GetSystemRhs(electrical_activation_numbering_subset));
+            double local_apex_elevation = (*it)->z();
+            
+            apex_elevation_ = MpiHappyHeart().AllReduce(local_apex_elevation,
+                                                        Wrappers::MpiNS::Op::Min);            
+            
+            if (MpiHappyHeart().IsRootProcessor())
+            {
+                std::cout << std::endl << "Apex elevation : ";
+                std::cout << std::scientific << std::setprecision(12) << GetApexElevation() << std::endl;
+            }
         }
         
         
@@ -377,6 +405,12 @@ namespace HappyHeart
         void VariationalFormulation::InitializeElectricalActivation()
         {
             ComputeElectricalActivationAtTime(GetTimeManager().GetTime());
+            
+            const auto& electrical_activation_numbering_subset = GetElectricalActivationNumberingSubset();
+            
+            GetNonCstSystemRhs(electrical_activation_numbering_subset).UpdateGhosts(__FILE__, __LINE__);
+            
+            GetNonCstSystemSolution(electrical_activation_numbering_subset).UpdateGhosts(__FILE__, __LINE__);
         }
         
         
@@ -450,50 +484,29 @@ namespace HappyHeart
         
         void VariationalFormulation::ComputeElectricalActivationAtTime(double time)
         {
-            const auto& electrical_activation_numbering_subset = GetElectricalActivationNumberingSubset();
-            
-            auto& u1 = GetNonCstSystemSolution(electrical_activation_numbering_subset);
-            
             const auto& find_coords_of_global_vector_operator = GetFindCoordsOfGlobalVectorOperator();
             
+            const auto coords_list_ptr = find_coords_of_global_vector_operator.GetCoordsList();
+            
+            double zi = 0.;
+            double apex_elevation = GetApexElevation();
             const double u_max = 15.;
             const double u_min = -5.;
             const double velocity = 0.5;
             const double long_bat = 0.8;
             const double tau_delay = -0.1;
             
+            const auto& electrical_activation_numbering_subset = GetElectricalActivationNumberingSubset();
+            
+            auto& u1 = GetNonCstSystemSolution(electrical_activation_numbering_subset);
+            
             Wrappers::Petsc::AccessVectorContent<Utilities::Access::read_and_write> u1_access(u1, __FILE__, __LINE__);
             
-            const unsigned int Nnode = u1_access.GetSize(__FILE__, __LINE__);
+            const unsigned int Nnode_processor = u1_access.GetSize(__FILE__, __LINE__);
             double t = 0.;
-            
-            const auto coords_list_ptr = find_coords_of_global_vector_operator.GetCoordsList();
-            double zi = 0.;
-
             double nagumo = 0.;
-            
-            double apex_elevation = 0.;
-            
-            
-//            const auto it = std::max_element(coords_list_ptr.cbegin(),
-//                                             coords_list_ptr.cend(),
-//                                             [](const auto& lhs, const auto& rhs)
-//                                               {
-//                                                   assert(!(!lhs));
-//                                                   assert(!(!rhs));
-//                                                   
-//                                                   return std::fabs(lhs->z()) < std::fabs(rhs->z());
-//                                               });
-//            
-//            double apex_elevation_2 = (*it)->z();
-//            
-            for (unsigned int i = 0; i < Nnode; ++i)
-            {
-                if (std::fabs(coords_list_ptr[i]->z()) > apex_elevation)
-                    apex_elevation = coords_list_ptr[i]->z();
-            }
 
-            for (unsigned int i = 0; i < Nnode; ++i)
+            for (unsigned int i = 0; i < Nnode_processor; ++i)
             {
                 zi = coords_list_ptr[i]->z();
                 
@@ -522,7 +535,7 @@ namespace HappyHeart
             UpdateVelocityBetweenTimeStep();
             ComputeGuessForNextTimeStep();
             
-            //++increment_file;
+            ++increment_file;
         }
         
         
@@ -538,7 +551,6 @@ namespace HappyHeart
             Wrappers::Petsc::AXPY(interpolation * GetTimeManager().GetTimeStep(), velocity_previous_time_iteration, system_solution, __FILE__, __LINE__);
             
             system_solution.UpdateGhosts(__FILE__, __LINE__);
-            
         }
         
         void VariationalFormulation::UpdateDisplacementAtNewtonIteration(Vec petsc_vec_displacement_at_newton_iteration)
@@ -706,7 +718,7 @@ namespace HappyHeart
                                   __FILE__, __LINE__);
             
 
-//            GetVectorStiffnessResidual().Print<MpiScale::processor_wise>(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/Active/RHS_Stiffness/rhs_cardiac_" + std::to_string(EnumUnderlyingType(run_case)) + "_" + std::to_string(increment_file) + "_" +  std::to_string(GetSnes().GetSnesIteration(__FILE__, __LINE__)) + ".hhdata", __FILE__, __LINE__);
+            //GetVectorStiffnessResidual().Print<MpiScale::processor_wise>(MpiHappyHeart(), "/Volumes/Data/gautier/HappyHeart/Computation/Cardiac/Active/RHS_Stiffness/rhs_cardiac_" + std::to_string(EnumUnderlyingType(run_case)) + "_" + std::to_string(increment_file) + "_" +  std::to_string(GetSnes().GetSnesIteration(__FILE__, __LINE__)) + ".hhdata", __FILE__, __LINE__);
             
             if (run_case == StaticOrDynamic::static_)
             {
@@ -879,8 +891,8 @@ namespace HappyHeart
                                                           evaluation_state_max
                                                           );
             
-//            if (NumericNS::AreEqual(variational_formulation.GetTimeManager().GetTime(), 3.e-03))
-//                ExceptionNS::PrintAndAbort(variational_formulation.MpiHappyHeart(), "STOP HERE");
+            if (NumericNS::AreEqual(variational_formulation.GetTimeManager().GetTime(), 3.e-03))
+                ExceptionNS::PrintAndAbort(variational_formulation.MpiHappyHeart(), "STOP HERE");
             
             return 0;
         }

Sources/ModelInstances/UnderDevelopment/CardiacMechanics/VariationalFormulation.hpp

@@ -84,6 +84,15 @@ namespace HappyHeart
             using ActiveStressPolicy =
                 GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensorNS::ActiveStressPolicyNS::AnalyticalPrestressFiber<EnumUnderlyingType(FiberIndex::fiber)>;
             
+            using HyperNone =
+                GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensorNS::HyperelasticityPolicyNS::None;
+            
+            using ViscoNone =
+                GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensorNS::ViscoelasticityPolicyNS::None;
+            
+            using ActiveNone =
+                GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensorNS::ActiveStressPolicyNS::None;
+            
             using StiffnessOperatorType =
                 GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensor<HyperelasticityLaw, ViscoelasticityPolicy, ActiveStressPolicy>;
             
@@ -292,6 +301,9 @@ namespace HappyHeart
             //! Verify if the tangent is quadratique.
             void VerifyTangent(GlobalVector& rhs);
             
+            //! Compute the apex elevation.
+            void ComputeApexElevation();
+            
         private:
             
             /*!
@@ -414,6 +426,10 @@ namespace HappyHeart
             
             const std::string& GetEnergyFile() const noexcept;
             
+        private:
+            
+            const double GetApexElevation() const noexcept;
+            
         private:
             
             /*!
@@ -524,6 +540,9 @@ namespace HappyHeart
             //! Indicates if the tangent was quadratique over the whole simulation.
             bool tangent_non_quadratique_ = false;
             
+        private:
+            
+            double apex_elevation_ = 0.;
             
             
             
@@ -532,15 +551,6 @@ namespace HappyHeart
         // #782 Here just for tests.
         private:
 
-            using HyperNone =
-                GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensorNS::HyperelasticityPolicyNS::None;
-            
-            using ViscoNone =
-                GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensorNS::ViscoelasticityPolicyNS::None;
-            
-            using ActiveNone =
-                GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensorNS::ActiveStressPolicyNS::None;
-                        
             using StiffnessOperatorHyper =
                 GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensor<HyperelasticityLaw, ViscoNone, ActiveNone>;
             

Sources/ModelInstances/UnderDevelopment/CardiacMechanics/VariationalFormulation.hxx

@@ -295,6 +295,12 @@ namespace HappyHeart
         {
             return energy_file_;
         }
+        
+        
+        inline const double VariationalFormulation::GetApexElevation() const noexcept
+        {
+            return apex_elevation_;
+        }
 
         
         // #782 Temporary force.

Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/SecondPiolaKirchhoffStressTensor/ActiveStressPolicy/AnalyticalPrestressFiber.hpp

@@ -156,6 +156,8 @@ namespace HappyHeart
                     
                     void SetDoUpdateSigmaC(const bool do_update);
                     
+                    const ScalarParameterAtQuadPt& GetSigmaC() const noexcept;
+                    
                 private:
                     
                     const TimeManager& GetTimeManager() const noexcept;
@@ -164,8 +166,6 @@ namespace HappyHeart
                     
                     ScalarParameterAtQuadPt& GetNonCstSigmaC() noexcept;
                     
-                    const ScalarParameterAtQuadPt& GetSigmaC() const noexcept;
-                    
                     void SetSigmaC(ScalarParameterAtQuadPt* sigma_c);
                     
                 private:

Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/SecondPiolaKirchhoffStressTensor/ActiveStressPolicy/AnalyticalPrestressFiber.hxx

@@ -80,6 +80,19 @@ namespace HappyHeart
                         // Fiber part of the operator.
                         const auto& tau_interpolate = fibers.GetValue(quad_pt, geom_elt);
                         
+                        if (geom_elt.GetIndex() == 193)
+                        {
+                            std::cout << "tau_interpolate      " << ' ' << std::scientific << std::setprecision(12) << tau_interpolate(0) << ' ' << tau_interpolate(1) << ' ' << tau_interpolate(2) << std::endl;
+                            
+                            int Nnode = geom_elt.Nvertex();
+                            for (int i = 0 ; i < Nnode ; ++i)
+                            {
+                                const auto& coords_in_geom_elt = geom_elt.GetCoord(i);
+                            
+                                std::cout << "coords_in_geom_elt     " << ' '<< coords_in_geom_elt << std::endl;
+                            }
+                        }
+                        
                         const int dimension = static_cast<int>(DimensionT);
                         
                         double norm = 0.;
@@ -140,7 +153,19 @@ namespace HappyHeart
                             new_sigma_c = GetNonCstSigmaC().GetValue(quad_pt, geom_elt);
                         }
                         
+                        if (geom_elt.GetIndex() == 193)
+                        {
+                            std::cout << "dW before      " << ' ' << std::scientific << std::setprecision(12) << dW(0) << ' ' << dW(1) << ' ' << dW(2) << ' ' << dW(3) << ' ' << dW(4) << ' ' << dW(5) << std::endl;
+                            std::cout << "new_sigma_c      " <<  ' ' << std::scientific << std::setprecision(12) << new_sigma_c <<  std::endl;
+                            std::cout << "tauxtau      " << ' ' << std::scientific << std::setprecision(12) << tauxtau(0) << ' ' << tauxtau(1) << ' ' << tauxtau(2) << ' ' << tauxtau(3) << ' ' << tauxtau(4) << ' ' << tauxtau(5) << std::endl;
+                        }
+                        
                         Add(new_sigma_c, tauxtau, dW);
+                        
+                        if (geom_elt.GetIndex() == 193)
+                        {
+                            std::cout << "dW after      " << ' ' << std::scientific << std::setprecision(12) << dW(0) << ' ' << dW(1) << ' ' << dW(2) << ' ' << dW(3) << ' ' << dW(4) << ' ' << dW(5) << std::endl << std::endl;
+                        }
                     }
                 }
                 

Sources/Parameters/Policy/AtDof.hxx

@@ -73,6 +73,42 @@ namespace HappyHeart
                 vector_with_ghost_content(vector_with_ghost, __FILE__, __LINE__);
 
                 const auto& basic_ref_elt = ref_felt_space.GetRefFElt(felt_space, unknown).GetBasicRefFElt();
+                
+                // #782
+                if (geom_elt.GetIndex() == 193)
+                {
+                    const auto Nnode = static_cast<unsigned int>(node_list.size());
+                    
+                    std::cout << " vector_with_ghost_content.GetValue(processor_wise_index)      " << std::endl;
+                    
+                    for (auto local_node_index = 0u; local_node_index < Nnode; ++local_node_index)
+                    {
+                        const auto& node_ptr = node_list[static_cast<std::size_t>(local_node_index)];
+                        assert(!(!node_ptr));
+                        const auto& node = *node_ptr;
+                        assert(node.DoBelongToNumberingSubset(numbering_subset));
+                        
+                        const auto& dof_list = node.GetDofList();
+                        
+                        const auto Ndof = static_cast<int>(dof_list.size());
+                        
+                        for (auto dof_index = 0; dof_index < Ndof; ++dof_index)
+                        {
+                            const auto& dof_ptr = dof_list[static_cast<std::size_t>(dof_index)];
+                            assert(!(!dof_ptr));
+                            const auto processor_wise_index = dof_ptr->GetProcessorWiseOrGhostIndex(numbering_subset);
+                            assert(dof_index < value.GetSize());
+                            
+                            std::cout << vector_with_ghost_content.GetValue(processor_wise_index) << ' ';
+                            
+                        }
+                        
+                        std::cout << std::endl;
+                        
+                    }
+                    
+                    std::cout << std::endl;
+                }
 
                 Private::AtDofNS::ComputeLocalValue(basic_ref_elt,
                                                     quad_pt,

Sources/ThirdParty/Wrappers/Mpi/MacroEncapsulation/Op.cpp

@@ -34,6 +34,8 @@ namespace HappyHeart
                         return MPI_SUM;
                     case Op::Max:
                         return MPI_MAX;
+                    case Op::Min:
+                        return MPI_MIN;
                     case Op::LogicalOr:
                         return MPI_LOR;
                 }

Sources/ThirdParty/Wrappers/Mpi/MacroEncapsulation/Op.hpp

@@ -39,6 +39,7 @@ namespace HappyHeart
             {
                 Sum,
                 Max,
+                Min,
                 LogicalOr
             };