From 42f80b8c720806620c96e5c86640580ecbbf17a4 Mon Sep 17 00:00:00 2001 From: Sebastien Gilles <sebastien.gilles@inria.fr> Date: Tue, 27 Oct 2015 12:06:45 +0100 Subject: [PATCH] #721 Elasticity: remove dof source and the template parameters for volumic and surfacic sources. --- HappyHeart.xcodeproj/project.pbxproj | 4 + .../Elasticity/ElasticityModel.hpp | 10 +- .../VariationalFormulationElasticity.cpp | 320 ++++++++++ .../VariationalFormulationElasticity.hpp | 35 +- .../VariationalFormulationElasticity.hxx | 589 ++---------------- 5 files changed, 389 insertions(+), 569 deletions(-) create mode 100644 Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.cpp diff --git a/HappyHeart.xcodeproj/project.pbxproj b/HappyHeart.xcodeproj/project.pbxproj index 947759780a..ab42cf1776 100644 --- a/HappyHeart.xcodeproj/project.pbxproj +++ b/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; diff --git a/Sources/ModelInstances/Elasticity/ElasticityModel.hpp b/Sources/ModelInstances/Elasticity/ElasticityModel.hpp index 1a4dabbe23..84b5b0dde2 100644 --- a/Sources/ModelInstances/Elasticity/ElasticityModel.hpp +++ b/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: diff --git a/Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.cpp b/Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.cpp new file mode 100644 index 0000000000..9b91713c24 --- /dev/null +++ b/Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.cpp @@ -0,0 +1,320 @@ +// +// 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 + diff --git a/Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.hpp b/Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.hpp index 3fa21e0521..acc0d2cd8e 100644 --- a/Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.hpp +++ b/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. diff --git a/Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.hxx b/Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.hxx index 11e943543b..63dc5eebb3 100644 --- a/Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.hxx +++ b/Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.hxx @@ -1,5 +1,5 @@ // -// VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT>.hxx +// VariationalFormulationElasticity.hxx // HappyHeart // // Created by Sebastien Gilles on 01/07/14. @@ -18,375 +18,8 @@ namespace HappyHeart { - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - template<class InputParameterDataT> - void VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> - ::SupplInit(const InputParameterDataT& 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))); - } - } - - - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - template<typename... DofSourcePolicyArgs> - VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> - ::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, - DofSourcePolicyArgs&&... dof_source_policy_args) - : parent(mpi, - time_manager, - god_of_dof, - std::move(boundary_condition_list)), - DofSourcePolicyT(std::forward<DofSourcePolicyArgs>(dof_source_policy_args)...), - 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."); - } - - - - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - void VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT>::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); - } - - - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - void VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT>::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(); - } - - - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - void VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT>::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__); - - DofSourcePolicyT::AddToRhs(rhs); - } - - - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - void VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT>::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); - } - - DofSourcePolicyT::AddToRhs(rhs); - } - - - - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - void VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> - ::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); - } - - - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - void VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT>::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__); - } - } - - - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - void VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT>::UpdateDisplacement() - { - GetNonCstVectorDisplacementPreviousTimeIteration().Copy(parent::GetSystemSolution(GetNumberingSubset()), - __FILE__, __LINE__); - } - - - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - void VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT>::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)); - } - } - } - - - - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - template<class InputParameterDataT> - void VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT>::DefineOperators(const InputParameterDataT& 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); - } - - - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline void VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline void VariationalFormulationElasticity ::UpdateDisplacementAndVelocity() { // Ordering of both calls is capital here! @@ -395,53 +28,33 @@ namespace HappyHeart } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline Wrappers::Petsc::Snes::SNESFunction VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline Wrappers::Petsc::Snes::SNESFunction VariationalFormulationElasticity ::ImplementSnesFunction() const noexcept { return nullptr; } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline Wrappers::Petsc::Snes::SNESJacobian VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline Wrappers::Petsc::Snes::SNESJacobian VariationalFormulationElasticity ::ImplementSnesJacobian() const noexcept { return nullptr; } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > + inline Wrappers::Petsc::Snes::SNESViewer - VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + VariationalFormulationElasticity ::ImplementSnesViewer() const noexcept { return nullptr; } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const GlobalVector& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const GlobalVector& VariationalFormulationElasticity ::GetVectorDisplacementPreviousTimeIteration() const noexcept { assert(!(!vector_displacement_previous_time_iteration_)); @@ -449,26 +62,16 @@ namespace HappyHeart } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline GlobalVector& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline GlobalVector& VariationalFormulationElasticity ::GetNonCstVectorDisplacementPreviousTimeIteration() noexcept { return const_cast<GlobalVector&>(GetVectorDisplacementPreviousTimeIteration()); } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const GlobalVector& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const GlobalVector& VariationalFormulationElasticity ::GetVectorVelocityPreviousTimeIteration() const noexcept { assert(!(!vector_velocity_previous_time_iteration_)); @@ -476,26 +79,16 @@ namespace HappyHeart } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline GlobalVector& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline GlobalVector& VariationalFormulationElasticity ::GetNonCstVectorVelocityPreviousTimeIteration() noexcept { return const_cast<GlobalVector&>(GetVectorVelocityPreviousTimeIteration()); } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const GlobalMatrix& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const GlobalMatrix& VariationalFormulationElasticity ::GetMatrixDisplacementPreviousTimeIteration() const noexcept { assert(!(!matrix_displacement_previous_time_iteration_)); @@ -503,26 +96,16 @@ namespace HappyHeart } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline GlobalMatrix& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline GlobalMatrix& VariationalFormulationElasticity ::GetNonCstMatrixDisplacementPreviousTimeIteration() noexcept { return const_cast<GlobalMatrix&>(GetMatrixDisplacementPreviousTimeIteration()); } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const GlobalMatrix& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const GlobalMatrix& VariationalFormulationElasticity ::GetMatrixVelocityPreviousTimeIteration() const noexcept { assert(!(!matrix_velocity_previous_time_iteration_)); @@ -530,26 +113,16 @@ namespace HappyHeart } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline GlobalMatrix& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline GlobalMatrix& VariationalFormulationElasticity ::GetNonCstMatrixVelocityPreviousTimeIteration() noexcept { return const_cast<GlobalMatrix&>(GetMatrixVelocityPreviousTimeIteration()); } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const GlobalMatrix& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const GlobalMatrix& VariationalFormulationElasticity ::GetMass() const noexcept { assert(!(!mass_)); @@ -557,26 +130,16 @@ namespace HappyHeart } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline GlobalMatrix& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline GlobalMatrix& VariationalFormulationElasticity ::GetNonCstMass() noexcept { return const_cast<GlobalMatrix&>(GetMass()); } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const GlobalMatrix& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const GlobalMatrix& VariationalFormulationElasticity ::GetStiffness() const noexcept { assert(!(!stiffness_)); @@ -584,27 +147,17 @@ namespace HappyHeart } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline GlobalMatrix& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline GlobalMatrix& VariationalFormulationElasticity ::GetNonCstStiffness() noexcept { return const_cast<GlobalMatrix&>(GetStiffness()); } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > + inline const GlobalVariationalOperatorNS::GradOnGradientBasedElasticityTensor& - VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + VariationalFormulationElasticity ::GetStiffnessOperator() const noexcept { assert(!(!stiffness_operator_)); @@ -612,13 +165,8 @@ namespace HappyHeart } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const GlobalVariationalOperatorNS::Mass& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const GlobalVariationalOperatorNS::Mass& VariationalFormulationElasticity ::GetMassOperator() const noexcept { assert(!(!mass_operator_)); @@ -626,26 +174,16 @@ namespace HappyHeart } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const NumberingSubset& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const NumberingSubset& VariationalFormulationElasticity ::GetNumberingSubset() const noexcept { return numbering_subset_; } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const ScalarParameter& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const ScalarParameter& VariationalFormulationElasticity ::GetVolumicMass() const noexcept { assert(!(!volumic_mass_)); @@ -653,13 +191,8 @@ namespace HappyHeart } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const ScalarParameter& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const ScalarParameter& VariationalFormulationElasticity ::GetYoungModulus() const noexcept { assert(!(!young_modulus_)); @@ -667,13 +200,8 @@ namespace HappyHeart } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const ScalarParameter& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const ScalarParameter& VariationalFormulationElasticity ::GetPoissonRatio() const noexcept { assert(!(!poisson_ratio_)); @@ -681,39 +209,24 @@ namespace HappyHeart } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const Unknown& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const Unknown& VariationalFormulationElasticity ::GetSolidDisplacement() const noexcept { return solid_displacement_; } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const FEltSpace& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const FEltSpace& VariationalFormulationElasticity ::GetMainFEltSpace() const noexcept { return main_felt_space_; } - template - < - unsigned int VolumicIndexT, - unsigned int SurfacicIndexT, - class DofSourcePolicyT - > - inline const FEltSpace& VariationalFormulationElasticity<VolumicIndexT, SurfacicIndexT, DofSourcePolicyT> + + inline const FEltSpace& VariationalFormulationElasticity ::GetNeumannFEltSpace() const noexcept { return neumann_felt_space_; -- GitLab