#820 Add midpoint position in implicit step fluid (at the moment we assume...

#820 Add midpoint position in implicit step fluid (at the moment we assume midpoint is used throughout the model; a static_assert checks it).

HappyHeart.xcodeproj/project.pbxproj

@@ -4398,8 +4398,6 @@
 		BE1EBB0F1CA0A59600EC0EAA /* Model.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; name = Model.cpp; path = Coloring/Model.cpp; sourceTree = "<group>"; };
 		BE1EBB101CA0A59600EC0EAA /* Model.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; name = Model.hpp; path = Coloring/Model.hpp; sourceTree = "<group>"; };
 		BE1EBB111CA0A59600EC0EAA /* Model.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; name = Model.hxx; path = Coloring/Model.hxx; sourceTree = "<group>"; };
-		BE2195771D40B28F00052F9C /* DeltaSolidVarf.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = DeltaSolidVarf.hpp; sourceTree = "<group>"; };
-		BE2195781D40B28F00052F9C /* DeltaSolidVarf.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = DeltaSolidVarf.hxx; sourceTree = "<group>"; };
 		BE2241BC19FA5FDE00B90563 /* InterfaceSpecialization.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = InterfaceSpecialization.cpp; sourceTree = "<group>"; };
 		BE2241BD19FA5FDE00B90563 /* InterfaceSpecialization.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = InterfaceSpecialization.hpp; sourceTree = "<group>"; };
 		BE2241BE19FA5FDE00B90563 /* InterfaceSpecialization.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = InterfaceSpecialization.hxx; sourceTree = "<group>"; };
@@ -9617,8 +9615,6 @@
 				BEB7B2A31D12EC9E00FC333C /* VariationalFormulationComputeMatrixBlock.hxx */,
 				BEAE11F11D07064A00967C19 /* VariationalFormulationAccessors.hxx */,
 				BEAE11F21D07064A00967C19 /* VariationalFormulationInit.hxx */,
-				BE2195771D40B28F00052F9C /* DeltaSolidVarf.hpp */,
-				BE2195781D40B28F00052F9C /* DeltaSolidVarf.hxx */,
 			);
 			path = InnerLoop;
 			sourceTree = "<group>";

Sources/ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/SolidDeltaResidual.hpp

@@ -131,13 +131,11 @@ namespace HappyHeart
                  * of the finite element space, assembling can only occur in a subset of the domain defined in the finite
                  * element space; if current \a domain is not a subset of finite element space one, assembling will occur
                  * upon the intersection of both.
-                 * \param[in] displacement_previous_iteration Vector that includes data from the previous iteration. (its nature
-                 * varies depending on the time scheme used).
                  *
                  */
                 template<class LinearAlgebraTupleT>
                 void Assemble(LinearAlgebraTupleT&& linear_algebra_tuple,
-                              const GlobalVector& displacement_previous_iteration,
+                              const GlobalVector& midpoint_position, // \todo #820 Assumes time scheme is midpoint!
                               const Domain& domain = Domain()) const;
                 
                 

Sources/ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/SolidDeltaResidual.hxx

@@ -75,7 +75,7 @@ namespace HappyHeart
                 ExtractLocalDofValues(local_felt_space,
                                       this->GetNthUnknown(),
                                       std::get<0>(additional_arguments),
-                                      local_operator.GetNonCstFormerLocalDisplacement());
+                                      local_operator.GetNonCstMidpointPosition()); // \todo #820 Assumes time scheme is midpoint...
             }
             
             

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx

@@ -197,7 +197,6 @@ namespace HappyHeart
                     auto& robin_vector = GetNonCstNonHomogeneousRobinVector();
                     robin_vector.ZeroEntries(__FILE__, __LINE__);
                     god_of_dof.template ApplyBoundaryCondition<BoundaryConditionMethod::pseudo_elimination>(non_homogeneous_bc, robin_vector);
-                    robin_vector.View(this->MpiHappyHeart(), __FILE__, __LINE__);
                     
                     {
                         const auto filename = parent::GetOutputDirectory(numbering_subset)
@@ -346,8 +345,24 @@ namespace HappyHeart
                                                 __FILE__, __LINE__,
                                                 PETSC_VIEWER_ASCII_MATLAB);
                     }
+                    
+                    {
+                        decltype(auto) midpoint_position = this->GetVariableHolder().GetMidpointPosition();
+                        const auto filename = parent::GetOutputDirectory(numbering_subset)
+                        + "/"
+                        + DifferentialPreffix(differential::yes)
+                        + "midpoint_position_" + GetIterationTag() + ".m";
+                        
+                        midpoint_position.View(parent::MpiHappyHeart(),
+                                                filename,
+                                                __FILE__, __LINE__,
+                                                PETSC_VIEWER_ASCII_MATLAB);
+                        
+                    }
 
                     {
+                        
+                        
                         auto& delta_residual = GetNonCstDeltaResidual();
                         GlobalVectorWithCoefficient with_coeff(delta_residual, 1.);
                         // ASSEMBLE SHOULD NOT BE CALLED DIRECTLY HERE! Use the free function in Private

Sources/ModelInstances/UnderDevelopment/Poromechanics/LocalVariationalOperatorInstances/SolidDeltaResidual.hpp

@@ -156,10 +156,10 @@ namespace HappyHeart
                 void ComputeEltArray();
                 
                 //! Constant accessor to the former local displacement required by ComputeEltArray().
-                const std::vector<double>& GetFormerLocalDisplacement() const noexcept;
+                const std::vector<double>& GetMidpointPosition() const noexcept;
                 
                 //! Non constant accessor to the former local displacement required by ComputeEltArray().
-                std::vector<double>& GetNonCstFormerLocalDisplacement() noexcept;
+                std::vector<double>& GetNonCstMidpointPosition() noexcept;
                 
             private:
                 

Sources/ModelInstances/UnderDevelopment/Poromechanics/LocalVariationalOperatorInstances/SolidDeltaResidual.hxx

@@ -138,7 +138,7 @@ namespace HappyHeart
             {
                 auto& elementary_data = parent::GetNonCstElementaryData();
                 
-                const auto& local_displacement = GetFormerLocalDisplacement();
+                const auto& local_displacement = GetMidpointPosition();
                 
                 auto& tangent_matrix = this->matrix_parent::template GetLocalMatrix<EnumUnderlyingType(LocalMatrixIndex::tangent_matrix)>();
                 const auto& rhs_part = this->vector_parent::template GetLocalVector<EnumUnderlyingType(LocalVectorIndex::rhs_part)>();
@@ -348,7 +348,7 @@ namespace HappyHeart
             template <class HyperelasticLawT>
             inline const std::vector<double>&
             SolidDeltaResidual<HyperelasticLawT>
-            ::GetFormerLocalDisplacement() const noexcept
+            ::GetMidpointPosition() const noexcept
             {
                 return former_local_displacement_;
             }
@@ -376,9 +376,9 @@ namespace HappyHeart
             template <class HyperelasticLawT>
             inline std::vector<double>&
             SolidDeltaResidual<HyperelasticLawT>
-            ::GetNonCstFormerLocalDisplacement() noexcept
+            ::GetNonCstMidpointPosition() noexcept
             {
-                return const_cast<std::vector<double>&>(GetFormerLocalDisplacement());
+                return const_cast<std::vector<double>&>(GetMidpointPosition());
             }
 
 

Sources/ModelInstances/UnderDevelopment/Poromechanics/ModelInitialize.hxx

@@ -48,10 +48,6 @@ namespace HappyHeart
             decltype(auto) solid_god_of_dof = parent::GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));
             
             decltype(auto) fluid_god_of_dof = parent::GetGodOfDof(EnumUnderlyingType(MeshIndex::fluid));
-            decltype(auto) solid_mesh = solid_god_of_dof.GetGeometricMeshRegion();
-            decltype(auto) felt_space = solid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::solid));
-            
-
             
       
             // \todo #820 Order below will have to be adapted probably: porosity that is used in almost all
@@ -74,6 +70,11 @@ namespace HappyHeart
 
   
             InitVariationalFormulations(input_parameter_data);
+            
+            static_assert(SolidVariationalFormulationPolicyT::GetTimeScheme() == HyperelasticityNS::TimeScheme::midpoint,
+                          "Currently differential part of dH assumes midpoint scheme (to be changed...)");
+            GetNonCstVariableHolder().SetMidpointPosition(&SolidVariationalFormulationPolicyT::GetVariationalFormulation().GetVectorsAndMatrices().GetMidpointPosition());
+            
             InitializeMassOperator();
             InitializeLoadOnSolidOperator();
             InitSolvers(input_parameter_data);

Sources/ModelInstances/UnderDevelopment/Poromechanics/SolidVariationalFormulationPolicy/Hyperelasticity/Policy.hpp

@@ -54,6 +54,10 @@ namespace HappyHeart
                 //! Alias to hyperelastic law type.
                 using hyperelastic_law_type = HyperelasticLawT;
                 
+                //! Returns time scheme enum
+                static constexpr HyperelasticityNS::TimeScheme GetTimeScheme();
+
+                
                 
             public:
                 

Sources/ModelInstances/UnderDevelopment/Poromechanics/SolidVariationalFormulationPolicy/Hyperelasticity/Policy.hxx

@@ -22,6 +22,18 @@ namespace HappyHeart
         
         namespace Policy
         {
+            
+            
+            template
+            <
+                class HyperelasticLawT,
+                HyperelasticityNS::TimeScheme TimeSchemeT
+            >
+            constexpr HyperelasticityNS::TimeScheme Hyperelasticity<HyperelasticLawT, TimeSchemeT>
+            ::GetTimeScheme()
+            {
+                return TimeSchemeT;
+            }
 
             
             template

Sources/ModelInstances/UnderDevelopment/Poromechanics/SolidVariationalFormulationPolicy/Hyperelasticity/VariationalFormulation.hpp

@@ -102,6 +102,7 @@ namespace HappyHeart
                     ::HappyHeart::GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensorNS::ActiveStressPolicyNS::None;
                 
                 using variable_holder_parent = ::HappyHeart::PoromechanicsNS::Crtp::VariableHolder<self>;
+             
                 
             private:
                 

Sources/ModelInstances/UnderDevelopment/Poromechanics/SolidVariationalFormulationPolicy/Hyperelasticity/VariationalFormulation.hxx

@@ -23,6 +23,7 @@ namespace HappyHeart
         namespace Policy
         {
             
+                       
             
             template
             <

Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.cpp

@@ -94,6 +94,13 @@ namespace HappyHeart
             
         }
         
+        
+        void VariableHolder::SetMidpointPosition(const GlobalVector* const midpoint_position)
+        {
+            assert(midpoint_position_ == nullptr && "Should only be assigned once!");
+            midpoint_position_ = midpoint_position;
+        }
+        
              
         
     } // namespace PoromechanicsNS

Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hpp

@@ -78,6 +78,13 @@ namespace HappyHeart
             
             ///@}
             
+            
+            //! TMP!
+            void SetMidpointPosition(const GlobalVector* const midpoint_position);
+            
+            
+            const GlobalVector& GetMidpointPosition() const noexcept;
+            
         public:
             
             //! Accessor to mass vector.
@@ -141,6 +148,9 @@ namespace HappyHeart
             
             //! Material parameters related to the solid.
             Solid::unique_ptr solid_ = nullptr;
+            
+            //! Midpoint position. // \todo #820 TMP: midpoint time scheme should not be the sole choice!
+            const GlobalVector* midpoint_position_ = nullptr;
 
             
         };

Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hxx

@@ -91,6 +91,14 @@ namespace HappyHeart
             assert(!(!solid_));
             return *solid_;
         }
+        
+        
+        inline const GlobalVector& VariableHolder::GetMidpointPosition() const noexcept
+        {
+            assert(!(!midpoint_position_) && "Make sure SetMidpointPosition() was correctly called!");
+            return *midpoint_position_;
+        }
+