#721 Elasticity: remove dof source and the template parameters for volumic and surfacic sources.

HappyHeart.xcodeproj/project.pbxproj

@@ -1040,6 +1040,7 @@
 		BEB917111AD3D37A0096A3D9 /* ExtendedUnknown.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BEB9170B1AD3D37A0096A3D9 /* ExtendedUnknown.hxx */; };
 		BEBC4FC51BA71048003633B3 /* LagrangianInterpolator.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BEBC4FC21BA71048003633B3 /* LagrangianInterpolator.hpp */; };
 		BEBC4FC61BA71048003633B3 /* LagrangianInterpolator.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BEBC4FC31BA71048003633B3 /* LagrangianInterpolator.hxx */; };
+		BEBEBBBB1BDF90180011D835 /* VariationalFormulationElasticity.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BEBEBBBA1BDF90180011D835 /* VariationalFormulationElasticity.cpp */; };
 		BEC157F01A249BB9007D20EB /* TupleHelper.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BEC157EF1A249BB9007D20EB /* TupleHelper.hpp */; };
 		BEC157F21A249C3D007D20EB /* Numeric.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BEC157F11A249C3D007D20EB /* Numeric.hpp */; };
 		BEC37E0717DDC0210021BFB7 /* Accelerate.framework in Frameworks */ = {isa = PBXBuildFile; fileRef = BEC37E0617DDC0210021BFB7 /* Accelerate.framework */; };
@@ -4258,6 +4259,7 @@
 		BEBEB22019C849C200E4EA1D /* AttributeProcessorHelper.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = AttributeProcessorHelper.hxx; sourceTree = "<group>"; };
 		BEBEB22119C849C200E4EA1D /* ComputeMatrixPattern.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; lineEnding = 0; path = ComputeMatrixPattern.cpp; sourceTree = "<group>"; xcLanguageSpecificationIdentifier = xcode.lang.cpp; };
 		BEBEB22219C849C200E4EA1D /* ComputeMatrixPattern.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; lineEnding = 0; path = ComputeMatrixPattern.hpp; sourceTree = "<group>"; xcLanguageSpecificationIdentifier = xcode.lang.cpp; };
+		BEBEBBBA1BDF90180011D835 /* VariationalFormulationElasticity.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = VariationalFormulationElasticity.cpp; sourceTree = "<group>"; };
 		BEC157EF1A249BB9007D20EB /* TupleHelper.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = TupleHelper.hpp; sourceTree = "<group>"; };
 		BEC157F11A249C3D007D20EB /* Numeric.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; name = Numeric.hpp; path = Sources/Utilities/Numeric/Numeric.hpp; sourceTree = SOURCE_ROOT; };
 		BEC1B0C519C3209400088A5D /* SConscript */ = {isa = PBXFileReference; explicitFileType = text.script.python; fileEncoding = 4; path = SConscript; sourceTree = "<group>"; };
@@ -6031,6 +6033,7 @@
 				BE63B4921A31C20E003A6523 /* ElasticityModel.cpp */,
 				BE63B4931A31C20E003A6523 /* ElasticityModel.hpp */,
 				BE63B4941A31C20E003A6523 /* ElasticityModel.hxx */,
+				BEBEBBBA1BDF90180011D835 /* VariationalFormulationElasticity.cpp */,
 				BE63B4991A31C20E003A6523 /* VariationalFormulationElasticity.hpp */,
 				BE63B49A1A31C20E003A6523 /* VariationalFormulationElasticity.hxx */,
 			);
@@ -9878,6 +9881,7 @@
 			isa = PBXSourcesBuildPhase;
 			buildActionMask = 2147483647;
 			files = (
+				BEBEBBBB1BDF90180011D835 /* VariationalFormulationElasticity.cpp in Sources */,
 				BE4478871AA740F800665010 /* main.cpp in Sources */,
 			);
 			runOnlyForDeploymentPostprocessing = 0;

Sources/ModelInstances/Elasticity/ElasticityModel.hpp

@@ -42,14 +42,8 @@ namespace HappyHeart
             friend parent;
             
             //! Variational formulation.
-            using variational_formulation_type =
-                VariationalFormulationElasticity
-                <
-                    EnumUnderlyingType(SourceIndex::volumic),
-                    EnumUnderlyingType(SourceIndex::surfacic)
-                >;
-            
-            
+            using variational_formulation_type = VariationalFormulationElasticity;
+
             
         public:
             

Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.cpp

+//
+//  VariationalFormulationElasticity.hxx
+//  HappyHeart
+//
+//  Created by Sebastien Gilles on 01/07/14.
+//  Copyright (c) 2014 Inria. All rights reserved.
+//
+
+#include "ModelInstances/Elasticity/VariationalFormulationElasticity.hpp"
+
+
+namespace HappyHeart
+{
+
+
+    namespace ElasticityNS
+    {
+        
+        
+        void VariationalFormulationElasticity
+        ::SupplInit(const InputParameterList& input_parameter_data)
+        {
+            const auto& god_of_dof = this->GetGodOfDof();
+            const auto& mesh = god_of_dof.GetGeometricMeshRegion();
+            
+            volumic_mass_ = InitScalarParameter<InputParameter::Solid::VolumicMass>("Volumic mass",
+                                                                                    mesh,
+                                                                                    input_parameter_data);
+            
+            if (!GetVolumicMass().IsConstant())
+                throw Exception("Current elastic model is restricted to a constant volumic mass!", __FILE__, __LINE__);
+            
+            young_modulus_ = InitScalarParameter<InputParameter::Solid::YoungModulus>("Young modulus",
+                                                                                      mesh,
+                                                                                      input_parameter_data);
+            
+            poisson_ratio_ = InitScalarParameter<InputParameter::Solid::PoissonRatio>("Poisson ratio",
+                                                                                      mesh,
+                                                                                      input_parameter_data);
+            
+            {
+                // \todo #583: Temporary to make code work; quite stupid nonetheless!
+                const auto degree_of_exactness = DEFAULT_DEGREE_OF_EXACTNESS;
+                auto&& default_rule = DetermineDefaultQuadratureRule(degree_of_exactness);
+                this->GetNonCstQuadratureRulePerTopology().emplace_back(std::move(default_rule));
+                assert(this->GetNonCstQuadratureRulePerTopology().size() == 1ul);
+            }
+            
+            
+            DefineOperators(input_parameter_data);
+            
+            {
+                // Assemble the stiffness matrix.
+                GlobalMatrixWithCoefficient matrix(GetNonCstStiffness(), 1.);
+                GetStiffnessOperator().Assemble(std::make_tuple(std::ref(matrix)));
+            }
+        }
+        
+        
+        
+        VariationalFormulationElasticity
+        ::VariationalFormulationElasticity(const Wrappers::Mpi& mpi,
+                                           const FEltSpace& main_felt_space,
+                                           const FEltSpace& neumann_felt_space,
+                                           const Unknown& solid_displacement,
+                                           const NumberingSubset& numbering_subset,
+                                           const TimeManager& time_manager,
+                                           const GodOfDof& god_of_dof,
+                                           DirichletBoundaryCondition::vector_shared_ptr&& boundary_condition_list)
+        : parent(mpi,
+                 time_manager,
+                 god_of_dof,
+                 std::move(boundary_condition_list)),
+        main_felt_space_(main_felt_space),
+        neumann_felt_space_(neumann_felt_space),
+        solid_displacement_(solid_displacement),
+        numbering_subset_(numbering_subset)
+        {
+            assert(time_manager.IsTimeStepConstant() && "Model's current implementation relies on this assumption.");
+        }
+        
+        
+        
+        
+        void VariationalFormulationElasticity::RunStaticCase()
+        {
+            AddSourcesToRhs();
+            
+            const auto& numbering_subset = GetNumberingSubset();
+            parent::GetNonCstSystemMatrix(numbering_subset, numbering_subset).Copy(GetStiffness(), __FILE__, __LINE__);
+            
+            parent::template ApplyEssentialBoundaryCondition<OperatorNS::Nature::nonlinear>(numbering_subset, numbering_subset);
+            parent::template SolveLinear<IsFactorized::no>(numbering_subset, numbering_subset);
+            parent::DebugPrintSolutionElasticWithOneUnknown(numbering_subset, 10);
+        }
+        
+        
+        
+        void VariationalFormulationElasticity::PrepareDynamicRuns()
+        {
+            // Assemble once and for all the system matrix in dynamic case; intermediate matrices used
+            // to compute rhs at each time iteration are also computed there.
+            ComputeDynamicMatrices();
+            UpdateDisplacement();
+        }
+        
+        
+        
+        void VariationalFormulationElasticity::ComputeDynamicSystemRhs()
+        {
+            auto& rhs = this->GetNonCstSystemRhs(GetNumberingSubset());
+            
+            // Compute the system RHS. The rhs is effectively zeroed through the first MatMult call.
+            const auto& displacement_previous_time_iteration_matrix = GetMatrixDisplacementPreviousTimeIteration();
+            const auto& velocity_previous_time_iteration_matrix = GetMatrixVelocityPreviousTimeIteration();
+            
+            const auto& displacement_previous_time_iteration_vector = GetVectorDisplacementPreviousTimeIteration();
+            const auto& velocity_previous_time_iteration_vector = GetVectorVelocityPreviousTimeIteration();
+            
+            Wrappers::Petsc::MatMult(displacement_previous_time_iteration_matrix, displacement_previous_time_iteration_vector, rhs, __FILE__, __LINE__);
+            Wrappers::Petsc::MatMultAdd(velocity_previous_time_iteration_matrix, velocity_previous_time_iteration_vector, rhs, rhs, __FILE__, __LINE__);
+        }
+        
+        
+        
+        void VariationalFormulationElasticity::AddSourcesToRhs()
+        {
+            const auto& numbering_subset = GetNumberingSubset();
+            auto& rhs = parent::GetNonCstSystemRhs(numbering_subset);
+            const auto time = parent::GetTimeManager().GetTime();
+            
+            if (this->template IsOperatorActivated<SourceType::volumic>())
+            {
+                GlobalVectorWithCoefficient force_vector(rhs, 1.);
+                this->template GetForceOperator<SourceType::volumic>().Assemble(std::make_tuple(std::ref(force_vector)),
+                                                                                time);
+            }
+            
+            if (this->template IsOperatorActivated<SourceType::surfacic>())
+            {
+                GlobalVectorWithCoefficient force_vector(rhs, 1.);
+                this->template GetForceOperator<SourceType::surfacic>().Assemble(std::make_tuple(std::ref(force_vector)),
+                                                                                 time);
+            }
+        }
+        
+        
+        void VariationalFormulationElasticity
+        ::AllocateMatricesAndVectors()
+        {
+            const auto& numbering_subset = GetNumberingSubset();
+            
+            parent::AllocateGlobalMatrix(numbering_subset, numbering_subset);
+            parent::AllocateGlobalVector(numbering_subset);
+            
+            const auto& system_matrix = parent::GetSystemMatrix(numbering_subset, numbering_subset);
+            const auto& system_rhs = parent::GetSystemRhs(numbering_subset);
+            
+            mass_ = std::make_unique<GlobalMatrix>(system_matrix);
+            stiffness_ = std::make_unique<GlobalMatrix>(system_matrix);
+            matrix_displacement_previous_time_iteration_ = std::make_unique<GlobalMatrix>(system_matrix);
+            matrix_velocity_previous_time_iteration_ = std::make_unique<GlobalMatrix>(system_matrix);
+            
+            vector_velocity_previous_time_iteration_ = std::make_unique<GlobalVector>(system_rhs);
+            vector_displacement_previous_time_iteration_ = std::make_unique<GlobalVector>(system_rhs);
+        }
+        
+        
+        
+        void VariationalFormulationElasticity::ComputeDynamicMatrices()
+        {
+            const auto& numbering_subset = GetNumberingSubset();
+            auto& system_matrix = this->GetNonCstSystemMatrix(numbering_subset, numbering_subset);
+            const auto& stiffness = GetStiffness();
+            
+            {
+                GlobalMatrixWithCoefficient mass(GetNonCstMass(), 1.);
+                GetMassOperator().Assemble(std::make_tuple(std::ref(mass)));
+            }
+            
+            const auto& mass = GetMass();
+            
+            {
+                // Compute the system matrix, which won't change afterwards!
+                system_matrix.Copy(stiffness, __FILE__, __LINE__);
+                system_matrix.Scale(0.5, __FILE__, __LINE__);
+
+                const auto coefficient =
+                    2. * GetVolumicMass().GetConstantValue() / Utilities::Square(parent::GetTimeManager().GetTimeStep());
+            
+                Wrappers::Petsc::AXPY(coefficient,
+                                      mass,
+                                      system_matrix,
+                                      __FILE__, __LINE__);
+            }
+            
+            {
+                // Displacement matrix.
+                auto& displacement_previous_time_iteration_matrix = GetNonCstMatrixDisplacementPreviousTimeIteration();
+                displacement_previous_time_iteration_matrix.Copy(mass, __FILE__, __LINE__);
+                
+                const auto coefficient =
+                2. * GetVolumicMass().GetConstantValue() / Utilities::Square(parent::GetTimeManager().GetTimeStep());
+                
+                displacement_previous_time_iteration_matrix.Scale(coefficient, __FILE__, __LINE__);
+                
+                Wrappers::Petsc::AXPY(-.5,
+                                      stiffness,
+                                      displacement_previous_time_iteration_matrix,
+                                      __FILE__, __LINE__);
+                
+            }
+            
+            {
+                // Velocity matrix.
+                auto& velocity_previous_time_iteration_matrix = GetNonCstMatrixVelocityPreviousTimeIteration();
+                velocity_previous_time_iteration_matrix.Copy(mass, __FILE__, __LINE__);
+                velocity_previous_time_iteration_matrix.Scale(2. * GetVolumicMass().GetConstantValue() / parent::GetTimeManager().GetTimeStep(),
+                                              __FILE__, __LINE__);
+            }
+        }
+        
+        
+        
+        void VariationalFormulationElasticity::UpdateDisplacement()
+        {
+            GetNonCstVectorDisplacementPreviousTimeIteration().Copy(parent::GetSystemSolution(GetNumberingSubset()),
+                                                      __FILE__, __LINE__);
+        }
+        
+        
+        
+        void VariationalFormulationElasticity::UpdateVelocity()
+        {
+            const auto& displacement_previous_time_iteration_vector = GetVectorDisplacementPreviousTimeIteration();
+            const auto& system_solution = parent::GetSystemSolution(GetNumberingSubset());
+            auto& velocity_previous_time_iteration_vector = GetNonCstVectorVelocityPreviousTimeIteration();
+            
+            assert(parent::GetTimeManager().GetStaticOrDynamic() == StaticOrDynamic::dynamic_);
+            
+            {
+                // Update first the velocity.
+                Wrappers::Petsc::AccessVectorContent<Utilities::Access::read_only> solution(system_solution, __FILE__, __LINE__);
+                Wrappers::Petsc::AccessVectorContent<Utilities::Access::read_only> displacement_prev(displacement_previous_time_iteration_vector, __FILE__, __LINE__);
+                Wrappers::Petsc::AccessVectorContent<Utilities::Access::read_and_write> velocity(velocity_previous_time_iteration_vector, __FILE__, __LINE__);
+                
+                const unsigned int size = velocity.GetSize(__FILE__, __LINE__);
+                assert(size == solution.GetSize(__FILE__, __LINE__));
+                assert(size == displacement_prev.GetSize(__FILE__, __LINE__));
+                
+                const double factor = 2. / parent::GetTimeManager().GetTimeStep();
+                
+                for (unsigned int i = 0u; i < size; ++i)
+                {
+                    velocity[i] *= -1.;
+                    velocity[i] += factor * (solution.GetValue(i) - displacement_prev.GetValue(i));
+                }
+            }
+        }
+
+
+        
+        
+        
+        void VariationalFormulationElasticity::DefineOperators(const InputParameterList& input_parameter_data)
+        {
+            const auto& god_of_dof = this->GetGodOfDof();
+            const auto& felt_space_highest_dimension = GetMainFEltSpace();
+            const auto& felt_space_neumann = GetNeumannFEltSpace();
+            
+            const auto mesh_dimension = god_of_dof.GetGeometricMeshRegion().GetDimension();
+            
+            const auto& displacement = GetSolidDisplacement();
+            
+            const auto& default_quadrature_rule_set = this->GetQuadratureRulePerTopologyType();
+            
+            this->template SetIfTaggedAsActivated<SourceType::volumic>("Volumic force",
+                                                                       input_parameter_data,
+                                                                       felt_space_highest_dimension,
+                                                                       displacement,
+                                                                       mesh_dimension,
+                                                                       default_quadrature_rule_set,
+                                                                       DoComputeProcessorWiseLocal2Global::no);
+            
+            this->template SetIfTaggedAsActivated<SourceType::surfacic>("Surfacic force",
+                                                                        input_parameter_data,
+                                                                        felt_space_neumann,
+                                                                        displacement,
+                                                                        mesh_dimension,
+                                                                        default_quadrature_rule_set,
+                                                                        DoComputeProcessorWiseLocal2Global::no);
+            
+            namespace GVO = GlobalVariationalOperatorNS;
+            
+            namespace IPL = Utilities::InputParameterListNS;
+            
+            stiffness_operator_ = std::make_unique<GVO
+            ::GradOnGradientBasedElasticityTensor>(felt_space_highest_dimension,
+                                                   displacement,
+                                                   mesh_dimension,
+                                                   default_quadrature_rule_set,
+                                                   DoComputeProcessorWiseLocal2Global::no,
+                                                   GetYoungModulus(),
+                                                   GetPoissonRatio(),
+                                                   ParameterNS::ReadGradientBasedElasticityTensorConfigurationFromFile(mesh_dimension,
+                                                                                                                       input_parameter_data));
+            
+            mass_operator_ = std::make_unique<GVO::Mass>(felt_space_highest_dimension,
+                                                         displacement,
+                                                         mesh_dimension,
+                                                         default_quadrature_rule_set,
+                                                         DoComputeProcessorWiseLocal2Global::no);
+        }
+
+      
+    } // namespace ElasticityNS
+    
+    
+} // namespace HappyHeart
+

Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.hpp

@@ -24,8 +24,8 @@
 
 # include "FormulationSolver/VariationalFormulation.hpp"
 # include "FormulationSolver/Crtp/VolumicAndSurfacicSource.hpp"
-# include "FormulationSolver/DofSourcePolicy/None.hpp"
-# include "FormulationSolver/DofSourcePolicy/DofSource.hpp"
+
+# include "ModelInstances/Elasticity/InputParameterList.hpp"
 
 
 namespace HappyHeart
@@ -35,27 +35,20 @@ namespace HappyHeart
     namespace ElasticityNS
     {
         
-
-        template
-        <
-            unsigned int VolumicIndexT,
-            unsigned int SurfacicIndexT,
-            class DofSourcePolicyT = ::HappyHeart::VariationalFormulationNS::DofSourcePolicyNS::None
-        >
+        
         class VariationalFormulationElasticity final
-        : public VariationalFormulation<VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT>>,
+        : public VariationalFormulation<VariationalFormulationElasticity>,
         public Crtp::VolumicAndSurfacicSource
         <
-            VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT>,
-            VolumicIndexT,
-            SurfacicIndexT
-        >,
-        public DofSourcePolicyT
+            VariationalFormulationElasticity,
+            EnumUnderlyingType(SourceIndex::volumic),
+            EnumUnderlyingType(SourceIndex::surfacic)
+        >        
         {
             
         private:
             
-            using self = VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT>;
+            using self = VariationalFormulationElasticity;
             
             using parent = VariationalFormulation<self>;
             
@@ -74,7 +67,6 @@ namespace HappyHeart
             ///@{
             
             //! Constructor.
-            template<typename... DofSourcePolicyArgs>
             explicit VariationalFormulationElasticity(const Wrappers::Mpi& mpi,
                                                       const FEltSpace& main_felt_space,
                                                       const FEltSpace& neumann_felt_space,
@@ -82,8 +74,7 @@ namespace HappyHeart
                                                       const NumberingSubset& numbering_subset,
                                                       const TimeManager& time_manager,
                                                       const GodOfDof& god_of_dof,
-                                                      DirichletBoundaryCondition::vector_shared_ptr&& boundary_condition_list,
-                                                      DofSourcePolicyArgs&&... dof_source_policy_args);
+                                                      DirichletBoundaryCondition::vector_shared_ptr&& boundary_condition_list);
             
             //! Destructor.
             ~VariationalFormulationElasticity() = default;
@@ -173,8 +164,7 @@ namespace HappyHeart
              *
              * \internal This method is called by base class method VariationalFormulation::Init().
              */
-            template<class InputParameterDataT>
-            void SupplInit(const InputParameterDataT& input_parameter_data);
+            void SupplInit(const InputParameterList& input_parameter_data);
             
             /*!
              * \brief Allocate the global matrices and vectors.
@@ -199,8 +189,7 @@ namespace HappyHeart
              * \param[in] input_parameter_data Object which hold the values of all the parameters defined in
              * the input file.
              */
-            template<class InputParameterDataT>
-            void DefineOperators(const InputParameterDataT& input_parameter_data);
+            void DefineOperators(const InputParameterList& input_parameter_data);
             
             
             /// \name Accessors to vectors and matrices specific to the elastic problem.