From 269b7cf183f8d0a6115c24edd0f7a311722fb81d Mon Sep 17 00:00:00 2001 From: Alex_Mueller Date: Tue, 6 Feb 2024 13:01:54 +0100 Subject: [PATCH] refactor member variables and template naming (#234) --- docs/website/03_contribution/howToEdit.md | 8 +- docs/website/doxygen/modules.hh | 72 ++-- ikarus/assembler/simpleassemblers.hh | 168 ++++---- ikarus/assembler/simpleassemblers.inl | 165 ++++---- ikarus/controlroutines/adaptivestepsizing.hh | 52 +-- ikarus/controlroutines/controlinfos.hh | 2 +- ikarus/controlroutines/loadcontrol.hh | 17 +- ikarus/controlroutines/loadcontrol.inl | 10 +- ikarus/controlroutines/pathfollowing.hh | 50 ++- ikarus/controlroutines/pathfollowing.inl | 14 +- .../controlroutines/pathfollowingfunctions.hh | 53 ++- ikarus/finiteelements/febases/autodifffe.hh | 65 ++-- ikarus/finiteelements/febases/powerbasisfe.hh | 50 +-- ikarus/finiteelements/febases/scalarfe.hh | 45 ++- ikarus/finiteelements/fehelper.hh | 20 +- ikarus/finiteelements/ferequirements.hh | 74 ++-- ikarus/finiteelements/fetraits.hh | 10 +- ikarus/finiteelements/mechanics/eas/eas2d.hh | 69 ++-- ikarus/finiteelements/mechanics/eas/eas3d.hh | 35 +- .../finiteelements/mechanics/easvariants.hh | 16 +- .../mechanics/enhancedassumedstrains.hh | 97 ++--- .../mechanics/kirchhoffloveshell.hh | 311 ++++++++------- .../finiteelements/mechanics/linearelastic.hh | 144 ++++--- ikarus/finiteelements/mechanics/loads.hh | 6 +- .../mechanics/loads/traction.hh | 97 +++-- .../finiteelements/mechanics/loads/volume.hh | 81 ++-- ikarus/finiteelements/mechanics/materials.hh | 6 +- .../mechanics/materials/interface.hh | 102 ++--- .../mechanics/materials/linearelasticity.hh | 78 ++-- .../mechanics/materials/materials.cpp | 4 +- .../mechanics/materials/neohooke.hh | 77 ++-- .../mechanics/materials/strainconversions.hh | 10 +- .../finiteelements/mechanics/materials/svk.hh | 107 ++--- .../mechanics/materials/tags.hh | 4 +- .../mechanics/materials/vanishingstress.hh | 224 +++++------ .../mechanics/membranestrains.hh | 41 +- .../mechanics/nonlinearelastic.hh | 196 +++++----- ikarus/finiteelements/physicshelper.hh | 175 ++++----- ikarus/io/resultevaluators.hh | 54 +-- ikarus/io/resultfunction.hh | 68 ++-- .../truncatedconjugategradient.hh | 71 ++-- ikarus/solver/linearsolver/linearsolver.cpp | 38 +- ikarus/solver/linearsolver/linearsolver.hh | 31 +- .../solver/nonlinearsolver/newtonraphson.hh | 102 +++-- ...wtonraphsonwithscalarsubsidiaryfunction.hh | 93 +++-- ikarus/solver/nonlinearsolver/trustregion.hh | 365 +++++++++--------- ikarus/utils/algorithms.hh | 206 +++++----- ikarus/utils/autodiffhelper.hh | 26 +- ikarus/utils/basis.hh | 84 ++-- ikarus/utils/concepts.hh | 288 +++++++------- ikarus/utils/defaultfunctions.hh | 44 +-- ikarus/utils/dirichletvalues.hh | 106 ++--- ikarus/utils/drawing/griddrawer.hh | 20 +- ikarus/utils/drawing/matplothelper.hh | 18 +- ikarus/utils/eigendunetransformations.hh | 78 ++-- ikarus/utils/eigensparseaddon.hh | 8 +- ikarus/utils/findlinesegment.hh | 10 +- ikarus/utils/functionhelper.hh | 6 +- ikarus/utils/functionsanitychecks.hh | 58 +-- ikarus/utils/init.hh | 36 +- ikarus/utils/linearalgebrahelper.hh | 318 +++++++-------- ikarus/utils/math.hh | 20 +- ikarus/utils/nonlinearoperator.hh | 166 ++++---- ikarus/utils/observer/controllogger.cpp | 10 +- ikarus/utils/observer/controllogger.hh | 34 +- ikarus/utils/observer/controlvtkwriter.hh | 73 ++-- ikarus/utils/observer/genericobserver.hh | 39 +- .../utils/observer/nonlinearsolverlogger.cpp | 18 +- .../utils/observer/nonlinearsolverlogger.hh | 26 +- ikarus/utils/observer/observer.hh | 177 ++++----- ikarus/utils/observer/observermessages.hh | 4 +- ikarus/utils/polyfit.hh | 10 +- ikarus/utils/pythonautodiffdefinitions.hh | 18 +- ikarus/utils/tensorproductquadrule.hh | 14 +- ikarus/utils/tensorutils.hh | 148 +++---- ikarus/utils/traits.hh | 110 +++--- tests/src/testfeelement.hh | 8 +- 77 files changed, 2858 insertions(+), 2900 deletions(-) diff --git a/docs/website/03_contribution/howToEdit.md b/docs/website/03_contribution/howToEdit.md index a4f466f92..43fff5334 100644 --- a/docs/website/03_contribution/howToEdit.md +++ b/docs/website/03_contribution/howToEdit.md @@ -63,10 +63,10 @@ double complicatedCalculation(double number, double anotherNumber) Look at the Markdown file (`03_contribution/howToEdit.md`) to see how tables, warnings and notes can be inserted. -| Grid Entity Interface || -| :------------ | :-----------: | -| `#!cpp GridViewType leafGridView()` | -| `#!cpp GridViewType levelGridView(int level)` | +| Grid Entity Interface | +|:------------------------------------------------:| +| `#!cpp GridViewType leafGridView()` | +| `#!cpp GridViewType levelGridView(int level)` | !!! warning "Insert a warning" Note that the **four** spaces at the beginning of this line are essential for the warning to be displayed diff --git a/docs/website/doxygen/modules.hh b/docs/website/doxygen/modules.hh index a6dda3744..b57bc1493 100644 --- a/docs/website/doxygen/modules.hh +++ b/docs/website/doxygen/modules.hh @@ -2,76 +2,76 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @defgroup group_main Ikarus Modules - * @brief All Ikarus modules + * \defgroup group_main Ikarus Modules + * \brief All Ikarus modules * * @{ */ -/** @defgroup assembler Assembler -@brief All finite element assembler +/** \defgroup assembler Assembler +\brief All finite element assembler */ -/** @defgroup controlroutines Control Routines -@brief Routines to follow a non-linear solution curve +/** \defgroup controlroutines Control Routines +\brief Routines to follow a non-linear solution curve */ -/** @defgroup finiteelements Finite Elements -@brief Local assemblers for linear algebra contributions to the global system +/** \defgroup finiteelements Finite Elements +\brief Local assemblers for linear algebra contributions to the global system * @{ */ -/** @defgroup mechanics Mechanics -@brief All finite elements in the context of mechanics +/** \defgroup mechanics Mechanics +\brief All finite elements in the context of mechanics */ -/** @defgroup Affordancestags FE Affordance Tags -@brief All finite element affordances +/** \defgroup Affordancestags FE Affordance Tags +\brief All finite element affordances */ -/** @defgroup FEParameterTags FE Parameter Tags -@brief All finite element parameter tags +/** \defgroup FEParameterTags FE Parameter Tags +\brief All finite element parameter tags */ /** @} */ -/** @defgroup materials Materials -@brief Materials for mechanical simulations +/** \defgroup materials Materials +\brief Materials for mechanical simulations * @{*/ -/** @defgroup Materialtags Material Tags -@brief All tags related to strain or stress types +/** \defgroup Materialtags Material Tags +\brief All tags related to strain or stress types */ /** @} */ -/** @defgroup io IO -@brief Functions for input and output +/** \defgroup io IO +\brief Functions for input and output * @{*/ -/** @defgroup resultevaluators Result Evaluator -@brief Computation of finite element results using `calculateAt` +/** \defgroup resultevaluators Result Evaluator +\brief Computation of finite element results using `calculateAt` */ /** @} */ -/** @defgroup pythonbindings Python bindings -@brief Header for Python Bindings +/** \defgroup pythonbindings Python bindings +\brief Header for Python Bindings */ -/** @defgroup solvers Solver -@brief Solvers for algebraic system of equations +/** \defgroup solvers Solver +\brief Solvers for algebraic system of equations */ -/** @defgroup observer Observer -@brief Solvers for algebraic system of equations +/** \defgroup observer Observer +\brief Solvers for algebraic system of equations */ -/** @defgroup utils Utilities -@brief Collection of several utilities +/** \defgroup utils Utilities +\brief Collection of several utilities */ -/** @defgroup tensor Tensor Utilities -@brief Collection of several utilities for dealing with Eigen tensors +/** \defgroup tensor Tensor Utilities +\brief Collection of several utilities for dealing with Eigen tensors */ -/** @defgroup algos Algorithms -@brief Stl-like algorithms for runtime and compile time +/** \defgroup algos Algorithms +\brief Stl-like algorithms for runtime and compile time */ -/** @defgroup traits Type traits -@brief Type traits for TMP +/** \defgroup traits Type traits +\brief Type traits for TMP */ /** @} */ diff --git a/ikarus/assembler/simpleassemblers.hh b/ikarus/assembler/simpleassemblers.hh index 19f390109..4a2fbacd9 100644 --- a/ikarus/assembler/simpleassemblers.hh +++ b/ikarus/assembler/simpleassemblers.hh @@ -24,54 +24,54 @@ namespace Ikarus { * \brief The FlatAssemblerBase takes care of common subtasks done by flat assemblers. * \ingroup assembler * - * \tparam FEContainer_ Type of the finite element container. - * \tparam DirichletValuesType_ Type of the Dirichlet values. + * \tparam FEC Type of the finite element container. + * \tparam DV Type of the Dirichlet values. */ -template +template class FlatAssemblerBase { public: - using FEContainerRaw = std::remove_cvref_t; ///< Type of the raw finite element container. + using FEContainerRaw = std::remove_cvref_t; ///< Type of the raw finite element container. using FERequirementType = typename FEContainerRaw::value_type::FERequirementType; ///< Type of the finite element requirement. using GlobalIndex = typename FEContainerRaw::value_type::GlobalIndex; ///< Type of the global index. - using Basis = typename DirichletValuesType_::Basis; ///< Type of the basis. + using Basis = typename DV::Basis; ///< Type of the basis. using GridView = typename Basis::GridView; ///< Type of the grid view. - using FEContainer = FEContainer_; ///< Type of the finite element container. + using FEContainer = FEC; ///< Type of the finite element container. using FEContainerType = std::conditional_t, const FEContainer, FEContainer>; ///< Type of the finite element container (reference or by value). - using DirichletValuesType = DirichletValuesType_; ///< Type of the Dirichlet values. + using DirichletValuesType = DV; ///< Type of the Dirichlet values. /** * \brief Constructor for FlatAssemblerBase. * * \param fes Finite element container. - * \param p_dirichletValues Reference to Dirichlet values. + * \param dirichletValues Reference to Dirichlet values. */ - FlatAssemblerBase(FEContainer&& fes, const DirichletValuesType& p_dirichletValues) - : feContainer{std::forward(fes)}, - dirichletValues{&p_dirichletValues} { - constraintsBelow_.reserve(dirichletValues->size()); + FlatAssemblerBase(FEContainer&& fes, const DirichletValuesType& dirichletValues) + : feContainer_{std::forward(fes)}, + dirichletValues_{&dirichletValues} { + constraintsBelow_.reserve(dirichletValues_->size()); size_t counter = 0; - for (auto iv : std::ranges::iota_view{decltype(dirichletValues->size())(0), dirichletValues->size()}) { + for (auto iv : std::ranges::iota_view{decltype(dirichletValues_->size())(0), dirichletValues_->size()}) { constraintsBelow_.emplace_back(counter); - if (dirichletValues->isConstrained(iv)) + if (dirichletValues_->isConstrained(iv)) ++counter; } - fixedDofs = dirichletValues->fixedDOFsize(); + fixedDofs_ = dirichletValues_->fixedDOFsize(); } /** * \brief Returns the size of the free degrees of freedom, which are not fixed by a Dirichlet boundary condition. * \return Size of the reduced degrees of freedom. */ - size_t reducedSize() { return dirichletValues->size() - fixedDofs; } + size_t reducedSize() { return dirichletValues_->size() - fixedDofs_; } /** * \brief Returns the size of nodes, i.e., the number of degrees of freedom. * \return Size of the degrees of freedom. */ - size_t size() { return dirichletValues->size(); } + size_t size() { return dirichletValues_->size(); } /** * \brief Creates the full-sized vector of size Dof and inserts the values of a reduced vector at the "free" @@ -86,7 +86,7 @@ public: * \brief Returns the container of finite elements. * \return Reference to the finite element container. */ - auto& finiteElements() const { return feContainer; } + auto& finiteElements() const { return feContainer_; } /** * \brief Returns the number of constraints below a given degrees of freedom index. @@ -102,7 +102,7 @@ public: * \param i Index of the degree of freedom. * \return True if the degree of freedom is fixed; false otherwise. */ - [[nodiscard]] bool isConstrained(size_t i) const { return dirichletValues->isConstrained(i); } + [[nodiscard]] bool isConstrained(size_t i) const { return dirichletValues_->isConstrained(i); } /** * \brief Coarse estimate of node connectivity, i.e., this relates to the bandwidth of a sparse matrix. @@ -111,33 +111,33 @@ public: * \return Size_t Coarse estimate of node connectivity. */ [[nodiscard]] size_t estimateOfConnectivity() const { - return dirichletValues->basis().gridView().size(GridView::dimension) * 8; + return dirichletValues_->basis().gridView().size(GridView::dimension) * 8; } private: - FEContainerType feContainer; - const DirichletValuesType* dirichletValues; + FEContainerType feContainer_; + const DirichletValuesType* dirichletValues_; std::vector constraintsBelow_{}; - size_t fixedDofs{}; + size_t fixedDofs_{}; }; #ifndef DOXYGEN template -FlatAssemblerBase(T&& fes, const DirichletValuesType& dirichletValues_) -> FlatAssemblerBase; +FlatAssemblerBase(T&& fes, const DirichletValuesType& dirichletValues) -> FlatAssemblerBase; #endif /** * \class ScalarAssembler * \brief ScalarAssembler assembles scalar quantities. * \ingroup assembler - * \tparam FEContainer_ Type of the finite element container. - * \tparam DirichletValuesType_ Type of the Dirichlet values. + * \tparam FEC Type of the finite element container. + * \tparam DV Type of the Dirichlet values. */ -template -class ScalarAssembler : public FlatAssemblerBase +template +class ScalarAssembler : public FlatAssemblerBase { - using FEContainerRaw = std::remove_cvref_t; ///< Type of the raw finite element container. - using Base = FlatAssemblerBase; + using FEContainerRaw = std::remove_cvref_t; ///< Type of the raw finite element container. + using Base = FlatAssemblerBase; public: using typename Base::Basis; @@ -150,10 +150,10 @@ public: * \brief Constructor for ScalarAssembler. * * \param fes Finite element container. - * \param dirichletValues_ Reference to Dirichlet values. + * \param dirichletValues Reference to Dirichlet values. */ - ScalarAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues_) - : FlatAssemblerBase(std::forward(fes), dirichletValues_) {} + ScalarAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues) + : FlatAssemblerBase(std::forward(fes), dirichletValues) {} /** * \brief Calculates the scalar quantity requested by feRequirements and returns a reference. @@ -171,33 +171,33 @@ private: * \return Reference to the calculated scalar quantity. */ double& getScalarImpl(const FERequirementType& feRequirements) { - scal = 0.0; + scal_ = 0.0; for (auto& fe : this->finiteElements()) { - scal += fe.calculateScalar(feRequirements); + scal_ += fe.calculateScalar(feRequirements); } - return scal; + return scal_; } - double scal{0.0}; + double scal_{0.0}; }; #ifndef DOXYGEN template -ScalarAssembler(T&& fes, const DirichletValuesType& dirichletValues_) -> ScalarAssembler; +ScalarAssembler(T&& fes, const DirichletValuesType& dirichletValues) -> ScalarAssembler; #endif /** * \class VectorFlatAssembler * \brief VectorFlatAssembler assembles vector quantities using a flat basis Indexing strategy. *\ingroup assembler - * \tparam FEContainer_ Type of the finite element container. - * \tparam DirichletValuesType_ Type of the Dirichlet values. + * \tparam FEC Type of the finite element container. + * \tparam DV Type of the Dirichlet values. */ -template -class VectorFlatAssembler : public ScalarAssembler +template +class VectorFlatAssembler : public ScalarAssembler { - using FEContainerRaw = std::remove_cvref_t; ///< Type of the raw finite element container. - using Base = ScalarAssembler; + using FEContainerRaw = std::remove_cvref_t; ///< Type of the raw finite element container. + using Base = ScalarAssembler; public: using typename Base::Basis; @@ -211,10 +211,10 @@ public: * \brief Constructor for VectorFlatAssembler. * * \param fes Finite element container. - * \param dirichletValues_ Reference to Dirichlet values. + * \param dirichletValues Reference to Dirichlet values. */ - VectorFlatAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues_) - : ScalarAssembler(std::forward(fes), dirichletValues_) {} + VectorFlatAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues) + : ScalarAssembler(std::forward(fes), dirichletValues) {} /** * \brief Calculates the vectorial quantity requested by feRequirements and returns a reference. @@ -253,15 +253,14 @@ private: Eigen::VectorXd& getReducedVectorImpl(const FERequirementType& feRequirements); - Eigen::VectorXd vecRaw{}; ///< Raw vector without changes for dirichlet degrees of freedom - Eigen::VectorXd vec{}; ///< Vector quantity. - Eigen::VectorXd vecRed{}; ///< Reduced vector quantity. + Eigen::VectorXd vecRaw_{}; ///< Raw vector without changes for dirichlet degrees of freedom + Eigen::VectorXd vec_{}; ///< Vector quantity. + Eigen::VectorXd vecRed_{}; ///< Reduced vector quantity. }; #ifndef DOXYGEN template -VectorFlatAssembler(T&& fes, const DirichletValuesType& dirichletValues_) - -> VectorFlatAssembler; +VectorFlatAssembler(T&& fes, const DirichletValuesType& dirichletValues) -> VectorFlatAssembler; #endif /** @@ -269,15 +268,15 @@ VectorFlatAssembler(T&& fes, const DirichletValuesType& dirichletValues_) * \brief SparseFlatAssembler assembles matrix quantities using a flat basis Indexing strategy. * The matrix is stored in a sparse matrix format. This format is exploited during the assembly process. * \ingroup assembler - * \tparam FEContainer_ Type of the finite element container. - * \tparam DirichletValuesType_ Type of the Dirichlet values. + * \tparam FEC Type of the finite element container. + * \tparam DV Type of the Dirichlet values. */ -template -class SparseFlatAssembler : public VectorFlatAssembler +template +class SparseFlatAssembler : public VectorFlatAssembler { public: - using FEContainerRaw = std::remove_cvref_t; ///< Type of the raw finite element container. - using Base = VectorFlatAssembler; + using FEContainerRaw = std::remove_cvref_t; ///< Type of the raw finite element container. + using Base = VectorFlatAssembler; using typename Base::Basis; using typename Base::DirichletValuesType; @@ -289,10 +288,10 @@ public: * \brief Constructor for SparseFlatAssembler. * * \param fes Finite element container. - * \param dirichletValues_ Reference to Dirichlet values. + * \param dirichletValues Reference to Dirichlet values. */ - SparseFlatAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues_) - : VectorFlatAssembler(std::forward(fes), dirichletValues_) {} + SparseFlatAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues) + : VectorFlatAssembler(std::forward(fes), dirichletValues) {} using GridView = typename Basis::GridView; ///< Type of the grid view. @@ -355,21 +354,20 @@ private: /** Pre-processes the reduced sparse matrix before assembly. */ void preProcessSparseMatrixReduced(Eigen::SparseMatrix& assemblyMat); - Eigen::SparseMatrix spMatRaw; ///< Raw sparse matrix without changes for dirichlet degrees of freedom. - Eigen::SparseMatrix spMat; ///< Sparse matrix. - Eigen::SparseMatrix spMatReduced; ///< Reduced sparse matrix. + Eigen::SparseMatrix spMatRaw_; ///< Raw sparse matrix without changes for dirichlet degrees of freedom. + Eigen::SparseMatrix spMat_; ///< Sparse matrix. + Eigen::SparseMatrix spMatReduced_; ///< Reduced sparse matrix. std::vector> - elementLinearIndices; ///< Vector storing indices of matrix entries in linear storage + elementLinearIndices_; ///< Vector storing indices of matrix entries in linear storage std::vector> - elementLinearReducedIndices; ///< Vector storing indices of matrix entries in linear storage - std::once_flag sparsePreProcessorRaw, sparsePreProcessor, - sparsePreProcessorReduced; ///< flags that store if the sparsity pattern construction happened + elementLinearReducedIndices_; ///< Vector storing indices of matrix entries in linear storage + std::once_flag sparsePreProcessorRaw_, sparsePreProcessor_, + sparsePreProcessorReduced_; ///< flags that store if the sparsity pattern construction happened }; #ifndef DOXYGEN -template -SparseFlatAssembler(T&& fes, const DirichletValuesType& dirichletValues_) - -> SparseFlatAssembler; +template +SparseFlatAssembler(FEC&& fes, const DV& dirichletValues) -> SparseFlatAssembler; #endif /** @@ -377,16 +375,16 @@ SparseFlatAssembler(T&& fes, const DirichletValuesType& dirichletValues_) * \brief DenseFlatAssembler assembles matrix quantities using a flat basis Indexing strategy. * The matrix is stored in a dense matrix format. This format is exploited during the assembly process. * \ingroup assembler - * \tparam FEContainer_ Type of the finite element container. - * \tparam DirichletValuesType_ Type of the Dirichlet values. + * \tparam FEC Type of the finite element container. + * \tparam DV Type of the Dirichlet values. * \note Requires Ikarus::Concepts::FlatIndexBasis. */ -template -class DenseFlatAssembler : public VectorFlatAssembler +template +class DenseFlatAssembler : public VectorFlatAssembler { public: - using FEContainerRaw = std::remove_cvref_t; ///< Type of the raw finite element container. - using Base = VectorFlatAssembler; ///< Type alias for the base class. + using FEContainerRaw = std::remove_cvref_t; ///< Type of the raw finite element container. + using Base = VectorFlatAssembler; ///< Type alias for the base class. using typename Base::Basis; ///< Type of the basis. using typename Base::DirichletValuesType; ///< Type of the Dirichlet values. @@ -398,10 +396,10 @@ public: * \brief Constructor for DenseFlatAssembler. * * \param fes Finite element container. - * \param dirichletValues_ Reference to Dirichlet values. + * \param dirichletValues Reference to Dirichlet values. */ - explicit DenseFlatAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues_) - : VectorFlatAssembler(std::forward(fes), dirichletValues_) {} + explicit DenseFlatAssembler(FEContainer&& fes, const DirichletValuesType& dirichletValues) + : VectorFlatAssembler(std::forward(fes), dirichletValues) {} /** * \brief Calculates the matrix quantity requested by feRequirements and returns a reference. @@ -439,15 +437,15 @@ private: Eigen::MatrixXd& getMatrixImpl(const FERequirementType& feRequirements); Eigen::MatrixXd& getReducedMatrixImpl(const FERequirementType& feRequirements); - Eigen::MatrixXd matRaw{}; ///< Raw dense matrix for assembly. - Eigen::MatrixXd mat{}; ///< Dense matrix quantity. - Eigen::MatrixXd matRed{}; ///< Reduced dense matrix quantity. + Eigen::MatrixXd matRaw_{}; ///< Raw dense matrix for assembly. + Eigen::MatrixXd mat_{}; ///< Dense matrix quantity. + Eigen::MatrixXd matRed_{}; ///< Reduced dense matrix quantity. }; #ifndef DOXYGEN // https://en.cppreference.com/w/cpp/language/class_template_argument_deduction -template -DenseFlatAssembler(T&& fes, const DirichletValuesType& dirichletValues_) -> DenseFlatAssembler; +template +DenseFlatAssembler(FEC&& fes, const DV& dirichletValues) -> DenseFlatAssembler; #endif } // namespace Ikarus diff --git a/ikarus/assembler/simpleassemblers.inl b/ikarus/assembler/simpleassemblers.inl index f3cb20962..cd8dd3d90 100644 --- a/ikarus/assembler/simpleassemblers.inl +++ b/ikarus/assembler/simpleassemblers.inl @@ -11,15 +11,14 @@ namespace Ikarus { -template -Eigen::VectorXd FlatAssemblerBase::createFullVector( - Eigen::Ref reducedVector) { +template +Eigen::VectorXd FlatAssemblerBase::createFullVector(Eigen::Ref reducedVector) { assert(reducedVector.size() == static_cast(this->reducedSize()) && "The reduced vector you passed has the wrong dimensions."); Eigen::Index reducedCounter = 0; Eigen::VectorXd fullVec(this->size()); for (Eigen::Index i = 0; i < fullVec.size(); ++i) { - if (dirichletValues->isConstrained(i)) { + if (dirichletValues_->isConstrained(i)) { ++reducedCounter; fullVec[i] = 0.0; continue; @@ -29,9 +28,9 @@ Eigen::VectorXd FlatAssemblerBase::createFullVector( return fullVec; } -template -void VectorFlatAssembler::assembleRawVectorImpl(const FERequirementType& feRequirements, - Eigen::VectorXd& assemblyVec) { +template +void VectorFlatAssembler::assembleRawVectorImpl(const FERequirementType& feRequirements, + Eigen::VectorXd& assemblyVec) { assemblyVec.setZero(this->size()); Eigen::VectorXd vecLocal; std::vector dofs; @@ -47,25 +46,24 @@ void VectorFlatAssembler::assembleRawVectorImpl(const FERequ } } -template -Eigen::VectorXd& VectorFlatAssembler::getRawVectorImpl(const FERequirementType& feRequirements) { - assembleRawVectorImpl(feRequirements, vecRaw); - return vecRaw; +template +Eigen::VectorXd& VectorFlatAssembler::getRawVectorImpl(const FERequirementType& feRequirements) { + assembleRawVectorImpl(feRequirements, vecRaw_); + return vecRaw_; } -template -Eigen::VectorXd& VectorFlatAssembler::getVectorImpl(const FERequirementType& feRequirements) { - assembleRawVectorImpl(feRequirements, vec); +template +Eigen::VectorXd& VectorFlatAssembler::getVectorImpl(const FERequirementType& feRequirements) { + assembleRawVectorImpl(feRequirements, vec_); for (auto i = 0U; i < this->size(); ++i) if (this->isConstrained(i)) - vec[i] = 0; - return vec; + vec_[i] = 0; + return vec_; } -template -Eigen::VectorXd& VectorFlatAssembler::getReducedVectorImpl( - const FERequirementType& feRequirements) { - vecRed.setZero(this->reducedSize()); +template +Eigen::VectorXd& VectorFlatAssembler::getReducedVectorImpl(const FERequirementType& feRequirements) { + vecRed_.setZero(this->reducedSize()); Eigen::VectorXd vecLocal; std::vector dofs; for (auto& fe : this->finiteElements()) { @@ -79,60 +77,58 @@ Eigen::VectorXd& VectorFlatAssembler::getReducedVectorImpl( ++localConstrainedCounter; continue; } else - vecRed(dofIndex[0] - this->constraintsBelow(dofIndex[0])) += vecLocal[localConstrainedCounter++]; + vecRed_(dofIndex[0] - this->constraintsBelow(dofIndex[0])) += vecLocal[localConstrainedCounter++]; } } - return vecRed; + return vecRed_; } -template -void SparseFlatAssembler::assembleRawMatrixImpl(const FERequirementType& feRequirements, - Eigen::SparseMatrix& assemblyMat) { +template +void SparseFlatAssembler::assembleRawMatrixImpl(const FERequirementType& feRequirements, + Eigen::SparseMatrix& assemblyMat) { assemblyMat.coeffs().setZero(); Eigen::MatrixXd A; for (size_t elementIndex = 0; const auto& fe : this->finiteElements()) { A.setZero(fe.size(), fe.size()); fe.calculateMatrix(feRequirements, A); - assert(static_cast(std::sqrt(elementLinearIndices[elementIndex].size())) == A.rows() && + assert(static_cast(std::sqrt(elementLinearIndices_[elementIndex].size())) == A.rows() && "The returned matrix has wrong rowSize!"); - assert(static_cast(std::sqrt(elementLinearIndices[elementIndex].size())) == A.cols() && + assert(static_cast(std::sqrt(elementLinearIndices_[elementIndex].size())) == A.cols() && "The returned matrix has wrong colSize!"); for (Eigen::Index linearIndex = 0; double matrixEntry : A.reshaped()) - assemblyMat.coeffs()(elementLinearIndices[elementIndex][linearIndex++]) += matrixEntry; + assemblyMat.coeffs()(elementLinearIndices_[elementIndex][linearIndex++]) += matrixEntry; ++elementIndex; } } -template -Eigen::SparseMatrix& SparseFlatAssembler::getRawMatrixImpl( - const FERequirementType& feRequirements) { - std::call_once(sparsePreProcessorRaw, [&]() { preProcessSparseMatrix(spMatRaw); }); - assembleRawMatrixImpl(feRequirements, spMatRaw); - return spMatRaw; +template +Eigen::SparseMatrix& SparseFlatAssembler::getRawMatrixImpl(const FERequirementType& feRequirements) { + std::call_once(sparsePreProcessorRaw_, [&]() { preProcessSparseMatrix(spMatRaw_); }); + assembleRawMatrixImpl(feRequirements, spMatRaw_); + return spMatRaw_; } -template -Eigen::SparseMatrix& SparseFlatAssembler::getMatrixImpl( - const FERequirementType& feRequirements) { - std::call_once(sparsePreProcessor, [&]() { preProcessSparseMatrix(spMat); }); - assembleRawMatrixImpl(feRequirements, spMat); +template +Eigen::SparseMatrix& SparseFlatAssembler::getMatrixImpl(const FERequirementType& feRequirements) { + std::call_once(sparsePreProcessor_, [&]() { preProcessSparseMatrix(spMat_); }); + assembleRawMatrixImpl(feRequirements, spMat_); for (auto i = 0U; i < this->size(); ++i) if (this->isConstrained(i)) - spMat.col(i) *= 0; + spMat_.col(i) *= 0; for (auto i = 0U; i < this->size(); ++i) if (this->isConstrained(i)) - spMat.row(i) *= 0; + spMat_.row(i) *= 0; for (auto i = 0U; i < this->size(); ++i) if (this->isConstrained(i)) - spMat.diagonal()[i] = 1; - return spMat; + spMat_.diagonal()[i] = 1; + return spMat_; } -template -Eigen::SparseMatrix& SparseFlatAssembler::getReducedMatrixImpl( +template +Eigen::SparseMatrix& SparseFlatAssembler::getReducedMatrixImpl( const FERequirementType& feRequirements) { - std::call_once(sparsePreProcessorReduced, [&]() { preProcessSparseMatrixReduced(spMatReduced); }); - spMatReduced.coeffs().setZero(); + std::call_once(sparsePreProcessorReduced_, [&]() { preProcessSparseMatrixReduced(spMatReduced_); }); + spMatReduced_.coeffs().setZero(); Eigen::MatrixXd A; std::vector dofs; for (size_t elementIndex = 0; const auto& fe : this->finiteElements()) { @@ -150,17 +146,17 @@ Eigen::SparseMatrix& SparseFlatAssembler::getReduced for (auto c = 0U; c < dofs.size(); ++c) { if (this->isConstrained(dofs[c][0])) continue; - spMatReduced.coeffs()(elementLinearReducedIndices[elementIndex][linearIndex++]) += A(r, c); + spMatReduced_.coeffs()(elementLinearReducedIndices_[elementIndex][linearIndex++]) += A(r, c); } } } ++elementIndex; } - return spMatReduced; + return spMatReduced_; } -template -void SparseFlatAssembler::createOccupationPattern(Eigen::SparseMatrix& assemblyMat) { +template +void SparseFlatAssembler::createOccupationPattern(Eigen::SparseMatrix& assemblyMat) { assemblyMat.resize(this->size(), this->size()); std::vector> vectorOfTriples; @@ -176,8 +172,8 @@ void SparseFlatAssembler::createOccupationPattern(Eigen::Spa assemblyMat.setFromTriplets(vectorOfTriples.begin(), vectorOfTriples.end()); } -template -void SparseFlatAssembler::createReducedOccupationPattern(Eigen::SparseMatrix& assemblyMat) { +template +void SparseFlatAssembler::createReducedOccupationPattern(Eigen::SparseMatrix& assemblyMat) { assemblyMat.resize(this->reducedSize(), this->reducedSize()); std::vector> vectorOfTriples; using std::size; @@ -203,55 +199,54 @@ void SparseFlatAssembler::createReducedOccupationPattern(Eig assemblyMat.setFromTriplets(vectorOfTriples.begin(), vectorOfTriples.end()); } -template -void SparseFlatAssembler::createLinearDOFsPerElement(Eigen::SparseMatrix& assemblyMat) { +template +void SparseFlatAssembler::createLinearDOFsPerElement(Eigen::SparseMatrix& assemblyMat) { std::vector dofs; for (auto&& fe : this->finiteElements()) { dofs.resize(0); fe.globalFlatIndices(dofs); - elementLinearIndices.emplace_back(Dune::power(dofs.size(), 2)); + elementLinearIndices_.emplace_back(Dune::power(dofs.size(), 2)); for (Eigen::Index linearIndexOfElement = 0; auto&& c : dofs) for (auto&& r : dofs) - elementLinearIndices.back()[linearIndexOfElement++] = assemblyMat.getLinearIndex(r[0], c[0]); + elementLinearIndices_.back()[linearIndexOfElement++] = assemblyMat.getLinearIndex(r[0], c[0]); } } -template -void SparseFlatAssembler::createLinearDOFsPerElementReduced( - Eigen::SparseMatrix& assemblyMat) { +template +void SparseFlatAssembler::createLinearDOFsPerElementReduced(Eigen::SparseMatrix& assemblyMat) { std::vector dofs; for (auto&& fe : this->finiteElements()) { dofs.resize(0); fe.globalFlatIndices(dofs); - elementLinearReducedIndices.emplace_back(); + elementLinearReducedIndices_.emplace_back(); for (auto r = 0U; r < dofs.size(); ++r) { if (this->isConstrained(dofs[r][0])) continue; for (auto c = 0U; c < dofs.size(); ++c) { if (this->isConstrained(dofs[c][0])) continue; - elementLinearReducedIndices.back().push_back(assemblyMat.getLinearIndex( + elementLinearReducedIndices_.back().push_back(assemblyMat.getLinearIndex( dofs[r][0] - this->constraintsBelow(dofs[r][0]), dofs[c][0] - this->constraintsBelow(dofs[c][0]))); } } } } -template -void SparseFlatAssembler::preProcessSparseMatrix(Eigen::SparseMatrix& assemblyMat) { +template +void SparseFlatAssembler::preProcessSparseMatrix(Eigen::SparseMatrix& assemblyMat) { createOccupationPattern(assemblyMat); createLinearDOFsPerElement(assemblyMat); } -template -void SparseFlatAssembler::preProcessSparseMatrixReduced(Eigen::SparseMatrix& assemblyMat) { +template +void SparseFlatAssembler::preProcessSparseMatrixReduced(Eigen::SparseMatrix& assemblyMat) { createReducedOccupationPattern(assemblyMat); createLinearDOFsPerElementReduced(assemblyMat); } -template -void DenseFlatAssembler::assembleRawMatrixImpl(const FERequirementType& feRequirements, - Eigen::MatrixXd& assemblyMat) { +template +void DenseFlatAssembler::assembleRawMatrixImpl(const FERequirementType& feRequirements, + Eigen::MatrixXd& assemblyMat) { assemblyMat.setZero(this->size(), this->size()); Eigen::MatrixXd matLocal; std::vector dofs; @@ -270,30 +265,30 @@ void DenseFlatAssembler::assembleRawMatrixImpl(const FERequi } } -template -Eigen::MatrixXd& DenseFlatAssembler::getRawMatrixImpl(const FERequirementType& feRequirements) { - assembleRawMatrixImpl(feRequirements, matRaw); - return matRaw; +template +Eigen::MatrixXd& DenseFlatAssembler::getRawMatrixImpl(const FERequirementType& feRequirements) { + assembleRawMatrixImpl(feRequirements, matRaw_); + return matRaw_; } -template -Eigen::MatrixXd& DenseFlatAssembler::getMatrixImpl(const FERequirementType& feRequirements) { - assembleRawMatrixImpl(feRequirements, mat); +template +Eigen::MatrixXd& DenseFlatAssembler::getMatrixImpl(const FERequirementType& feRequirements) { + assembleRawMatrixImpl(feRequirements, mat_); for (auto i = 0U; i < this->size(); ++i) if (this->isConstrained(i)) - mat.col(i).setZero(); + mat_.col(i).setZero(); for (auto i = 0U; i < this->size(); ++i) if (this->isConstrained(i)) - mat.row(i).setZero(); + mat_.row(i).setZero(); for (auto i = 0U; i < this->size(); ++i) if (this->isConstrained(i)) - mat(i, i) = 1; - return mat; + mat_(i, i) = 1; + return mat_; } -template -Eigen::MatrixXd& DenseFlatAssembler::getReducedMatrixImpl(const FERequirementType& feRequirements) { - matRed.setZero(this->reducedSize(), this->reducedSize()); +template +Eigen::MatrixXd& DenseFlatAssembler::getReducedMatrixImpl(const FERequirementType& feRequirements) { + matRed_.setZero(this->reducedSize(), this->reducedSize()); Eigen::MatrixXd matLocal; std::vector dofs; for (auto& fe : this->finiteElements()) { @@ -311,12 +306,12 @@ Eigen::MatrixXd& DenseFlatAssembler::getReducedMatrixImpl(co if (this->isConstrained(dofs[c][0])) { continue; } - matRed(dofs[r][0] - this->constraintsBelow(dofs[r][0]), dofs[c][0] - this->constraintsBelow(dofs[c][0])) += + matRed_(dofs[r][0] - this->constraintsBelow(dofs[r][0]), dofs[c][0] - this->constraintsBelow(dofs[c][0])) += matLocal(r, c); } } } } - return matRed; + return matRed_; } } // namespace Ikarus diff --git a/ikarus/controlroutines/adaptivestepsizing.hh b/ikarus/controlroutines/adaptivestepsizing.hh index 887685eb5..d7f9ab49d 100644 --- a/ikarus/controlroutines/adaptivestepsizing.hh +++ b/ikarus/controlroutines/adaptivestepsizing.hh @@ -2,11 +2,11 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file adaptivestepsizing.hh - * @brief Contains the AdaptiveStepSizing namespace with strategies for adaptive step sizing. + * \file adaptivestepsizing.hh + * \brief Contains the AdaptiveStepSizing namespace with strategies for adaptive step sizing. * * This file defines two strategies for adaptive step sizing: NoOp and IterationBased. - * @ingroup controlroutines + * \ingroup controlroutines */ #pragma once @@ -17,40 +17,41 @@ namespace Ikarus::AdaptiveStepSizing { /** - * @brief The NoOp strategy for adaptive step sizing. + * \brief The NoOp strategy for adaptive step sizing. * * This strategy does nothing and keeps the step size constant. */ struct NoOp { /** - * @brief Call operator for the NoOp strategy. + * \brief Call operator for the NoOp strategy. * * Does nothing. * - * @param solverInfo Information about the nonlinear solver. - * @param subsidiaryArgs Subsidiary arguments for adaptive step sizing. - * @param nonLinearOperator The nonlinear operator. + * \param solverInfo Information about the nonlinear solver. + * \param subsidiaryArgs Subsidiary arguments for adaptive step sizing. + * \param nonLinearOperator The nonlinear operator. + * \tparam NLO The nonlinear operator type. */ - template + template void operator()(const NonLinearSolverInformation& solverInfo, SubsidiaryArgs& subsidiaryArgs, - const NonLinearOperator& nonLinearOperator) {} + const NLO& nonLinearOperator) {} /** - * @brief Get the target iterations. - * @return Always returns 0, indicating no specific target. + * \brief Get the target iterations. + * \return Always returns 0, indicating no specific target. */ [[nodiscard]] int targetIterations() const { return 0; } /** - * @brief Set the target iterations. - * @param targetIterations The target iterations (ignored). + * \brief Set the target iterations. + * \param targetIterations The target iterations (ignored). */ void setTargetIterations([[maybe_unused]] int targetIterations) {} }; /** - * @brief The IterationBased strategy for adaptive step sizing. + * \brief The IterationBased strategy for adaptive step sizing. * * This strategy adjusts the step size based on the ratio of target iterations to previous iterations. * @@ -62,17 +63,18 @@ struct NoOp struct IterationBased { /** - * @brief Call operator of the IterationBased strategy. + * \brief Call operator of the IterationBased strategy. * * Adjusts the step size based on the ratio of target iterations to previous iterations. * - * @param solverInfo Information about the nonlinear solver. - * @param subsidiaryArgs Subsidiary arguments for adaptive step sizing. - * @tparam NonLinearOperator The nonlinear operator. + * \param solverInfo Information about the nonlinear solver. + * \param subsidiaryArgs Subsidiary arguments for adaptive step sizing. + * \param nonLinearOperator The nonlinear operator. + * \tparam NLO The nonlinear operator. */ - template + template void operator()(const NonLinearSolverInformation& solverInfo, SubsidiaryArgs& subsidiaryArgs, - const NonLinearOperator&) { + const NLO& nonLinearOperator) { if (subsidiaryArgs.currentStep == 0) return; if (targetIterations_ == 0) @@ -86,14 +88,14 @@ struct IterationBased } /** - * @brief Get the target iterations. - * @return The number of target iterations. + * \brief Get the target iterations. + * \return The number of target iterations. */ [[nodiscard]] int targetIterations() const { return targetIterations_; } /** - * @brief Set the target iterations. - * @param targetIterations The number of target iterations. + * \brief Set the target iterations. + * \param targetIterations The number of target iterations. */ void setTargetIterations(int targetIterations) { targetIterations_ = targetIterations; } diff --git a/ikarus/controlroutines/controlinfos.hh b/ikarus/controlroutines/controlinfos.hh index 4ce81ec78..d7167bf06 100644 --- a/ikarus/controlroutines/controlinfos.hh +++ b/ikarus/controlroutines/controlinfos.hh @@ -4,7 +4,7 @@ /** * \file * \brief Defines the ControlInformation structure for storing control results. - * @ingroup controlroutines + * \ingroup controlroutines */ #pragma once diff --git a/ikarus/controlroutines/loadcontrol.hh b/ikarus/controlroutines/loadcontrol.hh index fc8f542d6..ed7d82a66 100644 --- a/ikarus/controlroutines/loadcontrol.hh +++ b/ikarus/controlroutines/loadcontrol.hh @@ -24,9 +24,9 @@ namespace Ikarus { * This class represents the LoadControl control routine. It increments the last parameter of a nonlinear operator * and utilizes a nonlinear solver, such as Newton's method, to solve the resulting system at each step. * - * \tparam NonLinearSolver Type of the nonlinear solver used in the control routine. + * \tparam NLS Type of the nonlinear solver used in the control routine. */ -template +template class LoadControl : public IObservable { public: @@ -36,21 +36,20 @@ public: /** * \brief Constructor for LoadControl. * - * \param nonLinearSolver_ Shared pointer to the nonlinear solver. + * \param nonLinearSolver Shared pointer to the nonlinear solver. * \param loadSteps Number of load steps in the control routine. * \param tbeginEnd Array representing the range of load parameters [tbegin, tend]. */ - LoadControl(const std::shared_ptr& nonLinearSolver_, int loadSteps, - const std::array& tbeginEnd) - : nonLinearSolver{nonLinearSolver_}, + LoadControl(const std::shared_ptr& nonLinearSolver, int loadSteps, const std::array& tbeginEnd) + : nonLinearSolver_{nonLinearSolver}, loadSteps_{loadSteps}, parameterBegin_{tbeginEnd[0]}, parameterEnd_{tbeginEnd[1]}, stepSize_{(parameterEnd_ - parameterBegin_) / loadSteps_} { static_assert( requires { - nonLinearSolver->nonLinearOperator().lastParameter() = 0.0; - nonLinearSolver->nonLinearOperator().lastParameter() += 0.0; + nonLinearSolver_->nonLinearOperator().lastParameter() = 0.0; + nonLinearSolver_->nonLinearOperator().lastParameter() += 0.0; }, "The last parameter (load factor) must be assignable and incrementable with a double!"); } @@ -62,7 +61,7 @@ public: ControlInformation run(); private: - std::shared_ptr nonLinearSolver; + std::shared_ptr nonLinearSolver_; int loadSteps_; double parameterBegin_; double parameterEnd_; diff --git a/ikarus/controlroutines/loadcontrol.inl b/ikarus/controlroutines/loadcontrol.inl index 18802dcd3..63aa23658 100644 --- a/ikarus/controlroutines/loadcontrol.inl +++ b/ikarus/controlroutines/loadcontrol.inl @@ -10,16 +10,16 @@ #pragma once namespace Ikarus { -template -ControlInformation LoadControl::run() { +template +ControlInformation LoadControl::run() { ControlInformation info({false}); - auto& nonOp = nonLinearSolver->nonLinearOperator(); + auto& nonOp = nonLinearSolver_->nonLinearOperator(); this->notify(ControlMessages::CONTROL_STARTED, static_cast(this->name())); auto& loadParameter = nonOp.lastParameter(); loadParameter = 0.0; this->notify(ControlMessages::STEP_STARTED, 0, stepSize_); - auto solverInfo = nonLinearSolver->solve(); + auto solverInfo = nonLinearSolver_->solve(); info.solverInfos.push_back(solverInfo); info.totalIterations += solverInfo.iterations; if (not solverInfo.success) @@ -30,7 +30,7 @@ ControlInformation LoadControl::run() { for (int ls = 0; ls < loadSteps_; ++ls) { this->notify(ControlMessages::STEP_STARTED, ls, stepSize_); loadParameter += stepSize_; - solverInfo = nonLinearSolver->solve(); + solverInfo = nonLinearSolver_->solve(); info.solverInfos.push_back(solverInfo); info.totalIterations += solverInfo.iterations; if (not solverInfo.success) diff --git a/ikarus/controlroutines/pathfollowing.hh b/ikarus/controlroutines/pathfollowing.hh index 8aace93fe..1193c95db 100644 --- a/ikarus/controlroutines/pathfollowing.hh +++ b/ikarus/controlroutines/pathfollowing.hh @@ -27,19 +27,18 @@ namespace Impl { * This consteval function checks whether the provided template parameters meet the requirements * for path-following control routines. * - * \tparam NonLinearSolver Type of the nonlinear solver. - * \tparam PathFollowingType Type of the path-following strategy. - * \tparam AdaptiveStepSizing Type of the adaptive step sizing strategy. + * \tparam NLS Type of the nonlinear solver. + * \tparam PF Type of the path-following strategy. + * \tparam ASS Type of the adaptive step sizing strategy. * * \return True if the templates meet the requirements, false otherwise. */ - template + template consteval bool checkPathFollowingTemplates() { - return Concepts::PathFollowingStrategy and - Concepts::AdaptiveStepSizingStrategy> and - Concepts::NonLinearSolverCheckForPathFollowing; + return Concepts::PathFollowingStrategy and + Concepts::AdaptiveStepSizingStrategy> and + Concepts::NonLinearSolverCheckForPathFollowing; } } // namespace Impl @@ -69,13 +68,12 @@ namespace Impl { * Currently the following subsidiary functions are implemented \ref LoadControlSubsidiaryFunction, ArcLength and * DisplacementControl * \ingroup controlroutines - * \tparam NonLinearSolver Type of the nonlinear solver used in the control routine. - * \tparam PathFollowingType Type of the path-following strategy. - * \tparam AdaptiveStepSizing Type of the adaptive step sizing strategy. + * \tparam NLS Type of the nonlinear solver used in the control routine. + * \tparam PF Type of the path-following strategy. + * \tparam ASS Type of the adaptive step sizing strategy. */ -template -requires(Impl::checkPathFollowingTemplates()) +template +requires(Impl::checkPathFollowingTemplates()) class PathFollowing : public IObservable { public: @@ -85,19 +83,19 @@ public: /** * \brief Constructor for PathFollowing. * - * \param p_nonLinearSolver Shared pointer to the nonlinear solver. + * \param nonLinearSolver Shared pointer to the nonlinear solver. * \param steps Number of steps in the control routine. * \param stepSize Size of each step. - * \param p_pathFollowingType Type of the path-following function. - * \param p_adaptiveStepSizing Type of the adaptive step sizing strategy. + * \param pathFollowingType Type of the path-following function. + * \param adaptiveStepSizing Type of the adaptive step sizing strategy. */ - PathFollowing(const std::shared_ptr& p_nonLinearSolver, int steps, double stepSize, - PathFollowingType p_pathFollowingType = ArcLength{}, AdaptiveStepSizing p_adaptiveStepSizing = {}) - : nonLinearSolver{p_nonLinearSolver}, + PathFollowing(const std::shared_ptr& nonLinearSolver, int steps, double stepSize, + PF pathFollowingType = ArcLength{}, ASS adaptiveStepSizing = {}) + : nonLinearSolver_{nonLinearSolver}, steps_{steps}, stepSize_{stepSize}, - pathFollowingType_{p_pathFollowingType}, - adaptiveStepSizing{p_adaptiveStepSizing} {} + pathFollowingType_{pathFollowingType}, + adaptiveStepSizing_{adaptiveStepSizing} {} /** * \brief Executes the PathFollowing routine. @@ -107,11 +105,11 @@ public: ControlInformation run(); private: - std::shared_ptr nonLinearSolver; + std::shared_ptr nonLinearSolver_; int steps_; double stepSize_; - PathFollowingType pathFollowingType_; - AdaptiveStepSizing adaptiveStepSizing; + PF pathFollowingType_; + ASS adaptiveStepSizing_; }; } // namespace Ikarus diff --git a/ikarus/controlroutines/pathfollowing.inl b/ikarus/controlroutines/pathfollowing.inl index 675b07324..59ec5992f 100644 --- a/ikarus/controlroutines/pathfollowing.inl +++ b/ikarus/controlroutines/pathfollowing.inl @@ -21,11 +21,11 @@ namespace Ikarus { -template -requires(Impl::checkPathFollowingTemplates()) -ControlInformation PathFollowing::run() { +template +requires(Impl::checkPathFollowingTemplates()) +ControlInformation PathFollowing::run() { ControlInformation info; - auto& nonOp = nonLinearSolver->nonLinearOperator(); + auto& nonOp = nonLinearSolver_->nonLinearOperator(); this->notify(ControlMessages::CONTROL_STARTED, pathFollowingType_.name()); SubsidiaryArgs subsidiaryArgs; @@ -39,7 +39,7 @@ ControlInformation PathFollowingnotify(ControlMessages::STEP_STARTED, 0, subsidiaryArgs.stepSize); auto subsidiaryFunctionPtr = [&](auto&& args) { return pathFollowingType_(args); }; pathFollowingType_.initialPrediction(nonOp, subsidiaryArgs); - auto solverInfo = nonLinearSolver->solve(subsidiaryFunctionPtr, subsidiaryArgs); + auto solverInfo = nonLinearSolver_->solve(subsidiaryFunctionPtr, subsidiaryArgs); info.solverInfos.push_back(solverInfo); info.totalIterations += solverInfo.iterations; if (not solverInfo.success) @@ -51,13 +51,13 @@ ControlInformation PathFollowingnotify(ControlMessages::STEP_STARTED, subsidiaryArgs.currentStep, subsidiaryArgs.stepSize); pathFollowingType_.intermediatePrediction(nonOp, subsidiaryArgs); - solverInfo = nonLinearSolver->solve(subsidiaryFunctionPtr, subsidiaryArgs); + solverInfo = nonLinearSolver_->solve(subsidiaryFunctionPtr, subsidiaryArgs); info.solverInfos.push_back(solverInfo); info.totalIterations += solverInfo.iterations; diff --git a/ikarus/controlroutines/pathfollowingfunctions.hh b/ikarus/controlroutines/pathfollowingfunctions.hh index 4785d949b..e0ea5aae2 100644 --- a/ikarus/controlroutines/pathfollowingfunctions.hh +++ b/ikarus/controlroutines/pathfollowingfunctions.hh @@ -33,7 +33,7 @@ namespace Ikarus { * \brief Structure containing arguments for subsidiary functions. * * This structure holds various arguments used by subsidiary functions in control routines. - * @ingroup controlroutines + * \ingroup controlroutines */ struct SubsidiaryArgs { @@ -61,7 +61,7 @@ struct SubsidiaryArgs * factor increment. \f$\psi\f$ is the to-be-determined correction factor for the different dimensionalities between * \f$\mathrm{D}\mathbf{D}\f$ and \f$\mathrm{D} \lambda\f$. The scalar \f$\hat{s} \f$ defines the requested size of * the step. - * @ingroup controlroutines + * \ingroup controlroutines */ struct ArcLength { @@ -90,15 +90,15 @@ struct ArcLength * This method initializes the prediction step for the standard arc-length method it computes \f$\psi\f$ and * computes initial \f$\mathrm{D}\mathbf{D}\f$ and \f$\mathrm{D} \lambda\f$. * - * \tparam NonLinearOperator Type of the nonlinear operator. + * \tparam NLO Type of the nonlinear operator. * \param nonLinearOperator The nonlinear operator. * \param args The subsidiary function arguments. - * @ingroup controlroutines + * \ingroup controlroutines */ - template - void initialPrediction(NonLinearOperator& nonLinearOperator, SubsidiaryArgs& args) { + template + void initialPrediction(NLO& nonLinearOperator, SubsidiaryArgs& args) { SolverTypeTag solverTag; - using JacobianType = std::remove_cvref_t; + using JacobianType = std::remove_cvref_t; static_assert((traits::isSpecializationTypeAndNonTypes::value) or (traits::isSpecializationTypeNonTypeAndType::value), "Linear solver not implemented for the chosen derivative type of the non-linear operator"); @@ -114,9 +114,8 @@ struct ArcLength const auto& K = nonLinearOperator.derivative(); static constexpr bool isLinearSolver = - Ikarus::Concepts::LinearSolverCheck; + Ikarus::Concepts::LinearSolverCheck; static_assert(isLinearSolver, "Initial predictor step in the standard arc-length method doesn't have a linear solver"); @@ -142,12 +141,12 @@ struct ArcLength * * This method updates the prediction step for the standard arc-length method. * - * \tparam NonLinearOperator Type of the nonlinear operator. + * \tparam NLO Type of the nonlinear operator. * \param nonLinearOperator The nonlinear operator. * \param args The subsidiary function arguments. */ - template - void intermediatePrediction(NonLinearOperator& nonLinearOperator, SubsidiaryArgs& args) { + template + void intermediatePrediction(NLO& nonLinearOperator, SubsidiaryArgs& args) { nonLinearOperator.firstParameter() += args.DD; nonLinearOperator.lastParameter() += args.Dlambda; } @@ -169,7 +168,7 @@ private: * \mathrm{D} \lambda - \hat{s}, \f] * where \f$\mathrm{D}\mathbf{D}\f$ is the increment of the solution vector and \f$\mathrm{D} \lambda\f$ is the load * factor increment. The scalar \f$\hat{s} \f$ defines the requested size of the step. - * @ingroup controlroutines + * \ingroup controlroutines */ struct LoadControlSubsidiaryFunction { @@ -191,12 +190,12 @@ struct LoadControlSubsidiaryFunction * * This method initializes the prediction step for the load control method. * - * \tparam NonLinearOperator Type of the nonlinear operator. + * \tparam NLO Type of the nonlinear operator. * \param nonLinearOperator The nonlinear operator. * \param args The subsidiary function arguments. */ - template - void initialPrediction(NonLinearOperator& nonLinearOperator, SubsidiaryArgs& args) { + template + void initialPrediction(NLO& nonLinearOperator, SubsidiaryArgs& args) { args.Dlambda = args.stepSize; nonLinearOperator.lastParameter() = args.Dlambda; } @@ -206,12 +205,12 @@ struct LoadControlSubsidiaryFunction * * This method updates the prediction step for the load control method. * - * \tparam NonLinearOperator Type of the nonlinear operator. + * \tparam NLO Type of the nonlinear operator. * \param nonLinearOperator The nonlinear operator. * \param args The subsidiary function arguments. */ - template - void intermediatePrediction(NonLinearOperator& nonLinearOperator, SubsidiaryArgs& args) { + template + void intermediatePrediction(NLO& nonLinearOperator, SubsidiaryArgs& args) { nonLinearOperator.lastParameter() += args.Dlambda; } @@ -229,7 +228,7 @@ struct LoadControlSubsidiaryFunction * ||\mathrm{D}\mathbf{D}|| - \hat{s}, \f] * where \f$\mathrm{D}\mathbf{D}\f$ is the increment of the solution vector and \f$\mathrm{D} \lambda\f$ is the load * factor increment. The scalar \f$\hat{s} \f$ defines the requested size of the step. - * @ingroup controlroutines + * \ingroup controlroutines */ struct DisplacementControl { @@ -261,12 +260,12 @@ struct DisplacementControl * * This method initializes the prediction step for the displacement control method. * - * \tparam NonLinearOperator Type of the nonlinear operator. + * \tparam NLO Type of the nonlinear operator. * \param nonLinearOperator The nonlinear operator. * \param args The subsidiary function arguments. */ - template - void initialPrediction(NonLinearOperator& nonLinearOperator, SubsidiaryArgs& args) { + template + void initialPrediction(NLO& nonLinearOperator, SubsidiaryArgs& args) { args.DD(controlledIndices).array() = args.stepSize; nonLinearOperator.firstParameter() = args.DD; } @@ -276,12 +275,12 @@ struct DisplacementControl * * This method updates the prediction step for the displacement control method. * - * \tparam NonLinearOperator Type of the nonlinear operator. + * \tparam NLO Type of the nonlinear operator. * \param nonLinearOperator The nonlinear operator. * \param args The subsidiary function arguments. */ - template - void intermediatePrediction(NonLinearOperator& nonLinearOperator, SubsidiaryArgs& args) { + template + void intermediatePrediction(NLO& nonLinearOperator, SubsidiaryArgs& args) { nonLinearOperator.firstParameter() += args.DD; } diff --git a/ikarus/finiteelements/febases/autodifffe.hh b/ikarus/finiteelements/febases/autodifffe.hh index 6018974b3..c69e41ed2 100644 --- a/ikarus/finiteelements/febases/autodifffe.hh +++ b/ikarus/finiteelements/febases/autodifffe.hh @@ -2,8 +2,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file autoDiffFE.hh - * @brief Contains the AutoDiffFE class, an automatic differentiation wrapper for finite elements. + * \file autoDiffFE.hh + * \brief Contains the AutoDiffFE class, an automatic differentiation wrapper for finite elements. */ #pragma once @@ -18,31 +18,30 @@ namespace Ikarus { /** - * @brief AutoDiffFE class, an automatic differentiation wrapper for finite elements. + * \brief AutoDiffFE class, an automatic differentiation wrapper for finite elements. * - * @tparam RealFE_ The type of the original finite element, which does not implement the derivatives - * @tparam FERequirementType_ Type of the Finite Element Requirements. - * @tparam useEigenRef A boolean indicating whether to use Eigen::Ref for references types in calculateMatrix,... - * @tparam forceAutoDiff A boolean indicating whether to force the use of automatic differentiation, even when the + * \tparam FEImpl The type of the original finite element, which does not implement the derivatives + * \tparam FER Type of the Finite Element Requirements. + * \tparam useEigenRef A boolean indicating whether to use Eigen::Ref for references types in calculateMatrix,... + * \tparam forceAutoDiff A boolean indicating whether to force the use of automatic differentiation, even when the * real element implements the derivatives. */ -template , bool useEigenRef = false, - bool forceAutoDiff = false> -class AutoDiffFE : public RealFE_ +template , bool useEigenRef = false, bool forceAutoDiff = false> +class AutoDiffFE : public FEImpl { public: - using RealFE = RealFE_; ///< Type of the base finite element. - using Basis = typename RealFE::Basis; ///< Type of the basis. - using Traits = FETraits; ///< Type traits for local view. - using LocalView = typename Traits::LocalView; ///< Type of the local view. - using Element = typename Traits::Element; ///< Type of the element. + using RealFE = FEImpl; ///< Type of the base finite element. + using Basis = typename RealFE::Basis; ///< Type of the basis. + using Traits = FETraits; ///< Type traits for local view. + using LocalView = typename Traits::LocalView; ///< Type of the local view. + using Element = typename Traits::Element; ///< Type of the element. using FERequirementType = typename Traits::FERequirementType; ///< Type of the Finite Element Requirements. /** - * @brief Calculate the matrix associated with the finite element. + * \brief Calculate the matrix associated with the finite element. * - * @param req Finite Element Requirements. - * @param h Matrix to be calculated. + * \param req Finite Element Requirements. + * \param h Matrix to be calculated. */ void calculateMatrix(const FERequirementType& req, typename Traits::template MatrixType<> h) const { // real element implements calculateMatrix by itself, then we simply forward the call @@ -81,10 +80,10 @@ public: } /** - * @brief Calculate the vector associated with the finite element. + * \brief Calculate the vector associated with the finite element. * - * @param req Finite Element Requirements. - * @param g Vector to be calculated. + * \param req Finite Element Requirements. + * \param g Vector to be calculated. */ void calculateVector(const FERequirementType& req, typename Traits::template VectorType<> g) const { // real element implements calculateVector by itself, then we simply forward the call @@ -111,11 +110,11 @@ public: } /** - * @brief Calculate the local system associated with the finite element. + * \brief Calculate the local system associated with the finite element. * - * @param req Finite Element Requirements. - * @param h Matrix to be calculated. - * @param g Vector to be calculated. + * \param req Finite Element Requirements. + * \param h Matrix to be calculated. + * \param g Vector to be calculated. */ void calculateLocalSystem(const FERequirementType& req, typename Traits::template MatrixType<> h, typename Traits::template VectorType<> g) const { @@ -126,10 +125,10 @@ public: } /** - * @brief Calculate the scalar value associated with the finite element. + * \brief Calculate the scalar value associated with the finite element. * - * @param par Finite Element Requirements. - * @return The calculated scalar value. + * \param par Finite Element Requirements. + * \return The calculated scalar value. */ [[nodiscard]] double calculateScalar(const FERequirementType& par) const { // real element implements calculateScalar by itself, then we simply forward the call @@ -146,18 +145,18 @@ public: } /** - * @brief Get the reference to the base finite element. + * \brief Get the reference to the base finite element. * - * @return The reference to the base finite element. + * \return The reference to the base finite element. */ const RealFE& realFE() const { return *this; } /** - * @brief Constructor for the AutoDiffFE class. + * \brief Constructor for the AutoDiffFE class. * Forward the construction to the underlying element * - * @tparam Args Variadic template for constructor arguments. - * @param args Constructor arguments. + * \tparam Args Variadic template for constructor arguments. + * \param args Constructor arguments. */ template explicit AutoDiffFE(Args&&... args) diff --git a/ikarus/finiteelements/febases/powerbasisfe.hh b/ikarus/finiteelements/febases/powerbasisfe.hh index f674ce931..8b28c357e 100644 --- a/ikarus/finiteelements/febases/powerbasisfe.hh +++ b/ikarus/finiteelements/febases/powerbasisfe.hh @@ -2,8 +2,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file powerbasisfe.hh - * @brief Contains the PowerBasisFE class, which works with a power basis in FlatInterLeaved elements. + * \file powerbasisfe.hh + * \brief Contains the PowerBasisFE class, which works with a power basis in FlatInterLeaved elements. */ #pragma once @@ -13,52 +13,52 @@ namespace Ikarus { /** - * @brief PowerBasisFE class for working with a power basis in FlatInterLeaved elements. + * \brief PowerBasisFE class for working with a power basis in FlatInterLeaved elements. * * This class assumes the local indices in a FlatInterLeaved order. * - * @tparam Basis The type of the power basis. + * \tparam B The type of the power basis. */ -template +template class PowerBasisFE { public: - using Traits = FETraits; ///< Type of the traits. - using RootBasis = typename Traits::FlatBasis; ///< Type of the root basis. + using Traits = FETraits; ///< Type of the traits. + using FlatBasis = typename Traits::FlatBasis; ///< Type of the flat basis. using LocalView = typename Traits::LocalView; ///< Type of the local view. using GlobalIndex = typename Traits::GlobalIndex; ///< Type of the global index. using GridElement = typename Traits::Element; ///< Type of the grid element. /** - * @brief Constructor for the PowerBasisFE class. + * \brief Constructor for the PowerBasisFE class. * - * @param p_basis The power basis. - * @param element The local element. + * \param p_basis The power basis. + * \param element The local element. */ - explicit PowerBasisFE(const Basis& p_basis, const typename LocalView::Element& element) + explicit PowerBasisFE(const B& p_basis, const typename LocalView::Element& element) : localView_{p_basis.flat().localView()} { - static_assert(Ikarus::Concepts::PowerBasis, + static_assert(Ikarus::Concepts::PowerBasis, "You didn't pass a localview of a power basis to this method"); - static_assert(RootBasis::PreBasis::Node::degree() != 1, "The basis has only one children. Maybe use scalarFE.hh."); + static_assert(FlatBasis::PreBasis::Node::degree() != 1, "The basis has only one children. Maybe use scalarFE.hh."); localView_.bind(element); } - /** @brief Number of children in the powerBasis. */ - static constexpr int num_children = RootBasis::PreBasis::Node::degree(); + /** \brief Number of children in the powerBasis. */ + static constexpr int num_children = FlatBasis::PreBasis::Node::degree(); /** - * @brief Get the size of the local view. - * @return The size of the local view. + * \brief Get the size of the local view. + * \return The size of the local view. */ [[nodiscard]] constexpr size_t size() const { return localView_.size(); } /** - * @brief Get the global flat indices for the power basis. + * \brief Get the global flat indices for the power basis. * * The global indices are collected in a FlatInterLeaved order. I.e. x_0, y_0, z_0, ..., x_n, y_n, z_n * - * @param globalIndices Output vector to store global indices. + * \param globalIndices Output vector to store global indices. */ void globalFlatIndices(std::vector& globalIndices) const { globalIndices.clear(); @@ -72,20 +72,20 @@ public: } /** - * @brief Get the grid element associated with the local view. - * @return The grid element. + * \brief Get the grid element associated with the local view. + * \return The grid element. */ const GridElement& gridElement() const { return localView_.element(); } /** - * @brief Get the const reference to the local view. - * @return The const reference to the local view. + * \brief Get the const reference to the local view. + * \return The const reference to the local view. */ const LocalView& localView() const { return localView_; } /** - * @brief Get the reference to the local view. - * @return The reference to the local view. + * \brief Get the reference to the local view. + * \return The reference to the local view. */ LocalView& localView() { return localView_; } diff --git a/ikarus/finiteelements/febases/scalarfe.hh b/ikarus/finiteelements/febases/scalarfe.hh index 1e0c1cb27..9bddab76d 100644 --- a/ikarus/finiteelements/febases/scalarfe.hh +++ b/ikarus/finiteelements/febases/scalarfe.hh @@ -2,8 +2,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file scalarFE.hh - * @brief Contains the ScalarFieldFE class, a class for Single-DOF elements using a scalar basis. + * \file scalarFE.hh + * \brief Contains the ScalarFieldFE class, a class for Single-DOF elements using a scalar basis. */ #pragma once @@ -13,44 +13,43 @@ namespace Ikarus { /** - * @brief ScalarFieldFE class, a class for Single-DOF elements using a scalar basis. + * \brief ScalarFieldFE class, a class for Single-DOF elements using a scalar basis. * - * @tparam Basis_ Type of the scalar basis. + * \tparam B Type of the scalar basis. */ -template +template class ScalarFieldFE { public: - using Traits = FETraits; ///< Type of the traits. - using Basis = typename Traits::Basis; ///< Type of the basis. - using RootBasis = typename Traits::FlatBasis; ///< Type of the root basis. + using Traits = FETraits; ///< Type of the traits. + using FlatBasis = typename Traits::FlatBasis; ///< Type of the root basis. using LocalView = typename Traits::LocalView; ///< Type of the local view. using GlobalIndex = typename Traits::GlobalIndex; ///< Type of the global index. using GridElement = typename Traits::Element; ///< Type of the grid element. static constexpr int worlddim = Traits::worlddim; ///< Dimension of the world space. /** - * @brief Constructor for the ScalarFieldFE class. + * \brief Constructor for the ScalarFieldFE class. * - * @param basis The scalar basis. - * @param element The local element. + * \param basis The scalar basis. + * \param element The local element. */ - explicit ScalarFieldFE(const Basis& basis, const typename LocalView::Element& element) + explicit ScalarFieldFE(const B& basis, const typename LocalView::Element& element) : localView_{basis.flat().localView()} { - static_assert(RootBasis::PreBasis::Node::CHILDREN == 0, "This is no scalar basis!"); + static_assert(FlatBasis::PreBasis::Node::CHILDREN == 0, "This is no scalar basis!"); localView_.bind(element); } /** - * @brief Get the size of the finite element. + * \brief Get the size of the finite element. * - * @return The size of the finite element. + * \return The size of the finite element. */ [[nodiscard]] int size() const { return localView_.size(); } /** - * @brief Get the global flat indices of the finite element. + * \brief Get the global flat indices of the finite element. * - * @param globalIndices Vector to store global flat indices. + * \param globalIndices Vector to store global flat indices. */ void globalFlatIndices(std::vector& globalIndices) const { const auto& fe = localView_.tree().finiteElement(); @@ -59,23 +58,23 @@ public: } /** - * @brief Get the entity associated with the finite element. + * \brief Get the entity associated with the finite element. * - * @return The entity associated with the finite element. + * \return The entity associated with the finite element. */ const GridElement& gridElement() const { return localView_.element(); } /** - * @brief Get the local view of the finite element. + * \brief Get the local view of the finite element. * - * @return The local view of the finite element. + * \return The local view of the finite element. */ const LocalView& localView() const { return localView_; } /** - * @brief Get the local view of the finite element. + * \brief Get the local view of the finite element. * - * @return The local view of the finite element. + * \return The local view of the finite element. */ LocalView& localView() { return localView_; } diff --git a/ikarus/finiteelements/fehelper.hh b/ikarus/finiteelements/fehelper.hh index fc40e2486..05acbd2c1 100644 --- a/ikarus/finiteelements/fehelper.hh +++ b/ikarus/finiteelements/fehelper.hh @@ -10,24 +10,24 @@ namespace Ikarus::FEHelper { /** - * @brief Gets the local solution Dune block vector + * \brief Gets the local solution Dune block vector * - * @tparam Traits Type of the FE traits. - * @tparam ScalarType Scalar type for the local solution vector. + * \tparam Traits Type of the FE traits. + * \tparam ST Scalar type for the local solution vector. * - * @param x The global solution vector. - * @param localView Local view of the element. - * @param dx Optional global solution vector. + * \param x The global solution vector. + * \param localView Local view of the element. + * \param dx Optional global solution vector. * - * @return A Dune block vector representing the solution quantities at each node. + * \return A Dune block vector representing the solution quantities at each node. * */ -template +template auto localSolutionBlockVector(const typename Traits::FERequirementType::SolutionVectorTypeRaw& x, const typename Traits::LocalView& localView, - const std::optional>& dx = std::nullopt) { + const std::optional>& dx = std::nullopt) { constexpr int worldDim = Traits::worlddim; const auto& fe = localView.tree().child(0).finiteElement(); - Dune::BlockVector> localX(fe.size()); + Dune::BlockVector> localX(fe.size()); if (dx) for (auto i = 0U; i < localX.size(); ++i) for (auto j = 0U; j < worldDim; ++j) diff --git a/ikarus/finiteelements/ferequirements.hh b/ikarus/finiteelements/ferequirements.hh index 9ac043a95..6f0e3d662 100644 --- a/ikarus/finiteelements/ferequirements.hh +++ b/ikarus/finiteelements/ferequirements.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file ferequirements.hh - * @brief Definition of the LinearElastic class for finite element mechanics computations. - * @ingroup finiteelements + * \file ferequirements.hh + * \brief Definition of the LinearElastic class for finite element mechanics computations. + * \ingroup finiteelements */ #pragma once @@ -117,11 +117,11 @@ struct AffordanceCollectionImpl /** * \brief Concept to check if a given type is one of the predefined affordance enums or the AffordanceCollectionImpl. */ -template -concept FEAffordance = std::is_same_v, ScalarAffordances> or - std::is_same_v, VectorAffordances> or - std::is_same_v, MatrixAffordances> or - std::is_same_v, AffordanceCollectionImpl>; +template +concept FEAffordance = std::is_same_v, ScalarAffordances> or + std::is_same_v, VectorAffordances> or + std::is_same_v, MatrixAffordances> or + std::is_same_v, AffordanceCollectionImpl>; inline constexpr VectorAffordances forces = VectorAffordances::forces; @@ -164,19 +164,18 @@ namespace Impl { * and parameters needed. It provides methods to add affordances, insert parameters, and manage global solution * vectors. * - * \tparam SolutionVectorType_ Type of the solution vector, defaulting to std::reference_wrapper. - * \tparam ParameterType_ Type of the parameter, defaulting to std::reference_wrapper. + * \tparam SV Type of the solution vector, defaulting to std::reference_wrapper. + * \tparam PM Type of the parameter, defaulting to std::reference_wrapper. * */ -template , - typename ParameterType_ = std::reference_wrapper> +template , typename PM = std::reference_wrapper> class FERequirements { public: - using SolutionVectorType = SolutionVectorType_; - using SolutionVectorTypeRaw = typename Impl::DeduceRawVectorType>::Type; - using ParameterType = ParameterType_; - using ParameterTypeRaw = typename ParameterType_::type; + using SolutionVectorType = SV; + using SolutionVectorTypeRaw = typename Impl::DeduceRawVectorType>::Type; + using ParameterType = PM; + using ParameterTypeRaw = typename PM::type; /** * \brief Add an affordance to the requirements. @@ -190,13 +189,13 @@ public: template FERequirements& addAffordance(Affordance&& affordance) { if constexpr (std::is_same_v) - affordances.scalarAffordances = affordance; + affordances_.scalarAffordances = affordance; else if constexpr (std::is_same_v) - affordances.vectorAffordances = affordance; + affordances_.vectorAffordances = affordance; else if constexpr (std::is_same_v) - affordances.matrixAffordances = affordance; + affordances_.matrixAffordances = affordance; else if constexpr (std::is_same_v) - affordances = affordance; + affordances_ = affordance; return *this; } @@ -210,7 +209,7 @@ public: * \return Reference to the updated FERequirements instance. */ FERequirements& insertParameter(const FEParameter& key, ParameterTypeRaw& val) { - parameter.insert_or_assign(key, val); + parameter_.insert_or_assign(key, val); return *this; } @@ -224,7 +223,7 @@ public: * \return Reference to the updated FERequirements instance. */ FERequirements& insertGlobalSolution(const FESolutions& key, SolutionVectorTypeRaw& sol) { - sols.insert_or_assign(key, sol); + sols_.insert_or_assign(key, sol); return *this; } @@ -241,9 +240,9 @@ public: const SolutionVectorTypeRaw& getGlobalSolution(const FESolutions& key) const { try { if constexpr (std::is_same_v>) - return sols.at(key).get(); + return sols_.at(key).get(); else - return sols.at(key); + return sols_.at(key); } catch (std::out_of_range& oor) { DUNE_THROW(Dune::RangeError, std::string("Out of Range error: ") + std::string(oor.what()) + " in getGlobalSolution with key" + toString(key)); @@ -263,7 +262,7 @@ public: */ const ParameterTypeRaw& getParameter(FEParameter&& key) const { try { - return parameter.at(key).get(); + return parameter_.at(key).get(); } catch (std::out_of_range& oor) { DUNE_THROW(Dune::RangeError, std::string("Out of Range error: ") + std::string(oor.what()) + " in getParameter with key" + toString(key)); @@ -283,19 +282,19 @@ public: template bool hasAffordance(Affordance&& affordance) const { if constexpr (std::is_same_v) - return affordances.scalarAffordances == affordance; + return affordances_.scalarAffordances == affordance; else if constexpr (std::is_same_v) - return affordances.vectorAffordances == affordance; + return affordances_.vectorAffordances == affordance; else if constexpr (std::is_same_v) - return affordances.matrixAffordances == affordance; + return affordances_.matrixAffordances == affordance; else if constexpr (std::is_same_v) - return affordances == affordance; + return affordances_ == affordance; } private: - std::map sols; - std::map parameter; - AffordanceCollectionImpl affordances; + std::map sols_; + std::map parameter_; + AffordanceCollectionImpl affordances_; }; /** @@ -307,17 +306,16 @@ private: * and parameters needed. It provides methods to add affordances, insert parameters, and manage global solution * vectors. * - * \tparam SolutionVectorType_ Type of the solution vector, defaulting to std::reference_wrapper. - * \tparam ParameterType_ Type of the parameter, defaulting to std::reference_wrapper. + * \tparam SV Type of the solution vector, defaulting to std::reference_wrapper. + * \tparam PM Type of the parameter, defaulting to std::reference_wrapper. * */ -template , - typename ParameterType_ = std::reference_wrapper> +template , typename PM = std::reference_wrapper> class [[deprecated( "FErequirements is deprecaded and will be removed after v0.5. Use FERequirements instead.")]] FErequirements - : public FERequirements + : public FERequirements { - using Base = FERequirements; + using Base = FERequirements; public: FErequirements() = default; diff --git a/ikarus/finiteelements/fetraits.hh b/ikarus/finiteelements/fetraits.hh index c9052b7bc..8f45a4fcd 100644 --- a/ikarus/finiteelements/fetraits.hh +++ b/ikarus/finiteelements/fetraits.hh @@ -16,18 +16,18 @@ namespace Ikarus { /** * \brief Traits for handling finite elements. * - * \tparam Basis_ The basis type for the finite element. - * \tparam FERequirements_ The requirements for the finite element. + * \tparam B The basis type for the finite element. + * \tparam FER The requirements for the finite element. * \tparam useRef Boolean indicating whether to use Eigen::Ref for VectorType and MatrixType. */ -template , bool useRef = false> +template , bool useRef = false> struct FETraits { /** \brief Type of the basis of the finite element */ - using Basis = Basis_; + using Basis = B; /** \brief Type of the requirements for the finite element */ - using FERequirementType = FERequirements_; + using FERequirementType = FER; /** \brief Type of the flat basis */ using FlatBasis = typename Basis::FlatBasis; diff --git a/ikarus/finiteelements/mechanics/eas/eas2d.hh b/ikarus/finiteelements/mechanics/eas/eas2d.hh index 1995c70b8..8c0df0498 100644 --- a/ikarus/finiteelements/mechanics/eas/eas2d.hh +++ b/ikarus/finiteelements/mechanics/eas/eas2d.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file eas2d.hh - * @brief Definition of the types of EAS formulations for 2D elements. - * @ingroup eas + * \file eas2d.hh + * \brief Definition of the types of EAS formulations for 2D elements. + * \ingroup eas */ #pragma once @@ -14,15 +14,15 @@ namespace Ikarus::EAS { /** - * @brief Q1E4 structure for EAS with linear strains and 4 enhanced modes. + * \brief Q1E4 structure for EAS with linear strains and 4 enhanced modes. * * \details The Q1E4 structure represents an implementation of EAS for a * specific case where linear strains are considered, and the method includes * enhancement with 4 additional modes to improve the accuracy of finite element solutions. * - * @tparam Geometry The geometry type. + * \tparam GEO The geometry type. */ -template +template struct Q1E4 { static constexpr int strainSize = 3; @@ -30,9 +30,9 @@ struct Q1E4 using MType = Eigen::Matrix; Q1E4() = default; - explicit Q1E4(const Geometry& geometry_) - : geometry{std::make_shared(geometry_)}, - T0InverseTransformed{calcTransformationMatrix2D(geometry_)} {} + explicit Q1E4(const GEO& geometry) + : geometry_{std::make_shared(geometry)}, + T0InverseTransformed_{calcTransformationMatrix2D(geometry)} {} auto calcM(const Dune::FieldVector& quadPos) const { MType M; @@ -43,23 +43,24 @@ struct Q1E4 M(1, 1) = 2 * eta - 1.0; M(2, 2) = 2 * xi - 1.0; M(2, 3) = 2 * eta - 1.0; - const double detJ = geometry->integrationElement(quadPos); - M = T0InverseTransformed / detJ * M; + const double detJ = geometry_->integrationElement(quadPos); + M = T0InverseTransformed_ / detJ * M; return M; } - std::shared_ptr geometry; - Eigen::Matrix3d T0InverseTransformed; +private: + std::shared_ptr geometry_; + Eigen::Matrix3d T0InverseTransformed_; }; /** - * @brief Structure representing EAS for Q1 with 5 enhanced strains. + * \brief Structure representing EAS for Q1 with 5 enhanced strains. * * This structure defines the EAS for Q1 elements with 5 enhanced strains. * - * @tparam Geometry The geometry type. + * \tparam GEO The geometry type. */ -template +template struct Q1E5 { static constexpr int strainSize = 3; @@ -67,9 +68,9 @@ struct Q1E5 using MType = Eigen::Matrix; Q1E5() = default; - explicit Q1E5(const Geometry& geometry_) - : geometry{std::make_shared(geometry_)}, - T0InverseTransformed{calcTransformationMatrix2D(geometry_)} {} + explicit Q1E5(const GEO& geometry) + : geometry_{std::make_shared(geometry)}, + T0InverseTransformed_{calcTransformationMatrix2D(geometry)} {} auto calcM(const Dune::FieldVector& quadPos) const { MType M; @@ -81,23 +82,24 @@ struct Q1E5 M(2, 2) = 2 * xi - 1.0; M(2, 3) = 2 * eta - 1.0; M(2, 4) = (2 * xi - 1.0) * (2 * eta - 1.0); - const double detJ = geometry->integrationElement(quadPos); - M = T0InverseTransformed / detJ * M; + const double detJ = geometry_->integrationElement(quadPos); + M = T0InverseTransformed_ / detJ * M; return M; } - std::shared_ptr geometry; - Eigen::Matrix3d T0InverseTransformed; +private: + std::shared_ptr geometry_; + Eigen::Matrix3d T0InverseTransformed_; }; /** - * @brief Structure representing EAS for Q1 with 7 enhanced strains. + * \brief Structure representing EAS for Q1 with 7 enhanced strains. * * This structure defines the EAS for Q1 elements with 7 enhanced strains. * - * @tparam Geometry The geometry type. + * \tparam GEO The geometry type. */ -template +template struct Q1E7 { static constexpr int strainSize = 3; @@ -105,9 +107,9 @@ struct Q1E7 using MType = Eigen::Matrix; Q1E7() = default; - explicit Q1E7(const Geometry& geometry_) - : geometry{std::make_shared(geometry_)}, - T0InverseTransformed{calcTransformationMatrix2D(geometry_)} {} + explicit Q1E7(const GEO& geometry) + : geometry_{std::make_shared(geometry)}, + T0InverseTransformed_{calcTransformationMatrix2D(geometry)} {} auto calcM(const Dune::FieldVector& quadPos) const { MType M; @@ -121,12 +123,13 @@ struct Q1E7 M(0, 4) = (2 * xi - 1.0) * (2 * eta - 1.0); M(1, 5) = (2 * xi - 1.0) * (2 * eta - 1.0); M(2, 6) = (2 * xi - 1.0) * (2 * eta - 1.0); - const double detJ = geometry->integrationElement(quadPos); - M = T0InverseTransformed / detJ * M; + const double detJ = geometry_->integrationElement(quadPos); + M = T0InverseTransformed_ / detJ * M; return M; } - std::shared_ptr geometry; - Eigen::Matrix3d T0InverseTransformed; +private: + std::shared_ptr geometry_; + Eigen::Matrix3d T0InverseTransformed_; }; } // namespace Ikarus::EAS diff --git a/ikarus/finiteelements/mechanics/eas/eas3d.hh b/ikarus/finiteelements/mechanics/eas/eas3d.hh index eae1d930b..93d7c990b 100644 --- a/ikarus/finiteelements/mechanics/eas/eas3d.hh +++ b/ikarus/finiteelements/mechanics/eas/eas3d.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file eas3d.hh - * @brief Definition of the types of EAS formulations for 3D elements. - * @ingroup mechanics + * \file eas3d.hh + * \brief Definition of the types of EAS formulations for 3D elements. + * \ingroup mechanics */ #pragma once @@ -14,13 +14,13 @@ namespace Ikarus::EAS { /** - * @brief Structure representing EAS for H1 with 9 enhanced strains. + * \brief Structure representing EAS for H1 with 9 enhanced strains. * * This structure defines the EAS for H1 elements with 9 enhanced strains. * - * @tparam Geometry The geometry type. + * \tparam GEO The geometry type. */ -template +template struct H1E9 { static constexpr int strainSize = 6; @@ -28,9 +28,9 @@ struct H1E9 using MType = Eigen::Matrix; H1E9() = default; - explicit H1E9(const Geometry& geometry_) - : geometry{std::make_shared(geometry_)}, - T0InverseTransformed{calcTransformationMatrix3D(geometry_)} {} + explicit H1E9(const GEO& geometry) + : geometry_{std::make_shared(geometry)}, + T0InverseTransformed_{calcTransformationMatrix3D(geometry)} {} auto calcM(const Dune::FieldVector& quadPos) const { MType M; @@ -47,20 +47,22 @@ struct H1E9 M(4, 6) = 2 * zeta - 1.0; M(5, 7) = 2 * eta - 1.0; M(5, 8) = 2 * zeta - 1.0; - const double detJ = geometry->integrationElement(quadPos); - M = T0InverseTransformed / detJ * M; + const double detJ = geometry_->integrationElement(quadPos); + M = T0InverseTransformed_ / detJ * M; return M; } - std::shared_ptr geometry; - Eigen::Matrix T0InverseTransformed; + +private: + std::shared_ptr geometry_; + Eigen::Matrix T0InverseTransformed_; }; /** - * @brief Structure representing EAS for H1 with 21 enhanced strains. + * \brief Structure representing EAS for H1 with 21 enhanced strains. * * This structure defines the EAS for H1 elements with 21 enhanced strains. * - * @tparam Geometry The geometry type. + * \tparam Geometry The geometry type. */ template struct H1E21 @@ -109,7 +111,8 @@ struct H1E21 return M; } +private: std::shared_ptr geometry; - Eigen::Matrix T0InverseTransformed; + Eigen::Matrix T0InverseTransformed; }; } // namespace Ikarus::EAS diff --git a/ikarus/finiteelements/mechanics/easvariants.hh b/ikarus/finiteelements/mechanics/easvariants.hh index 650919258..793bc7b10 100644 --- a/ikarus/finiteelements/mechanics/easvariants.hh +++ b/ikarus/finiteelements/mechanics/easvariants.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file easvariants.hh - * @brief Definition of the EAS variants. - * @ingroup mechanics + * \file easvariants.hh + * \brief Definition of the EAS variants. + * \ingroup mechanics */ #pragma once @@ -20,19 +20,19 @@ namespace Ikarus::EAS { * \details EASVariants holds the different variants of EAS formulations for * both 2D and 3D elements. * - * @tparam Geometry The geometry type. + * \tparam GEO The geometry type. */ -template +template struct Variants { - static constexpr int worldDim = Geometry::coorddimension; + static constexpr int worldDim = GEO::coorddimension; static_assert((worldDim == 2) or (worldDim == 3), "EAS variants are only available 2D and 3D elements."); /** \brief Variant type holding different EAS formulations for 2D elements. */ - using EAS2D = std::variant, Q1E5, Q1E7>; + using EAS2D = std::variant, Q1E5, Q1E7>; /** \brief Variant type holding different EAS formulations for 3D elements. */ - using EAS3D = std::variant, H1E21>; + using EAS3D = std::variant, H1E21>; /** \brief Type of the EAS variant depending on the grid dimension. */ using type = std::conditional_t; diff --git a/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh b/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh index d0b5ea8b7..f3d1f2dac 100644 --- a/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh +++ b/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file enhancedassumedstrains.hh - * @brief Definition of the EAS class. - * @ingroup mechanics + * \file enhancedassumedstrains.hh + * \brief Definition of the EAS class. + * \ingroup mechanics */ #pragma once @@ -21,32 +21,33 @@ namespace Ikarus { /** - * @brief Wrapper class for using Enhanced Assumed Strains (EAS) with displacement based elements. + * \brief Wrapper class for using Enhanced Assumed Strains (EAS) with displacement based elements. * - * @ingroup mechanics + * \ingroup mechanics * * This class extends a displacement-based element to support Enhanced Assumed Strains. * - * @tparam DisplacementBasedElement The base displacement-based element type. + * \tparam DFE The base displacement-based element type. */ -template -class EnhancedAssumedStrains : public DisplacementBasedElement +template +class EnhancedAssumedStrains : public DFE { public: - using FERequirementType = typename DisplacementBasedElement::FERequirementType; - using LocalView = typename DisplacementBasedElement::LocalView; - using Geometry = typename DisplacementBasedElement::Geometry; - using GridView = typename DisplacementBasedElement::GridView; - using Traits = typename DisplacementBasedElement::Traits; + using DisplacementBasedElement = DFE; + using FERequirementType = typename DisplacementBasedElement::FERequirementType; + using LocalView = typename DisplacementBasedElement::LocalView; + using Geometry = typename DisplacementBasedElement::Geometry; + using GridView = typename DisplacementBasedElement::GridView; + using Traits = typename DisplacementBasedElement::Traits; using DisplacementBasedElement::localView; /** -* @brief Constructor for Enhanced Assumed Strains elements. +* \brief Constructor for Enhanced Assumed Strains elements. * \details Disabling this forwarding constructor if the argument provided is EnhancedAssumedStrains itself, to forward the // calls to the implicit copy constructor -* @tparam Args Variadic template for constructor arguments. -* @param args Constructor arguments forwarded to the base class. +* \tparam Args Variadic template for constructor arguments. +* \param args Constructor arguments forwarded to the base class. */ template requires( @@ -55,12 +56,12 @@ forward the : DisplacementBasedElement(std::forward(args)...) {} /** - * @brief Calculates a scalar quantity for the element. + * \brief Calculates a scalar quantity for the element. * * This function calculates a scalar quantity for the element based on the FERequirementType. * - * @param par The FERequirementType object. - * @return Computed scalar quantity. + * \param par The FERequirementType object. + * \return Computed scalar quantity. */ double calculateScalar(const FERequirementType& par) const { if (isDisplacementBased()) @@ -70,35 +71,35 @@ forward the } /** - * @brief Checks if the element is displacement-based and the EAS is turned off. + * \brief Checks if the element is displacement-based and the EAS is turned off. * - * @return True if the element is displacement-based, false otherwise. + * \return True if the element is displacement-based, false otherwise. */ bool isDisplacementBased() const { return std::holds_alternative(easVariant_); } /** - * @brief Calculates vectorial quantities for the element. + * \brief Calculates vectorial quantities for the element. * * This function calculates the vectorial quantities for the element based on the FERequirementType. * - * @param par The FERequirementType object. - * @param force Vector to store the calculated forces. + * \param par The FERequirementType object. + * \param force Vector to store the calculated forces. */ inline void calculateVector(const FERequirementType& par, typename Traits::template VectorType<> force) const { calculateVectorImpl(par, force); } /** - * @brief Gets the variant representing the type of Enhanced Assumed Strains (EAS). + * \brief Gets the variant representing the type of Enhanced Assumed Strains (EAS). * - * @return Const reference to the EAS variant. + * \return Const reference to the EAS variant. */ const auto& easVariant() const { return easVariant_; } /** - * @brief Gets the number of EAS parameters based on the current EAS type. + * \brief Gets the number of EAS parameters based on the current EAS type. * - * @return Number of EAS parameters. + * \return Number of EAS parameters. */ auto getNumberOfEASParameters() const { return std::visit( @@ -112,12 +113,12 @@ forward the } /** - * @brief Calculates the matrix for the element. + * \brief Calculates the matrix for the element. * * This function calculates the matrix for the element based on the FERequirementType. * - * @param par The FERequirementType object. - * @param K Matrix to store the calculated stiffness. + * \param par The FERequirementType object. + * \param K Matrix to store the calculated stiffness. */ void calculateMatrix(const FERequirementType& par, typename Traits::template MatrixType<> K) const { using namespace Dune::DerivativeDirections; @@ -136,9 +137,9 @@ forward the if constexpr (not std::is_same_v) { constexpr int enhancedStrainSize = EAST::enhancedStrainSize; Eigen::Matrix D; - calculateDAndLMatrix(easFunction, par, D, L); + calculateDAndLMatrix(easFunction, par, D, L_); - K.template triangularView() -= L.transpose() * D.inverse() * L; + K.template triangularView() -= L_.transpose() * D.inverse() * L_; K.template triangularView() = K.transpose(); } }, @@ -146,18 +147,18 @@ forward the } /** - * @brief Calculates a requested result at a specific local position using the Enhanced Assumed Strains (EAS) + * \brief Calculates a requested result at a specific local position using the Enhanced Assumed Strains (EAS) * method. * * This function calculates the results at the specified local coordinates . * It takes into account the displacement-based element calculations and, if applicable, incorporates the EAS method * for enhanced accuracy. * - * @param req The result requirements. - * @param local The local coordinates at which results are to be calculated. - * @return calculated result + * \param req The result requirements. + * \param local The local coordinates at which results are to be calculated. + * \return calculated result * - * @tparam resType The type representing the requested result. + * \tparam resType The type representing the requested result. */ template auto calculateAt(const FERequirementType& req, const Dune::FieldVector& local) const { @@ -182,8 +183,8 @@ forward the if constexpr (not std::is_same_v) { constexpr int enhancedStrainSize = EAST::enhancedStrainSize; Eigen::Matrix D; - calculateDAndLMatrix(easFunction, req, D, L); - const auto alpha = (-D.inverse() * L * disp).eval(); + calculateDAndLMatrix(easFunction, req, D, L_); + const auto alpha = (-D.inverse() * L_ * disp).eval(); const auto M = easFunction.calcM(local); const auto CEval = C(local); auto easStress = (CEval * M * alpha).eval(); @@ -197,18 +198,18 @@ forward the } /** - * @brief Sets the EAS type for 2D elements. + * \brief Sets the EAS type for 2D elements. * - * @param numberOfEASParameters_ The number of EAS parameters + * \param numberOfEASParameters The number of EAS parameters */ - void setEASType(int numberOfEASParameters_) { + void setEASType(int numberOfEASParameters) { const auto& numberOfNodes = DisplacementBasedElement::numberOfNodes(); if (not((numberOfNodes == 4 and Traits::mydim == 2) or (numberOfNodes == 8 and Traits::mydim == 3)) and (not isDisplacementBased())) DUNE_THROW(Dune::NotImplemented, "EAS only supported for Q1 or H1 elements"); if constexpr (Traits::mydim == 2) { - switch (numberOfEASParameters_) { + switch (numberOfEASParameters) { case 0: easVariant_ = std::monostate(); break; @@ -226,7 +227,7 @@ forward the break; } } else if constexpr (Traits::mydim == 3) { - switch (numberOfEASParameters_) { + switch (numberOfEASParameters) { case 0: easVariant_ = std::monostate(); break; @@ -276,9 +277,9 @@ protected: if constexpr (not std::is_same_v) { constexpr int enhancedStrainSize = EAST::enhancedStrainSize; Eigen::Matrix D; - calculateDAndLMatrix(easFunction, par, D, L); + calculateDAndLMatrix(easFunction, par, D, L_); - const auto alpha = (-D.inverse() * L * disp).eval(); + const auto alpha = (-D.inverse() * L_ * disp).eval(); for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) { const auto M = easFunction.calcM(gp.position()); @@ -299,7 +300,7 @@ protected: private: using EASVariant = EAS::Variants::type; EASVariant easVariant_; - mutable Eigen::MatrixXd L; + mutable Eigen::MatrixXd L_; template void calculateDAndLMatrix(const auto& easFunction, const auto& par, Eigen::Matrix& DMat, diff --git a/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh b/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh index 1dcf62a0f..fcd8be992 100644 --- a/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh +++ b/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh @@ -2,8 +2,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file kirchhoffloveshell.hh - * @brief Definition of the KirchhoffLoveShell class for Kirchhoff-Love shell elements in Ikarus. + * \file kirchhoffloveshell.hh + * \brief Definition of the KirchhoffLoveShell class for Kirchhoff-Love shell elements in Ikarus. */ #pragma once @@ -25,25 +25,23 @@ namespace Ikarus { /** - * @brief Kirchhoff-Love shell finite element class. + * \brief Kirchhoff-Love shell finite element class. * * This class represents Kirchhoff-Love shell finite elements. * - * @ingroup mechanics + * \ingroup mechanics * - * @tparam Basis_ The basis type for the finite element. - * @tparam FERequirements_ The type representing the requirements for finite element calculations. - * @tparam useEigenRef A boolean indicating whether to use Eigen references for efficiency. + * \tparam B The basis type for the finite element. + * \tparam FER The type representing the requirements for finite element calculations. + * \tparam useEigenRef A boolean indicating whether to use Eigen references for efficiency. */ -template , bool useEigenRef = false> -class KirchhoffLoveShell : public PowerBasisFE, - public Volume, - FETraits>, - public Traction, - FETraits> +template , bool useEigenRef = false> +class KirchhoffLoveShell : public PowerBasisFE, + public Volume, FETraits>, + public Traction, FETraits> { public: - using Traits = FETraits; + using Traits = FETraits; using Basis = typename Traits::Basis; using FlatBasis = typename Traits::FlatBasis; using FERequirementType = typename Traits::FERequirementType; @@ -52,8 +50,8 @@ public: using GridView = typename Traits::GridView; using Element = typename Traits::Element; using BasePowerFE = PowerBasisFE; // Handles globalIndices function - using VolumeType = Volume, Traits>; - using TractionType = Traction, Traits>; + using VolumeType = Volume; + using TractionType = Traction; using LocalBasisType = decltype(std::declval().tree().child(0).finiteElement().localBasis()); static constexpr int myDim = Traits::mydim; @@ -68,84 +66,84 @@ public: * membrane strain, bending strain, Jacobian matrices of deformed and reference geometries, Hessian matrices of * deformed and reference geometries, the normal vector, and the normalized normal vector. * - * \tparam ScalarType The scalar type for the matrix and vector elements. + * \tparam ST The scalar type for the matrix and vector elements. */ - template + template struct KinematicVariables { - Eigen::Matrix C; ///< material tangent - Eigen::Vector3 epsV; ///< membrane strain in Voigt notation - Eigen::Vector3 kappaV; ///< bending strain in Voigt notation - Eigen::Matrix j; ///< Jacobian of the deformed geometry - Eigen::Matrix J; ///< Jacobian of the reference geometry - Eigen::Matrix3 h; ///< Hessian of the deformed geometry - Eigen::Matrix3 H; ///< Hessian of the reference geometry - Eigen::Vector3 a3N; ///< Normal vector of the deformed geometry - Eigen::Vector3 a3; ///< normalized normal vector of the deformed geometry + Eigen::Matrix C; ///< material tangent + Eigen::Vector3 epsV; ///< membrane strain in Voigt notation + Eigen::Vector3 kappaV; ///< bending strain in Voigt notation + Eigen::Matrix j; ///< Jacobian of the deformed geometry + Eigen::Matrix J; ///< Jacobian of the reference geometry + Eigen::Matrix3 h; ///< Hessian of the deformed geometry + Eigen::Matrix3 H; ///< Hessian of the reference geometry + Eigen::Vector3 a3N; ///< Normal vector of the deformed geometry + Eigen::Vector3 a3; ///< normalized normal vector of the deformed geometry }; /** - * @brief Constructor for the KirchhoffLoveShell class. + * \brief Constructor for the KirchhoffLoveShell class. * * Initializes the KirchhoffLoveShell instance with the given parameters. * - * @tparam VolumeLoad The type representing the volume load function. - * @tparam NeumannBoundaryLoad The type representing the Neumann boundary load function. - * @param globalBasis The global basis for the finite element. - * @param element The local element to bind. - * @param emod Young's modulus of the material. - * @param nu Poisson's ratio of the material. - * @param thickness Thickness of the shell. - * @param p_volumeLoad The volume load function (optional, default is utils::LoadDefault). - * @param p_neumannBoundary The Neumann boundary patch (optional, default is nullptr). - * @param p_neumannBoundaryLoad The Neumann boundary load function (optional, default is LoadDefault). + * \tparam VolumeLoad The type representing the volume load function. + * \tparam NeumannBoundaryLoad The type representing the Neumann boundary load function. + * \param globalBasis The global basis for the finite element. + * \param element The local element to bind. + * \param emod Young's modulus of the material. + * \param nu Poisson's ratio of the material. + * \param thickness Thickness of the shell. + * \param volumeLoad The volume load function (optional, default is utils::LoadDefault). + * \param neumannBoundary The Neumann boundary patch (optional, default is nullptr). + * \param neumannBoundaryLoad The Neumann boundary load function (optional, default is LoadDefault). */ template KirchhoffLoveShell(const Basis& globalBasis, const typename LocalView::Element& element, double emod, double nu, - double thickness, VolumeLoad p_volumeLoad = {}, - const BoundaryPatch* p_neumannBoundary = nullptr, - NeumannBoundaryLoad p_neumannBoundaryLoad = {}) + double thickness, VolumeLoad volumeLoad = {}, + const BoundaryPatch* neumannBoundary = nullptr, + NeumannBoundaryLoad neumannBoundaryLoad = {}) : BasePowerFE(globalBasis, element), - VolumeType(p_volumeLoad), - TractionType(p_neumannBoundary, p_neumannBoundaryLoad), + VolumeType(volumeLoad), + TractionType(neumannBoundary, neumannBoundaryLoad), emod_{emod}, nu_{nu}, thickness_{thickness} { this->localView().bind(element); - auto& first_child = this->localView().tree().child(0); - const auto& fe = first_child.finiteElement(); - geo_ = std::make_shared(this->localView().element().geometry()); - numberOfNodes_ = fe.size(); - order_ = 2 * (fe.localBasis().order()); - localBasis = Dune::CachedLocalBasis(fe.localBasis()); + auto& firstChild = this->localView().tree().child(0); + const auto& fe = firstChild.finiteElement(); + geo_ = std::make_shared(this->localView().element().geometry()); + numberOfNodes_ = fe.size(); + order_ = 2 * (fe.localBasis().order()); + localBasis_ = Dune::CachedLocalBasis(fe.localBasis()); if constexpr (requires { this->localView().element().impl().getQuadratureRule(order_); }) if (this->localView().element().impl().isTrimmed()) - localBasis.bind(this->localView().element().impl().getQuadratureRule(order_), Dune::bindDerivatives(0, 1, 2)); + localBasis_.bind(this->localView().element().impl().getQuadratureRule(order_), Dune::bindDerivatives(0, 1, 2)); else - localBasis.bind(Dune::QuadratureRules::rule(this->localView().element().type(), order_), - Dune::bindDerivatives(0, 1, 2)); + localBasis_.bind(Dune::QuadratureRules::rule(this->localView().element().type(), order_), + Dune::bindDerivatives(0, 1, 2)); else - localBasis.bind(Dune::QuadratureRules::rule(this->localView().element().type(), order_), - Dune::bindDerivatives(0, 1, 2)); + localBasis_.bind(Dune::QuadratureRules::rule(this->localView().element().type(), order_), + Dune::bindDerivatives(0, 1, 2)); } public: /** - * @brief Get the displacement function and nodal displacements. + * \brief Get the displacement function and nodal displacements. * * Retrieves the displacement function and nodal displacements based on the given FERequirements. * - * @tparam ScalarType The scalar type used for calculations. - * @param par The FERequirements. - * @param dx Optional additional displacement vector. - * @return The displacement function. + * \tparam ST The scalar type used for calculations. + * \param par The FERequirements. + * \param dx Optional additional displacement vector. + * \return The displacement function. */ - template + template auto displacementFunction(const FERequirementType& par, - const std::optional>& dx = std::nullopt) const { + const std::optional>& dx = std::nullopt) const { const auto& d = par.getGlobalSolution(Ikarus::FESolutions::displacement); auto disp = Ikarus::FEHelper::localSolutionBlockVector(d, this->localView(), dx); - Dune::StandardLocalFunction uFunction(localBasis, disp, geo_); + Dune::StandardLocalFunction uFunction(localBasis_, disp, geo_); return uFunction; } @@ -154,44 +152,44 @@ public: [[nodiscard]] int order() const { return order_; } /** - * @brief Calculate the scalar value. + * \brief Calculate the scalar value. * * Calculates the scalar value based on the given FERequirements. * - * @param req The FERequirements. - * @return The calculated scalar value. + * \param req The FERequirements. + * \return The calculated scalar value. */ double calculateScalar(const FERequirementType& req) const { return calculateScalarImpl(req); } /** - * @brief Calculate the vector associated with the given FERequirementType. + * \brief Calculate the vector associated with the given FERequirementType. * - * @tparam ScalarType The scalar type for the calculation. - * @param req The FERequirementType object specifying the requirements for the calculation. - * @param force The vector to store the calculated result. + * \tparam ST The scalar type for the calculation. + * \param req The FERequirementType object specifying the requirements for the calculation. + * \param force The vector to store the calculated result. */ void calculateVector(const FERequirementType& req, typename Traits::template VectorType<> force) const { calculateVectorImpl(req, force); } /** - * @brief Calculate the matrix associated with the given FERequirementType. + * \brief Calculate the matrix associated with the given FERequirementType. * - * @tparam ScalarType The scalar type for the calculation. - * @param req The FERequirementType object specifying the requirements for the calculation. - * @param K The matrix to store the calculated result. + * \tparam ST The scalar type for the calculation. + * \param req The FERequirementType object specifying the requirements for the calculation. + * \param K The matrix to store the calculated result. */ void calculateMatrix(const FERequirementType& req, typename Traits::template MatrixType<> K) const { calculateMatrixImpl(req, K); } /** - * @brief Calculates a requested result at a specific local position. + * \brief Calculates a requested result at a specific local position. * - * @param req The FERequirementType object holding the global solution. - * @param local Local position vector. - * @return calculated result + * \param req The FERequirementType object holding the global solution. + * \param local Local position vector. + * \return calculated result * - * @tparam resType The type representing the requested result. + * \tparam resType The type representing the requested result. */ template auto calculateAt([[maybe_unused]] const FERequirementType& req, @@ -201,8 +199,8 @@ public: private: std::shared_ptr geo_; - Dune::CachedLocalBasis> localBasis; - DefaultMembraneStrain membraneStrain; + Dune::CachedLocalBasis> localBasis_; + DefaultMembraneStrain membraneStrain_; double emod_; double nu_; double thickness_; @@ -231,9 +229,9 @@ protected: */ auto computeMaterialAndStrains(const Dune::FieldVector& gpPos, int gpIndex, const Geometry& geo, const auto& uFunction) const { - using ScalarType = typename std::remove_cvref_t::ctype; + using ST = typename std::remove_cvref_t::ctype; - KinematicVariables kin; + KinematicVariables kin; using namespace Dune; using namespace Dune::DerivativeDirections; const auto [X, Jd, Hd] = geo_->impl().zeroFirstAndSecondDerivativeOfPosition(gpPos); @@ -247,27 +245,27 @@ protected: const Eigen::Matrix GInv = G.inverse(); kin.C = materialTangent(GInv); - kin.epsV = membraneStrain.value(gpPos, geo, uFunction); + kin.epsV = membraneStrain_.value(gpPos, geo, uFunction); - const auto& Ndd = localBasis.evaluateSecondDerivatives(gpIndex); - const auto uasMatrix = Dune::viewAsEigenMatrixAsDynFixed(uFunction.coefficientsRef()); - const auto hessianu = Ndd.transpose().template cast() * uasMatrix; - kin.h = kin.H + hessianu; - const Eigen::Matrix gradu = toEigen(uFunction.evaluateDerivative( + const auto& Ndd = localBasis_.evaluateSecondDerivatives(gpIndex); + const auto uasMatrix = Dune::viewAsEigenMatrixAsDynFixed(uFunction.coefficientsRef()); + const auto hessianu = Ndd.transpose().template cast() * uasMatrix; + kin.h = kin.H + hessianu; + const Eigen::Matrix gradu = toEigen(uFunction.evaluateDerivative( gpIndex, Dune::wrt(spatialAll), Dune::on(Dune::DerivativeDirections::referenceElement))); - kin.j = kin.J + gradu.transpose(); - kin.a3N = (kin.j.row(0).cross(kin.j.row(1))); - kin.a3 = kin.a3N.normalized(); - Eigen::Vector bV = kin.h * kin.a3; + kin.j = kin.J + gradu.transpose(); + kin.a3N = (kin.j.row(0).cross(kin.j.row(1))); + kin.a3 = kin.a3N.normalized(); + Eigen::Vector bV = kin.h * kin.a3; bV(2) *= 2; // Voigt notation requires the two here const auto BV = toVoigt(toEigen(geo_->impl().secondFundamentalForm(gpPos))); kin.kappaV = BV - bV; return kin; } - template - void calculateMatrixImpl(const FERequirementType& par, typename Traits::template MatrixType K, - const std::optional>& dx = std::nullopt) const { + template + void calculateMatrixImpl(const FERequirementType& par, typename Traits::template MatrixType K, + const std::optional>& dx = std::nullopt) const { K.setZero(); using namespace Dune::DerivativeDirections; using namespace Dune; @@ -278,27 +276,27 @@ protected: const auto intElement = geo_->integrationElement(gp.position()) * gp.weight(); const auto [C, epsV, kappaV, jE, J, h, H, a3N, a3] = computeMaterialAndStrains(gp.position(), gpIndex, *geo_, uFunction); - const Eigen::Vector membraneForces = thickness_ * C * epsV; - const Eigen::Vector moments = Dune::power(thickness_, 3) / 12.0 * C * kappaV; + const Eigen::Vector membraneForces = thickness_ * C * epsV; + const Eigen::Vector moments = Dune::power(thickness_, 3) / 12.0 * C * kappaV; - const auto& Nd = localBasis.evaluateJacobian(gpIndex); - const auto& Ndd = localBasis.evaluateSecondDerivatives(gpIndex); + const auto& Nd = localBasis_.evaluateJacobian(gpIndex); + const auto& Ndd = localBasis_.evaluateSecondDerivatives(gpIndex); for (size_t i = 0; i < numberOfNodes_; ++i) { - Eigen::Matrix bopIMembrane = - membraneStrain.derivative(gp.position(), jE, Nd, *geo_, uFunction, localBasis, i); + Eigen::Matrix bopIMembrane = + membraneStrain_.derivative(gp.position(), jE, Nd, *geo_, uFunction, localBasis_, i); - Eigen::Matrix bopIBending = bopBending(jE, h, Nd, Ndd, i, a3N, a3); + Eigen::Matrix bopIBending = bopBending(jE, h, Nd, Ndd, i, a3N, a3); for (size_t j = i; j < numberOfNodes_; ++j) { auto KBlock = K.template block(worldDim * i, worldDim * j); - Eigen::Matrix bopJMembrane = - membraneStrain.derivative(gp.position(), jE, Nd, *geo_, uFunction, localBasis, j); - Eigen::Matrix bopJBending = bopBending(jE, h, Nd, Ndd, j, a3N, a3); + Eigen::Matrix bopJMembrane = + membraneStrain_.derivative(gp.position(), jE, Nd, *geo_, uFunction, localBasis_, j); + Eigen::Matrix bopJBending = bopBending(jE, h, Nd, Ndd, j, a3N, a3); KBlock += thickness_ * bopIMembrane.transpose() * C * bopJMembrane * intElement; KBlock += Dune::power(thickness_, 3) / 12.0 * bopIBending.transpose() * C * bopJBending * intElement; - Eigen::Matrix kgMembraneIJ = - membraneStrain.secondDerivative(gp.position(), Nd, *geo_, uFunction, localBasis, membraneForces, i, j); - Eigen::Matrix kgBendingIJ = kgBending(jE, h, Nd, Ndd, a3N, a3, moments, i, j); + Eigen::Matrix kgMembraneIJ = + membraneStrain_.secondDerivative(gp.position(), Nd, *geo_, uFunction, localBasis_, membraneForces, i, j); + Eigen::Matrix kgBendingIJ = kgBending(jE, h, Nd, Ndd, a3N, a3, moments, i, j); KBlock += kgMembraneIJ * intElement; KBlock += kgBendingIJ * intElement; } @@ -311,9 +309,9 @@ protected: TractionType::calculateMatrixImpl(par, K, dx); } - template - void calculateVectorImpl(const FERequirementType& par, typename Traits::template VectorType force, - const std::optional>& dx = std::nullopt) const { + template + void calculateVectorImpl(const FERequirementType& par, typename Traits::template VectorType force, + const std::optional>& dx = std::nullopt) const { force.setZero(); using namespace Dune::DerivativeDirections; using namespace Dune; @@ -324,15 +322,15 @@ protected: for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { const auto [C, epsV, kappaV, jE, J, h, H, a3N, a3] = computeMaterialAndStrains(gp.position(), gpIndex, *geo_, uFunction); - const Eigen::Vector membraneForces = thickness_ * C * epsV; - const Eigen::Vector moments = Dune::power(thickness_, 3) / 12.0 * C * kappaV; + const Eigen::Vector membraneForces = thickness_ * C * epsV; + const Eigen::Vector moments = Dune::power(thickness_, 3) / 12.0 * C * kappaV; - const auto& Nd = localBasis.evaluateJacobian(gpIndex); - const auto& Ndd = localBasis.evaluateSecondDerivatives(gpIndex); + const auto& Nd = localBasis_.evaluateJacobian(gpIndex); + const auto& Ndd = localBasis_.evaluateSecondDerivatives(gpIndex); for (size_t i = 0; i < numberOfNodes_; ++i) { - Eigen::Matrix bopIMembrane = - membraneStrain.derivative(gp.position(), jE, Nd, *geo_, uFunction, localBasis, i); - Eigen::Matrix bopIBending = bopBending(jE, h, Nd, Ndd, i, a3N, a3); + Eigen::Matrix bopIMembrane = + membraneStrain_.derivative(gp.position(), jE, Nd, *geo_, uFunction, localBasis_, i); + Eigen::Matrix bopIBending = bopBending(jE, h, Nd, Ndd, i, a3N, a3); force.template segment<3>(3 * i) += bopIMembrane.transpose() * membraneForces * geo_->integrationElement(gp.position()) * gp.weight(); force.template segment<3>(3 * i) += @@ -347,22 +345,21 @@ protected: TractionType::calculateVectorImpl(par, force, dx); } - template + template auto calculateScalarImpl(const FERequirementType& par, - const std::optional>& dx = std::nullopt) const - -> ScalarType { + const std::optional>& dx = std::nullopt) const -> ST { using namespace Dune::DerivativeDirections; using namespace Dune; const auto uFunction = displacementFunction(par, dx); const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); - ScalarType energy = 0.0; + ST energy = 0.0; for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { const auto [C, epsV, kappaV, j, J, h, H, a3N, a3] = computeMaterialAndStrains(gp.position(), gpIndex, *geo_, uFunction); - const ScalarType membraneEnergy = 0.5 * thickness_ * epsV.dot(C * epsV); - const ScalarType bendingEnergy = 0.5 * Dune::power(thickness_, 3) / 12.0 * kappaV.dot(C * kappaV); + const ST membraneEnergy = 0.5 * thickness_ * epsV.dot(C * epsV); + const ST bendingEnergy = 0.5 * Dune::power(thickness_, 3) / 12.0 * kappaV.dot(C * kappaV); energy += (membraneEnergy + bendingEnergy) * geo_->integrationElement(gp.position()) * gp.weight(); } @@ -374,12 +371,11 @@ protected: return energy; } - template - Eigen::Matrix kgBending(const Eigen::Matrix& jcur, - const Eigen::Matrix3& h, const auto& dN, const auto& ddN, - const Eigen::Vector3& a3N, const Eigen::Vector3& a3, - const Eigen::Vector3& S, int I, int J) const { - Eigen::Matrix kg; + template + Eigen::Matrix kgBending(const Eigen::Matrix& jcur, const Eigen::Matrix3& h, const auto& dN, + const auto& ddN, const Eigen::Vector3& a3N, const Eigen::Vector3& a3, + const Eigen::Vector3& S, int I, int J) const { + Eigen::Matrix kg; kg.setZero(); const auto& dN1i = dN(I, 0); @@ -387,33 +383,33 @@ protected: const auto& dN2i = dN(I, 1); const auto& dN2j = dN(J, 1); - const Eigen::Matrix P = + const Eigen::Matrix P = 1.0 / a3N.norm() * (Eigen::Matrix::Identity() - a3 * a3.transpose()); const auto a1dxI = Eigen::Matrix::Identity() * dN1i; // the auto here allows the exploitation of the identity matrices, // due to Eigen's expression templates - const auto a2dxI = Eigen::Matrix::Identity() * dN2i; - const auto a1dxJ = Eigen::Matrix::Identity() * dN1j; - const auto a2dxJ = Eigen::Matrix::Identity() * dN2j; - const auto a1 = jcur.row(0); - const auto a2 = jcur.row(1); - const Eigen::Matrix a3NdI = a1dxI.colwise().cross(a2) - a2dxI.colwise().cross(a1); - const Eigen::Matrix a3NdJ = a1dxJ.colwise().cross(a2) - a2dxJ.colwise().cross(a1); - Eigen::Matrix a3dI = P * a3NdI; - Eigen::Matrix a3dJ = P * a3NdJ; + const auto a2dxI = Eigen::Matrix::Identity() * dN2i; + const auto a1dxJ = Eigen::Matrix::Identity() * dN1j; + const auto a2dxJ = Eigen::Matrix::Identity() * dN2j; + const auto a1 = jcur.row(0); + const auto a2 = jcur.row(1); + const Eigen::Matrix a3NdI = a1dxI.colwise().cross(a2) - a2dxI.colwise().cross(a1); + const Eigen::Matrix a3NdJ = a1dxJ.colwise().cross(a2) - a2dxJ.colwise().cross(a1); + Eigen::Matrix a3dI = P * a3NdI; + Eigen::Matrix a3dJ = P * a3NdJ; for (int i = 0; i < 3; ++i) { - const auto a_albe = h.row(i).transpose(); - const auto& ddNI = ddN(I, i); - const auto& ddNJ = ddN(J, i); - Eigen::Vector3 vecd = P * a_albe; + const auto a_albe = h.row(i).transpose(); + const auto& ddNI = ddN(I, i); + const auto& ddNJ = ddN(J, i); + Eigen::Vector3 vecd = P * a_albe; - Eigen::Matrix a3Ndd = + Eigen::Matrix a3Ndd = 1.0 / a3N.squaredNorm() * ((3 * a3 * a3.transpose() - Eigen::Matrix::Identity()) * (a3.dot(a_albe)) - a_albe * a3.transpose() - a3 * a_albe.transpose()); - Eigen::Matrix secondDerivativeDirectorIJ = skew(((dN2i * dN1j - dN1i * dN2j) * vecd).eval()); + Eigen::Matrix secondDerivativeDirectorIJ = skew(((dN2i * dN1j - dN1i * dN2j) * vecd).eval()); kg -= (a3NdI.transpose() * a3Ndd * a3NdJ + secondDerivativeDirectorIJ + (ddNI * a3dJ + ddNJ * a3dI.transpose())) * S[i] * (i == 2 ? 2 : 1); } @@ -421,24 +417,23 @@ protected: return kg; } - template - Eigen::Matrix bopBending(const Eigen::Matrix& jcur, - const Eigen::Matrix3& h, const auto& dN, const auto& ddN, - const int node, const Eigen::Vector3& a3N, - const Eigen::Vector3& a3) const { - const Eigen::Matrix a1dxI = + template + Eigen::Matrix bopBending(const Eigen::Matrix& jcur, const Eigen::Matrix3& h, const auto& dN, + const auto& ddN, const int node, const Eigen::Vector3& a3N, + const Eigen::Vector3& a3) const { + const Eigen::Matrix a1dxI = Eigen::Matrix::Identity() * dN(node, 0); // this should be double // but the cross-product below complains otherwise - const Eigen::Matrix a2dxI = Eigen::Matrix::Identity() * dN(node, 1); - const auto a1 = jcur.row(0); - const auto a2 = jcur.row(1); - const Eigen::Matrix a3NdI = + const Eigen::Matrix a2dxI = Eigen::Matrix::Identity() * dN(node, 1); + const auto a1 = jcur.row(0); + const auto a2 = jcur.row(1); + const Eigen::Matrix a3NdI = a1dxI.colwise().cross(a2) - a2dxI.colwise().cross(a1); // minus needed since order has // to be swapped to get column-wise cross product working - const Eigen::Matrix a3d1 = + const Eigen::Matrix a3d1 = 1.0 / a3N.norm() * (Eigen::Matrix::Identity() - a3 * a3.transpose()) * a3NdI; - Eigen::Matrix bop = -(h * a3d1 + (a3 * ddN.row(node)).transpose()); + Eigen::Matrix bop = -(h * a3d1 + (a3 * ddN.row(node)).transpose()); bop.row(2) *= 2; return bop; diff --git a/ikarus/finiteelements/mechanics/linearelastic.hh b/ikarus/finiteelements/mechanics/linearelastic.hh index bd0b1e6f5..8ba60a696 100644 --- a/ikarus/finiteelements/mechanics/linearelastic.hh +++ b/ikarus/finiteelements/mechanics/linearelastic.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file linearelastic.hh - * @brief Definition of the LinearElastic class for finite element mechanics computations. - * @ingroup mechanics + * \file linearelastic.hh + * \brief Definition of the LinearElastic class for finite element mechanics computations. + * \ingroup mechanics */ #pragma once @@ -34,23 +34,21 @@ namespace Ikarus { /** - * @brief LinearElastic class represents a linear elastic finite element. + * \brief LinearElastic class represents a linear elastic finite element. * - * @ingroup mechanics + * \ingroup mechanics * - * @tparam Basis_ The basis type for the finite element. - * @tparam FERequirements_ The requirements for the finite element. - * @tparam useEigenRef A boolean flag indicating whether to use Eigen references. + * \tparam B The basis type for the finite element. + * \tparam FER The requirements for the finite element. + * \tparam useEigenRef A boolean flag indicating whether to use Eigen references. */ -template , bool useEigenRef = false> -class LinearElastic : public PowerBasisFE, - public Volume, - FETraits>, - public Traction, - FETraits> +template , bool useEigenRef = false> +class LinearElastic : public PowerBasisFE, + public Volume, FETraits>, + public Traction, FETraits> { public: - using Traits = FETraits; + using Traits = FETraits; using Basis = typename Traits::Basis; using FlatBasis = typename Traits::FlatBasis; using FERequirementType = typename Traits::FERequirementType; @@ -59,73 +57,73 @@ public: using GridView = typename Traits::GridView; using Element = typename Traits::Element; using BaseDisp = PowerBasisFE; // Handles globalIndices function - using VolumeType = Volume, Traits>; - using TractionType = Traction, Traits>; + using VolumeType = Volume; + using TractionType = Traction; static constexpr int myDim = Traits::mydim; using LocalBasisType = decltype(std::declval().tree().child(0).finiteElement().localBasis()); /** - * @brief Constructor for the LinearElastic class. + * \brief Constructor for the LinearElastic class. * - * @tparam VolumeLoad The type for the volume load function. - * @tparam NeumannBoundaryLoad The type for the Neumann boundary load function. - * @param globalBasis The global basis for the finite element. - * @param element The element for which the finite element is constructed. - * @param emod Young's modulus. - * @param nu Poisson's ratio. - * @param p_volumeLoad Volume load function (default is LoadDefault). - * @param p_neumannBoundary Neumann boundary patch (default is nullptr). - * @param p_neumannBoundaryLoad Neumann boundary load function (default is LoadDefault). + * \tparam VolumeLoad The type for the volume load function. + * \tparam NeumannBoundaryLoad The type for the Neumann boundary load function. + * \param globalBasis The global basis for the finite element. + * \param element The element for which the finite element is constructed. + * \param emod Young's modulus. + * \param nu Poisson's ratio. + * \param volumeLoad Volume load function (default is LoadDefault). + * \param neumannBoundary Neumann boundary patch (default is nullptr). + * \param neumannBoundaryLoad Neumann boundary load function (default is LoadDefault). */ template LinearElastic(const Basis& globalBasis, const typename LocalView::Element& element, double emod, double nu, - VolumeLoad p_volumeLoad = {}, const BoundaryPatch* p_neumannBoundary = nullptr, - NeumannBoundaryLoad p_neumannBoundaryLoad = {}) + VolumeLoad volumeLoad = {}, const BoundaryPatch* neumannBoundary = nullptr, + NeumannBoundaryLoad neumannBoundaryLoad = {}) : BaseDisp(globalBasis, element), - VolumeType(p_volumeLoad), - TractionType(p_neumannBoundary, p_neumannBoundaryLoad), + VolumeType(volumeLoad), + TractionType(neumannBoundary, neumannBoundaryLoad), emod_{emod}, nu_{nu} { this->localView().bind(element); - auto& first_child = this->localView().tree().child(0); - const auto& fe = first_child.finiteElement(); - geo_ = std::make_shared(this->localView().element().geometry()); - numberOfNodes_ = fe.size(); - order_ = 2 * (fe.localBasis().order()); - localBasis = Dune::CachedLocalBasis(fe.localBasis()); + auto& firstChild = this->localView().tree().child(0); + const auto& fe = firstChild.finiteElement(); + geo_ = std::make_shared(this->localView().element().geometry()); + numberOfNodes_ = fe.size(); + order_ = 2 * (fe.localBasis().order()); + localBasis_ = Dune::CachedLocalBasis(fe.localBasis()); if constexpr (requires { this->localView().element().impl().getQuadratureRule(order_); }) if (this->localView().element().impl().isTrimmed()) - localBasis.bind(this->localView().element().impl().getQuadratureRule(order_), Dune::bindDerivatives(0, 1)); + localBasis_.bind(this->localView().element().impl().getQuadratureRule(order_), Dune::bindDerivatives(0, 1)); else - localBasis.bind(Dune::QuadratureRules::rule(this->localView().element().type(), order_), - Dune::bindDerivatives(0, 1)); + localBasis_.bind(Dune::QuadratureRules::rule(this->localView().element().type(), order_), + Dune::bindDerivatives(0, 1)); else - localBasis.bind(Dune::QuadratureRules::rule(this->localView().element().type(), order_), - Dune::bindDerivatives(0, 1)); + localBasis_.bind(Dune::QuadratureRules::rule(this->localView().element().type(), order_), + Dune::bindDerivatives(0, 1)); } /** - * @brief Gets the displacement function for the given FERequirementType and optional displacement vector. + * \brief Gets the displacement function for the given FERequirementType and optional displacement vector. * - * @tparam ScalarType The scalar type for the displacement vector. - * @param par The FERequirementType object. - * @param dx Optional displacement vector. - * @return The displacement function. + * \tparam ScalarType The scalar type for the displacement vector. + * \param par The FERequirementType object. + * \param dx Optional displacement vector. + * \return The displacement function. */ template auto displacementFunction(const FERequirementType& par, const std::optional>& dx = std::nullopt) const { const auto& d = par.getGlobalSolution(Ikarus::FESolutions::displacement); auto disp = Ikarus::FEHelper::localSolutionBlockVector(d, this->localView(), dx); - Dune::StandardLocalFunction uFunction(localBasis, disp, geo_); + Dune::StandardLocalFunction uFunction(localBasis_, disp, geo_); return uFunction; } /** - * @brief Gets the strain function for the given FERequirementType and optional displacement vector. + * \brief Gets the strain function for the given FERequirementType and optional displacement vector. * - * @tparam ScalarType The scalar type for the strain vector. - * @param par The FERequirementType object. - * @param dx Optional displacement vector. - * @return The strain function. + * \tparam ScalarType The scalar type for the strain vector. + * \param par The FERequirementType object. + * \param dx Optional displacement vector. + * \return The strain function. */ template auto strainFunction(const FERequirementType& par, @@ -134,9 +132,9 @@ public: } /** - * @brief Gets the material tangent matrix for the linear elastic material. + * \brief Gets the material tangent matrix for the linear elastic material. * - * @return The material tangent matrix. + * \return The material tangent matrix. */ auto materialTangent() const { if constexpr (myDim == 2) @@ -146,10 +144,10 @@ public: } /** - * @brief Gets the material tangent function for the given FERequirementType. + * \brief Gets the material tangent function for the given FERequirementType. * - * @param par The FERequirementType object. - * @return The material tangent function. + * \param par The FERequirementType object. + * \return The material tangent function. */ auto materialTangentFunction([[maybe_unused]] const FERequirementType& par) const { return [&]([[maybe_unused]] auto gp) { return materialTangent(); }; @@ -160,28 +158,28 @@ public: [[nodiscard]] int order() const { return order_; } /** - * @brief Calculates the scalar energy for the given FERequirementType. + * \brief Calculates the scalar energy for the given FERequirementType. * - * @param par The FERequirementType object. - * @return The scalar energy. + * \param par The FERequirementType object. + * \return The scalar energy. */ inline double calculateScalar(const FERequirementType& par) const { return calculateScalarImpl(par); } /** - * @brief Calculates the vector force for the given FERequirementType. + * \brief Calculates the vector force for the given FERequirementType. * - * @param par The FERequirementType object. - * @param force Vector to store the calculated force. + * \param par The FERequirementType object. + * \param force Vector to store the calculated force. */ inline void calculateVector(const FERequirementType& par, typename Traits::template VectorType<> force) const { calculateVectorImpl(par, force); } /** - * @brief Calculates the matrix stiffness for the given FERequirementType. + * \brief Calculates the matrix stiffness for the given FERequirementType. * - * @param par The FERequirementType object. - * @param K Matrix to store the calculated stiffness. + * \param par The FERequirementType object. + * \param K Matrix to store the calculated stiffness. */ void calculateMatrix(const FERequirementType& par, typename Traits::template MatrixType<> K) const { const auto eps = strainFunction(par); @@ -206,13 +204,13 @@ public: } /** - * @brief Calculates a requested result at a specific local position. + * \brief Calculates a requested result at a specific local position. * - * @param req The FERequirementType object holding the global solution. - * @param local Local position vector. - * @return calculated result + * \param req The FERequirementType object holding the global solution. + * \param local Local position vector. + * \return calculated result * - * @tparam resType The type representing the requested result. + * \tparam resType The type representing the requested result. */ template auto calculateAt(const FERequirementType& req, const Dune::FieldVector& local) const { @@ -229,7 +227,7 @@ public: private: std::shared_ptr geo_; - Dune::CachedLocalBasis> localBasis; + Dune::CachedLocalBasis> localBasis_; double emod_; double nu_; size_t numberOfNodes_{0}; diff --git a/ikarus/finiteelements/mechanics/loads.hh b/ikarus/finiteelements/mechanics/loads.hh index 181526d82..b1b3d9d99 100644 --- a/ikarus/finiteelements/mechanics/loads.hh +++ b/ikarus/finiteelements/mechanics/loads.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file loads.hh - * @brief Header file for types of loads in Ikarus finite element mechanics. - * @ingroup mechanics + * \file loads.hh + * \brief Header file for types of loads in Ikarus finite element mechanics. + * \ingroup mechanics */ #pragma once diff --git a/ikarus/finiteelements/mechanics/loads/traction.hh b/ikarus/finiteelements/mechanics/loads/traction.hh index 961eb4800..d8f89846e 100644 --- a/ikarus/finiteelements/mechanics/loads/traction.hh +++ b/ikarus/finiteelements/mechanics/loads/traction.hh @@ -11,14 +11,14 @@ namespace Ikarus { /** - * @brief Traction class represents distributed traction load that can be applied. + * \brief Traction class represents distributed traction load that can be applied. * - * @ingroup mechanics + * \ingroup mechanics * - * @tparam DisplacementBasedElement The type of the displacement-based finite element. - * @tparam Traits Type of traits for handling finite elements. + * \tparam DFE The type of the displacement-based finite element. + * \tparam Traits Type of traits for handling finite elements. */ -template +template class Traction { public: @@ -29,67 +29,66 @@ public: static constexpr int worldDim = Traits::worlddim; /** - * @brief Constructor for the Loads class. + * \brief Constructor for the Loads class. * - * @tparam NeumannBoundaryLoad The type for the Neumann boundary load function. - * @param p_neumannBoundary Neumann boundary patch. - * @param p_neumannBoundaryLoad Neumann boundary load function. + * \tparam NBL The type for the Neumann boundary load function. + * \param neumannBoundary Neumann boundary patch. + * \param neumannBoundaryLoad Neumann boundary load function. */ - template - explicit Traction(const BoundaryPatch* p_neumannBoundary, NeumannBoundaryLoad p_neumannBoundaryLoad) - : neumannBoundary{p_neumannBoundary} { - if constexpr (!std::is_same_v) - neumannBoundaryLoad = p_neumannBoundaryLoad; + template + explicit Traction(const BoundaryPatch* neumannBoundary, NBL neumannBoundaryLoad) + : neumannBoundary_{neumannBoundary} { + if constexpr (!std::is_same_v) + neumannBoundaryLoad_ = neumannBoundaryLoad; - assert(((not p_neumannBoundary and not neumannBoundaryLoad) or (p_neumannBoundary and neumannBoundaryLoad)) && + assert(((not neumannBoundary and not neumannBoundaryLoad_) or (neumannBoundary and neumannBoundaryLoad_)) && "If you pass a Neumann boundary you should also pass the function for the Neumann load!"); } /** - * @brief Calculate the scalar value. + * \brief Calculate the scalar value. * * Calculates the scalar value based on the given FERequirements. * - * @param req The FERequirements. - * @return The calculated scalar value. + * \param req The FERequirements. + * \return The calculated scalar value. */ double calculateScalar(const FERequirementType& req) const { return calculateScalarImpl(req); } /** - * @brief Calculate the vector associated with the given FERequirementType. + * \brief Calculate the vector associated with the given FERequirementType. * - * @tparam ScalarType The scalar type for the calculation. - * @param req The FERequirementType object specifying the requirements for the calculation. - * @param force The vector to store the calculated result. + * \tparam ScalarType The scalar type for the calculation. + * \param req The FERequirementType object specifying the requirements for the calculation. + * \param force The vector to store the calculated result. */ void calculateVector(const FERequirementType& req, typename Traits::template VectorType<> force) const { calculateVectorImpl(req, force); } /** - * @brief Calculate the matrix associated with the given FERequirementType. + * \brief Calculate the matrix associated with the given FERequirementType. * - * @tparam ScalarType The scalar type for the calculation. - * @param req The FERequirementType object specifying the requirements for the calculation. - * @param K The matrix to store the calculated result. + * \tparam ScalarType The scalar type for the calculation. + * \param req The FERequirementType object specifying the requirements for the calculation. + * \param K The matrix to store the calculated result. */ void calculateMatrix(const FERequirementType& req, typename Traits::template MatrixType<> K) const { calculateMatrixImpl(req, K); } protected: - template + template auto calculateScalarImpl(const FERequirementType& par, - const std::optional>& dx = std::nullopt) const - -> ScalarType { - if (not neumannBoundary and not neumannBoundaryLoad) + const std::optional>& dx = std::nullopt) const -> ST { + if (not neumannBoundary_ and not neumannBoundaryLoad_) return 0.0; - ScalarType energy = 0.0; + ST energy = 0.0; const auto uFunction = dbElement().displacementFunction(par, dx); const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); auto& element = dbElement().localView().element(); - for (auto&& intersection : intersections(neumannBoundary->gridView(), element)) { - if (not neumannBoundary->contains(intersection)) + for (auto&& intersection : intersections(neumannBoundary_->gridView(), element)) { + if (not neumannBoundary_->contains(intersection)) continue; const auto& quadLine = Dune::QuadratureRules::rule(intersection.type(), dbElement().order()); @@ -104,7 +103,7 @@ protected: const auto uVal = uFunction.evaluate(quadPos); // Value of the Neumann data at the current position - auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda); + auto neumannValue = neumannBoundaryLoad_(intersection.geometry().global(curQuad.position()), lambda); energy -= neumannValue.dot(uVal) * curQuad.weight() * integrationElement; } @@ -112,10 +111,10 @@ protected: return energy; } - template - void calculateVectorImpl(const FERequirementType& par, typename Traits::template VectorType force, - const std::optional> dx = std::nullopt) const { - if (not neumannBoundary and not neumannBoundaryLoad) + template + void calculateVectorImpl(const FERequirementType& par, typename Traits::template VectorType force, + const std::optional> dx = std::nullopt) const { + if (not neumannBoundary_ and not neumannBoundaryLoad_) return; using namespace Dune::DerivativeDirections; using namespace Dune; @@ -123,8 +122,8 @@ protected: const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); auto& element = dbElement().localView().element(); - for (auto&& intersection : intersections(neumannBoundary->gridView(), element)) { - if (not neumannBoundary->contains(intersection)) + for (auto&& intersection : intersections(neumannBoundary_->gridView(), element)) { + if (not neumannBoundary_->contains(intersection)) continue; /// Integration rule along the boundary @@ -139,29 +138,27 @@ protected: const auto udCi = uFunction.evaluateDerivative(quadPos, wrt(coeff(i))); /// Value of the Neumann data at the current position - auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda); + auto neumannValue = neumannBoundaryLoad_(intersection.geometry().global(curQuad.position()), lambda); force.template segment(worldDim * i) -= udCi * neumannValue * curQuad.weight() * intElement; } } } } - template + template void calculateMatrixImpl(const FERequirementType& par, typename Traits::template MatrixType<> K, - const std::optional>& dx = std::nullopt) const {} + const std::optional>& dx = std::nullopt) const {} private: std::function(const Dune::FieldVector&, const double&)> - neumannBoundaryLoad; - const BoundaryPatch* neumannBoundary; + neumannBoundaryLoad_; + const BoundaryPatch* neumannBoundary_; /** - * @brief Const accessor to the underlying displacement-based finite element (CRTP). + * \brief Const accessor to the underlying displacement-based finite element (CRTP). * - * @return Const reference to the underlying displacement-based finite element. + * \return Const reference to the underlying displacement-based finite element. */ - constexpr const DisplacementBasedElement& dbElement() const { - return static_cast(*this); - } + constexpr const DFE& dbElement() const { return static_cast(*this); } }; } // namespace Ikarus diff --git a/ikarus/finiteelements/mechanics/loads/volume.hh b/ikarus/finiteelements/mechanics/loads/volume.hh index b8bd49fe0..85b9eb380 100644 --- a/ikarus/finiteelements/mechanics/loads/volume.hh +++ b/ikarus/finiteelements/mechanics/loads/volume.hh @@ -9,14 +9,14 @@ namespace Ikarus { /** - * @brief Volume class represents distributed volume load that can be applied. + * \brief Volume class represents distributed volume load that can be applied. * - * @ingroup mechanics + * \ingroup mechanics * - * @tparam DisplacementBasedElement The type of the displacement-based finite element. - * @tparam Traits Type of traits for handling finite elements. + * \tparam DFE The type of the displacement-based finite element. + * \tparam Traits Type of traits for handling finite elements. */ -template +template class Volume { public: @@ -25,73 +25,72 @@ public: static constexpr int worldDim = Traits::worlddim; /** - * @brief Constructor for the Loads class. + * \brief Constructor for the Loads class. * - * @tparam VolumeLoad The type for the volume load function. + * \tparam VL The type for the volume load function. * - * @param p_volumeLoad Volume load function. + * \param volumeLoad Volume load function. */ - template - explicit Volume(VolumeLoad p_volumeLoad = {}) { - if constexpr (!std::is_same_v) - volumeLoad = p_volumeLoad; + template + explicit Volume(VL volumeLoad = {}) { + if constexpr (!std::is_same_v) + volumeLoad_ = volumeLoad; } /** - * @brief Calculate the scalar value. + * \brief Calculate the scalar value. * * Calculates the scalar value based on the given FERequirements. * - * @param req The FERequirements. - * @return The calculated scalar value. + * \param req The FERequirements. + * \return The calculated scalar value. */ double calculateScalar(const FERequirementType& req) const { return calculateScalarImpl(req); } /** - * @brief Calculate the vector associated with the given FERequirementType. + * \brief Calculate the vector associated with the given FERequirementType. * - * @tparam ScalarType The scalar type for the calculation. - * @param req The FERequirementType object specifying the requirements for the calculation. - * @param force The vector to store the calculated result. + * \tparam ScalarType The scalar type for the calculation. + * \param req The FERequirementType object specifying the requirements for the calculation. + * \param force The vector to store the calculated result. */ void calculateVector(const FERequirementType& req, typename Traits::template VectorType<> force) const { calculateVectorImpl(req, force); } /** - * @brief Calculate the matrix associated with the given FERequirementType. + * \brief Calculate the matrix associated with the given FERequirementType. * - * @tparam ScalarType The scalar type for the calculation. - * @param req The FERequirementType object specifying the requirements for the calculation. - * @param K The matrix to store the calculated result. + * \tparam ScalarType The scalar type for the calculation. + * \param req The FERequirementType object specifying the requirements for the calculation. + * \param K The matrix to store the calculated result. */ void calculateMatrix(const FERequirementType& req, typename Traits::template MatrixType<> K) const { calculateMatrixImpl(req, K); } protected: - template + template auto calculateScalarImpl(const FERequirementType& par, - const std::optional>& dx = std::nullopt) const - -> ScalarType { - if (not volumeLoad) + const std::optional>& dx = std::nullopt) const -> ST { + if (not volumeLoad_) return 0.0; - ScalarType energy = 0.0; + ST energy = 0.0; const auto uFunction = dbElement().displacementFunction(par, dx); const auto geo = dbElement().geometry(); const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { const auto uVal = uFunction.evaluate(gpIndex); - Eigen::Vector fext = volumeLoad(geo.global(gp.position()), lambda); + Eigen::Vector fext = volumeLoad_(geo.global(gp.position()), lambda); energy -= uVal.dot(fext) * geo.integrationElement(gp.position()) * gp.weight(); } return energy; } - template - void calculateVectorImpl(const FERequirementType& par, typename Traits::template VectorType force, - const std::optional>& dx = std::nullopt) const { - if (not volumeLoad) + template + void calculateVectorImpl(const FERequirementType& par, typename Traits::template VectorType force, + const std::optional>& dx = std::nullopt) const { + if (not volumeLoad_) return; using namespace Dune::DerivativeDirections; using namespace Dune; @@ -100,7 +99,7 @@ protected: const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { - const Eigen::Vector fext = volumeLoad(geo.global(gp.position()), lambda); + const Eigen::Vector fext = volumeLoad_(geo.global(gp.position()), lambda); const double intElement = geo.integrationElement(gp.position()) * gp.weight(); for (size_t i = 0; i < dbElement().numberOfNodes(); ++i) { const auto udCi = uFunction.evaluateDerivative(gpIndex, wrt(coeff(i))); @@ -109,20 +108,18 @@ protected: } } - template + template void calculateMatrixImpl(const FERequirementType& par, typename Traits::template MatrixType<> K, - const std::optional>& dx = std::nullopt) const {} + const std::optional>& dx = std::nullopt) const {} private: - std::function(const Dune::FieldVector&, const double&)> volumeLoad; + std::function(const Dune::FieldVector&, const double&)> volumeLoad_; /** - * @brief Const accessor to the underlying displacement-based finite element (CRTP). + * \brief Const accessor to the underlying displacement-based finite element (CRTP). * - * @return Const reference to the underlying displacement-based finite element. + * \return Const reference to the underlying displacement-based finite element. */ - constexpr const DisplacementBasedElement& dbElement() const { - return static_cast(*this); - } + constexpr const DFE& dbElement() const { return static_cast(*this); } }; } // namespace Ikarus diff --git a/ikarus/finiteelements/mechanics/materials.hh b/ikarus/finiteelements/mechanics/materials.hh index 6cd010c4e..cbfd3ac2b 100644 --- a/ikarus/finiteelements/mechanics/materials.hh +++ b/ikarus/finiteelements/mechanics/materials.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file materials.hh - * @brief Header file for material models in Ikarus finite element mechanics. - * @ingroup mechanics + * \file materials.hh + * \brief Header file for material models in Ikarus finite element mechanics. + * \ingroup mechanics */ #pragma once diff --git a/ikarus/finiteelements/mechanics/materials/interface.hh b/ikarus/finiteelements/mechanics/materials/interface.hh index d35fe08fd..5707d8108 100644 --- a/ikarus/finiteelements/mechanics/materials/interface.hh +++ b/ikarus/finiteelements/mechanics/materials/interface.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file interface.hh - * @brief Contains the Material interface class and related template functions for material properties. - * @ingroup materials + * \file interface.hh + * \brief Contains the Material interface class and related template functions for material properties. + * \ingroup materials */ #pragma once @@ -19,43 +19,43 @@ namespace Ikarus { #ifndef DOXYGEN -template +template struct Material; -template +template struct VanishingStress; #endif /** - * @brief Template function for checking if the strain size is correct. + * \brief Template function for checking if the strain size is correct. * * The given strain quantity has to be a Eigen::Vector6 or a Eigen::Matrix3 * - * @tparam Material Type of the material. - * @tparam Strains Type of the strains. + * \tparam MAT Type of the material. + * \tparam S Type of the strains. */ -template +template consteval bool hasCorrectSize() { - if constexpr (Concepts::EigenVector6 or Concepts::EigenMatrix33) + if constexpr (Concepts::EigenVector6 or Concepts::EigenMatrix33) return true; - if constexpr (Material::isReduced and Concepts::EigenVector) { - return Strains::RowsAtCompileTime == Material::freeStrains; + if constexpr (MAT::isReduced and Concepts::EigenVector) { + return S::RowsAtCompileTime == MAT::freeStrains; } else return false; } /** - * @brief Template concept for ensuring correct strain size. + * \brief Template concept for ensuring correct strain size. * - * @tparam Material Type of the material. - * @tparam Strains Type of the strains. + * \tparam MAT Type of the material. + * \tparam S Type of the strains. */ -template -concept CorrectStrainSize = hasCorrectSize(); +template +concept CorrectStrainSize = hasCorrectSize(); /** - * @brief Interface classf or materials. - * @ingroup materials + * \brief Interface classf or materials. + * \ingroup materials * \details Consider a hyper elastic material with the free Helmholtz energy * \f[\require{cases}\psi: \begin{cases}\mathbb{R}^{3\times 3} &\rightarrow \mathbb{R} \\ \BC * &\mapsto \psi(\BC) \end{cases}.\f] @@ -70,48 +70,48 @@ concept CorrectStrainSize = hasCorrectSize(); * tensor](https://en.wikipedia.org/wiki/Finite_strain_theory#Cauchy_strain_tensor_(right_Cauchy%E2%80%93Green_deformation_tensor)), * the [deformation gradient](https://en.wikipedia.org/wiki/Finite_strain_theory#Deformation_gradient_tensor) * \f$\mathbf{F}\f$ or linear strains. The current supported tags are given by Ikarus::StrainTags. - * @tparam MaterialImpl Type of the underlying material implementation. + * \tparam MI Type of the underlying material implementation. */ -template +template struct Material { - using MaterialImplType = MaterialImpl; ///< Type of material implementation + using MaterialImpl = MI; ///< Type of material implementation /** - * @brief Static constant for determining if the material has vanishing stress components (is reduced). + * \brief Static constant for determining if the material has vanishing stress components (is reduced). */ static constexpr bool isReduced = traits::isSpecializationNonTypeAndTypes::value; /** - * @brief Const accessor to the underlying material (CRTP). + * \brief Const accessor to the underlying material (CRTP). * - * @return Const reference to the underlying material. + * \return Const reference to the underlying material. */ constexpr const MaterialImpl& impl() const { return static_cast(*this); } /** - * @brief Accessor to the underlying material (CRTP). + * \brief Accessor to the underlying material (CRTP). * - * @return Reference to the underlying material. + * \return Reference to the underlying material. */ constexpr MaterialImpl& impl() { return static_cast(*this); } /** - * @brief Get the name of the implemented material. + * \brief Get the name of the implemented material. * - * @return Name of the material. + * \return Name of the material. */ [[nodiscard]] constexpr std::string name() const { return impl().nameImpl(); } /** - * @brief Return the stored potential energy of the material. + * \brief Return the stored potential energy of the material. * *\details This function return the free Helmholtz energy of the material * - * @tparam tag Strain tag indicating which strain tensor components are passed. - * @tparam Derived The underlying Eigen type. - * @param Eraw The strain tensor components passed in Voigt notation or matrix notation. - * @return Scalar return of stored energy. + * \tparam tag Strain tag indicating which strain tensor components are passed. + * \tparam Derived The underlying Eigen type. + * \param Eraw The strain tensor components passed in Voigt notation or matrix notation. + * \return Scalar return of stored energy. */ template requires CorrectStrainSize @@ -129,13 +129,13 @@ struct Material } /** - * @brief Get the stresses of the material. + * \brief Get the stresses of the material. * - * @tparam tag Strain tag indicating which strain tensor components are passed. - * @tparam voigt Boolean indicating whether to return Voigt-shaped result. - * @tparam Derived The underlying Eigen type. - * @param Eraw The strain tensor components passed in Voigt notation or matrix notation. - * @return Vectorial or Matrix return of stresses. + * \tparam tag Strain tag indicating which strain tensor components are passed. + * \tparam voigt Boolean indicating whether to return Voigt-shaped result. + * \tparam Derived The underlying Eigen type. + * \param Eraw The strain tensor components passed in Voigt notation or matrix notation. + * \return Vectorial or Matrix return of stresses. */ template requires CorrectStrainSize @@ -150,13 +150,13 @@ struct Material } /** - * @brief Get the tangentModuli of the material. + * \brief Get the tangentModuli of the material. * - * @tparam tag Strain tag indicating which strain tensor components are passed. - * @tparam voigt Boolean indicating whether to return Voigt-shaped result. - * @tparam Derived The underlying Eigen type. - * @param Eraw The strain tensor components passed in Voigt notation or matrix notation. - * @return Tangent moduli in Voigt notation or as fourth-order tensor. + * \tparam tag Strain tag indicating which strain tensor components are passed. + * \tparam voigt Boolean indicating whether to return Voigt-shaped result. + * \tparam Derived The underlying Eigen type. + * \param Eraw The strain tensor components passed in Voigt notation or matrix notation. + * \return Tangent moduli in Voigt notation or as fourth-order tensor. */ template requires CorrectStrainSize @@ -171,16 +171,16 @@ struct Material } /** - * @brief Rebind material to a different scalar type. + * \brief Rebind material to a different scalar type. * * Useful for using automatic differentiation. * - * @tparam ScalarTypeOther The scalar type to rebind to. - * @return Rebound material. + * \tparam STO The scalar type to rebind to. + * \return Rebound material. */ - template + template auto rebind() const { - return impl().template rebind(); + return impl().template rebind(); } private: diff --git a/ikarus/finiteelements/mechanics/materials/linearelasticity.hh b/ikarus/finiteelements/mechanics/materials/linearelasticity.hh index c41a4b6d8..35837c2f5 100644 --- a/ikarus/finiteelements/mechanics/materials/linearelasticity.hh +++ b/ikarus/finiteelements/mechanics/materials/linearelasticity.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file linearelasticity.hh - * @brief Contains the LinearElasticityT class implementation. - * @ingroup materials + * \file linearelasticity.hh + * \brief Contains the LinearElasticityT class implementation. + * \ingroup materials */ #pragma once @@ -12,13 +12,12 @@ #include "svk.hh" #include -#include namespace Ikarus { /** - * @brief Implementation of the Linear Elasticity material model. - * @ingroup materials + * \brief Implementation of the Linear Elasticity material model. + * \ingroup materials * The energy is computed as * \f[ \psi(\Bvep) = \frac{\lambda}{2} (\tr \Bvep)^2 +\mu \tr (\Bvep^2) ,\f] * where \f$ \Bvep \f$ denotes the linear strain tensor. @@ -30,23 +29,23 @@ namespace Ikarus { * \f[ \BBC(\Bvep) =\fracpt{^2\psi(\Bvep)}{\Bvep^2} = \lambda \tr \Bvep \CI +2 \mu \CI^{\mathrm{sym}},\f] * where \f$ \CI_{IJKL} = \de_{IJ}\de_{KL}\f$ and \f$ \CI_{IJKL}^\mathrm{sym} = \frac{1}{2}(\de_{IK}\de_{JL}+ * \de_{IL}\de_{JK})\f$. - * @tparam ScalarType_ The scalar type used in the material. + * \tparam ST The scalar type used in the material. */ -template -struct LinearElasticityT : public Material> +template +struct LinearElasticityT : Material> { [[nodiscard]] constexpr std::string nameImpl() const noexcept { return "LinearElasticity"; } - using ScalarType = ScalarType_; + using ScalarType = ST; using Base = StVenantKirchhoffT; /** - * @brief Constructor for LinearElasticityT. + * \brief Constructor for LinearElasticityT. * - * @param mpt The LamesFirstParameterAndShearModulus object. + * \param mpt The LamesFirstParameterAndShearModulus object. */ explicit LinearElasticityT(const LamesFirstParameterAndShearModulus& mpt) - : svk{mpt} {} + : svk_{mpt} {} using field_type = ScalarType; static constexpr int worldDimension = 3; @@ -69,62 +68,63 @@ struct LinearElasticityT : public Material> static constexpr double derivativeFactor = 1; /** - * @brief Calculates the stored energy in the material. + * \brief Calculates the stored energy in the material. * - * @tparam Derived The underlying Eigen type. - * @param E The strain tensor components. - * @return Scalar return of stored energy. + * \tparam Derived The underlying Eigen type. + * \param E The strain tensor components. + * \return Scalar return of stored energy. */ template ScalarType storedEnergyImpl(const Eigen::MatrixBase& E) const { - return svk.template storedEnergyImpl(E); + return svk_.template storedEnergyImpl(E); } /** - * @brief Calculates the stresses in the material. + * \brief Calculates the stresses in the material. * - * @tparam voigt Boolean indicating whether to return Voigt-shaped result. - * @tparam Derived The underlying Eigen type. - * @param E The strain tensor components. - * @return Matrix valued return of stresses. + * \tparam voigt Boolean indicating whether to return Voigt-shaped result. + * \tparam Derived The underlying Eigen type. + * \param E The strain tensor components. + * \return Matrix valued return of stresses. */ template auto stressesImpl(const Eigen::MatrixBase& E) const { - return svk.template stressesImpl(E); + return svk_.template stressesImpl(E); } /** - * @brief Calculates the tangent moduli in the material. + * \brief Calculates the tangent moduli in the material. * - * @tparam voigt Boolean indicating whether to return Voigt-shaped result. - * @tparam Derived The underlying Eigen type. - * @param E The strain tensor components. - * @return Tangent moduli as fourth-order tensor. + * \tparam voigt Boolean indicating whether to return Voigt-shaped result. + * \tparam Derived The underlying Eigen type. + * \param E The strain tensor components. + * \return Tangent moduli as fourth-order tensor. */ template auto tangentModuliImpl(const Eigen::MatrixBase& E) const { - return svk.template tangentModuliImpl(E); + return svk_.template tangentModuliImpl(E); } /** - * @brief Rebind material to a different scalar type. + * \brief Rebind material to a different scalar type. * - * @tparam ScalarTypeOther The scalar type to rebind to. - * @return Rebound material. + * \tparam STO The scalar type to rebind to. + * \return Rebound material. */ - template + template auto rebind() const { - LinearElasticityT leRebound; - leRebound.svk = svk.template rebind(); + LinearElasticityT leRebound; + leRebound.svk = svk_.template rebind(); return leRebound; } - StVenantKirchhoffT svk; +private: + StVenantKirchhoffT svk_; }; /** - * @brief Convenience typedef for LinearElasticity with double as ScalarType. + * \brief Convenience typedef for LinearElasticity with double as ScalarType. */ -typedef LinearElasticityT LinearElasticity; +using LinearElasticity = LinearElasticityT; } // namespace Ikarus diff --git a/ikarus/finiteelements/mechanics/materials/materials.cpp b/ikarus/finiteelements/mechanics/materials/materials.cpp index 7be454036..6bb3c9001 100644 --- a/ikarus/finiteelements/mechanics/materials/materials.cpp +++ b/ikarus/finiteelements/mechanics/materials/materials.cpp @@ -2,8 +2,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file materials.cpp - * @brief Explicit instantiations for material templates. + * \file materials.cpp + * \brief Explicit instantiations for material templates. */ #include diff --git a/ikarus/finiteelements/mechanics/materials/neohooke.hh b/ikarus/finiteelements/mechanics/materials/neohooke.hh index 239fcfb3b..2e0d80149 100644 --- a/ikarus/finiteelements/mechanics/materials/neohooke.hh +++ b/ikarus/finiteelements/mechanics/materials/neohooke.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file neohooke.hh - * @brief Implementation of the NeoHooke material model. - * @ingroup materials + * \file neohooke.hh + * \brief Implementation of the NeoHooke material model. + * \ingroup materials */ #pragma once @@ -15,8 +15,8 @@ namespace Ikarus { /** - * @brief Implementation of the Neo-Hookean material model. -* @ingroup materials + * \brief Implementation of the Neo-Hookean material model. +* \ingroup materials * The energy is computed as * \f[ \psi(\BC) = \frac{\mu}{2} (\tr \BC-3- 2 \log \sqrt{\det \BC}) + \frac{\lambda}{2} (\log \sqrt{\det \BC})^2 ,\f] * where \f$ \BC \f$ denotes the right Cauchy-Green strain tensor. @@ -30,21 +30,21 @@ namespace Ikarus { * where \f$ \CI_{IJKL} = \frac{1}{2}({(\BC^{-1})}^{IK}{(\BC^{-1})}^{JL}+{(\BC^{-1})}^{IL} {(\BC^{-1})}^{JK}).\f$ * * \remark See \cite bonet2008nonlinear, Section 6.4.3 for a discussion of this material - * @tparam ScalarType_ The scalar type for the strains and stresses,.... + * \tparam ST The scalar type for the strains and stresses,.... */ -template -struct NeoHookeT : public Material> +template +struct NeoHookeT : public Material> { [[nodiscard]] constexpr std::string nameImpl() const noexcept { return "NeoHooke"; } /** - * @brief Constructor for NeoHookeT. - * @param mpt The Lame's parameters (first parameter and shear modulus). + * \brief Constructor for NeoHookeT. + * \param mpt The Lame's parameters (first parameter and shear modulus). */ explicit NeoHookeT(const LamesFirstParameterAndShearModulus& mpt) - : lambdaAndmu{mpt} {} + : lambdaAndmu_{mpt} {} - using ScalarType = ScalarType_; + using ScalarType = ST; static constexpr int worldDimension = 3; using StrainMatrix = Eigen::Matrix; using StressMatrix = StrainMatrix; @@ -60,10 +60,10 @@ struct NeoHookeT : public Material> static constexpr double derivativeFactor = 2; /** - * @brief Computes the stored energy in the Neo-Hookean material model. - * @tparam Derived The derived type of the input matrix. - * @param C The right Cauchy-Green tensor. - * @return ScalarType The stored energy. + * \brief Computes the stored energy in the Neo-Hookean material model. + * \tparam Derived The derived type of the input matrix. + * \param C The right Cauchy-Green tensor. + * \return ScalarType The stored energy. */ template ScalarType storedEnergyImpl(const Eigen::MatrixBase& C) const noexcept { @@ -71,18 +71,18 @@ struct NeoHookeT : public Material> if constexpr (!Concepts::EigenVector) { const auto traceC = C.trace(); const auto logdetF = log(sqrt(C.determinant())); - return lambdaAndmu.mu / 2.0 * (traceC - 3 - 2 * logdetF) + lambdaAndmu.lambda / 2.0 * logdetF * logdetF; + return lambdaAndmu_.mu / 2.0 * (traceC - 3 - 2 * logdetF) + lambdaAndmu_.lambda / 2.0 * logdetF * logdetF; } else static_assert(!Concepts::EigenVector, "NeoHooke energy can only be called with a matrix and not a vector in Voigt notation"); } /** - * @brief Computes the stresses in the Neo-Hookean material model. - * @tparam voigt A boolean indicating whether to return stresses in Voigt notation. - * @tparam Derived The derived type of the input matrix. - * @param C The right Cauchy-Green tensor. - * @return StressMatrix The stresses. + * \brief Computes the stresses in the Neo-Hookean material model. + * \tparam voigt A boolean indicating whether to return stresses in Voigt notation. + * \tparam Derived The derived type of the input matrix. + * \param C The right Cauchy-Green tensor. + * \return StressMatrix The stresses. */ template auto stressesImpl(const Eigen::MatrixBase& C) const { @@ -91,7 +91,7 @@ struct NeoHookeT : public Material> if constexpr (!Concepts::EigenVector) { const auto logdetF = log(sqrt(C.determinant())); const auto invC = C.inverse().eval(); - return (lambdaAndmu.mu * (StrainMatrix::Identity() - invC) + lambdaAndmu.lambda * logdetF * invC).eval(); + return (lambdaAndmu_.mu * (StrainMatrix::Identity() - invC) + lambdaAndmu_.lambda * logdetF * invC).eval(); } else static_assert(!Concepts::EigenVector, "NeoHooke can only be called with a matrix and not a vector in Voigt notation"); @@ -100,11 +100,11 @@ struct NeoHookeT : public Material> } /** - * @brief Computes the tangent moduli in the Neo-Hookean material model. - * @tparam voigt A boolean indicating whether to return tangent moduli in Voigt notation. - * @tparam Derived The derived type of the input matrix. - * @param C The right Cauchy-Green tensor. - * @return Eigen::TensorFixedSize> The tangent moduli. + * \brief Computes the tangent moduli in the Neo-Hookean material model. + * \tparam voigt A boolean indicating whether to return tangent moduli in Voigt notation. + * \tparam Derived The derived type of the input matrix. + * \param C The right Cauchy-Green tensor. + * \return Eigen::TensorFixedSize> The tangent moduli. */ template auto tangentModuliImpl(const Eigen::MatrixBase& C) const { @@ -115,8 +115,8 @@ struct NeoHookeT : public Material> const auto CTinv = tensorView(invC, std::array({3, 3})); static_assert(Eigen::TensorFixedSize>::NumIndices == 2); Eigen::TensorFixedSize> moduli = - (lambdaAndmu.lambda * dyadic(CTinv, CTinv) + - 2 * (lambdaAndmu.mu - lambdaAndmu.lambda * logdetF) * symTwoSlots(fourthOrderIKJL(invC, invC), {2, 3})) + (lambdaAndmu_.lambda * dyadic(CTinv, CTinv) + + 2 * (lambdaAndmu_.mu - lambdaAndmu_.lambda * logdetF) * symTwoSlots(fourthOrderIKJL(invC, invC), {2, 3})) .eval(); return moduli; } else @@ -124,21 +124,22 @@ struct NeoHookeT : public Material> } /** - * @brief Rebinds the material to a different scalar type. - * @tparam ScalarTypeOther The target scalar type. - * @return NeoHookeT The rebound NeoHooke material. + * \brief Rebinds the material to a different scalar type. + * \tparam STO The target scalar type. + * \return NeoHookeT The rebound NeoHooke material. */ - template + template auto rebind() const { - return NeoHookeT(lambdaAndmu); + return NeoHookeT(lambdaAndmu_); } - LamesFirstParameterAndShearModulus lambdaAndmu; +private: + LamesFirstParameterAndShearModulus lambdaAndmu_; }; /** - * @brief Alias for NeoHookeT with double as the default scalar type. + * \brief Alias for NeoHookeT with double as the default scalar type. */ -typedef NeoHookeT NeoHooke; +using NeoHooke = NeoHookeT; } // namespace Ikarus diff --git a/ikarus/finiteelements/mechanics/materials/strainconversions.hh b/ikarus/finiteelements/mechanics/materials/strainconversions.hh index 7ea884d41..491179ea1 100644 --- a/ikarus/finiteelements/mechanics/materials/strainconversions.hh +++ b/ikarus/finiteelements/mechanics/materials/strainconversions.hh @@ -8,7 +8,7 @@ * This file provides implementation for strain-related functions, including the creation and transformation of * different strain tensors. * - * @ingroup materials + * \ingroup materials */ #pragma once @@ -36,7 +36,7 @@ namespace Impl { * * This function creates Green-Lagrangian strains based on the input strain matrix. * What to do is decided by the provided strain tag - * @ingroup materials + * \ingroup materials * \tparam tag Type of the strain tag. * \tparam Derived Type of the Eigen matrix. * \param eMB Eigen matrix representing the input strain. @@ -61,7 +61,7 @@ auto createGreenLagrangianStrains(const Eigen::MatrixBase& eMB) { * * This function creates deformation gradient based on the input strain matrix. * What to do is decided by the provided strain tag - * @ingroup materials + * \ingroup materials * \tparam tag Type of the strain tag. * \tparam Derived Type of the Eigen matrix. * \param eMB Eigen matrix representing the input strain. @@ -92,7 +92,7 @@ decltype(auto) createDeformationGradient(const Eigen::MatrixBase& eMB) * * This function creates Right Cauchy-Green tensor based on the input strain matrix. * What to do is decided by the provided strain tag - * @ingroup materials + * \ingroup materials * \tparam tag Type of the strain tag. * \tparam Derived Type of the Eigen matrix. * \param eMB Eigen matrix representing the input strain. @@ -121,7 +121,7 @@ decltype(auto) createRightCauchyGreen(const Eigen::MatrixBase& eMB) { * \brief Transform strain from one type to another. * * This function transforms one strain component matrix from one type to another, based on the provided strain tags - * @ingroup materials + * \ingroup materials * \tparam from Type of the source strain tag. * \tparam to Type of the target strain tag. * \tparam Derived Type of the Eigen matrix. diff --git a/ikarus/finiteelements/mechanics/materials/svk.hh b/ikarus/finiteelements/mechanics/materials/svk.hh index 9ae956fb8..084e5b6c9 100644 --- a/ikarus/finiteelements/mechanics/materials/svk.hh +++ b/ikarus/finiteelements/mechanics/materials/svk.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file svk.hh - * @brief Implementation of the Saint Venant-Kirchhoff material model. - * @ingroup materials + * \file svk.hh + * \brief Implementation of the Saint Venant-Kirchhoff material model. + * \ingroup materials */ // SPDX-License-Identifier: LGPL-3.0-or-later @@ -17,8 +17,8 @@ namespace Ikarus { /** - * @brief Implementation of the Saint Venant-Kirchhoff material model. - * @ingroup materials + * \brief Implementation of the Saint Venant-Kirchhoff material model. + * \ingroup materials * The energy is computed as * \f[ \psi(\BE) = \frac{\lambda}{2} (\tr \BE)^2 +\mu \tr (\BE^2) ,\f] * where \f$ \BE \f$ denotes the Green-Lagrangian strain. @@ -30,21 +30,21 @@ namespace Ikarus { * \f[ \BBC(\BE) =\fracpt{^2\psi(\BE)}{\BE^2} = \lambda \tr \BE \CI +2 \mu \CI^{\mathrm{sym}},\f] * where \f$ \CI_{IJKL} = \de_{IJ}\de_{KL}\f$ and \f$ \CI_{IJKL}^\mathrm{sym} = \frac{1}{2}(\de_{IK}\de_{JL}+ * \de_{IL}\de_{JK})\f$. - * @tparam ScalarType_ The scalar type used in the material. + * \tparam ST The scalar type used in the material. */ -template -struct StVenantKirchhoffT : public Material> +template +struct StVenantKirchhoffT : public Material> { [[nodiscard]] constexpr std::string nameImpl() const { return "StVenantKirchhoff"; } /** - * @brief Constructor for StVenantKirchhoffT. - * @param mpt The material parameters (Lamé's first parameter and shear modulus). + * \brief Constructor for StVenantKirchhoffT. + * \param mpt The material parameters (Lamé's first parameter and shear modulus). */ explicit StVenantKirchhoffT(const LamesFirstParameterAndShearModulus& mpt) - : materialParameter{mpt} {} + : materialParameter_{mpt} {} - using ScalarType = ScalarType_; + using ScalarType = ST; static constexpr int worldDimension = 3; using StrainMatrix = Eigen::Matrix; using StressMatrix = StrainMatrix; @@ -64,10 +64,10 @@ struct StVenantKirchhoffT : public Material> static constexpr double derivativeFactor = 1; /** - * @brief Computes the stored energy in the Saint Venant-Kirchhoff material model. - * @tparam Derived The derived type of the input matrix. - * @param E The Green-Lagrangian strain. - * @return ScalarType The stored energy. + * \brief Computes the stored energy in the Saint Venant-Kirchhoff material model. + * \tparam Derived The derived type of the input matrix. + * \param E The Green-Lagrangian strain. + * \return ScalarType The stored energy. */ template ScalarType storedEnergyImpl(const Eigen::MatrixBase& E) const { @@ -76,19 +76,19 @@ struct StVenantKirchhoffT : public Material> const ScalarType traceE = E.template segment<3>(0).sum(); const ScalarType squaredNorm = E.template segment<3>(0).squaredNorm() + E.template segment<3>(3).squaredNorm() / ScalarType(2.0); - return materialParameter.lambda / ScalarType(2.0) * traceE * traceE + materialParameter.mu * squaredNorm; + return materialParameter_.lambda / ScalarType(2.0) * traceE * traceE + materialParameter_.mu * squaredNorm; } else { const auto traceE = E.trace(); - return materialParameter.lambda / ScalarType(2.0) * traceE * traceE + materialParameter.mu * E.squaredNorm(); + return materialParameter_.lambda / ScalarType(2.0) * traceE * traceE + materialParameter_.mu * E.squaredNorm(); } } /** - * @brief Computes the stresses in the Saint Venant-Kirchhoff material model. - * @tparam voigt A boolean indicating whether to return stresses in Voigt notation. - * @tparam Derived The derived type of the input matrix. - * @param E The Green-Lagrangian strain. - * @return StressMatrix The stresses. + * \brief Computes the stresses in the Saint Venant-Kirchhoff material model. + * \tparam voigt A boolean indicating whether to return stresses in Voigt notation. + * \tparam Derived The derived type of the input matrix. + * \param E The Green-Lagrangian strain. + * \return StressMatrix The stresses. */ template auto stressesImpl(const Eigen::MatrixBase& E) const { @@ -100,14 +100,14 @@ struct StVenantKirchhoffT : public Material> Eigen::Matrix S; const ScalarType traceE = Ed.template segment<3>(0).sum(); S.diagonal().array() = - materialParameter.lambda * traceE + 2 * materialParameter.mu * Ed.template segment<3>(0).array(); - S(1, 2) = S(2, 1) = materialParameter.mu * Ed(3); - S(0, 2) = S(2, 0) = materialParameter.mu * Ed(4); - S(0, 1) = S(1, 0) = materialParameter.mu * Ed(5); // no two since E voigt has 2* on off-diagonal terms + materialParameter_.lambda * traceE + 2 * materialParameter_.mu * Ed.template segment<3>(0).array(); + S(1, 2) = S(2, 1) = materialParameter_.mu * Ed(3); + S(0, 2) = S(2, 0) = materialParameter_.mu * Ed(4); + S(0, 1) = S(1, 0) = materialParameter_.mu * Ed(5); // no two since E voigt has 2* on off-diagonal terms return S; } else { static_assert(Concepts::EigenMatrix33); - return (materialParameter.lambda * Ed.trace() * StrainMatrix::Identity() + 2 * materialParameter.mu * Ed) + return (materialParameter_.lambda * Ed.trace() * StrainMatrix::Identity() + 2 * materialParameter_.mu * Ed) .eval(); } } else { @@ -115,64 +115,65 @@ struct StVenantKirchhoffT : public Material> static_assert(Concepts::EigenVector6); Eigen::Matrix S; const ScalarType traceE = Ed.template segment<3>(0).sum(); - S.template segment<3>(0).array() = traceE * materialParameter.lambda; - S.template segment<3>(0) += materialParameter.mu * 2 * Ed.template segment<3>(0); + S.template segment<3>(0).array() = traceE * materialParameter_.lambda; + S.template segment<3>(0) += materialParameter_.mu * 2 * Ed.template segment<3>(0); S.template segment<3>(3) = - materialParameter.mu * Ed.template segment<3>(3); // no two since E voigt has 2* on off-diagonal terms + materialParameter_.mu * Ed.template segment<3>(3); // no two since E voigt has 2* on off-diagonal terms return S; } else { Eigen::Matrix S; - S.template segment<3>(0).array() = Ed.trace() * materialParameter.lambda; - S.template segment<3>(0) += 2 * materialParameter.mu * Ed.diagonal(); - S(3) = 2 * materialParameter.mu * Ed(1, 2); - S(4) = 2 * materialParameter.mu * Ed(0, 2); - S(5) = 2 * materialParameter.mu * Ed(0, 1); + S.template segment<3>(0).array() = Ed.trace() * materialParameter_.lambda; + S.template segment<3>(0) += 2 * materialParameter_.mu * Ed.diagonal(); + S(3) = 2 * materialParameter_.mu * Ed(1, 2); + S(4) = 2 * materialParameter_.mu * Ed(0, 2); + S(5) = 2 * materialParameter_.mu * Ed(0, 1); return S; } } } /** - * @brief Computes the tangent moduli in the Saint Venant-Kirchhoff material model. - * @tparam voigt A boolean indicating whether to return tangent moduli in Voigt notation. - * @tparam Derived The derived type of the input matrix. - * @param E The Green-Lagrangian strain. - * @return Eigen::TensorFixedSize> The tangent moduli. + * \brief Computes the tangent moduli in the Saint Venant-Kirchhoff material model. + * \tparam voigt A boolean indicating whether to return tangent moduli in Voigt notation. + * \tparam Derived The derived type of the input matrix. + * \param E The Green-Lagrangian strain. + * \return Eigen::TensorFixedSize> The tangent moduli. */ template auto tangentModuliImpl([[maybe_unused]] const Eigen::MatrixBase& E) const { static_assert(Concepts::EigenMatrixOrVoigtNotation3); if constexpr (!voigt) { Eigen::TensorFixedSize> moduli; - moduli = - materialParameter.lambda * identityFourthOrder() + 2 * materialParameter.mu * symmetricIdentityFourthOrder(); + moduli = materialParameter_.lambda * identityFourthOrder() + + 2 * materialParameter_.mu * symmetricIdentityFourthOrder(); return moduli; } else { Eigen::Matrix moduli; moduli.setZero(); - moduli.template block<3, 3>(0, 0).array() = materialParameter.lambda; - moduli.template block<3, 3>(0, 0).diagonal().array() += 2 * materialParameter.mu; - moduli.template block<3, 3>(3, 3).diagonal().array() = materialParameter.mu; + moduli.template block<3, 3>(0, 0).array() = materialParameter_.lambda; + moduli.template block<3, 3>(0, 0).diagonal().array() += 2 * materialParameter_.mu; + moduli.template block<3, 3>(3, 3).diagonal().array() = materialParameter_.mu; return moduli; } } /** - * @brief Rebinds the material to a different scalar type. - * @tparam ScalarTypeOther The target scalar type. - * @return StVenantKirchhoffT The rebound StVenantKirchhoff material. + * \brief Rebinds the material to a different scalar type. + * \tparam ScalarTypeOther The target scalar type. + * \return StVenantKirchhoffT The rebound StVenantKirchhoff material. */ template auto rebind() const { - return StVenantKirchhoffT(materialParameter); + return StVenantKirchhoffT(materialParameter_); } - LamesFirstParameterAndShearModulus materialParameter; +private: + LamesFirstParameterAndShearModulus materialParameter_; }; /** - * @brief Alias for StVenantKirchhoffT with double as the default scalar type. + * \brief Alias for StVenantKirchhoffT with double as the default scalar type. */ -typedef StVenantKirchhoffT StVenantKirchhoff; +using StVenantKirchhoff = StVenantKirchhoffT; } // namespace Ikarus diff --git a/ikarus/finiteelements/mechanics/materials/tags.hh b/ikarus/finiteelements/mechanics/materials/tags.hh index 5a785064b..d2e1bd31e 100644 --- a/ikarus/finiteelements/mechanics/materials/tags.hh +++ b/ikarus/finiteelements/mechanics/materials/tags.hh @@ -2,8 +2,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file tags.hh - * @brief Definition of several material related enums + * \file tags.hh + * \brief Definition of several material related enums */ #pragma once diff --git a/ikarus/finiteelements/mechanics/materials/vanishingstress.hh b/ikarus/finiteelements/mechanics/materials/vanishingstress.hh index 310fedf4d..18394d6bb 100644 --- a/ikarus/finiteelements/mechanics/materials/vanishingstress.hh +++ b/ikarus/finiteelements/mechanics/materials/vanishingstress.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file vanishingstress.hh - * @brief Defines the VanishingStress material model and related functions. - * @ingroup materials + * \file vanishingstress.hh + * \brief Defines the VanishingStress material model and related functions. + * \ingroup materials */ // SPDX-License-Identifier: LGPL-3.0-or-later @@ -20,7 +20,7 @@ namespace Ikarus { namespace Impl { /** - * @brief Represents a pair of stress matrix indices (row and column). + * \brief Represents a pair of stress matrix indices (row and column). */ struct StressIndexPair { @@ -29,10 +29,10 @@ namespace Impl { }; /** - * @brief Helper function to create an array of free Voigt indices. - * @tparam size The size of the fixed pairs array. - * @param fixed An array of StressIndexPair representing fixed indices. - * @return std::array The array of free Voigt indices. + * \brief Helper function to create an array of free Voigt indices. + * \tparam size The size of the fixed pairs array. + * \param fixed An array of StressIndexPair representing fixed indices. + * \return std::array The array of free Voigt indices. */ template consteval auto createfreeVoigtIndices(const std::array& fixed) { @@ -46,10 +46,10 @@ namespace Impl { } /** - * @brief Helper function to create an array of fixed Voigt indices. - * @tparam size The size of the fixed pairs array. - * @param fixed An array of StressIndexPair representing fixed indices. - * @return std::array The array of fixed Voigt indices. + * \brief Helper function to create an array of fixed Voigt indices. + * \tparam size The size of the fixed pairs array. + * \param fixed An array of StressIndexPair representing fixed indices. + * \return std::array The array of fixed Voigt indices. */ template consteval auto createFixedVoigtIndices(const std::array& fixed) { @@ -60,10 +60,10 @@ namespace Impl { } /** - * @brief Helper function to count the number of diagonal indices in the fixed pairs array. - * @tparam size The size of the fixed pairs array. - * @param fixed An array of StressIndexPair representing fixed indices. - * @return constexpr size_t The number of diagonal indices. + * \brief Helper function to count the number of diagonal indices in the fixed pairs array. + * \tparam size The size of the fixed pairs array. + * \param fixed An array of StressIndexPair representing fixed indices. + * \return constexpr size_t The number of diagonal indices. */ template constexpr size_t countDiagonalIndices(const std::array& fixed) { @@ -78,74 +78,74 @@ namespace Impl { } // namespace Impl /** - * @brief VanishingStress material model that enforces stress components to be zero. - * @ingroup materials - * @tparam stressIndexPair An array of StressIndexPair representing fixed stress components. - * @tparam MaterialImpl The underlying material model. + * \brief VanishingStress material model that enforces stress components to be zero. + * \ingroup materials + * \tparam stressIndexPair An array of StressIndexPair representing fixed stress components. + * \tparam MI The underlying material model. */ -template -struct VanishingStress : public Material> +template +struct VanishingStress : public Material> { /** - * @brief Constructor for VanishingStress. - * @param mat The underlying material model. - * @param p_tol Tolerance for stress reduction. + * \brief Constructor for VanishingStress. + * \param mat The underlying material model. + * \param tol Tolerance for stress reduction. */ - explicit VanishingStress(MaterialImpl mat, typename MaterialImpl::ScalarType p_tol = 1e-12) - : matImpl{mat}, - tol{p_tol} {} + explicit VanishingStress(MI mat, typename MI::ScalarType tol = 1e-12) + : matImpl_{mat}, + tol_{tol} {} - using Underlying = MaterialImpl; ///< The underlying material type. + using Underlying = MI; ///< The underlying material type. static constexpr auto fixedPairs = stressIndexPair; ///< Array of fixed stress components. static constexpr auto freeVoigtIndices = createfreeVoigtIndices(fixedPairs); ///< Free Voigt indices. static constexpr auto fixedVoigtIndices = createFixedVoigtIndices(fixedPairs); ///< Fixed Voigt indices. static constexpr auto fixedDiagonalVoigtIndicesSize = - countDiagonalIndices(fixedPairs); ///< Number of fixed diagonal indices. - static constexpr auto freeStrains = freeVoigtIndices.size(); ///< Number of free strains. - using ScalarType = typename MaterialImpl::ScalarType; ///< Scalar type. + countDiagonalIndices(fixedPairs); ///< Number of fixed diagonal indices. + static constexpr auto freeStrains = freeVoigtIndices.size(); ///< Number of free strains. + using ScalarType = typename Underlying::ScalarType; ///< Scalar type. [[nodiscard]] constexpr std::string nameImpl() const noexcept { - auto matName = matImpl.name() + "_Vanishing("; + auto matName = matImpl_.name() + "_Vanishing("; for (auto p : fixedPairs) matName += "(" + std::to_string(p.row) + std::to_string(p.col) + ")"; matName += ")"; return matName; } - static constexpr auto strainTag = MaterialImpl::strainTag; ///< Strain tag. - static constexpr auto stressTag = MaterialImpl::stressTag; ///< Stress tag. - static constexpr auto tangentModuliTag = MaterialImpl::tangentModuliTag; ///< Tangent moduli tag. - static constexpr bool energyAcceptsVoigt = MaterialImpl::energyAcceptsVoigt; ///< Energy accepts Voigt notation. - static constexpr bool stressToVoigt = true; ///< Stress to Voigt notation. - static constexpr bool stressAcceptsVoigt = true; ///< Stress accepts Voigt notation. - static constexpr bool moduliToVoigt = true; ///< Moduli to Voigt notation. - static constexpr bool moduliAcceptsVoigt = true; ///< Moduli accepts Voigt notation. - static constexpr double derivativeFactor = 1; ///< Derivative factor. + static constexpr auto strainTag = Underlying::strainTag; ///< Strain tag. + static constexpr auto stressTag = Underlying::stressTag; ///< Stress tag. + static constexpr auto tangentModuliTag = Underlying::tangentModuliTag; ///< Tangent moduli tag. + static constexpr bool energyAcceptsVoigt = Underlying::energyAcceptsVoigt; ///< Energy accepts Voigt notation. + static constexpr bool stressToVoigt = true; ///< Stress to Voigt notation. + static constexpr bool stressAcceptsVoigt = true; ///< Stress accepts Voigt notation. + static constexpr bool moduliToVoigt = true; ///< Moduli to Voigt notation. + static constexpr bool moduliAcceptsVoigt = true; ///< Moduli accepts Voigt notation. + static constexpr double derivativeFactor = 1; ///< Derivative factor. /** - * @brief Computes the stored energy for the VanishingStress material. - * @tparam Derived The derived type of the input matrix. - * @param E The Green-Lagrangian strain. - * @return ScalarType The stored energy. + * \brief Computes the stored energy for the VanishingStress material. + * \tparam Derived The derived type of the input matrix. + * \param E The Green-Lagrangian strain. + * \return ScalarType The stored energy. */ template ScalarType storedEnergyImpl(const Eigen::MatrixBase& E) const { const auto [nonOp, Esol] = reduceStress(E); - return matImpl.storedEnergyImpl(Esol); + return matImpl_.storedEnergyImpl(Esol); } /** - * @brief Computes the stresses for the VanishingStress material. - * @tparam voigt A boolean indicating whether to return stresses in Voigt notation. - * @tparam Derived The derived type of the input matrix. - * @param E The Green-Lagrangian strain. - * @return StressMatrix The stresses. + * \brief Computes the stresses for the VanishingStress material. + * \tparam voigt A boolean indicating whether to return stresses in Voigt notation. + * \tparam Derived The derived type of the input matrix. + * \param E The Green-Lagrangian strain. + * \return StressMatrix The stresses. */ template auto stressesImpl(const Eigen::MatrixBase& E) const { const auto [nonOp, Esol] = reduceStress(E); - auto stressesRed = matImpl.template stresses(Esol); + auto stressesRed = matImpl_.template stresses(Esol); if constexpr (voigt) { return removeCol(stressesRed, fixedVoigtIndices); @@ -154,16 +154,16 @@ struct VanishingStress : public Material auto tangentModuliImpl(const Eigen::MatrixBase& E) const { const auto [nonOp, Esol] = reduceStress(E); - auto C = matImpl.template tangentModuli(Esol); + auto C = matImpl_.template tangentModuli(Esol); if constexpr (voigt) return staticCondensation(C, fixedVoigtIndices); else @@ -171,22 +171,22 @@ struct VanishingStress : public Material auto rebind() const { - auto reboundMatImpl = matImpl.template rebind(); - return VanishingStress(reboundMatImpl, tol); + auto reboundMatImpl = matImpl_.template rebind(); + return VanishingStress(reboundMatImpl, tol_); } private: /** - * @brief Converts the input strain matrix to the appropriate form for stress reduction. - * @tparam Derived The derived type of the input matrix. - * @param E The input strain matrix. - * @return decltype(auto) The converted strain matrix. + * \brief Converts the input strain matrix to the appropriate form for stress reduction. + * \tparam Derived The derived type of the input matrix. + * \param E The input strain matrix. + * \return decltype(auto) The converted strain matrix. */ template decltype(auto) maybeFromVoigt(const Eigen::MatrixBase& E) const { @@ -197,9 +197,9 @@ private: } /** - * @brief Initializes unknown strains based on fixed indices. - * @tparam Derived The derived type of the input matrix. - * @param E The input strain matrix. + * \brief Initializes unknown strains based on fixed indices. + * \tparam Derived The derived type of the input matrix. + * \param E The input strain matrix. */ template void initUnknownStrains(Eigen::MatrixBase& E) const { @@ -216,14 +216,14 @@ private: } /** - * @brief Reduces stress components to satisfy the vanishing stress condition. - * @tparam Derived The derived type of the input matrix. - * @param p_Eraw The input strain matrix. - * @return std::pair The stress reduction result. + * \brief Reduces stress components to satisfy the vanishing stress condition. + * \tparam Derived The derived type of the input matrix. + * \param Eraw The input strain matrix. + * \return std::pair The stress reduction result. */ template - auto reduceStress(const Eigen::MatrixBase& p_Eraw) const { - auto E = maybeFromVoigt(p_Eraw); + auto reduceStress(const Eigen::MatrixBase& Eraw) const { + auto E = maybeFromVoigt(Eraw); initUnknownStrains(E); std::array fixedDiagonalVoigtIndices; @@ -234,25 +234,25 @@ private: } auto f = [&](auto&) { - auto S = matImpl.template stresses(E); + auto S = matImpl_.template stresses(E); return S(fixedDiagonalVoigtIndices).eval(); }; auto df = [&](auto&) { - auto moduli = (matImpl.template tangentModuli(E)).eval(); - return (moduli(fixedDiagonalVoigtIndices, fixedDiagonalVoigtIndices) / MaterialImpl::derivativeFactor).eval(); + auto moduli = (matImpl_.template tangentModuli(E)).eval(); + return (moduli(fixedDiagonalVoigtIndices, fixedDiagonalVoigtIndices) / Underlying::derivativeFactor).eval(); }; auto Er = E(fixedDiagonalVoigtIndices, fixedDiagonalVoigtIndices).eval().template cast(); auto nonOp = Ikarus::NonLinearOperator(functions(f, df), parameter(Er)); auto nr = Ikarus::makeNewtonRaphson( nonOp, [&](auto& r, auto& A) { return (A.inverse() * r).eval(); }, - [&](auto& /* Ex33 */, auto& Ecomps) { + [&](auto& /* Ex33 */, auto& ecomps) { for (int ri = 0; auto i : fixedDiagonalVoigtIndices) { auto indexPair = fromVoigt(i); - E(indexPair[0], indexPair[1]) += Ecomps(ri++); + E(indexPair[0], indexPair[1]) += ecomps(ri++); } }); - nr->setup({.tol = tol, .maxIter = 100}); + nr->setup({.tol = tol_, .maxIter = 100}); if (!static_cast(nr->solve())) DUNE_THROW(Dune::MathError, "The stress reduction of material " << nameImpl() << " was unsuccessful\n" << "The strains are\n" @@ -261,17 +261,17 @@ private: return std::make_pair(nonOp, E); } - MaterialImpl matImpl; ///< The underlying material model. - double tol{}; ///< Tolerance for stress reduction. + Underlying matImpl_; ///< The underlying material model. + double tol_{}; ///< Tolerance for stress reduction. }; /** - * @brief Factory function to create a VanishingStress material with specified stress indices. - * @tparam stressIndexPair The array of StressIndexPair representing fixed stress components. - * @tparam MaterialImpl The underlying material model. - * @param mat The underlying material model. - * @param p_tol Tolerance for stress reduction. - * @return VanishingStress The created VanishingStress material. + * \brief Factory function to create a VanishingStress material with specified stress indices. + * \tparam stressIndexPair The array of StressIndexPair representing fixed stress components. + * \tparam MaterialImpl The underlying material model. + * \param mat The underlying material model. + * \param p_tol Tolerance for stress reduction. + * \return VanishingStress The created VanishingStress material. */ template auto makeVanishingStress(MaterialImpl mat, typename MaterialImpl::ScalarType p_tol = 1e-12) { @@ -279,41 +279,41 @@ auto makeVanishingStress(MaterialImpl mat, typename MaterialImpl::ScalarType p_t } /** - * @brief Factory function to create a VanishingStress material for plane stress conditions. - * @tparam MaterialImpl The underlying material model. - * @param mat The underlying material model. - * @param p_tol Tolerance for stress reduction. - * @return VanishingStress The created VanishingStress material for plane stress. + * \brief Factory function to create a VanishingStress material for plane stress conditions. + * \tparam MaterialImpl The underlying material model. + * \param mat The underlying material model. + * \param tol Tolerance for stress reduction. + * \return VanishingStress The created VanishingStress material for plane stress. */ template -auto planeStress(const MaterialImpl& mat, typename MaterialImpl::ScalarType p_tol = 1e-8) { +auto planeStress(const MaterialImpl& mat, typename MaterialImpl::ScalarType tol = 1e-8) { return makeVanishingStress( - mat, p_tol); + mat, tol); } /** - * @brief Factory function to create a VanishingStress material for a shell material with zero normal stress + * \brief Factory function to create a VanishingStress material for a shell material with zero normal stress * condition. - * @tparam MaterialImpl The underlying material model. - * @param mat The underlying material model. - * @param p_tol Tolerance for stress reduction. - * @return VanishingStress The created VanishingStress material for plane stress. + * \tparam MaterialImpl The underlying material model. + * \param mat The underlying material model. + * \param tol Tolerance for stress reduction. + * \return VanishingStress The created VanishingStress material for plane stress. */ template -auto shellMaterial(const MaterialImpl& mat, typename MaterialImpl::ScalarType p_tol = 1e-8) { - return makeVanishingStress(mat, p_tol); +auto shellMaterial(const MaterialImpl& mat, typename MaterialImpl::ScalarType tol = 1e-8) { + return makeVanishingStress(mat, tol); } /** - * @brief Factory function to create a VanishingStress material for a beam material with two zero normal stress + * \brief Factory function to create a VanishingStress material for a beam material with two zero normal stress * condition. - * @tparam MaterialImpl The underlying material model. - * @param mat The underlying material model. - * @param p_tol Tolerance for stress reduction. - * @return VanishingStress The created VanishingStress material for plane stress. + * \tparam MaterialImpl The underlying material model. + * \param mat The underlying material model. + * \param p_tol Tolerance for stress reduction. + * \return VanishingStress The created VanishingStress material for plane stress. */ template -auto beamMaterial(const MaterialImpl& mat, typename MaterialImpl::ScalarType p_tol = 1e-8) { - return makeVanishingStress(mat, p_tol); +auto beamMaterial(const MaterialImpl& mat, typename MaterialImpl::ScalarType tol = 1e-8) { + return makeVanishingStress(mat, tol); } } // namespace Ikarus diff --git a/ikarus/finiteelements/mechanics/membranestrains.hh b/ikarus/finiteelements/mechanics/membranestrains.hh index 21101c43a..5d873533a 100644 --- a/ikarus/finiteelements/mechanics/membranestrains.hh +++ b/ikarus/finiteelements/mechanics/membranestrains.hh @@ -23,12 +23,12 @@ struct DefaultMembraneStrain * \param geo The geometry object providing the transposed Jacobian. * \param uFunction The function representing the displacement field. * - * \tparam Geometry The type of the geometry object. + * \tparam GEO The type of the geometry object. * * \return The strain vector at the given integration point. */ - template - auto value(const Dune::FieldVector& gpPos, const Geometry& geo, const auto& uFunction) const + template + auto value(const Dune::FieldVector& gpPos, const GEO& geo, const auto& uFunction) const -> Eigen::Vector3::ctype> { using ScalarType = typename std::remove_cvref_t::ctype; Eigen::Vector3 epsV; @@ -56,16 +56,15 @@ struct DefaultMembraneStrain * \param localBasis The local basis object. * \param node The FE node index. * - * \tparam Geometry The type of the geometry object. - * \tparam ScalarType The scalar type for the matrix elements. + * \tparam GEO The type of the geometry object. + * \tparam ST The scalar type for the matrix elements. * * \return The strain-displacement matrix for the given node and integration point. */ - template - auto derivative(const Dune::FieldVector& gpPos, const Eigen::Matrix& jcur, - const auto& dNAtGp, const Geometry& geo, const auto& uFunction, const auto& localBasis, - const int node) const { - Eigen::Matrix bop; + template + auto derivative(const Dune::FieldVector& gpPos, const Eigen::Matrix& jcur, const auto& dNAtGp, + const GEO& geo, const auto& uFunction, const auto& localBasis, const int node) const { + Eigen::Matrix bop; bop.row(0) = jcur.row(0) * dNAtGp(node, 0); bop.row(1) = jcur.row(1) * dNAtGp(node, 1); bop.row(2) = jcur.row(0) * dNAtGp(node, 1) + jcur.row(1) * dNAtGp(node, 0); @@ -87,21 +86,21 @@ struct DefaultMembraneStrain * \param I The index of the first node. * \param J The index of the second node. * - * \tparam Geometry The type of the geometry object. - * \tparam ScalarType The scalar type for the matrix elements. + * \tparam GEO The type of the geometry object. + * \tparam ST The scalar type for the matrix elements. * * \return The second derivative of the membrane strain. */ - template - auto secondDerivative(const Dune::FieldVector& gpPos, const auto& dNAtGp, const Geometry& geo, - const auto& uFunction, const auto& localBasis, const Eigen::Vector3& S, int I, + template + auto secondDerivative(const Dune::FieldVector& gpPos, const auto& dNAtGp, const GEO& geo, + const auto& uFunction, const auto& localBasis, const Eigen::Vector3& S, int I, int J) const { - const auto& dN1i = dNAtGp(I, 0); - const auto& dN1j = dNAtGp(J, 0); - const auto& dN2i = dNAtGp(I, 1); - const auto& dN2j = dNAtGp(J, 1); - const ScalarType NS = dN1i * dN1j * S[0] + dN2i * dN2j * S[1] + (dN1i * dN2j + dN2i * dN1j) * S[2]; - Eigen::Matrix kg = Eigen::Matrix::Identity() * NS; + const auto& dN1i = dNAtGp(I, 0); + const auto& dN1j = dNAtGp(J, 0); + const auto& dN2i = dNAtGp(I, 1); + const auto& dN2j = dNAtGp(J, 1); + const ST NS = dN1i * dN1j * S[0] + dN2i * dN2j * S[1] + (dN1i * dN2j + dN2i * dN1j) * S[2]; + Eigen::Matrix kg = Eigen::Matrix::Identity() * NS; return kg; } }; diff --git a/ikarus/finiteelements/mechanics/nonlinearelastic.hh b/ikarus/finiteelements/mechanics/nonlinearelastic.hh index c976ebfc2..3e4cf9b10 100644 --- a/ikarus/finiteelements/mechanics/nonlinearelastic.hh +++ b/ikarus/finiteelements/mechanics/nonlinearelastic.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file nonlinearelastic.hh - * @brief Definition of the NonLinearElastic class for finite element mechanics computations. - * @ingroup mechanics + * \file nonlinearelastic.hh + * \brief Definition of the NonLinearElastic class for finite element mechanics computations. + * \ingroup mechanics */ #pragma once @@ -32,24 +32,22 @@ namespace Ikarus { /** - * @brief NonLinearElastic class represents a non-linear elastic finite element. + * \brief NonLinearElastic class represents a non-linear elastic finite element. * - * @ingroup mechanics + * \ingroup mechanics * - * @tparam Basis_ The basis type for the finite element. - * @tparam Material_ The material type for the finite element. - * @tparam FERequirements_ The requirements for the finite element. - * @tparam useEigenRef A boolean flag indicating whether to use Eigen references. + * \tparam B The basis type for the finite element. + * \tparam MAT The material type for the finite element. + * \tparam FER The requirements for the finite element. + * \tparam useEigenRef A boolean flag indicating whether to use Eigen references. */ -template , bool useEigenRef = false> -class NonLinearElastic : public PowerBasisFE, - public Volume, - FETraits>, - public Traction, - FETraits> +template , bool useEigenRef = false> +class NonLinearElastic : public PowerBasisFE, + public Volume, FETraits>, + public Traction, FETraits> { public: - using Traits = FETraits; + using Traits = FETraits; using Basis = typename Traits::Basis; using FlatBasis = typename Traits::FlatBasis; using FERequirementType = typename Traits::FERequirementType; @@ -58,76 +56,76 @@ public: using GridView = typename Traits::GridView; using Element = typename Traits::Element; using BasePowerFE = PowerBasisFE; // Handles globalIndices function - using Material = Material_; - using VolumeType = Volume, Traits>; - using TractionType = Traction, Traits>; + using Material = MAT; + using VolumeType = Volume; + using TractionType = Traction; using LocalBasisType = decltype(std::declval().tree().child(0).finiteElement().localBasis()); static constexpr int myDim = Traits::mydim; static constexpr auto strainType = StrainTags::greenLagrangian; /** - * @brief Constructor for the NonLinearElastic class. + * \brief Constructor for the NonLinearElastic class. * - * @tparam VolumeLoad The type for the volume load function. - * @tparam NeumannBoundaryLoad The type for the Neumann boundary load function. - * @param globalBasis The global basis for the finite element. - * @param element The element for which the finite element is constructed. - * @param p_mat The material for the non-linear elastic element. - * @param p_volumeLoad Volume load function (default is LoadDefault). - * @param p_neumannBoundary Neumann boundary patch (default is nullptr). - * @param p_neumannBoundaryLoad Neumann boundary load function (default is LoadDefault). + * \tparam VolumeLoad The type for the volume load function. + * \tparam NeumannBoundaryLoad The type for the Neumann boundary load function. + * \param globalBasis The global basis for the finite element. + * \param element The element for which the finite element is constructed. + * \param mat The material for the non-linear elastic element. + * \param volumeLoad Volume load function (default is LoadDefault). + * \param neumannBoundary Neumann boundary patch (default is nullptr). + * \param neumannBoundaryLoad Neumann boundary load function (default is LoadDefault). */ template - NonLinearElastic(const Basis& globalBasis, const typename LocalView::Element& element, const Material& p_mat, - VolumeLoad p_volumeLoad = {}, const BoundaryPatch* p_neumannBoundary = nullptr, - NeumannBoundaryLoad p_neumannBoundaryLoad = {}) + NonLinearElastic(const Basis& globalBasis, const typename LocalView::Element& element, const Material& mat, + VolumeLoad volumeLoad = {}, const BoundaryPatch* neumannBoundary = nullptr, + NeumannBoundaryLoad neumannBoundaryLoad = {}) : BasePowerFE(globalBasis, element), - VolumeType(p_volumeLoad), - TractionType(p_neumannBoundary, p_neumannBoundaryLoad), - mat{p_mat} { + VolumeType(volumeLoad), + TractionType(neumannBoundary, neumannBoundaryLoad), + mat_{mat} { this->localView().bind(element); - auto& first_child = this->localView().tree().child(0); - const auto& fe = first_child.finiteElement(); - geo_ = std::make_shared(this->localView().element().geometry()); - numberOfNodes_ = fe.size(); - order_ = 2 * (fe.localBasis().order()); - localBasis = Dune::CachedLocalBasis(fe.localBasis()); + auto& firstChild = this->localView().tree().child(0); + const auto& fe = firstChild.finiteElement(); + geo_ = std::make_shared(this->localView().element().geometry()); + numberOfNodes_ = fe.size(); + order_ = 2 * (fe.localBasis().order()); + localBasis_ = Dune::CachedLocalBasis(fe.localBasis()); if constexpr (requires { this->localView().element().impl().getQuadratureRule(order_); }) if (this->localView().element().impl().isTrimmed()) - localBasis.bind(this->localView().element().impl().getQuadratureRule(order_), Dune::bindDerivatives(0, 1)); + localBasis_.bind(this->localView().element().impl().getQuadratureRule(order_), Dune::bindDerivatives(0, 1)); else - localBasis.bind(Dune::QuadratureRules::rule(this->localView().element().type(), order_), - Dune::bindDerivatives(0, 1)); + localBasis_.bind(Dune::QuadratureRules::rule(this->localView().element().type(), order_), + Dune::bindDerivatives(0, 1)); else - localBasis.bind(Dune::QuadratureRules::rule(this->localView().element().type(), order_), - Dune::bindDerivatives(0, 1)); + localBasis_.bind(Dune::QuadratureRules::rule(this->localView().element().type(), order_), + Dune::bindDerivatives(0, 1)); } /** - * @brief Get the displacement function for the given FERequirementType. + * \brief Get the displacement function for the given FERequirementType. * - * @tparam ScalarType The scalar type for the displacement function. - * @param par The FERequirementType object. - * @param dx Optional displacement vector. - * @return A StandardLocalFunction representing the displacement function. + * \tparam ScalarType The scalar type for the displacement function. + * \param par The FERequirementType object. + * \param dx Optional displacement vector. + * \return A StandardLocalFunction representing the displacement function. */ template auto displacementFunction(const FERequirementType& par, const std::optional>& dx = std::nullopt) const { const auto& d = par.getGlobalSolution(Ikarus::FESolutions::displacement); auto disp = Ikarus::FEHelper::localSolutionBlockVector(d, this->localView(), dx); - Dune::StandardLocalFunction uFunction(localBasis, disp, geo_); + Dune::StandardLocalFunction uFunction(localBasis_, disp, geo_); return uFunction; } /** - * @brief The strain function for the given FERequirementType. + * \brief The strain function for the given FERequirementType. * - * @tparam ScalarType The scalar type for the strain function. - * @param par The FERequirementType object. - * @param dx Optional displacement vector. - * @return The strain function calculated using greenLagrangeStrains. + * \tparam ScalarType The scalar type for the strain function. + * \param par The FERequirementType object. + * \param dx Optional displacement vector. + * \return The strain function calculated using greenLagrangeStrains. */ template inline auto strainFunction(const FERequirementType& par, @@ -136,57 +134,57 @@ public: } /** - * @brief Get the material tangent for the given strain. + * \brief Get the material tangent for the given strain. * - * @tparam ScalarType The scalar type for the material and strain. - * @tparam strainDim The dimension of the strain vector. - * @tparam voigt Flag indicating whether to use Voigt notation. - * @param strain The strain vector. - * @return The material tangent calculated using the material's tangentModuli function. + * \tparam ScalarType The scalar type for the material and strain. + * \tparam strainDim The dimension of the strain vector. + * \tparam voigt Flag indicating whether to use Voigt notation. + * \param strain The strain vector. + * \return The material tangent calculated using the material's tangentModuli function. */ template auto materialTangent(const Eigen::Vector& strain) const { if constexpr (std::is_same_v) - return mat.template tangentModuli(strain); + return mat_.template tangentModuli(strain); else { - decltype(auto) matAD = mat.template rebind(); + decltype(auto) matAD = mat_.template rebind(); return matAD.template tangentModuli(strain); } } /** - * @brief Get the internal energy for the given strain. + * \brief Get the internal energy for the given strain. * - * @tparam ScalarType The scalar type for the material and strain. - * @tparam strainDim The dimension of the strain vector. - * @param strain The strain vector. - * @return The internal energy calculated using the material's storedEnergy function. + * \tparam ScalarType The scalar type for the material and strain. + * \tparam strainDim The dimension of the strain vector. + * \param strain The strain vector. + * \return The internal energy calculated using the material's storedEnergy function. */ template auto getInternalEnergy(const Eigen::Vector& strain) const { if constexpr (std::is_same_v) - return mat.template storedEnergy(strain); + return mat_.template storedEnergy(strain); else { - decltype(auto) matAD = mat.template rebind(); + decltype(auto) matAD = mat_.template rebind(); return matAD.template storedEnergy(strain); } } /** - * @brief Get the stress for the given strain. + * \brief Get the stress for the given strain. * - * @tparam ScalarType The scalar type for the material and strain. - * @tparam strainDim The dimension of the strain vector. - * @tparam voigt A boolean indicating whether to use the Voigt notation for stress. - * @param strain The strain vector. - * @return The stress vector calculated using the material's stresses function. + * \tparam ScalarType The scalar type for the material and strain. + * \tparam strainDim The dimension of the strain vector. + * \tparam voigt A boolean indicating whether to use the Voigt notation for stress. + * \param strain The strain vector. + * \return The stress vector calculated using the material's stresses function. */ template auto getStress(const Eigen::Vector& strain) const { if constexpr (std::is_same_v) - return mat.template stresses(strain); + return mat_.template stresses(strain); else { - decltype(auto) matAD = mat.template rebind(); + decltype(auto) matAD = mat_.template rebind(); return matAD.template stresses(strain); } } @@ -196,31 +194,31 @@ public: [[nodiscard]] int order() const { return order_; } /** - * @brief Calculate the scalar value associated with the given FERequirementType. + * \brief Calculate the scalar value associated with the given FERequirementType. * - * @tparam ScalarType The scalar type for the calculation. - * @param par The FERequirementType object specifying the requirements for the calculation. - * @return The calculated scalar value. + * \tparam ScalarType The scalar type for the calculation. + * \param par The FERequirementType object specifying the requirements for the calculation. + * \return The calculated scalar value. */ double calculateScalar(const FERequirementType& par) const { return calculateScalarImpl(par); } /** - * @brief Calculate the vector associated with the given FERequirementType. + * \brief Calculate the vector associated with the given FERequirementType. * - * @tparam ScalarType The scalar type for the calculation. - * @param par The FERequirementType object specifying the requirements for the calculation. - * @param force The vector to store the calculated result. + * \tparam ScalarType The scalar type for the calculation. + * \param par The FERequirementType object specifying the requirements for the calculation. + * \param force The vector to store the calculated result. */ void calculateVector(const FERequirementType& par, typename Traits::template VectorType<> force) const { calculateVectorImpl(par, force); } /** - * @brief Calculate the matrix associated with the given FERequirementType. + * \brief Calculate the matrix associated with the given FERequirementType. * - * @tparam ScalarType The scalar type for the calculation. - * @param par The FERequirementType object specifying the requirements for the calculation. - * @param K The matrix to store the calculated result. + * \tparam ScalarType The scalar type for the calculation. + * \param par The FERequirementType object specifying the requirements for the calculation. + * \param K The matrix to store the calculated result. */ void calculateMatrix(const FERequirementType& par, typename Traits::template MatrixType<> K) const { using namespace Dune::DerivativeDirections; @@ -247,13 +245,13 @@ public: } /** - * @brief Calculates a requested result at a specific local position. + * \brief Calculates a requested result at a specific local position. * - * @param req The FERequirementType object holding the global solution. - * @param local Local position vector. - * @return calculated result + * \param req The FERequirementType object holding the global solution. + * \param local Local position vector. + * \return calculated result * - * @tparam resType The type representing the requested result. + * \tparam resType The type representing the requested result. */ template auto calculateAt(const FERequirementType& req, const Dune::FieldVector& local) const { @@ -266,14 +264,14 @@ public: const auto E = (0.5 * (H.transpose() + H + H.transpose() * H)).eval(); const auto EVoigt = toVoigt(E); - return mat.template stresses(EVoigt); + return mat_.template stresses(EVoigt); } } private: std::shared_ptr geo_; - Dune::CachedLocalBasis> localBasis; - Material mat; + Dune::CachedLocalBasis> localBasis_; + Material mat_; size_t numberOfNodes_{0}; int order_{}; diff --git a/ikarus/finiteelements/physicshelper.hh b/ikarus/finiteelements/physicshelper.hh index ecf3b3b38..a28380f56 100644 --- a/ikarus/finiteelements/physicshelper.hh +++ b/ikarus/finiteelements/physicshelper.hh @@ -19,12 +19,8 @@ namespace Ikarus { * \param E Young's modulus. * \param nu Poisson's ratio. * \return 3x3 material tangent matrix. - * \deprecated Use the material library Ikarus::LinearElasticity. */ -[[deprecated( - "These are hard-coded function you should use the material library Ikarus::LinearElasticity")]] inline Eigen:: - Matrix3d - planeStressLinearElasticMaterialTangent(double E, double nu) { +inline Eigen::Matrix3d planeStressLinearElasticMaterialTangent(double E, double nu) { Eigen::Matrix3d C; C.setZero(); C(0, 0) = C(1, 1) = 1; @@ -40,12 +36,8 @@ namespace Ikarus { * \param E Young's modulus. * \param nu Poisson's ratio. * \return 6x6 material tangent matrix. - * \deprecated Use the material library Ikarus::LinearElasticity. */ -[[deprecated( - "These are hard-coded function you should use the material library: Ikarus::LinearElasticity")]] inline Eigen:: - Matrix - linearElasticMaterialTangent3D(double E, double nu) { +inline Eigen::Matrix linearElasticMaterialTangent3D(double E, double nu) { Eigen::Matrix C; C.setZero(); C(0, 0) = C(1, 1) = C(2, 2) = 1 - nu; @@ -101,15 +93,13 @@ struct LamesFirstParameterAndShearModulus /** * \brief Concept for checking if a type is a valid material parameter tuple. * - * \tparam MaterialParameter Type to check. + * \tparam MP Type to check. */ -template -concept MaterialParameterTuple = std::is_same_v or - std::is_same_v or - std::is_same_v or - std::is_same_v or - std::is_same_v or - std::is_same_v; +template +concept MPTuple = + std::is_same_v or std::is_same_v or + std::is_same_v or std::is_same_v or + std::is_same_v or std::is_same_v; /** * \brief Conversion utility for Lame's constants. @@ -125,17 +115,17 @@ struct ConvertLameConstants !std::is_same_v) { if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& nu = vp.nu; + const auto& E = vp_.emodul; + const auto& nu = vp_.nu; return Dune::FloatCmp::eq(nu, 0.5) ? std::numeric_limits::infinity() : E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)); } else if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& mu = vp.mu; + const auto& E = vp_.emodul; + const auto& mu = vp_.mu; return mu * (E - 2.0 * mu) / (3.0 * mu - E); } else if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& K = vp.K; + const auto& E = vp_.emodul; + const auto& K = vp_.K; return 3.0 * K * (3.0 * K - E) / (9.0 * K - E); } else assert(false && "Your LameParameter request is not implemented"); @@ -146,20 +136,20 @@ struct ConvertLameConstants !std::is_same_v) { if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& nu = vp.nu; + const auto& E = vp_.emodul; + const auto& nu = vp_.nu; return E / (3.0 * (1.0 - 2.0 * nu)); } else if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& mu = vp.mu; + const auto& E = vp_.emodul; + const auto& mu = vp_.mu; return E * mu / (3.0 * (3.0 * mu - E)); } else if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& lambda = vp.lambda; - return (E + 3.0 * lambda + calcR(vp)) / 6.0; + const auto& E = vp_.emodul; + const auto& lambda = vp_.lambda; + return (E + 3.0 * lambda + calcR(vp_)) / 6.0; } else if constexpr (std::is_same_v) { - const auto& lambda = vp.lambda; - const auto& mu = vp.mu; + const auto& lambda = vp_.lambda; + const auto& mu = vp_.mu; return lambda + 2.0 * mu / 3.0; } else assert(false && "Your LameParameter request is not implemented"); @@ -170,20 +160,20 @@ struct ConvertLameConstants !std::is_same_v) { if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& nu = vp.nu; + const auto& E = vp_.emodul; + const auto& nu = vp_.nu; return E / (2.0 * (1.0 + nu)); } else if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& K = vp.K; + const auto& E = vp_.emodul; + const auto& K = vp_.K; return 3.0 * K * E / (9.0 * K - E); } else if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& lambda = vp.lambda; - return (E - 3.0 * lambda + calcR(vp)) / 4.0; + const auto& E = vp_.emodul; + const auto& lambda = vp_.lambda; + return (E - 3.0 * lambda + calcR(vp_)) / 4.0; } else if constexpr (std::is_same_v) { - const auto& K = vp.K; - const auto& lambda = vp.lambda; + const auto& K = vp_.K; + const auto& lambda = vp_.lambda; return 3.0 * (K - lambda) / 2.0; } else assert(false && "Your LameParameter request is not implemented"); @@ -191,28 +181,28 @@ struct ConvertLameConstants constexpr double toPWaveModulus() { if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& nu = vp.nu; + const auto& E = vp_.emodul; + const auto& nu = vp_.nu; return E * (1.0 - nu) / ((1.0 + nu) * (1.0 - 2.0 * nu)); } else if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& mu = vp.mu; + const auto& E = vp_.emodul; + const auto& mu = vp_.mu; return mu * (4.0 * mu - E) / (3.0 * mu - E); } else if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& K = vp.K; + const auto& E = vp_.emodul; + const auto& K = vp_.K; return 3.0 * K * (3.0 * K + E) / (9.0 * K - E); } else if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& lambda = vp.lambda; - return (E - lambda + calcR(vp)) / 2.0; + const auto& E = vp_.emodul; + const auto& lambda = vp_.lambda; + return (E - lambda + calcR(vp_)) / 2.0; } else if constexpr (std::is_same_v) { - const auto& K = vp.K; - const auto& lambda = vp.lambda; + const auto& K = vp_.K; + const auto& lambda = vp_.lambda; return 3.0 * K - 2.0 * lambda; } else if constexpr (std::is_same_v) { - const auto& lambda = vp.lambda; - const auto& mu = vp.mu; + const auto& lambda = vp_.lambda; + const auto& mu = vp_.mu; return lambda + 2.0 * mu; } else assert(false && "Your LameParameter request is not implemented"); @@ -222,24 +212,24 @@ struct ConvertLameConstants requires(!std::is_same_v) { if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& mu = vp.mu; + const auto& E = vp_.emodul; + const auto& mu = vp_.mu; return E / (2.0 * mu) - 1.0; } else if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& K = vp.K; + const auto& E = vp_.emodul; + const auto& K = vp_.K; return (3.0 * K - E) / (6.0 * K); } else if constexpr (std::is_same_v) { - const auto& E = vp.emodul; - const auto& lambda = vp.lambda; - return 2.0 * lambda / (E + lambda + calcR(vp)); + const auto& E = vp_.emodul; + const auto& lambda = vp_.lambda; + return 2.0 * lambda / (E + lambda + calcR(vp_)); } else if constexpr (std::is_same_v) { - const auto& K = vp.K; - const auto& lambda = vp.lambda; + const auto& K = vp_.K; + const auto& lambda = vp_.lambda; return lambda / (3 * K - lambda); } else if constexpr (std::is_same_v) { - const auto& lambda = vp.lambda; - const auto& mu = vp.mu; + const auto& lambda = vp_.lambda; + const auto& mu = vp_.mu; return lambda / (2.0 * (lambda + mu)); } else assert(false && "Your LameParameter request is not implemented"); @@ -252,10 +242,10 @@ struct ConvertLameConstants !std::is_same_v) { if constexpr (std::is_same_v) { - return 9.0 * vp.K * (vp.K - vp.lambda) / (3.0 * vp.K - vp.lambda); + return 9.0 * vp_.K * (vp_.K - vp_.lambda) / (3.0 * vp_.K - vp_.lambda); } else if constexpr (std::is_same_v) { - const auto& lambda = vp.lambda; - const auto& mu = vp.mu; + const auto& lambda = vp_.lambda; + const auto& mu = vp_.mu; return mu * (3.0 * lambda + 2.0 * mu) / (lambda + mu); } else assert(false && "Your LameParameter request is not implemented"); @@ -263,48 +253,49 @@ struct ConvertLameConstants private: friend ConvertLameConstants convertLameConstants( - const YoungsModulusAndPoissonsRatio& p_vp); + const YoungsModulusAndPoissonsRatio& valuePair); friend ConvertLameConstants convertLameConstants( - const YoungsModulusAndShearModulus& p_vp); + const YoungsModulusAndShearModulus& valuePair); friend ConvertLameConstants convertLameConstants( - const YoungsModulusAndBulkModulus& p_vp); + const YoungsModulusAndBulkModulus& valuePair); friend ConvertLameConstants convertLameConstants( - const LamesFirstParameterAndShearModulus& p_vp); + const LamesFirstParameterAndShearModulus& valuePair); friend ConvertLameConstants convertLameConstants( - const BulkModulusAndLamesFirstParameter& p_vp); - ConvertLameConstants(ValuePair&& p_vp) - : vp(p_vp) {} - ConvertLameConstants(const ValuePair& p_vp) - : vp(p_vp) {} + const BulkModulusAndLamesFirstParameter& valuePair); + ConvertLameConstants(ValuePair&& valuePair) + : vp_(valuePair) {} + ConvertLameConstants(const ValuePair& valuePair) + : vp_(valuePair) {} - double calcR(const YoungsModulusAndLamesFirstParameter& vp_) { - const auto& E = vp_.emodul; - const auto& lambda = vp_.lambda; + double calcR(const YoungsModulusAndLamesFirstParameter& valuePair) { + const auto& E = valuePair.emodul; + const auto& lambda = valuePair.lambda; return std::sqrt(E * E + 9 * lambda * lambda + 2 * E * lambda); } - ValuePair vp; + ValuePair vp_; }; inline ConvertLameConstants convertLameConstants( - const YoungsModulusAndPoissonsRatio& p_vp) { - return {p_vp}; + const YoungsModulusAndPoissonsRatio& valuePair) { + return {valuePair}; } inline ConvertLameConstants convertLameConstants( - const YoungsModulusAndShearModulus& p_vp) { - return {p_vp}; + const YoungsModulusAndShearModulus& valuePair) { + return {valuePair}; } -inline ConvertLameConstants convertLameConstants(const YoungsModulusAndBulkModulus& p_vp) { - return {p_vp}; +inline ConvertLameConstants convertLameConstants( + const YoungsModulusAndBulkModulus& valuePair) { + return {valuePair}; } inline ConvertLameConstants convertLameConstants( - const LamesFirstParameterAndShearModulus& p_vp) { - return {p_vp}; + const LamesFirstParameterAndShearModulus& valuePair) { + return {valuePair}; } inline ConvertLameConstants convertLameConstants( - const BulkModulusAndLamesFirstParameter& p_vp) { - return {p_vp}; + const BulkModulusAndLamesFirstParameter& valuePair) { + return {valuePair}; } /** diff --git a/ikarus/io/resultevaluators.hh b/ikarus/io/resultevaluators.hh index 59135068a..169a1e833 100644 --- a/ikarus/io/resultevaluators.hh +++ b/ikarus/io/resultevaluators.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file resultevaluators.hh - * @brief Ikarus Result Evaluators for Stress Analysis - * @ingroup resultevaluators + * \file resultevaluators.hh + * \brief Ikarus Result Evaluators for Stress Analysis + * \ingroup resultevaluators * */ @@ -23,20 +23,20 @@ class FieldVector; namespace Ikarus::ResultEvaluators { /** - * @brief Struct for calculating von Mises stress - * @ingroup resultevaluators - * @details The VonMises struct provides a function call operator to calculate von Mises stress. - * @tparam dim dimension of stress state + * \brief Struct for calculating von Mises stress + * \ingroup resultevaluators + * \details The VonMises struct provides a function call operator to calculate von Mises stress. + * \tparam dim dimension of stress state */ template requires(dim == 2 or dim == 3) struct VonMises { /** - * @brief Calculate the result quantity (von Mises stress) - * @param resultArray EigenMatrix containing the stress state - * @param comp component of result (not used here) - * @return von Mises stress + * \brief Calculate the result quantity (von Mises stress) + * \param resultArray EigenMatrix containing the stress state + * \param comp component of result (not used here) + * \return von Mises stress */ double operator()(const auto& resultArray, [[maybe_unused]] const int comp) const { if constexpr (dim == 2) { @@ -59,31 +59,31 @@ struct VonMises } /** - * @brief Get the name of the result type (VonMises) - * @return String representing the name + * \brief Get the name of the result type (VonMises) + * \return String representing the name */ static std::string name() { return "VonMises"; } /** - * @brief Get the number of components in the result (always 1 for VonMises) - * @return Number of components + * \brief Get the number of components in the result (always 1 for VonMises) + * \return Number of components */ static int ncomps() { return 1; } }; /** - * @brief Struct for calculating principal stresses - * @ingroup resultevaluators - * @details The PrincipalStress struct provides a function call operator to calculate principal stresses. - * @remark Only 2D stresses are supported + * \brief Struct for calculating principal stresses + * \ingroup resultevaluators + * \details The PrincipalStress struct provides a function call operator to calculate principal stresses. + * \remark Only 2D stresses are supported */ struct PrincipalStress { /** - * @brief Calculate the result quantity (principal stress) - * @param resultArray EigenMatrix containing the stress state - * @param comp component of result - * @return principal stress + * \brief Calculate the result quantity (principal stress) + * \param resultArray EigenMatrix containing the stress state + * \param comp component of result + * \return principal stress */ double operator()(const auto& resultArray, const int comp) const { const auto s_x = resultArray(0, 0); @@ -97,14 +97,14 @@ struct PrincipalStress } /** - * @brief Get the name of the result type (PrincipalStress) - * @return String representing the name + * \brief Get the name of the result type (PrincipalStress) + * \return String representing the name */ static std::string name() { return "PrincipalStress"; } /** - * @brief Get the number of components in the result (always 2 for PrincipalStress) - * @return Number of components + * \brief Get the number of components in the result (always 2 for PrincipalStress) + * \return Number of components */ static int ncomps() { return 2; } }; diff --git a/ikarus/io/resultfunction.hh b/ikarus/io/resultfunction.hh index f13445575..7804b2482 100644 --- a/ikarus/io/resultfunction.hh +++ b/ikarus/io/resultfunction.hh @@ -2,9 +2,9 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file resultfunction.hh - * @brief Ikarus Result Evaluators for Stress Analysis - * @ingroup io + * \file resultfunction.hh + * \brief Ikarus Result Evaluators for Stress Analysis + * \ingroup io * */ @@ -27,54 +27,54 @@ namespace Impl { } // namespace Impl /** - * @brief Wrapper to evaluate results for a vtkwriter. - * @details + * \brief Wrapper to evaluate results for a vtkwriter. + * \details * Usage: - * @code - * auto resultFunction = std::make_shared>(&fes, feReq); + * \code + * auto resultFunction = std::make_shared>(&fes, feReq); * * vtkWriter.addPointData(Dune::Vtk::Function( resultFunction)); * // or with Dunes native Vtk * vtkWriter.addVertexData(resultFunction); - * @endcode - * @ingroup io - * @tparam ElementType_ Type of the finite element - * @tparam resType requested result type - * @tparam UserFunction Type of the user-defined function for custom result evaluation (default is + * \endcode + * \ingroup io + * \tparam FE Type of the finite element + * \tparam resType requested result type + * \tparam UserFunction Type of the user-defined function for custom result evaluation (default is DefaultUserFunction) */ -template -class ResultFunction : public Dune::VTKFunction +template +class ResultFunction : public Dune::VTKFunction { public: - using ElementType = ElementType_; - using FERequirementType = typename ElementType::FERequirementType; - using GridView = typename ElementType::GridView; + using FiniteElement = FE; + using FERequirementType = typename FiniteElement::FERequirementType; + using GridView = typename FiniteElement::GridView; using ctype = typename GridView::ctype; constexpr static int griddim = GridView::dimension; using Entity = typename GridView::template Codim<0>::Entity; /** - * @brief Evaluate the component at a given entity and local coordinates. + * \brief Evaluate the component at a given entity and local coordinates. * * This function is required by the Dune::VTKFunction interface. * - * @param comp Stress component index - * @param e Entity on which to evaluate the stress - * @param local Local coordinates within the entity - * @return Stress component value + * \param comp Stress component index + * \param e Entity on which to evaluate the stress + * \param local Local coordinates within the entity + * \return Stress component value */ double evaluate(int comp, const Entity& e, const Dune::FieldVector& local) const override { - auto index = gridView.indexSet().index(e); + auto index = gridView_.indexSet().index(e); return evaluateComponent(index, local, comp); } /** - * @brief Get the number of components. + * \brief Get the number of components. * * This function is required by the Dune::VTKFunction interface. * - * @return Number of stress components + * \return Number of stress components */ [[nodiscard]] int ncomps() const override { if constexpr (std::is_same_v) { @@ -88,11 +88,11 @@ public: } /** - * @brief Get the name of the result type. + * \brief Get the name of the result type. * * This function is required by the Dune::VTKFunction interface. * - * @return String representing the name of the result type + * \return String representing the name of the result type */ [[nodiscard]] constexpr std::string name() const override { if constexpr (std::is_same_v) @@ -102,15 +102,15 @@ public: } /** - * @brief Constructor for ResultFunction. + * \brief Constructor for ResultFunction. * * Constructs a ResultFunction object with given finite elements, ferequirements * - * @param fes Pointer to a vector of finite elements - * @param req FERequirements for evaluation + * \param fes Pointer to a vector of finite elements + * \param req FERequirements for evaluation */ - ResultFunction(std::vector* fes, const FERequirementType& req) - : gridView{fes->at(0).localView().globalBasis().gridView()}, + ResultFunction(std::vector* fes, const FERequirementType& req) + : gridView_{fes->at(0).localView().globalBasis().gridView()}, feRequirements_{req}, fes_{fes}, userFunction_{UserFunction{}} {} @@ -125,9 +125,9 @@ private: return result(comp); } - GridView gridView; + GridView gridView_; FERequirementType feRequirements_; - std::vector* fes_; + std::vector* fes_; [[no_unique_address]] std::string name_{}; UserFunction userFunction_; }; diff --git a/ikarus/linearalgebra/truncatedconjugategradient.hh b/ikarus/linearalgebra/truncatedconjugategradient.hh index 8d1e5cbeb..28b59560e 100644 --- a/ikarus/linearalgebra/truncatedconjugategradient.hh +++ b/ikarus/linearalgebra/truncatedconjugategradient.hh @@ -6,8 +6,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later /** - * @file truncatedconjugategradient.hh - * @brief Definition of TruncatedConjugateGradient class for solving linear systems using truncated conjugate gradient + * \file truncatedconjugategradient.hh + * \brief Definition of TruncatedConjugateGradient class for solving linear systems using truncated conjugate gradient * method. * * This file defines the TruncatedConjugateGradient class, which is an iterative solver for solving linear systems @@ -52,18 +52,18 @@ namespace internal { /** * @internal - * @brief Low-level truncated conjugate gradient algorithm. - * @tparam MatrixType Type of the matrix A. - * @tparam Rhs Type of the right-hand side vector b. - * @tparam Dest Type of the solution vector x. - * @tparam Preconditioner Type of the preconditioner. - * @param mat The matrix A. - * @param rhs The right-hand side vector b. - * @param x On input and initial solution, on output the computed solution. - * @param precond A preconditioner being able to efficiently solve for an approximation of Ax=b (regardless of b). - * @param iters On input the max number of iterations, on output the number of performed iterations. - * @param tol_error On input the tolerance error, on output an estimation of the relative error. - * @param _info Information about the truncated conjugate gradient algorithm. + * \brief Low-level truncated conjugate gradient algorithm. + * \tparam MatrixType Type of the matrix A. + * \tparam Rhs Type of the right-hand side vector b. + * \tparam Dest Type of the solution vector x. + * \tparam Preconditioner Type of the preconditioner. + * \param mat The matrix A. + * \param rhs The right-hand side vector b. + * \param x On input and initial solution, on output the computed solution. + * \param precond A preconditioner being able to efficiently solve for an approximation of Ax=b (regardless of b). + * \param iters On input the max number of iterations, on output the number of performed iterations. + * \param tol_error On input the tolerance error, on output an estimation of the relative error. + * \param _info Information about the truncated conjugate gradient algorithm. */ template void truncated_conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, const Preconditioner& precond, @@ -186,20 +186,19 @@ namespace internal { /** * @class TruncatedConjugateGradient - * @brief Iterative solver for solving linear systems using the truncated conjugate gradient method. - * @tparam MatrixType_ Type of the matrix A. - * @tparam UpLo_ Type of the triangular part of the matrix (Lower or Upper or both). - * @tparam Preconditioner_ Type of the preconditioner. + * \brief Iterative solver for solving linear systems using the truncated conjugate gradient method. + * \tparam M Type of the matrix A. + * \tparam upLo Type of the triangular part of the matrix (Lower or Upper or both). + * \tparam PC Type of the preconditioner. */ -template -class TruncatedConjugateGradient - : public IterativeSolverBase > +template +class TruncatedConjugateGradient : public IterativeSolverBase > { public: typedef IterativeSolverBase Base; TruncatedConjugateGradient(TruncatedConjugateGradient&& other) noexcept : Base(std::move(other)), - algInfo{other.algInfo} {} + algInfo_{other.algInfo_} {} private: using Base::m_error; @@ -207,31 +206,31 @@ private: using Base::m_isInitialized; using Base::m_iterations; using Base::matrix; - mutable TCGInfo algInfo; + mutable TCGInfo algInfo_; public: - typedef MatrixType_ MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef Preconditioner_ Preconditioner; + using MatrixType = M; + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename MatrixType::RealScalar; + using Preconditioner = PC; enum { - UpLo = UpLo_ + UpLo = upLo }; public: /** - * @brief Get information about the truncated conjugate gradient algorithm. - * @return Information about the algorithm. + * \brief Get information about the truncated conjugate gradient algorithm. + * \return Information about the algorithm. */ - TCGInfo getInfo() { return algInfo; } + TCGInfo getInfo() { return algInfo_; } /** - * @brief Set information about the truncated conjugate gradient algorithm. - * @param _alginfo Information about the algorithm. + * \brief Set information about the truncated conjugate gradient algorithm. + * \param _alginfo Information about the algorithm. */ - void setInfo(TCGInfo _alginfo) { this->algInfo = _alginfo; } + void setInfo(TCGInfo alginfo) { this->algInfo_ = alginfo; } /** Default constructor. */ TruncatedConjugateGradient() : Base() {} @@ -275,11 +274,9 @@ public: RowMajorWrapper row_mat(matrix()); internal::truncated_conjugate_gradient(SelfAdjointWrapper(row_mat), b, x, Base::m_preconditioner, m_iterations, - m_error, algInfo); + m_error, algInfo_); m_info = m_error <= Base::m_tolerance ? Success : NoConvergence; } - -protected: }; } // end namespace Eigen diff --git a/ikarus/solver/linearsolver/linearsolver.cpp b/ikarus/solver/linearsolver/linearsolver.cpp index b515cc450..fe0e5ba49 100644 --- a/ikarus/solver/linearsolver/linearsolver.cpp +++ b/ikarus/solver/linearsolver/linearsolver.cpp @@ -15,65 +15,65 @@ namespace Ikarus { template -LinearSolverTemplate::LinearSolverTemplate(const SolverTypeTag& p_solverTypeTag) - : solverTypeTag{p_solverTypeTag} { +LinearSolverTemplate::LinearSolverTemplate(const SolverTypeTag& solverTypeTag) + : solverTypeTag_{solverTypeTag} { using namespace Eigen; switch (solverTypeTag) { case SolverTypeTag::si_ConjugateGradient: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::si_LeastSquaresConjugateGradient: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::si_BiCGSTAB: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::sd_SimplicialLLT: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::sd_SimplicialLDLT: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::sd_SparseLU: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::sd_SparseQR: - solverimpl = std::make_unique>>>(); + solverimpl_ = std::make_unique>>>(); break; case SolverTypeTag::sd_CholmodSupernodalLLT: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::sd_UmfPackLU: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::sd_SuperLU: DUNE_THROW(Dune::NotImplemented, "Not implemented yet."); break; // Dense Solver case SolverTypeTag::d_PartialPivLU: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::d_FullPivLU: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::d_HouseholderQR: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::d_ColPivHouseholderQR: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::d_FullPivHouseholderQR: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::d_CompleteOrthogonalDecomposition: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::d_LLT: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::d_LDLT: - solverimpl = std::make_unique>>(); + solverimpl_ = std::make_unique>>(); break; case SolverTypeTag::none: break; diff --git a/ikarus/solver/linearsolver/linearsolver.hh b/ikarus/solver/linearsolver/linearsolver.hh index 7b40bfd7b..5ab882b3a 100644 --- a/ikarus/solver/linearsolver/linearsolver.hh +++ b/ikarus/solver/linearsolver/linearsolver.hh @@ -60,21 +60,22 @@ enum class MatrixTypeTag /** * \class LinearSolverTemplate * \brief A type-erased class which wraps most of the linear solvers available in Eigen. - * \tparam ScalarType The scalar type of the linear system (default: double). + * \tparam ST The scalar type of the linear system (default: double). * \ingroup solvers */ -template +template class LinearSolverTemplate { public: + using ScalarType = ST; using SparseMatrixType = Eigen::SparseMatrix; using DenseMatrixType = Eigen::MatrixX; /** * \brief Constructor for LinearSolverTemplate. - * \param p_solverTypeTag The solver type tag representing the type of the linear solver. + * \param solverTypeTag The solver type tag representing the type of the linear solver. */ - explicit LinearSolverTemplate(const SolverTypeTag& p_solverTypeTag); + explicit LinearSolverTemplate(const SolverTypeTag& solverTypeTag); /** * \brief Destructor for LinearSolverTemplate. @@ -95,7 +96,7 @@ public: * \brief Copy constructor. * \param other The LinearSolverTemplate to copy. */ - LinearSolverTemplate(const LinearSolverTemplate& other) { *this = LinearSolverTemplate(other.solverTypeTag); } + LinearSolverTemplate(const LinearSolverTemplate& other) { *this = LinearSolverTemplate(other.solverTypeTag_); } /** * \brief Move constructor. * \param other The LinearSolverTemplate to move. @@ -174,8 +175,8 @@ private: Solver solver; }; - std::unique_ptr solverimpl; - SolverTypeTag solverTypeTag{SolverTypeTag::none}; + std::unique_ptr solverimpl_; + SolverTypeTag solverTypeTag_{SolverTypeTag::none}; public: /** @@ -187,7 +188,7 @@ public: template requires std::is_same_v || std::is_same_v inline LinearSolverTemplate& compute(const MatrixType& A) { - solverimpl->compute(A); + solverimpl_->compute(A); return *this; } @@ -199,7 +200,7 @@ public: template requires std::is_same_v || std::is_same_v inline void analyzePattern(const MatrixType& A) { - solverimpl->analyzePattern(A); + solverimpl_->analyzePattern(A); } /** @@ -210,7 +211,7 @@ public: template requires std::is_same_v || std::is_same_v inline void factorize(const MatrixType& A) { - solverimpl->factorize(A); + solverimpl_->factorize(A); } /** @@ -218,29 +219,29 @@ public: * \param x The solution vector. * \param b The right-hand side vector. */ - void solve(Eigen::VectorX& x, const Eigen::VectorX& b) { solverimpl->solve(x, b); } + void solve(Eigen::VectorX& x, const Eigen::VectorX& b) { solverimpl_->solve(x, b); } /** * \brief Solve the linear system for a `n` times `3` matrix. * \param x The solution matrix. * \param b The right-hand side matrix. */ - void solve(Eigen::MatrixX3& x, const Eigen::MatrixX3& b) { solverimpl->solve(x, b); } + void solve(Eigen::MatrixX3& x, const Eigen::MatrixX3& b) { solverimpl_->solve(x, b); } /** * \brief Solve the linear system for a `n` times `2` matrix. * \param x The solution matrix. * \param b The right-hand side matrix. */ - void solve(Eigen::MatrixX2& x, const Eigen::MatrixX2& b) { solverimpl->solve(x, b); } + void solve(Eigen::MatrixX2& x, const Eigen::MatrixX2& b) { solverimpl_->solve(x, b); } /** * \brief Solve the linear system for a `n` times `n` matrix. * \param x The solution matrix. * \param b The right-hand side matrix. */ - void solve(Eigen::MatrixX& x, const Eigen::MatrixX& b) { solverimpl->solve(x, b); } + void solve(Eigen::MatrixX& x, const Eigen::MatrixX& b) { solverimpl_->solve(x, b); } }; -typedef LinearSolverTemplate LinearSolver; +using LinearSolver = LinearSolverTemplate; extern template class LinearSolverTemplate; } // namespace Ikarus diff --git a/ikarus/solver/nonlinearsolver/newtonraphson.hh b/ikarus/solver/nonlinearsolver/newtonraphson.hh index f8342a545..b05f46ffc 100644 --- a/ikarus/solver/nonlinearsolver/newtonraphson.hh +++ b/ikarus/solver/nonlinearsolver/newtonraphson.hh @@ -31,48 +31,46 @@ struct NewtonRaphsonSettings /** * \class NewtonRaphson * \brief Implementation of the Newton-Raphson method for solving nonlinear equations. - * \tparam NonLinearOperatorImpl Type of the nonlinear operator to solve. - * \tparam LinearSolver Type of the linear solver used internally (default is SolverDefault). - * \tparam UpdateFunctionTypeImpl Type of the update function (default is UpdateDefault). + * \tparam NLO Type of the nonlinear operator to solve. + * \tparam LS Type of the linear solver used internally (default is SolverDefault). + * \tparam UF Type of the update function (default is UpdateDefault). * \relates makeNewtonRaphson * \ingroup solvers */ -template +template class NewtonRaphson : public IObservable { public: ///< Compile-time boolean indicating if the linear solver satisfies the non-linear solver concept static constexpr bool isLinearSolver = - Ikarus::Concepts::LinearSolverCheck; + Ikarus::Concepts::LinearSolverCheck; ///< Type representing the parameter vector of the nonlinear operator. - using ValueType = typename NonLinearOperatorImpl::template ParameterValue<0>; + using ValueType = typename NLO::template ParameterValue<0>; - using UpdateFunctionType = UpdateFunctionTypeImpl; ///< Type representing the update function. - using NonLinearOperator = NonLinearOperatorImpl; ///< Type of the non-linear operator + using UpdateFunction = UF; ///< Type representing the update function. + using NonLinearOperator = NLO; ///< Type of the non-linear operator /** * \brief Constructor for NewtonRaphson. - * \param p_nonLinearOperator Nonlinear operator to solve. - * \param p_linearSolver Linear solver used internally (default is SolverDefault). - * \param p_updateFunction Update function (default is UpdateDefault). + * \param nonLinearOperator Nonlinear operator to solve. + * \param linearSolver Linear solver used internally (default is SolverDefault). + * \param updateFunction Update function (default is UpdateDefault). */ - explicit NewtonRaphson(const NonLinearOperatorImpl& p_nonLinearOperator, LinearSolver&& p_linearSolver = {}, - UpdateFunctionTypeImpl p_updateFunction = {}) - : nonLinearOperator_{p_nonLinearOperator}, - linearSolver{std::move(p_linearSolver)}, - updateFunction{p_updateFunction} { - if constexpr (std::is_same_v) - correction.setZero(nonLinearOperator().value().size()); + explicit NewtonRaphson(const NonLinearOperator& nonLinearOperator, LS&& linearSolver = {}, + UpdateFunction updateFunction = {}) + : nonLinearOperator_{nonLinearOperator}, + linearSolver_{std::move(linearSolver)}, + updateFunction_{updateFunction} { + if constexpr (std::is_same_v) + correction_.setZero(this->nonLinearOperator().value().size()); } /** * \brief Set up the solver with the given settings. - * \param p_settings Newton-Raphson settings. + * \param settings Newton-Raphson settings. */ - void setup(const NewtonRaphsonSettings& p_settings) { settings = p_settings; } + void setup(const NewtonRaphsonSettings& settings) { settings_ = settings; } #ifndef DOXYGEN struct NoPredictor @@ -81,22 +79,22 @@ public: #endif /** * \brief Solve the nonlinear system. - * \param dx_predictor Predictor for the solution increment (default is NoPredictor). + * \param dxPredictor Predictor for the solution increment (default is NoPredictor). * \return Information about the solution process. */ template requires std::is_same_v || - std::is_convertible_v> + std::is_convertible_v> [[nodiscard( "The solve method returns information of the solution process. You should store this information and check if " "it was successful")]] Ikarus::NonLinearSolverInformation - solve(const SolutionType& dx_predictor = NoPredictor{}) { + solve(const SolutionType& dxPredictor = NoPredictor{}) { this->notify(NonLinearSolverMessages::INIT); Ikarus::NonLinearSolverInformation solverInformation; solverInformation.success = true; auto& x = nonLinearOperator().firstParameter(); if constexpr (not std::is_same_v) - updateFunction(x, dx_predictor); + updateFunction_(x, dxPredictor); nonLinearOperator().updateAll(); const auto& rx = nonLinearOperator().value(); const auto& Ax = nonLinearOperator().derivative(); @@ -104,18 +102,18 @@ public: decltype(rNorm) dNorm; int iter{0}; if constexpr (isLinearSolver) - linearSolver.analyzePattern(Ax); - while (rNorm > settings.tol && iter < settings.maxIter) { + linearSolver_.analyzePattern(Ax); + while (rNorm > settings_.tol && iter < settings_.maxIter) { this->notify(NonLinearSolverMessages::ITERATION_STARTED); if constexpr (isLinearSolver) { - linearSolver.factorize(Ax); - linearSolver.solve(correction, -rx); - dNorm = correction.norm(); - updateFunction(x, correction); + linearSolver_.factorize(Ax); + linearSolver_.solve(correction_, -rx); + dNorm = correction_.norm(); + updateFunction_(x, correction_); } else { - correction = -linearSolver(rx, Ax); - dNorm = norm(correction); - updateFunction(x, correction); + correction_ = -linearSolver_(rx, Ax); + dNorm = norm(correction_); + updateFunction_(x, correction_); } this->notify(NonLinearSolverMessages::CORRECTIONNORM_UPDATED, static_cast(dNorm)); this->notify(NonLinearSolverMessages::SOLUTION_CHANGED); @@ -125,7 +123,7 @@ public: this->notify(NonLinearSolverMessages::ITERATION_ENDED); ++iter; } - if (iter == settings.maxIter) + if (iter == settings_.maxIter) solverInformation.success = false; solverInformation.iterations = iter; solverInformation.residualNorm = static_cast(rNorm); @@ -142,29 +140,27 @@ public: auto& nonLinearOperator() { return nonLinearOperator_; } private: - NonLinearOperatorImpl nonLinearOperator_; - typename NonLinearOperatorImpl::ValueType correction; - LinearSolver linearSolver; - UpdateFunctionType updateFunction; - NewtonRaphsonSettings settings; + NonLinearOperator nonLinearOperator_; + typename NonLinearOperator::ValueType correction_; + LS linearSolver_; + UpdateFunction updateFunction_; + NewtonRaphsonSettings settings_; }; /** * \brief Function to create a NewtonRaphson solver instance. - * \tparam NonLinearOperatorImpl Type of the nonlinear operator to solve. - * \tparam LinearSolver Type of the linear solver used internally (default is SolverDefault). - * \tparam UpdateFunctionType Type of the update function (default is UpdateDefault). - * \param p_nonLinearOperator Nonlinear operator to solve. - * \param p_linearSolver Linear solver used internally (default is SolverDefault). - * \param p_updateFunction Update function (default is UpdateDefault). + * \tparam NLO Type of the nonlinear operator to solve. + * \tparam LS Type of the linear solver used internally (default is SolverDefault). + * \tparam UF Type of the update function (default is UpdateDefault). + * \param nonLinearOperator Nonlinear operator to solve. + * \param linearSolver Linear solver used internally (default is SolverDefault). + * \param updateFunction Update function (default is UpdateDefault). * \return Shared pointer to the NewtonRaphson solver instance. */ -template -auto makeNewtonRaphson(const NonLinearOperatorImpl& p_nonLinearOperator, LinearSolver&& p_linearSolver = {}, - UpdateFunctionType&& p_updateFunction = {}) { - return std::make_shared>( - p_nonLinearOperator, std::forward(p_linearSolver), std::move(p_updateFunction)); +template +auto makeNewtonRaphson(const NLO& nonLinearOperator, LS&& linearSolver = {}, UF&& updateFunction = {}) { + return std::make_shared>(nonLinearOperator, std::forward(linearSolver), + std::move(updateFunction)); } } // namespace Ikarus diff --git a/ikarus/solver/nonlinearsolver/newtonraphsonwithscalarsubsidiaryfunction.hh b/ikarus/solver/nonlinearsolver/newtonraphsonwithscalarsubsidiaryfunction.hh index c0533951f..8d038bbd3 100644 --- a/ikarus/solver/nonlinearsolver/newtonraphsonwithscalarsubsidiaryfunction.hh +++ b/ikarus/solver/nonlinearsolver/newtonraphsonwithscalarsubsidiaryfunction.hh @@ -14,8 +14,6 @@ #include #include #include -#include -#include #include #include @@ -36,46 +34,43 @@ struct NewtonRaphsonWithSubsidiaryFunctionSettings * This class provides a Newton-Raphson solver for solving nonlinear systems with a subsidiary function. * It uses a linear solver to handle the linear system arising in each iteration. * - * \tparam NonLinearOperatorImpl Type of the nonlinear operator. - * \tparam LinearSolver Type of the linear solver used internally (default is SolverDefault). - * \tparam UpdateFunctionTypeImpl Type of the update function (default is UpdateDefault). + * \tparam NLO Type of the nonlinear operator. + * \tparam LS Type of the linear solver used internally (default is SolverDefault). + * \tparam UF Type of the update function (default is UpdateDefault). */ -template +template class NewtonRaphsonWithSubsidiaryFunction : public IObservable { public: ///< Compile-time boolean indicating if the linear solver satisfies the non-linear solver concept static constexpr bool isLinearSolver = - Ikarus::Concepts::LinearSolverCheck; + Ikarus::Concepts::LinearSolverCheck; ///< Type representing the parameter vector of the nonlinear operator. - using ValueType = typename NonLinearOperatorImpl::template ParameterValue<0>; + using ValueType = typename NLO::template ParameterValue<0>; ///< Type representing the update function. - using UpdateFunctionType = UpdateFunctionTypeImpl; - using NonLinearOperator = NonLinearOperatorImpl; ///< Type of the non-linear operator + using UpdateFunctionType = UF; + using NonLinearOperator = NLO; ///< Type of the non-linear operator /** * \brief Constructor for NewtonRaphsonWithSubsidiaryFunction. * - * \param p_nonLinearOperator Nonlinear operator to solve. - * \param p_linearSolver Linear solver used internally (default is SolverDefault). - * \param p_updateFunction Update function (default is UpdateDefault). + * \param nonLinearOperator Nonlinear operator to solve. + * \param linearSolver Linear solver used internally (default is SolverDefault). + * \param updateFunction Update function (default is UpdateDefault). */ - explicit NewtonRaphsonWithSubsidiaryFunction(const NonLinearOperatorImpl& p_nonLinearOperator, - LinearSolver&& p_linearSolver = {}, - UpdateFunctionType p_updateFunction = {}) - : nonLinearOperator_{p_nonLinearOperator}, - linearSolver{std::move(p_linearSolver)}, - updateFunction{p_updateFunction} {} + explicit NewtonRaphsonWithSubsidiaryFunction(const NLO& nonLinearOperator, LS&& linearSolver = {}, + UpdateFunctionType updateFunction = {}) + : nonLinearOperator_{nonLinearOperator}, + linearSolver_{std::move(linearSolver)}, + updateFunction_{updateFunction} {} /** * \brief Setup the Newton-Raphson solver with subsidiary function. * * \param p_settings Settings for the solver. */ - void setup(const NewtonRaphsonWithSubsidiaryFunctionSettings& p_settings) { settings = p_settings; } + void setup(const NewtonRaphsonWithSubsidiaryFunctionSettings& settings) { settings_ = settings; } #ifndef DOXYGEN struct NoPredictor @@ -90,17 +85,17 @@ public: * \tparam SubsidiaryType Type of the subsidiary function. * \param subsidiaryFunction Subsidiary function to be solved. * \param subsidiaryArgs Additional arguments for the subsidiary function. - * \param dx_predictor Predictor for the solution increment (default is NoPredictor). + * \param dxPredictor Predictor for the solution increment (default is NoPredictor). * \return Information about the solution process. */ template requires std::is_same_v || - std::is_convertible_v> + std::is_convertible_v> [[nodiscard( "The solve method returns information of the solution process. You should store this information and check if " "it was successful")]] NonLinearSolverInformation solve(SubsidiaryType& subsidiaryFunction, SubsidiaryArgs& subsidiaryArgs, - const SolutionType& dx_predictor = NoPredictor{}) { + const SolutionType& dxPredictor = NoPredictor{}) { this->notify(NonLinearSolverMessages::INIT); /// Initializations @@ -108,7 +103,7 @@ public: solverInformation.success = true; auto& x = nonLinearOperator().firstParameter(); // x = D (Displacements) if constexpr (not std::is_same_v) - updateFunction(x, dx_predictor); + updateFunction_(x, dxPredictor); auto& lambda = nonLinearOperator().lastParameter(); /// Determine Fext0 @@ -140,10 +135,10 @@ public: decltype(rNorm) dNorm; int iter{0}; if constexpr (isLinearSolver) - linearSolver.analyzePattern(Ax); + linearSolver_.analyzePattern(Ax); /// Iterative solving scheme - while (rNorm > settings.tol && iter < settings.maxIter) { + while (rNorm > settings_.tol && iter < settings_.maxIter) { this->notify(NonLinearSolverMessages::ITERATION_STARTED); /// Two-step solving procedure @@ -152,8 +147,8 @@ public: sol2d.resize(rx.rows(), 2); if constexpr (isLinearSolver) { - linearSolver.factorize(Ax); - linearSolver.solve(sol2d, residual2d); + linearSolver_.factorize(Ax); + linearSolver_.solve(sol2d, residual2d); } else { sol2d = linearSolver(residual2d, Ax); } @@ -164,8 +159,8 @@ public: (subsidiaryArgs.dfdDD.dot(sol2d.col(1)) + subsidiaryArgs.dfdDlambda); deltaD = sol2d.col(0) + deltalambda * sol2d.col(1); - updateFunction(x, deltaD); - updateFunction(subsidiaryArgs.DD, deltaD); + updateFunction_(x, deltaD); + updateFunction_(subsidiaryArgs.DD, deltaD); lambda += deltalambda; subsidiaryArgs.Dlambda += deltalambda; @@ -182,7 +177,7 @@ public: ++iter; } - if (iter == settings.maxIter) + if (iter == settings_.maxIter) solverInformation.success = false; solverInformation.iterations = iter; solverInformation.residualNorm = rNorm; @@ -200,28 +195,26 @@ public: auto& nonLinearOperator() { return nonLinearOperator_; } private: - NonLinearOperatorImpl nonLinearOperator_; - LinearSolver linearSolver; - UpdateFunctionType updateFunction; - NewtonRaphsonWithSubsidiaryFunctionSettings settings; + NLO nonLinearOperator_; + LinearSolver linearSolver_; + UpdateFunctionType updateFunction_; + NewtonRaphsonWithSubsidiaryFunctionSettings settings_; }; /** * \brief Function to create a NewtonRaphson with subsidiary function solver instance. - * \tparam NonLinearOperatorImpl Type of the nonlinear operator to solve. - * \tparam LinearSolver Type of the linear solver used internally (default is SolverDefault). - * \tparam UpdateFunctionType Type of the update function (default is UpdateDefault). - * \param p_nonLinearOperator Nonlinear operator to solve. - * \param p_linearSolver Linear solver used internally (default is SolverDefault). - * \param p_updateFunction Update function (default is UpdateDefault). + * \tparam NLO Type of the nonlinear operator to solve. + * \tparam LS Type of the linear solver used internally (default is SolverDefault). + * \tparam UF Type of the update function (default is UpdateDefault). + * \param nonLinearOperator Nonlinear operator to solve. + * \param linearSolver Linear solver used internally (default is SolverDefault). + * \param updateFunction Update function (default is UpdateDefault). * \return Shared pointer to the NewtonRaphson solver instance. */ -template -auto makeNewtonRaphsonWithSubsidiaryFunction(const NonLinearOperatorImpl& p_nonLinearOperator, - LinearSolver&& p_linearSolver = {}, - UpdateFunctionType&& p_updateFunction = {}) { - return std::make_shared>( - p_nonLinearOperator, std::forward(p_linearSolver), std::move(p_updateFunction)); +template +auto makeNewtonRaphsonWithSubsidiaryFunction(const NLO& nonLinearOperator, LS&& linearSolver = {}, + UF&& updateFunction = {}) { + return std::make_shared>( + nonLinearOperator, std::forward(linearSolver), std::move(updateFunction)); } } // namespace Ikarus diff --git a/ikarus/solver/nonlinearsolver/trustregion.hh b/ikarus/solver/nonlinearsolver/trustregion.hh index 08b13116e..4a759e37c 100644 --- a/ikarus/solver/nonlinearsolver/trustregion.hh +++ b/ikarus/solver/nonlinearsolver/trustregion.hh @@ -113,51 +113,49 @@ struct Stats * This code is heavily inspired by the trust-region implementation of Manopt. * \ingroup solvers -* \tparam NonLinearOperatorImpl Type of the nonlinear operator to solve. +* \tparam NLO Type of the nonlinear operator to solve. * \tparam preConditioner Type of preconditioner to use (default is IncompleteCholesky). -* \tparam UpdateFunctionTypeImpl Type of the update function +* \tparam UF Type of the update function */ -template +template class TrustRegion : public IObservable { public: - using ValueType = typename NonLinearOperatorImpl::template ParameterValue<0>; ///< Type of the parameter vector of - ///< the nonlinear operator - using CorrectionType = typename NonLinearOperatorImpl::DerivativeType; ///< Type of the correction of x += deltaX. - using UpdateFunctionType = UpdateFunctionTypeImpl; ///< Type of the update function. + using ValueType = typename NLO::template ParameterValue<0>; ///< Type of the parameter vector of + ///< the nonlinear operator + using CorrectionType = typename NLO::DerivativeType; ///< Type of the correction of x += deltaX. + using UpdateFunction = UF; ///< Type of the update function. - using NonLinearOperator = NonLinearOperatorImpl; ///< Type of the non-linear operator + using NonLinearOperator = NLO; ///< Type of the non-linear operator - using ScalarType = - std::remove_cvref_t>; ///< Type of the scalar - ///< cost + using ScalarType = std::remove_cvref_t>; ///< Type of the scalar + ///< cost - using MatrixType = - std::remove_cvref_t>; ///< Type of the Hessian + using MatrixType = std::remove_cvref_t>; ///< Type of the Hessian /** * \brief Constructs a TrustRegion solver instance. - * \param p_nonLinearOperator Nonlinear operator to solve. - * \param p_updateFunction Update function + * \param nonLinearOperator Nonlinear operator to solve. + * \param updateFunction Update function */ - explicit TrustRegion(const NonLinearOperatorImpl& p_nonLinearOperator, UpdateFunctionTypeImpl p_updateFunction = {}) - : nonLinearOperator_{p_nonLinearOperator}, - updateFunction{p_updateFunction}, - xOld{nonLinearOperator().firstParameter()} { - eta.setZero(gradient().size()); - Heta.setZero(gradient().size()); + explicit TrustRegion(const NLO& nonLinearOperator, UF updateFunction = {}) + : nonLinearOperator_{nonLinearOperator}, + updateFunction_{updateFunction}, + xOld_{this->nonLinearOperator().firstParameter()} { + eta_.setZero(gradient().size()); + Heta_.setZero(gradient().size()); } /** * \brief Sets up the TrustRegion solver with the provided settings and checks feasibility. * \param p_settings TrustRegionSettings containing the solver configuration. */ - void setup(const TrustRegionSettings& p_settings) { - options = p_settings; - assert(options.rho_prime < 0.25 && "options.rho_prime must be strictly smaller than 1/4."); - assert(options.Delta_bar > 0 && "options.Delta_bar must be positive."); - assert(options.Delta0 > 0 && options.Delta0 < options.Delta_bar && + void setup(const TrustRegionSettings& settings) { + options_ = settings; + assert(options_.rho_prime < 0.25 && "options.rho_prime must be strictly smaller than 1/4."); + assert(options_.Delta_bar > 0 && "options.Delta_bar must be positive."); + assert(options_.Delta0 > 0 && options_.Delta0 < options_.Delta_bar && "options.Delta0 must be positive and smaller than Delta_bar."); } @@ -169,190 +167,190 @@ public: /** * \brief Solves the nonlinear optimization problem using the TrustRegion algorithm. * \tparam SolutionType Type of the solution predictor (default is NoPredictor). - * \param dx_predictor Solution predictor. + * \param dxPredictor Solution predictor. * \return NonLinearSolverInformation containing information about the solver result. */ template requires std::is_same_v || std::is_convertible_v - NonLinearSolverInformation solve(const SolutionType& dx_predictor = NoPredictor{}) { + NonLinearSolverInformation solve(const SolutionType& dxPredictor = NoPredictor{}) { this->notify(NonLinearSolverMessages::INIT); - stats = Stats{}; - info = AlgoInfo{}; + stats_ = Stats{}; + info_ = AlgoInfo{}; NonLinearSolverInformation solverInformation; nonLinearOperator().updateAll(); - stats.energy = energy(); - auto& x = nonLinearOperator().firstParameter(); - xOld = x; - stats.gradNorm = norm(gradient()); + stats_.energy = energy(); + auto& x = nonLinearOperator().firstParameter(); + xOld_ = x; + stats_.gradNorm = norm(gradient()); if constexpr (not std::is_same_v) - updateFunction(x, dx_predictor); - truncatedConjugateGradient.analyzePattern(hessian()); + updateFunction(x, dxPredictor); + truncatedConjugateGradient_.analyzePattern(hessian()); - innerInfo.Delta = options.Delta0; + innerInfo_.Delta = options_.Delta0; spdlog::info( " | iter | inner_i | rho | energy | energy_p | energy_inc | norm(g) | Delta | norm(corr) | " "InnerBreakReason"); spdlog::info("{:-^143}", "-"); while (not stoppingCriterion()) { this->notify(NonLinearSolverMessages::ITERATION_STARTED); - if (options.useRand) { - if (stats.outerIter == 0) { - eta.setRandom(); - while (eta.dot(hessian() * eta) > innerInfo.Delta * innerInfo.Delta) - eta *= eps; // eps is sqrt(sqrt(maschine-precision)) + if (options_.useRand) { + if (stats_.outerIter == 0) { + eta_.setRandom(); + while (eta_.dot(hessian() * eta_) > innerInfo_.Delta * innerInfo_.Delta) + eta_ *= eps_; // eps is sqrt(sqrt(maschine-precision)) } else - eta.setZero(); + eta_.setZero(); } else - eta.setZero(); + eta_.setZero(); solveInnerProblem(); - stats.innerIterSum += innerInfo.numInnerIter; - - info.stopReasonString = tcg_stop_reason[static_cast(innerInfo.stop_tCG)]; - Heta = hessian() * eta; - if (options.useRand and stats.outerIter == 0) { - info.used_cauchy = false; - info.randomPredictionString = " Used Random correction predictor"; - info.cauchystr = " "; - double tau_c; + stats_.innerIterSum += innerInfo_.numInnerIter; + + info_.stopReasonString = tcg_stop_reason_[static_cast(innerInfo_.stop_tCG)]; + Heta_ = hessian() * eta_; + if (options_.useRand and stats_.outerIter == 0) { + info_.used_cauchy = false; + info_.randomPredictionString = " Used Random correction predictor"; + info_.cauchystr = " "; + double tauC; // Check the curvature const Eigen::VectorXd Hg = hessian() * gradient(); const auto g_Hg = (gradient().dot(Hg)); if (g_Hg <= 0) - tau_c = 1; + tauC = 1; else - tau_c = std::min(Dune::power(stats.gradNorm, 3) / (innerInfo.Delta * g_Hg), 1.0); + tauC = std::min(Dune::power(stats_.gradNorm, 3) / (innerInfo_.Delta * g_Hg), 1.0); // generate the Cauchy point. - const Eigen::VectorXd eta_c = -tau_c * innerInfo.Delta / stats.gradNorm * gradient(); - const Eigen::VectorXd Heta_c = -tau_c * innerInfo.Delta / stats.gradNorm * Hg; + const Eigen::VectorXd etaC = -tauC * innerInfo_.Delta / stats_.gradNorm * gradient(); + const Eigen::VectorXd HetaC = -tauC * innerInfo_.Delta / stats_.gradNorm * Hg; - const double mdle = stats.energy + gradient().dot(eta) + .5 * Heta.dot(eta); - const double mdlec = stats.energy + gradient().dot(eta_c) + .5 * Heta_c.dot(eta_c); - if (mdlec < mdle && stats.outerIter == 0) { - eta = eta_c; - Heta = Heta_c; - info.used_cauchy = true; + const double mdle = stats_.energy + gradient().dot(eta_) + .5 * Heta_.dot(eta_); + const double mdlec = stats_.energy + gradient().dot(etaC) + .5 * HetaC.dot(etaC); + if (mdlec < mdle && stats_.outerIter == 0) { + eta_ = etaC; + Heta_ = HetaC; + info_.used_cauchy = true; - info.cauchystr = " USED CAUCHY POINT!"; + info_.cauchystr = " USED CAUCHY POINT!"; } } else - info.cauchystr = " "; + info_.cauchystr = " "; - stats.etaNorm = eta.norm(); + stats_.etaNorm = eta_.norm(); - updateFunction(x, eta); + updateFunction_(x, eta_); // Calculate energy of our proposed update step nonLinearOperator().template update<0>(); - stats.energyProposal = energy(); + stats_.energyProposal = energy(); // Will we accept the proposal or not? // Check the performance of the quadratic model against the actual energy. - auto rhonum = stats.energy - stats.energyProposal; - auto rhoden = -eta.dot(gradient() + 0.5 * Heta); + auto rhonum = stats_.energy - stats_.energyProposal; + auto rhoden = -eta_.dot(gradient() + 0.5 * Heta_); /* Close to convergence the proposed energy and the real energy almost coincide. * Therefore, the performance check of our model becomes ill-conditioned * The regularisation fixes this */ - const auto rho_reg = std::max(1.0, abs(stats.energy)) * eps * options.rho_reg; - rhonum = rhonum + rho_reg; - rhoden = rhoden + rho_reg; + const auto rhoReg = std::max(1.0, abs(stats_.energy)) * eps_ * options_.rho_reg; + rhonum = rhonum + rhoReg; + rhoden = rhoden + rhoReg; - const bool model_decreased = rhoden > 0.0; + const bool modelDecreased = rhoden > 0.0; - if (!model_decreased) - info.stopReasonString.append(", model did not decrease"); + if (!modelDecreased) + info_.stopReasonString.append(", model did not decrease"); - stats.rho = rhonum / rhoden; - stats.rho = stats.rho < 0.0 ? -1.0 : stats.rho; // move rho to the domain [-1.0,inf] + stats_.rho = rhonum / rhoden; + stats_.rho = stats_.rho < 0.0 ? -1.0 : stats_.rho; // move rho to the domain [-1.0,inf] - info.trstr = " "; + info_.trstr = " "; // measure if energy decreased - const bool energyDecreased = Dune::FloatCmp::ge(stats.energy - stats.energyProposal, -1e-12); + const bool energyDecreased = Dune::FloatCmp::ge(stats_.energy - stats_.energyProposal, -1e-12); // If the model behaves badly or if the energy increased we reduce the trust region size - if (stats.rho < 1e-4 || not model_decreased || std::isnan(stats.rho) || not energyDecreased) { - info.trstr = "TR-"; - innerInfo.Delta /= 4.0; - info.consecutive_TRplus = 0; - info.consecutive_TRminus++; - if (info.consecutive_TRminus >= 5 && options.verbosity >= 1) { - info.consecutive_TRminus = -std::numeric_limits::infinity(); + if (stats_.rho < 1e-4 || not modelDecreased || std::isnan(stats_.rho) || not energyDecreased) { + info_.trstr = "TR-"; + innerInfo_.Delta /= 4.0; + info_.consecutive_TRplus = 0; + info_.consecutive_TRminus++; + if (info_.consecutive_TRminus >= 5 && options_.verbosity >= 1) { + info_.consecutive_TRminus = -std::numeric_limits::infinity(); spdlog::info(" +++ Detected many consecutive TR- (radius decreases)."); spdlog::info(" +++ Consider decreasing options.Delta_bar by an order of magnitude."); } - } else if (stats.rho > 0.99 && (innerInfo.stop_tCG == Eigen::TCGStopReason::negativeCurvature || - innerInfo.stop_tCG == Eigen::TCGStopReason::exceededTrustRegion)) { - info.trstr = "TR+"; - innerInfo.Delta = std::min(3.5 * innerInfo.Delta, options.Delta_bar); - info.consecutive_TRminus = 0; - info.consecutive_TRplus++; - if (info.consecutive_TRplus >= 5 && options.verbosity >= 1) { - info.consecutive_TRplus = -std::numeric_limits::infinity(); + } else if (stats_.rho > 0.99 && (innerInfo_.stop_tCG == Eigen::TCGStopReason::negativeCurvature || + innerInfo_.stop_tCG == Eigen::TCGStopReason::exceededTrustRegion)) { + info_.trstr = "TR+"; + innerInfo_.Delta = std::min(3.5 * innerInfo_.Delta, options_.Delta_bar); + info_.consecutive_TRminus = 0; + info_.consecutive_TRplus++; + if (info_.consecutive_TRplus >= 5 && options_.verbosity >= 1) { + info_.consecutive_TRplus = -std::numeric_limits::infinity(); spdlog::info(" +++ Detected many consecutive TR+ (radius increases)"); spdlog::info(" +++ Consider increasing options.Delta_bar by an order of magnitude"); } else { - info.consecutive_TRplus = 0; - info.consecutive_TRminus = 0; + info_.consecutive_TRplus = 0; + info_.consecutive_TRminus = 0; } } - if (model_decreased && stats.rho > options.rho_prime && energyDecreased) { - if (stats.energyProposal > stats.energy) + if (modelDecreased && stats_.rho > options_.rho_prime && energyDecreased) { + if (stats_.energyProposal > stats_.energy) spdlog::info( "Energy function increased by {} (step size: {}). Since this is small we accept the step and hope for " "convergence of the gradient norm.", - stats.energyProposal - stats.energy, stats.etaNorm); + stats_.energyProposal - stats_.energy, stats_.etaNorm); - info.acceptProposal = true; - info.accstr = "acc"; - info.consecutive_Rejected = 0; + info_.acceptProposal = true; + info_.accstr = "acc"; + info_.consecutive_Rejected = 0; } else { - info.acceptProposal = false; - info.accstr = "REJ"; + info_.acceptProposal = false; + info_.accstr = "REJ"; - if (info.consecutive_Rejected >= 5) - innerInfo.Delta /= 2; + if (info_.consecutive_Rejected >= 5) + innerInfo_.Delta /= 2; else - innerInfo.Delta = std::min(innerInfo.Delta, stats.etaNorm / 2.0); - ++info.consecutive_Rejected; + innerInfo_.Delta = std::min(innerInfo_.Delta, stats_.etaNorm / 2.0); + ++info_.consecutive_Rejected; } - stats.outerIter++; + stats_.outerIter++; - if (options.verbosity == 1) + if (options_.verbosity == 1) logState(); - info.randomPredictionString = ""; + info_.randomPredictionString = ""; - if (info.acceptProposal) { - stats.energy = stats.energyProposal; + if (info_.acceptProposal) { + stats_.energy = stats_.energyProposal; nonLinearOperator_.updateAll(); - xOld = x; - this->notify(NonLinearSolverMessages::CORRECTIONNORM_UPDATED, stats.etaNorm); - this->notify(NonLinearSolverMessages::RESIDUALNORM_UPDATED, stats.gradNorm); + xOld_ = x; + this->notify(NonLinearSolverMessages::CORRECTIONNORM_UPDATED, stats_.etaNorm); + this->notify(NonLinearSolverMessages::RESIDUALNORM_UPDATED, stats_.gradNorm); this->notify(NonLinearSolverMessages::SOLUTION_CHANGED); } else { - x = xOld; - eta.setZero(); + x = xOld_; + eta_.setZero(); } nonLinearOperator_.updateAll(); - stats.gradNorm = gradient().norm(); + stats_.gradNorm = gradient().norm(); this->notify(NonLinearSolverMessages::ITERATION_ENDED); } - spdlog::info("{}", info.reasonString); - spdlog::info("Total iterations: {} Total CG Iterations: {}", stats.outerIter, stats.innerIterSum); + spdlog::info("{}", info_.reasonString); + spdlog::info("Total iterations: {} Total CG Iterations: {}", stats_.outerIter, stats_.innerIterSum); solverInformation.success = - (info.stop == StopReason::correctionNormTolReached) or (info.stop == StopReason::gradientNormTolReached); + (info_.stop == StopReason::correctionNormTolReached) or (info_.stop == StopReason::gradientNormTolReached); - solverInformation.iterations = stats.outerIter; - solverInformation.residualNorm = stats.gradNorm; + solverInformation.iterations = stats_.outerIter; + solverInformation.residualNorm = stats_.gradNorm; if (solverInformation.success) this->notify(NonLinearSolverMessages::FINISHED_SUCESSFULLY, solverInformation.iterations); return solverInformation; @@ -368,15 +366,15 @@ private: spdlog::info( "{:>3s} {:>3s} {:>6d} {:>9d} {:>6.2f} {:>9.2e} {:>9.2e} {:>11.2e} {:>9.2e} {:>9.2e} {:>11.2e} " "{:<73}", - info.accstr, info.trstr, stats.outerIter, innerInfo.numInnerIter, stats.rho, stats.energy, stats.energyProposal, - stats.energyProposal - stats.energy, stats.gradNorm, innerInfo.Delta, stats.etaNorm, - info.stopReasonString + info.cauchystr + info.randomPredictionString); + info_.accstr, info_.trstr, stats_.outerIter, innerInfo_.numInnerIter, stats_.rho, stats_.energy, + stats_.energyProposal, stats_.energyProposal - stats_.energy, stats_.gradNorm, innerInfo_.Delta, stats_.etaNorm, + info_.stopReasonString + info_.cauchystr + info_.randomPredictionString); } void logFinalState() { spdlog::info("{:>3s} {:>3s} {:>6d} {:>9d} {: ^6} {: ^9} {: ^9} {: ^11} {:>9.2e} {: ^9} {: ^11} {:<73}", - info.accstr, info.trstr, stats.outerIter, innerInfo.numInnerIter, " ", " ", " ", " ", stats.gradNorm, - " ", " ", info.stopReasonString + info.cauchystr + info.randomPredictionString); + info_.accstr, info_.trstr, stats_.outerIter, innerInfo_.numInnerIter, " ", " ", " ", " ", + stats_.gradNorm, " ", " ", info_.stopReasonString + info_.cauchystr + info_.randomPredictionString); } inline const auto& energy() { return nonLinearOperator().value(); } @@ -386,82 +384,82 @@ private: bool stoppingCriterion() { std::ostringstream stream; /** Gradient correction tolerance reached */ - if (stats.gradNorm < options.grad_tol && stats.outerIter != 0) { + if (stats_.gradNorm < options_.grad_tol && stats_.outerIter != 0) { logFinalState(); spdlog::info("CONVERGENCE: Energy: {:1.16e} norm(gradient): {:1.16e}", nonLinearOperator().value(), - stats.gradNorm); - stream << "Gradient norm tolerance reached; options.tolerance = " << options.grad_tol; + stats_.gradNorm); + stream << "Gradient norm tolerance reached; options.tolerance = " << options_.grad_tol; - info.reasonString = stream.str(); + info_.reasonString = stream.str(); - info.stop = StopReason::gradientNormTolReached; + info_.stop = StopReason::gradientNormTolReached; return true; - } else if (stats.etaNorm < options.corr_tol && stats.outerIter != 0) { + } else if (stats_.etaNorm < options_.corr_tol && stats_.outerIter != 0) { logFinalState(); spdlog::info("CONVERGENCE: Energy: {:1.16e} norm(correction): {:1.16e}", nonLinearOperator().value(), - stats.etaNorm); - stream << "Displacement norm tolerance reached; = " << options.corr_tol << "." << std::endl; + stats_.etaNorm); + stream << "Displacement norm tolerance reached; = " << options_.corr_tol << "." << std::endl; - info.reasonString = stream.str(); - info.stop = StopReason::correctionNormTolReached; + info_.reasonString = stream.str(); + info_.stop = StopReason::correctionNormTolReached; return true; } /** Maximum Time reached */ - if (stats.time >= options.maxtime) { + if (stats_.time >= options_.maxtime) { logFinalState(); - stream << "Max time exceeded; options.maxtime = " << options.maxtime << "."; - info.reasonString = stream.str(); - info.stop = StopReason::maximumTimeReached; + stream << "Max time exceeded; options.maxtime = " << options_.maxtime << "."; + info_.reasonString = stream.str(); + info_.stop = StopReason::maximumTimeReached; return true; } /** Maximum Iterations reached */ - if (stats.outerIter >= options.maxiter) { + if (stats_.outerIter >= options_.maxiter) { logFinalState(); - stream << "Max iteration count reached; options.maxiter = " << options.maxiter << "."; - info.reasonString = stream.str(); - info.stop = StopReason::maximumIterationsReached; + stream << "Max iteration count reached; options.maxiter = " << options_.maxiter << "."; + info_.reasonString = stream.str(); + info_.stop = StopReason::maximumIterationsReached; return true; } return false; } void solveInnerProblem() { - truncatedConjugateGradient.setInfo(innerInfo); + truncatedConjugateGradient_.setInfo(innerInfo_); int attempts = 0; - truncatedConjugateGradient.factorize(hessian()); + truncatedConjugateGradient_.factorize(hessian()); // If the preconditioner is IncompleteCholesky the factorization may fail if we have negative diagonal entries and // the initial shift is too small. Therefore, if the factorization fails we increase the initial shift by a factor // of 5. if constexpr (preConditioner == PreConditioner::IncompleteCholesky) { - while (truncatedConjugateGradient.info() != Eigen::Success) { - choleskyInitialShift *= 5; - truncatedConjugateGradient.preconditioner().setInitialShift(choleskyInitialShift); - truncatedConjugateGradient.factorize(hessian()); + while (truncatedConjugateGradient_.info() != Eigen::Success) { + choleskyInitialShift_ *= 5; + truncatedConjugateGradient_.preconditioner().setInitialShift(choleskyInitialShift_); + truncatedConjugateGradient_.factorize(hessian()); if (attempts > 5) DUNE_THROW(Dune::MathError, "Factorization of preconditioner failed!"); ++attempts; } - if (truncatedConjugateGradient.info() == Eigen::Success) - choleskyInitialShift = 1e-3; + if (truncatedConjugateGradient_.info() == Eigen::Success) + choleskyInitialShift_ = 1e-3; } - eta = truncatedConjugateGradient.solveWithGuess(-gradient(), eta); - innerInfo = truncatedConjugateGradient.getInfo(); + eta_ = truncatedConjugateGradient_.solveWithGuess(-gradient(), eta_); + innerInfo_ = truncatedConjugateGradient_.getInfo(); } - NonLinearOperatorImpl nonLinearOperator_; - UpdateFunctionType updateFunction; - typename NonLinearOperatorImpl::template ParameterValue<0> xOld; - CorrectionType eta; - CorrectionType Heta; - TrustRegionSettings options; - AlgoInfo info; - double choleskyInitialShift = 1e-3; - Eigen::TCGInfo innerInfo; - Stats stats; - static constexpr double eps = 0.0001220703125; // 0.0001220703125 is sqrt(sqrt(maschine-precision)) - std::array tcg_stop_reason{ + NLO nonLinearOperator_; + UpdateFunction updateFunction_; + typename NLO::template ParameterValue<0> xOld_; + CorrectionType eta_; + CorrectionType Heta_; + TrustRegionSettings options_; + AlgoInfo info_; + double choleskyInitialShift_ = 1e-3; + Eigen::TCGInfo innerInfo_; + Stats stats_; + static constexpr double eps_ = 0.0001220703125; // 0.0001220703125 is sqrt(sqrt(maschine-precision)) + std::array tcg_stop_reason_{ {"negative curvature", "exceeded trust region", "reached target residual-kappa (linear)", "reached target residual-theta (superlinear)", "maximum inner iterations", "model increased"} }; @@ -472,24 +470,23 @@ private: typename Eigen::DiagonalPreconditioner, typename Eigen::IncompleteCholesky>>; Eigen::TruncatedConjugateGradient - truncatedConjugateGradient; + truncatedConjugateGradient_; }; /** - * @brief Creates an instance of the TrustRegion solver. + * \brief Creates an instance of the TrustRegion solver. * - * @tparam NonLinearOperatorImpl Type of the nonlinear operator to solve. - * @tparam preConditioner Type of the preconditioner used internally (default is IncompleteCholesky). - * @tparam UpdateFunctionType Type of the update function (default is UpdateDefault). - * @param p_nonLinearOperator Nonlinear operator to solve. - * @param p_updateFunction Update function (default is UpdateDefault). - * @return Shared pointer to the TrustRegion solver instance. + * \tparam NLO Type of the nonlinear operator to solve. + * \tparam preConditioner Type of the preconditioner used internally (default is IncompleteCholesky). + * \tparam UF Type of the update function (default is UpdateDefault). + * \param nonLinearOperator Nonlinear operator to solve. + * \param updateFunction Update function (default is UpdateDefault). + * \return Shared pointer to the TrustRegion solver instance. */ -template -auto makeTrustRegion(const NonLinearOperatorImpl& p_nonLinearOperator, UpdateFunctionType&& p_updateFunction = {}) { - return std::make_shared>(p_nonLinearOperator, - p_updateFunction); +template +auto makeTrustRegion(const NLO& nonLinearOperator, UF&& updateFunction = {}) { + return std::make_shared>(nonLinearOperator, updateFunction); } } // namespace Ikarus diff --git a/ikarus/utils/algorithms.hh b/ikarus/utils/algorithms.hh index c529fce9e..65a253117 100644 --- a/ikarus/utils/algorithms.hh +++ b/ikarus/utils/algorithms.hh @@ -17,13 +17,13 @@ namespace Ikarus::utils { /** - * @brief Sorts and removes duplicate elements from a random access range. + * \brief Sorts and removes duplicate elements from a random access range. * \ingroup algos * *\details * This function sorts the elements of the given random access range and removes duplicate elements, * leaving only unique elements in the range. * - * @param r The random access range to be modified. + * \param r The random access range to be modified. */ void makeUniqueAndSort(std::ranges::random_access_range auto& r) { sort(r.begin(), r.end()); @@ -31,65 +31,65 @@ void makeUniqueAndSort(std::ranges::random_access_range auto& r) { } /** - * @brief Appends a value to the range if it is not already present. + * \brief Appends a value to the range if it is not already present. *\details * This function appends a value to the given random access range only if the value is not already present in the *range. It returns the index of the value in the range, whether it was added or already existed. \ingroup algos - * @tparam Value The type of the value to be appended. - * @param r The random access range to be modified. - * @param v The value to be appended. - * @return The index of the value in the range. + * \tparam T The type of the value to be appended. + * \param r The random access range to be modified. + * \param v The value to be appended. + * \return The index of the value in the range. */ -template -auto appendUnique(std::ranges::random_access_range auto& r, Value&& v) { +template +auto appendUnique(std::ranges::random_access_range auto& r, T&& v) { static_assert(std::is_same_v>); const auto it = find(begin(r), end(r), v); size_t index = std::distance(begin(r), it); if (it == end(r)) - r.push_back(std::forward(v)); + r.push_back(std::forward(v)); return index; } /** - * @brief Prints the contents of a container to the specified output stream. + * \brief Prints the contents of a container to the specified output stream. *\details * This function prints the contents of the given container to the specified output stream. * \ingroup algos - * @tparam Container The type of the container to be printed. - * @param c The container whose contents will be printed. - * @param os The output stream where the contents will be printed. Default is std::cout. + * \tparam C The type of the container to be printed. + * \param c The container whose contents will be printed. + * \param os The output stream where the contents will be printed. Default is std::cout. */ -template -void printContent(Container&& c, std::ostream& os = std::cout) { +template +void printContent(C&& c, std::ostream& os = std::cout) { std::ranges::for_each(c, [&os](auto&& var) { os << var << '\n'; }); } /** - * @brief Transforms a value range to a pointer range. + * \brief Transforms a value range to a pointer range. * * \ingroup algos - * @tparam Container The type of the container containing values. - * @param cont The container whose values will be transformed to pointers. - * @return A subrange containing pointers to the values. + * \tparam C The type of the container containing values. + * \param cont The container whose values will be transformed to pointers. + * \return A subrange containing pointers to the values. */ -template -auto transformValueRangeToPointerRange(Container& cont) { +template +auto transformValueRangeToPointerRange(C& cont) { auto transformValueToPointer = [](auto&& obj) { return &obj; }; return (std::ranges::subrange(cont.begin(), cont.end()) | std::views::transform(transformValueToPointer)); } /** - * @brief Transforms a pointer range to a reference range. + * \brief Transforms a pointer range to a reference range. *\details * This function transforms a range of pointers to a range of references. * \ingroup algos - * @tparam Container The type of the container containing pointers. - * @param cont The container whose pointers will be transformed to references. - * @return A subrange containing references to the pointed objects. + * \tparam C The type of the container containing pointers. + * \param cont The container whose pointers will be transformed to references. + * \return A subrange containing references to the pointed objects. */ -template -auto transformPointerRangeToReferenceRange(Container& cont) { +template +auto transformPointerRangeToReferenceRange(C& cont) { auto transformValueToPointer = [](auto&& obj) -> auto& { return *obj; }; return (std::ranges::subrange(cont.begin(), cont.end()) | std::views::transform(transformValueToPointer)); } @@ -202,13 +202,13 @@ namespace Impl { } // namespace Impl /** - * @brief Finds the index of the first element in the tuple satisfying a predicate. + * \brief Finds the index of the first element in the tuple satisfying a predicate. * - * @tparam Tuple Type of the input tuple. - * @tparam Predicate Type of the predicate function. - * @param tuple Input tuple. - * @param pred Predicate function to check each element. - * @return Index of the first element satisfying the predicate. If no element satisfies the + * \tparam Tuple Type of the input tuple. + * \tparam Predicate Type of the predicate function. + * \param tuple Input tuple. + * \param pred Predicate function to check each element. + * \return Index of the first element satisfying the predicate. If no element satisfies the * predicate, it returns the size of the tuple. * * \details This function takes a tuple and a predicate function and finds the index of the first element in the @@ -233,13 +233,13 @@ constexpr size_t find_if(Tuple&& tuple, Predicate pred) { } /** - * @brief Checks if none of the elements in the tuple satisfy a given predicate. + * \brief Checks if none of the elements in the tuple satisfy a given predicate. * - * @tparam Tuple Type of the input tuple. - * @tparam Predicate Type of the predicate function. - * @param tuple Input tuple. - * @param pred Predicate function to check each element. - * @return bool True if none of the elements satisfy the predicate, false otherwise. + * \tparam Tuple Type of the input tuple. + * \tparam Predicate Type of the predicate function. + * \param tuple Input tuple. + * \param pred Predicate function to check each element. + * \return bool True if none of the elements satisfy the predicate, false otherwise. * * * \ingroup algos @@ -250,13 +250,13 @@ bool none_of(Tuple&& tuple, Predicate pred) { } /** - * @brief Checks if any of the elements in the tuple satisfy a given predicate. + * \brief Checks if any of the elements in the tuple satisfy a given predicate. * - * @tparam Tuple Type of the input tuple. - * @tparam Predicate Type of the predicate function. - * @param tuple Input tuple. - * @param pred Predicate function to check each element. - * @return bool True if any of the elements satisfy the predicate, false otherwise. + * \tparam Tuple Type of the input tuple. + * \tparam Predicate Type of the predicate function. + * \param tuple Input tuple. + * \param pred Predicate function to check each element. + * \return bool True if any of the elements satisfy the predicate, false otherwise. * * * \ingroup algos @@ -267,13 +267,13 @@ bool any_of(Tuple&& tuple, Predicate pred) { } /** - * @brief Filters the elements of a tuple based on a given predicate. + * \brief Filters the elements of a tuple based on a given predicate. * - * @tparam Tuple Type of the input tuple. - * @tparam Predicate Type of the predicate function. - * @param tuple Input tuple. - * @param pred Predicate function to filter the elements. - * @return auto Tuple containing elements that satisfy the predicate. + * \tparam Tuple Type of the input tuple. + * \tparam Predicate Type of the predicate function. + * \param tuple Input tuple. + * \param pred Predicate function to filter the elements. + * \return auto Tuple containing elements that satisfy the predicate. * * \details This function applies the given predicate to each element of the tuple. It constructs a new tuple * containing only those elements for which the predicate returns true. The resulting tuple is returned. @@ -290,11 +290,11 @@ auto filter(Tuple&& tuple, Predicate pred) { } /** - * @brief Creates a tuple with unique types from the given tuple. + * \brief Creates a tuple with unique types from the given tuple. * - * @tparam Types Variadic template parameters representing types. - * @param tuple Input tuple. - * @return auto Tuple with unique types. + * \tparam Types Variadic template parameters representing types. + * \param tuple Input tuple. + * \return auto Tuple with unique types. * * \details This function takes a tuple and returns a new tuple containing only unique types from the input tuple. * @@ -306,13 +306,13 @@ constexpr auto unique([[maybe_unused]] std::tuple&& tuple) { } /** - * @brief Counts the number of elements in the tuple satisfying the given predicate. + * \brief Counts the number of elements in the tuple satisfying the given predicate. * - * @tparam Tuple Type of the tuple. - * @tparam Predicate Predicate function determining whether an element satisfies the condition. - * @param tuple Input tuple. - * @param pred Predicate function. - * @return constexpr size_t Number of elements satisfying the predicate. + * \tparam Tuple Type of the tuple. + * \tparam Predicate Predicate function determining whether an element satisfies the condition. + * \param tuple Input tuple. + * \param pred Predicate function. + * \return constexpr size_t Number of elements satisfying the predicate. * * \details This function counts the number of elements in the tuple that satisfy the given predicate. * @@ -329,11 +329,11 @@ constexpr size_t count_if(Tuple&& tuple, Predicate pred) { } /** - * @brief Finds the index of the first element in the tuple that is a specialization of the given template type. + * \brief Finds the index of the first element in the tuple that is a specialization of the given template type. * - * @tparam Type Template type to search for. - * @tparam Tuple Type of the tuple. - * @return Index of the first specialization in the tuple + * \tparam Type Template type to search for. + * \tparam Tuple Type of the tuple. + * \return Index of the first specialization in the tuple * * \details This function finds the index of the first element in the tuple that is a specialization of the given * template type. It returns the size of the tuple if the template type is not found @@ -347,12 +347,12 @@ constexpr int findTypeSpecialization() { } /** - * @brief Gets the specialization of the given template type from the tuple. + * \brief Gets the specialization of the given template type from the tuple. * - * @tparam Type Template type to search for. - * @tparam Tuple Type of the tuple. - * @param tuple The tuple containing elements. - * @return The specialization element + * \tparam Type Template type to search for. + * \tparam Tuple Type of the tuple. + * \param tuple The tuple containing elements. + * \return The specialization element * * \details This function retrieves the specialization of a template type from the tuple. * @@ -367,11 +367,11 @@ auto getSpecialization(Tuple&& tuple) { } /** - * @brief Checks if a tuple has a specialization of a template type. + * \brief Checks if a tuple has a specialization of a template type. * - * @tparam Type Template type to check for. - * @tparam Tuple Type of the tuple. - * @return true if the tuple has the specialization; otherwise, false. + * \tparam Type Template type to check for. + * \tparam Tuple Type of the tuple. + * \return true if the tuple has the specialization; otherwise, false. * * \details This function checks if a tuple has a specialization of a template type. * It uses `find_if` to search for the type and returns true if the index is less than the tuple size; otherwise, @@ -387,11 +387,11 @@ constexpr bool hasTypeSpecialization() { } /** - * @brief Counts the occurrences of a specialization of a template type in a tuple. + * \brief Counts the occurrences of a specialization of a template type in a tuple. * - * @tparam Type Template type to count occurrences for. - * @tparam Tuple Type of the tuple. - * @return The count of occurrences of the specialization. + * \tparam Type Template type to count occurrences for. + * \tparam Tuple Type of the tuple. + * \return The count of occurrences of the specialization. * * \details This function counts the occurrences of a specialization of a template type in a tuple. * @@ -404,10 +404,10 @@ constexpr bool countTypeSpecialization() { } /** - * @brief Variable template for counting the occurrences of a specialization of a template type in a tuple. + * \brief Variable template for counting the occurrences of a specialization of a template type in a tuple. * - * @tparam Type Template type to count occurrences for. - * @tparam Tuple Type of the tuple. + * \tparam Type Template type to count occurrences for. + * \tparam Tuple Type of the tuple. * * \details This variable template provides a compile-time constant for the count of occurrences * of a specialization of a template type in a tuple. @@ -418,12 +418,12 @@ template