From 83860fd41d01a11a20fa265e15f23e29235a8768 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Mon, 5 Sep 2016 19:28:37 +0200 Subject: [PATCH 1/7] [PL] Clarify KelvinV/M types in BMatrixPolicy. Especially point out the difference between the KelvinVector/MatrixType from ProcessLib. --- ProcessLib/Deformation/BMatrixPolicy.h | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/ProcessLib/Deformation/BMatrixPolicy.h b/ProcessLib/Deformation/BMatrixPolicy.h index eaee10bc612..612c834eeaa 100644 --- a/ProcessLib/Deformation/BMatrixPolicy.h +++ b/ProcessLib/Deformation/BMatrixPolicy.h @@ -39,13 +39,15 @@ struct KelvinVectorDimensions<3> /// Kelvin vector type for given displacement dimension. /// \note The Eigen vector is always a fixed size vector in contrast to the -/// BMatrixPolicyType::StressVectorType. +/// BMatrixPolicyType::KelvinVectorType. template using KelvinVectorType = Eigen::Matrix::value, 1, Eigen::ColMajor>; /// Kelvin matrix type for given displacement dimension. +/// \note The Eigen matrix is always a fixed size matrix in contrast to the +/// BMatrixPolicyType::KelvinMatrixType. template using KelvinMatrixType = Eigen::Matrix::value, @@ -80,11 +82,13 @@ class BMatrixPolicyType /// Rhs residual using NodalForceVectorType = VectorType<_number_of_dof>; - /// Sigma - using StressVectorType = VectorType<_kelvin_vector_size>; + /// This type can be different (fixed vs. dynamic size) from the + /// ProcessLib::KelvinVectorType. + using KelvinVectorType = VectorType<_kelvin_vector_size>; - /// C - using ModulusMatrixType = + /// This type can be different (fixed vs. dynamic size) from the + /// ProcessLib::KelvinMatrixType. + using KelvinMatrixType = MatrixType<_kelvin_vector_size, _kelvin_vector_size>; using BMatrixType = MatrixType<_kelvin_vector_size, _number_of_dof>; From 77b246ad387a8ddee575a00d4512b6d2519a775c Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Wed, 7 Sep 2016 15:46:57 +0200 Subject: [PATCH 2/7] [PL] Add an explicit break in a switch-case. --- ProcessLib/Deformation/LinearBMatrix.h | 1 + 1 file changed, 1 insertion(+) diff --git a/ProcessLib/Deformation/LinearBMatrix.h b/ProcessLib/Deformation/LinearBMatrix.h index 093c222a1ad..67a62e072fe 100644 --- a/ProcessLib/Deformation/LinearBMatrix.h +++ b/ProcessLib/Deformation/LinearBMatrix.h @@ -47,6 +47,7 @@ void computeBMatrix(DNDX_Type const& dNdx, BMatrixType& b_matrix) b_matrix(3, NPOINTS + i) = dNdx(0, i) / std::sqrt(2); b_matrix(0, i) = dNdx(0, i); } + break; default: break; } From b0413d8390de19ff5d8bbcd08fd56eec2bf4aa43 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Mon, 4 Jul 2016 16:22:58 +0200 Subject: [PATCH 3/7] [Mat] Add MechanicsBase for solid models. --- MaterialLib/SolidModels/MechanicsBase.h | 97 +++++++++++++++++++++++++ 1 file changed, 97 insertions(+) create mode 100644 MaterialLib/SolidModels/MechanicsBase.h diff --git a/MaterialLib/SolidModels/MechanicsBase.h b/MaterialLib/SolidModels/MechanicsBase.h new file mode 100644 index 00000000000..e4e891099b6 --- /dev/null +++ b/MaterialLib/SolidModels/MechanicsBase.h @@ -0,0 +1,97 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef MATERIALLIB_SOLIDMODELS_MECHANICSBASE_H_ +#define MATERIALLIB_SOLIDMODELS_MECHANICSBASE_H_ + +#include + +#include "ProcessLib/Deformation/BMatrixPolicy.h" +#include "ProcessLib/Parameter/Parameter.h" + +namespace MeshLib +{ +class Element; +} + +namespace Solids +{ +template +struct MechanicsBase +{ + /// The MaterialStateVariables may store material model specific state + /// (other than sigma and eps), which are usually material history + /// dependent. The objects are stored by the user (usually in assembly per + /// integration point) and are created via \ref + /// createMaterialStateVariables(). + struct MaterialStateVariables + { + virtual ~MaterialStateVariables() = default; + virtual void pushBackState() = 0; + }; + + using KelvinVector = ProcessLib::KelvinVectorType; + using KelvinMatrix = ProcessLib::KelvinMatrixType; + + void computeConstitutiveRelation( + double const t, + ProcessLib::SpatialPosition const& x, + double const dt, + Eigen::Matrix const& eps_prev, + Eigen::Matrix const& eps, + Eigen::Matrix const& sigma_prev, + Eigen::Matrix& sigma, + Eigen::Matrix& + C, + MaterialStateVariables& material_state_variables) + { + // TODO Avoid copies of data: + // Using MatrixBase not possible because template functions + // cannot be virtual. Maybe there is a workaround for this. Using + // Map> makes the interface (for the material model + // implementation) unnecessary difficult. + KelvinVector const eps_prev_{eps_prev}; + KelvinVector const eps_{eps}; + KelvinVector const sigma_prev_{sigma_prev}; + KelvinVector sigma_{sigma}; + KelvinMatrix C_{C}; + + computeConstitutiveRelation(t, + x, + dt, + eps_prev_, + eps_, + sigma_prev_, + sigma_, + C_, + material_state_variables); + + sigma = sigma_; + C = C_; + } + + virtual void computeConstitutiveRelation( + double const t, + ProcessLib::SpatialPosition const& x, + double const dt, + KelvinVector const& eps_prev, + KelvinVector const& eps, + KelvinVector const& sigma_prev, + KelvinVector& sigma, + KelvinMatrix& C, + MaterialStateVariables& material_state_variables) = 0; + + virtual std::unique_ptr + createMaterialStateVariables() = 0; + + virtual ~MechanicsBase() = default; +}; +} + +#endif // MATERIALLIB_SOLIDMODELS_MECHANICSBASE_H_ From 5d8c7fb7a2910fe3029f7c11c6b3303a1c90de4e Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Mon, 4 Jul 2016 16:24:06 +0200 Subject: [PATCH 4/7] [Mat] LinearElasticIsotropic model implementation. --- .../t_poissons_ratio.md | 1 + .../t_youngs_modulus.md | 1 + .../CreateLinearElasticIsotropic.h | 49 ++++++++ .../SolidModels/LinearElasticIsotropic.cpp | 19 +++ .../SolidModels/LinearElasticIsotropic.h | 112 ++++++++++++++++++ 5 files changed, 182 insertions(+) create mode 100644 Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/t_poissons_ratio.md create mode 100644 Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/t_youngs_modulus.md create mode 100644 MaterialLib/SolidModels/CreateLinearElasticIsotropic.h create mode 100644 MaterialLib/SolidModels/LinearElasticIsotropic.cpp create mode 100644 MaterialLib/SolidModels/LinearElasticIsotropic.h diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/t_poissons_ratio.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/t_poissons_ratio.md new file mode 100644 index 00000000000..ca4a1dc843a --- /dev/null +++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/t_poissons_ratio.md @@ -0,0 +1 @@ +Poisson's ratio of a linear elastic material. diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/t_youngs_modulus.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/t_youngs_modulus.md new file mode 100644 index 00000000000..1917432bca6 --- /dev/null +++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/t_youngs_modulus.md @@ -0,0 +1 @@ +Young's modulus of a linear elastic material. diff --git a/MaterialLib/SolidModels/CreateLinearElasticIsotropic.h b/MaterialLib/SolidModels/CreateLinearElasticIsotropic.h new file mode 100644 index 00000000000..ef2e76e29a9 --- /dev/null +++ b/MaterialLib/SolidModels/CreateLinearElasticIsotropic.h @@ -0,0 +1,49 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef MATERIALLIB_SOLIDMODELS_CREATELINEARELASTICISOTROPIC_H_ +#define MATERIALLIB_SOLIDMODELS_CREATELINEARELASTICISOTROPIC_H_ + +#include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter +#include "LinearElasticIsotropic.h" + +namespace Solids +{ + +template +std::unique_ptr> +createLinearElasticIsotropic( + std::vector> const& parameters, + BaseLib::ConfigTree const& config) +{ + config.checkConfigParameter("type", "LinearElasticIsotropic"); + DBUG("Create LinearElasticIsotropic material"); + + // Youngs modulus + auto& youngs_modulus = ProcessLib::findParameter( + config, "youngs_modulus", parameters, 1); + + DBUG("Use '%s' as youngs_modulus parameter.", youngs_modulus.name.c_str()); + + // Poissons ratio + auto& poissons_ratio = ProcessLib::findParameter( + config, "poissons_ratio", parameters, 1); + + DBUG("Use '%s' as poissons_ratio parameter.", poissons_ratio.name.c_str()); + + typename LinearElasticIsotropic::MaterialProperties mp{ + youngs_modulus, poissons_ratio}; + + return std::unique_ptr>{ + new LinearElasticIsotropic{mp}}; +} + +} // namespace Solids + +#endif // MATERIALLIB_SOLIDMODELS_CREATELINEARELASTICISOTROPIC_H_ diff --git a/MaterialLib/SolidModels/LinearElasticIsotropic.cpp b/MaterialLib/SolidModels/LinearElasticIsotropic.cpp new file mode 100644 index 00000000000..72b6e11281e --- /dev/null +++ b/MaterialLib/SolidModels/LinearElasticIsotropic.cpp @@ -0,0 +1,19 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include "LinearElasticIsotropic.h" + +namespace Solids +{ + +template class LinearElasticIsotropic<2>; +template class LinearElasticIsotropic<3>; + + +} // namespace Solids diff --git a/MaterialLib/SolidModels/LinearElasticIsotropic.h b/MaterialLib/SolidModels/LinearElasticIsotropic.h new file mode 100644 index 00000000000..bfa4dccc9de --- /dev/null +++ b/MaterialLib/SolidModels/LinearElasticIsotropic.h @@ -0,0 +1,112 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef MATERIALLIB_SOLIDMODELS_LINEARELASTICISOTROPIC_H_ +#define MATERIALLIB_SOLIDMODELS_LINEARELASTICISOTROPIC_H_ + +#include "MechanicsBase.h" + +namespace Solids +{ +template +class LinearElasticIsotropic final : public MechanicsBase +{ +public: + /// Variables specific to the material model + class MaterialProperties + { + using P = ProcessLib::Parameter; + using X = ProcessLib::SpatialPosition; + + public: + MaterialProperties(P const& youngs_modulus, P const& poissons_ratio) + : _youngs_modulus(youngs_modulus), _poissons_ratio(poissons_ratio) + { + } + + /// Lamé's first parameter. + double lambda(double const t, X const& x) const + { + return _youngs_modulus(t, x)[0] * _poissons_ratio(t, x)[0] / + (1 + _poissons_ratio(t, x)[0]) / + (1 - 2 * _poissons_ratio(t, x)[0]); + } + + /// Lamé's second parameter, the shear modulus. + double mu(double const t, X const& x) const + { + return _youngs_modulus(t, x)[0] / + (2 * (1 + _poissons_ratio(t, x)[0])); + } + + private: + P const& _youngs_modulus; + P const& _poissons_ratio; + }; + + struct MaterialStateVariables + : public MechanicsBase::MaterialStateVariables + { + void pushBackState() {} + }; + + std::unique_ptr< + typename MechanicsBase::MaterialStateVariables> + createMaterialStateVariables() override + { + return std::unique_ptr< + typename MechanicsBase::MaterialStateVariables>{ + new MaterialStateVariables}; + } + +public: + static int const KelvinVectorSize = + ProcessLib::KelvinVectorDimensions::value; + using KelvinVector = ProcessLib::KelvinVectorType; + using KelvinMatrix = ProcessLib::KelvinMatrixType; + + explicit LinearElasticIsotropic( + MaterialProperties const& material_properties) + : _mp(material_properties) + { + } + + void computeConstitutiveRelation( + double const t, + ProcessLib::SpatialPosition const& x, + double const /*dt*/, + KelvinVector const& eps_prev, + KelvinVector const& eps, + KelvinVector const& sigma_prev, + KelvinVector& sigma, + KelvinMatrix& C, + typename MechanicsBase::MaterialStateVariables& + /*material_state_variables*/) override + { + C.setZero(); + + C.template topLeftCorner<3, 3>().setConstant(_mp.lambda(t, x)); + C.noalias() += 2 * _mp.mu(t, x) * KelvinMatrix::Identity(); + + sigma.noalias() = sigma_prev + C * (eps - eps_prev); + } + +private: + MaterialProperties _mp; +}; + +} // namespace Solids + +namespace Solids +{ +extern template class LinearElasticIsotropic<2>; +extern template class LinearElasticIsotropic<3>; +} // namespace Solids + +#endif // MATERIALLIB_SOLIDMODELS_LINEARELASTICISOTROPIC_H_ From 92aaf548b28ddee5e1553541624c46d40f6925c1 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Mon, 4 Jul 2016 16:25:55 +0200 Subject: [PATCH 5/7] [PL] SmallDeformation process/fem implementation. --- .../SMALL_DEFORMATION/c_SMALL_DEFORMATION.md | 1 + .../c_LinearElasticIsotropic.md | 1 + .../i_constitutive_relation.md | 0 .../process_variables/i_process_variables.md | 0 .../process_variables/t_process_variable.md | 1 + .../process/SMALL_DEFORMATION/t_dimension.md | 1 + MaterialLib/SolidModels/MechanicsBase.h | 19 +- ProcessLib/CMakeLists.txt | 2 + .../SmallDeformation/CreateLocalAssemblers.h | 109 ++++++ .../CreateSmallDeformationProcess.h | 105 ++++++ .../SmallDeformation/LocalDataInitializer.h | 291 +++++++++++++++ .../SmallDeformation/SmallDeformationFEM.h | 338 ++++++++++++++++++ .../SmallDeformationProcess-fwd.h | 18 + .../SmallDeformationProcess.cpp | 22 ++ .../SmallDeformationProcess.h | 158 ++++++++ .../SmallDeformationProcessData.h | 53 +++ 16 files changed, 1116 insertions(+), 3 deletions(-) create mode 100644 Documentation/ProjectFile/process/SMALL_DEFORMATION/c_SMALL_DEFORMATION.md create mode 100644 Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/c_LinearElasticIsotropic.md create mode 100644 Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/i_constitutive_relation.md create mode 100644 Documentation/ProjectFile/process/SMALL_DEFORMATION/process_variables/i_process_variables.md create mode 100644 Documentation/ProjectFile/process/SMALL_DEFORMATION/process_variables/t_process_variable.md create mode 100644 Documentation/ProjectFile/process/SMALL_DEFORMATION/t_dimension.md create mode 100644 ProcessLib/SmallDeformation/CreateLocalAssemblers.h create mode 100644 ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h create mode 100644 ProcessLib/SmallDeformation/LocalDataInitializer.h create mode 100644 ProcessLib/SmallDeformation/SmallDeformationFEM.h create mode 100644 ProcessLib/SmallDeformation/SmallDeformationProcess-fwd.h create mode 100644 ProcessLib/SmallDeformation/SmallDeformationProcess.cpp create mode 100644 ProcessLib/SmallDeformation/SmallDeformationProcess.h create mode 100644 ProcessLib/SmallDeformation/SmallDeformationProcessData.h diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/c_SMALL_DEFORMATION.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/c_SMALL_DEFORMATION.md new file mode 100644 index 00000000000..7ad85966769 --- /dev/null +++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/c_SMALL_DEFORMATION.md @@ -0,0 +1 @@ +Deformation process with linear kinematics. diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/c_LinearElasticIsotropic.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/c_LinearElasticIsotropic.md new file mode 100644 index 00000000000..73a3d9d4c22 --- /dev/null +++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/LinearElasticIsotropic/c_LinearElasticIsotropic.md @@ -0,0 +1 @@ +Isotropic linear elastic material model. diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/i_constitutive_relation.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/constitutive_relation/i_constitutive_relation.md new file mode 100644 index 00000000000..e69de29bb2d diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/process_variables/i_process_variables.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/process_variables/i_process_variables.md new file mode 100644 index 00000000000..e69de29bb2d diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/process_variables/t_process_variable.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/process_variables/t_process_variable.md new file mode 100644 index 00000000000..d55b23cfca7 --- /dev/null +++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/process_variables/t_process_variable.md @@ -0,0 +1 @@ +The displacement vector field. diff --git a/Documentation/ProjectFile/process/SMALL_DEFORMATION/t_dimension.md b/Documentation/ProjectFile/process/SMALL_DEFORMATION/t_dimension.md new file mode 100644 index 00000000000..fe666ef74ad --- /dev/null +++ b/Documentation/ProjectFile/process/SMALL_DEFORMATION/t_dimension.md @@ -0,0 +1 @@ +Displacement vector dimension. The displacement dimension must be equal to the mesh dimension. diff --git a/MaterialLib/SolidModels/MechanicsBase.h b/MaterialLib/SolidModels/MechanicsBase.h index e4e891099b6..29c39fb72d2 100644 --- a/MaterialLib/SolidModels/MechanicsBase.h +++ b/MaterialLib/SolidModels/MechanicsBase.h @@ -22,6 +22,11 @@ class Element; namespace Solids { +/// Interface for mechanical solid material models. Provides updates of the +/// stress for a given current state and also a tangent at that position. If the +/// implemented material model stores an internal state, the nested +/// MaterialStateVariables class should be used; it's only responsibility is to +/// provide state's push back possibility. template struct MechanicsBase { @@ -36,9 +41,16 @@ struct MechanicsBase virtual void pushBackState() = 0; }; + /// Polymorphic creator for MaterialStateVariables objects specific for a + /// material model. + virtual std::unique_ptr + createMaterialStateVariables() = 0; + using KelvinVector = ProcessLib::KelvinVectorType; using KelvinMatrix = ProcessLib::KelvinMatrixType; + /// Dynamic size Kelvin vector and matrix wrapper for the polymorphic + /// constitutive relation compute function. void computeConstitutiveRelation( double const t, ProcessLib::SpatialPosition const& x, @@ -76,6 +88,10 @@ struct MechanicsBase C = C_; } + /// Computation of the constitutive relation for specific material model. + /// This should be implemented in the derived model. Fixed Kelvin vector and + /// matrix size version; for dynamic size arguments there is an overloaded + /// wrapper function. virtual void computeConstitutiveRelation( double const t, ProcessLib::SpatialPosition const& x, @@ -87,9 +103,6 @@ struct MechanicsBase KelvinMatrix& C, MaterialStateVariables& material_state_variables) = 0; - virtual std::unique_ptr - createMaterialStateVariables() = 0; - virtual ~MechanicsBase() = default; }; } diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt index ebc1893a65a..188c9850630 100644 --- a/ProcessLib/CMakeLists.txt +++ b/ProcessLib/CMakeLists.txt @@ -12,6 +12,8 @@ APPEND_SOURCE_FILES(SOURCES Parameter) add_subdirectory(GroundwaterFlow) APPEND_SOURCE_FILES(SOURCES GroundwaterFlow) +APPEND_SOURCE_FILES(SOURCES SmallDeformation) + add_subdirectory(TES) APPEND_SOURCE_FILES(SOURCES TES) diff --git a/ProcessLib/SmallDeformation/CreateLocalAssemblers.h b/ProcessLib/SmallDeformation/CreateLocalAssemblers.h new file mode 100644 index 00000000000..aaa61b6933b --- /dev/null +++ b/ProcessLib/SmallDeformation/CreateLocalAssemblers.h @@ -0,0 +1,109 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + */ + +#ifndef PROCESSLIB_SMALLDEFORMATION_CREATE_LOCAL_ASSEMBLERS_H_ +#define PROCESSLIB_SMALLDEFORMATION_CREATE_LOCAL_ASSEMBLERS_H_ + +#include + +#include + +#include "NumLib/DOF/LocalToGlobalIndexMap.h" + +#include "LocalDataInitializer.h" + +namespace ProcessLib +{ +namespace SmallDeformation +{ +namespace detail +{ +template + class LocalAssemblerImplementation, + typename LocalAssemblerInterface, typename... ExtraCtorArgs> +void createLocalAssemblers( + NumLib::LocalToGlobalIndexMap const& dof_table, + std::vector const& mesh_elements, + unsigned const integration_order, + std::vector>& local_assemblers, + ExtraCtorArgs&&... extra_ctor_args) +{ + // Shape matrices initializer + using LocalDataInitializer = + LocalDataInitializer; + + DBUG("Create local assemblers."); + // Populate the vector of local assemblers. + local_assemblers.resize(mesh_elements.size()); + + LocalDataInitializer initializer(dof_table); + + DBUG("Calling local assembler builder for all mesh elements."); + GlobalExecutor::transformDereferenced( + initializer, + mesh_elements, + local_assemblers, + integration_order, + std::forward(extra_ctor_args)...); +} + +} // namespace detail + +/*! Creates local assemblers for each element of the given \c mesh. + * + * \tparam LocalAssemblerImplementation the individual local assembler type + * \tparam LocalAssemblerInterface the general local assembler interface + * \tparam ExtraCtorArgs types of additional constructor arguments. + * Those arguments will be passed to the constructor of + * \c LocalAssemblerImplementation. + * + * The first two template parameters cannot be deduced from the arguments. + * Therefore they always have to be provided manually. + */ +template + class LocalAssemblerImplementation, + typename LocalAssemblerInterface, typename... ExtraCtorArgs> +void createLocalAssemblers( + const unsigned dimension, + std::vector const& mesh_elements, + NumLib::LocalToGlobalIndexMap const& dof_table, + unsigned const integration_order, + std::vector>& local_assemblers, + ExtraCtorArgs&&... extra_ctor_args) +{ + DBUG("Create local assemblers."); + + switch (dimension) + { + case 2: + detail::createLocalAssemblers<2, DisplacementDim, + LocalAssemblerImplementation>( + dof_table, mesh_elements, integration_order, local_assemblers, + std::forward(extra_ctor_args)...); + break; + case 3: + detail::createLocalAssemblers<3, DisplacementDim, + LocalAssemblerImplementation>( + dof_table, mesh_elements, integration_order, local_assemblers, + std::forward(extra_ctor_args)...); + break; + default: + OGS_FATAL( + "Meshes with dimension different than two and three are not " + "supported."); + } +} +} // SmallDeformation + +} // ProcessLib + +#endif // PROCESSLIB_SMALLDEFORMATION_CREATE_LOCAL_ASSEMBLERS_H_ diff --git a/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h b/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h new file mode 100644 index 00000000000..6fdb4236529 --- /dev/null +++ b/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h @@ -0,0 +1,105 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef PROCESS_LIB_CREATESMALLDEFORMATIONPROCESS_H_ +#define PROCESS_LIB_CREATESMALLDEFORMATIONPROCESS_H_ + +#include + +#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h" +#include "ProcessLib/Utils/ParseSecondaryVariables.h" + +#include "SmallDeformationProcess.h" +#include "SmallDeformationProcessData.h" + +namespace ProcessLib +{ +namespace SmallDeformation +{ +template +class SmallDeformationProcess; + +extern template class SmallDeformationProcess<2>; +extern template class SmallDeformationProcess<3>; + +template +std::unique_ptr> +createSmallDeformationProcess( + MeshLib::Mesh& mesh, + std::vector const& variables, + std::vector> const& parameters, + BaseLib::ConfigTree const& config) +{ + //! \ogs_file_param{process__type} + config.checkConfigParameter("type", "SMALL_DEFORMATION"); + DBUG("Create SmallDeformationProcess."); + + // Process variable. + auto process_variables = findProcessVariables( + variables, config, + {//! \ogs_file_param_special{process__SMALL_DEFORMATION__process_variables__process_variable} + "process_variable"}); + + DBUG("Associate displacement with process variable \'%s\'.", + process_variables.back().get().getName().c_str()); + + if (process_variables.back().get().getNumberOfComponents() != + DisplacementDim) + { + OGS_FATAL( + "Number of components of the process variable '%s' is different " + "from the displacement dimension: got %d, expected %d", + process_variables.back().get().getName().c_str(), + process_variables.back().get().getNumberOfComponents(), + DisplacementDim); + } + + // Constitutive relation. + // read type; + auto const constitutive_relation_config = + //! \ogs_file_param{process__SMALL_DEFORMATION__constitutive_relation} + config.getConfigSubtree("constitutive_relation"); + + auto const type = + constitutive_relation_config.peekConfigParameter("type"); + + std::unique_ptr> material = nullptr; + if (type == "LinearElasticIsotropic") + { + material = Solids::createLinearElasticIsotropic( + parameters, constitutive_relation_config); + } + else + { + OGS_FATAL( + "Cannot construct constitutive relation of given type \'%s\'.", + type.c_str()); + } + + SmallDeformationProcessData process_data{ + std::move(material)}; + + SecondaryVariableCollection secondary_variables; + + NumLib::NamedFunctionCaller named_function_caller( + {"SmallDeformation_displacement"}); + + ProcessLib::parseSecondaryVariables(config, secondary_variables, + named_function_caller); + + return std::unique_ptr>{ + new SmallDeformationProcess{ + mesh, parameters, std::move(process_variables), + std::move(process_data), std::move(secondary_variables), + std::move(named_function_caller)}}; +} +} // namespace SmallDeformation +} // namespace ProcessLib + +#endif // PROCESS_LIB_CREATESMALLDEFORMATIONPROCESS_H_ diff --git a/ProcessLib/SmallDeformation/LocalDataInitializer.h b/ProcessLib/SmallDeformation/LocalDataInitializer.h new file mode 100644 index 00000000000..59363f8e5c3 --- /dev/null +++ b/ProcessLib/SmallDeformation/LocalDataInitializer.h @@ -0,0 +1,291 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef PROCESSLIB_SMALLDEFORMATION_LOCALDATAINITIALIZER_H_ +#define PROCESSLIB_SMALLDEFORMATION_LOCALDATAINITIALIZER_H_ + +#include +#include +#include +#include +#include + +#include "MeshLib/Elements/Elements.h" +#include "NumLib/DOF/LocalToGlobalIndexMap.h" +#include "NumLib/Fem/Integration/GaussIntegrationPolicy.h" + + +#ifndef OGS_MAX_ELEMENT_DIM +static_assert(false, "The macro OGS_MAX_ELEMENT_DIM is undefined."); +#endif + +#ifndef OGS_MAX_ELEMENT_ORDER +static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined."); +#endif + +// The following macros decide which element types will be compiled, i.e. +// which element types will be available for use in simulations. + +#ifdef OGS_ENABLE_ELEMENT_SIMPLEX +#define ENABLED_ELEMENT_TYPE_SIMPLEX 1u +#else +#define ENABLED_ELEMENT_TYPE_SIMPLEX 0u +#endif + +#ifdef OGS_ENABLE_ELEMENT_CUBOID +#define ENABLED_ELEMENT_TYPE_CUBOID 1u << 1 +#else +#define ENABLED_ELEMENT_TYPE_CUBOID 0u +#endif + +#ifdef OGS_ENABLE_ELEMENT_PRISM +#define ENABLED_ELEMENT_TYPE_PRISM 1u << 2 +#else +#define ENABLED_ELEMENT_TYPE_PRISM 0u +#endif + +#ifdef OGS_ENABLE_ELEMENT_PYRAMID +#define ENABLED_ELEMENT_TYPE_PYRAMID 1u << 3 +#else +#define ENABLED_ELEMENT_TYPE_PYRAMID 0u +#endif + +// Dependent element types. +// Faces of tets, pyramids and prisms are triangles +#define ENABLED_ELEMENT_TYPE_TRI ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_PYRAMID) \ + |(ENABLED_ELEMENT_TYPE_PRISM)) +// Faces of hexes, pyramids and prisms are quads +#define ENABLED_ELEMENT_TYPE_QUAD ((ENABLED_ELEMENT_TYPE_CUBOID) | (ENABLED_ELEMENT_TYPE_PYRAMID) \ + |(ENABLED_ELEMENT_TYPE_PRISM)) + +// All enabled element types +#define OGS_ENABLED_ELEMENTS \ + ( (ENABLED_ELEMENT_TYPE_SIMPLEX) \ + | (ENABLED_ELEMENT_TYPE_CUBOID ) \ + | (ENABLED_ELEMENT_TYPE_PYRAMID) \ + | (ENABLED_ELEMENT_TYPE_PRISM ) ) + + +// Include only what is needed (Well, the conditions are not sharp). +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 +#include "NumLib/Fem/ShapeFunction/ShapeTet10.h" +#include "NumLib/Fem/ShapeFunction/ShapeTet4.h" +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 +#include "NumLib/Fem/ShapeFunction/ShapeTri3.h" +#include "NumLib/Fem/ShapeFunction/ShapeTri6.h" +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 +#include "NumLib/Fem/ShapeFunction/ShapeHex20.h" +#include "NumLib/Fem/ShapeFunction/ShapeHex8.h" +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 +#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h" +#include "NumLib/Fem/ShapeFunction/ShapeQuad8.h" +#include "NumLib/Fem/ShapeFunction/ShapeQuad9.h" +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 +#include "NumLib/Fem/ShapeFunction/ShapePrism15.h" +#include "NumLib/Fem/ShapeFunction/ShapePrism6.h" +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 +#include "NumLib/Fem/ShapeFunction/ShapePyra13.h" +#include "NumLib/Fem/ShapeFunction/ShapePyra5.h" +#endif + +namespace ProcessLib +{ + +/// The LocalDataInitializer is a functor creating a local assembler data with +/// corresponding to the mesh element type shape functions and calling +/// initialization of the new local assembler data. +/// For example for MeshLib::Quad a local assembler data with template argument +/// NumLib::ShapeQuad4 is created. +template < + typename LocalAssemblerInterface, + template class LocalAssemblerData, + unsigned GlobalDim, int DisplacementDim, typename... ConstructorArgs> +class LocalDataInitializer final +{ +public: + using LADataIntfPtr = std::unique_ptr; + + explicit LocalDataInitializer( + NumLib::LocalToGlobalIndexMap const& dof_table) + : _dof_table(dof_table) + { + // /// Quads and Hexahedra /////////////////////////////////// + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 \ + && OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 1 + _builder[std::type_index(typeid(MeshLib::Quad))] = + makeLocalAssemblerBuilder(); +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 \ + && OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1 + _builder[std::type_index(typeid(MeshLib::Hex))] = + makeLocalAssemblerBuilder(); +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 \ + && OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2 + _builder[std::type_index(typeid(MeshLib::Quad8))] = + makeLocalAssemblerBuilder(); + _builder[std::type_index(typeid(MeshLib::Quad9))] = + makeLocalAssemblerBuilder(); +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 \ + && OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2 + _builder[std::type_index(typeid(MeshLib::Hex20))] = + makeLocalAssemblerBuilder(); +#endif + + + // /// Simplices //////////////////////////////////////////////// + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 \ + && OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 1 + _builder[std::type_index(typeid(MeshLib::Tri))] = + makeLocalAssemblerBuilder(); +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 \ + && OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1 + _builder[std::type_index(typeid(MeshLib::Tet))] = + makeLocalAssemblerBuilder(); +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 \ + && OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2 + _builder[std::type_index(typeid(MeshLib::Tri6))] = + makeLocalAssemblerBuilder(); +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 \ + && OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2 + _builder[std::type_index(typeid(MeshLib::Tet10))] = + makeLocalAssemblerBuilder(); +#endif + + + // /// Prisms //////////////////////////////////////////////////// + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 \ + && OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1 + _builder[std::type_index(typeid(MeshLib::Prism))] = + makeLocalAssemblerBuilder(); +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 \ + && OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2 + _builder[std::type_index(typeid(MeshLib::Prism15))] = + makeLocalAssemblerBuilder(); +#endif + + // /// Pyramids ////////////////////////////////////////////////// + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 \ + && OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1 + _builder[std::type_index(typeid(MeshLib::Pyramid))] = + makeLocalAssemblerBuilder(); +#endif + +#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 \ + && OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2 + _builder[std::type_index(typeid(MeshLib::Pyramid13))] = + makeLocalAssemblerBuilder(); +#endif + } + + /// Sets the provided \c data_ptr to the newly created local assembler data. + /// + /// \attention + /// The index \c id is not necessarily the mesh item's id. Especially when + /// having multiple meshes it will differ from the latter. + void operator()( + std::size_t const id, + MeshLib::Element const& mesh_item, + LADataIntfPtr& data_ptr, + unsigned const integration_order, + ConstructorArgs&&... args) const + { + auto const type_idx = std::type_index(typeid(mesh_item)); + auto const it = _builder.find(type_idx); + + if (it != _builder.end()) { + auto const num_local_dof = _dof_table.getNumberOfElementDOF(id); + data_ptr = it->second( + mesh_item, num_local_dof, integration_order, + std::forward(args)...); + } else { + OGS_FATAL("You are trying to build a local assembler for an unknown mesh element type (%s)." + " Maybe you have disabled this mesh element type in your build configuration.", + type_idx.name()); + } + } + +private: + using LADataBuilder = std::function; + + template + using IntegrationMethod = typename NumLib::GaussIntegrationPolicy< + typename ShapeFunction::MeshElement>::IntegrationMethod; + + template + using LAData = + LocalAssemblerData, + GlobalDim, DisplacementDim>; + + /// Generates a function that creates a new LocalAssembler of type LAData + template + static LADataBuilder makeLocalAssemblerBuilder() + { + return [](MeshLib::Element const& e, + std::size_t const local_matrix_size, + unsigned const integration_order, + ConstructorArgs&&... args) + { + return LADataIntfPtr{ + new LAData{ + e, local_matrix_size, integration_order, + std::forward(args)... + }}; + }; + } + + /// Mapping of element types to local assembler constructors. + std::unordered_map _builder; + + NumLib::LocalToGlobalIndexMap const& _dof_table; +}; + +} // namespace ProcessLib + + +#undef ENABLED_ELEMENT_TYPE_SIMPLEX +#undef ENABLED_ELEMENT_TYPE_CUBOID +#undef ENABLED_ELEMENT_TYPE_PYRAMID +#undef ENABLED_ELEMENT_TYPE_PRISM +#undef ENABLED_ELEMENT_TYPE_TRI +#undef ENABLED_ELEMENT_TYPE_QUAD +#undef OGS_ENABLED_ELEMENTS + +#endif // PROCESSLIB_SMALLDEFORMATION_LOCALDATAINITIALIZER_H_ diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h new file mode 100644 index 00000000000..2c43afe7b3e --- /dev/null +++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h @@ -0,0 +1,338 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef PROCESS_LIB_SMALLDEFORMATION_FEM_H_ +#define PROCESS_LIB_SMALLDEFORMATION_FEM_H_ + +#include +#include + +#include "MaterialLib/SolidModels/LinearElasticIsotropic.h" +#include "NumLib/Extrapolation/ExtrapolatableElement.h" +#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h" +#include "NumLib/Fem/ShapeMatrixPolicy.h" +#include "ProcessLib/Deformation/BMatrixPolicy.h" +#include "ProcessLib/Deformation/LinearBMatrix.h" +#include "ProcessLib/LocalAssemblerInterface.h" +#include "ProcessLib/LocalAssemblerTraits.h" +#include "ProcessLib/Parameter/Parameter.h" +#include "ProcessLib/Utils/InitShapeMatrices.h" + +#include "SmallDeformationProcessData.h" + +template +struct IntegrationPointData final +{ + explicit IntegrationPointData( + Solids::MechanicsBase& solid_material) + : _solid_material(solid_material), + _material_state_variables( + _solid_material.createMaterialStateVariables()) + { + } + +#if defined(_MSC_VER) && _MSC_VER < 1900 + // The default generated move-ctor is correctly generated for other + // compilers. + explicit IntegrationPointData(IntegrationPointData&& other) + : _b_matrices(std::move(other._b_matrices)), + _sigma(std::move(other._sigma)), + _sigma_prev(std::move(other._sigma_prev)), + _eps(std::move(other._eps)), + _eps_prev(std::move(other._eps_prev)), + _solid_material(other._solid_material), + _material_state_variables(std::move(other._material_state_variables)), + _C(std::move(other._C)), + _detJ(std::move(other._detJ)) + { + } +#endif // _MSC_VER + + typename BMatricesType::BMatrixType _b_matrices; + typename BMatricesType::KelvinVectorType _sigma, _sigma_prev; + typename BMatricesType::KelvinVectorType _eps, _eps_prev; + + Solids::MechanicsBase& _solid_material; + std::unique_ptr< + typename Solids::MechanicsBase::MaterialStateVariables> + _material_state_variables; + + typename BMatricesType::KelvinMatrixType _C; + double _detJ; + + void pushBackState() + { + _eps_prev = _eps; + _sigma_prev = _sigma; + _material_state_variables->pushBackState(); + } +}; + +namespace ProcessLib +{ +namespace SmallDeformation +{ + +/// Used by for extrapolation of the integration point values. It is ordered +/// (and stored) by integration points. +template +struct SecondaryData +{ + std::vector N; + std::vector _sigmaXX; + std::vector _sigmaYY; + std::vector _sigmaXY; +}; + +struct SmallDeformationLocalAssemblerInterface + : public ProcessLib::LocalAssemblerInterface, + public NumLib::ExtrapolatableElement +{ + virtual std::vector const& getIntPtSigmaXX( + std::vector& /*cache*/) const = 0; + + virtual std::vector const& getIntPtSigmaYY( + std::vector& /*cache*/) const = 0; + + virtual std::vector const& getIntPtSigmaXY( + std::vector& /*cache*/) const = 0; +}; + +template +class SmallDeformationLocalAssembler + : public SmallDeformationLocalAssemblerInterface +{ +public: + using ShapeMatricesType = + ShapeMatrixPolicyType; + using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType; + using NodalVectorType = typename ShapeMatricesType::NodalVectorType; + using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices; + using BMatricesType = BMatrixPolicyType; + + using BMatrixType = typename BMatricesType::BMatrixType; + using StiffnessMatrixType = typename BMatricesType::StiffnessMatrixType; + using NodalForceVectorType = typename BMatricesType::NodalForceVectorType; + using NodalDisplacementVectorType = + typename BMatricesType::NodalForceVectorType; + + SmallDeformationLocalAssembler(SmallDeformationLocalAssembler const&) = + delete; + SmallDeformationLocalAssembler(SmallDeformationLocalAssembler&&) = delete; + + SmallDeformationLocalAssembler( + MeshLib::Element const& e, + std::size_t const local_matrix_size, + unsigned const integration_order, + SmallDeformationProcessData& process_data) + : _process_data(process_data), + _localA(local_matrix_size, local_matrix_size), + _localRhs(local_matrix_size), + _integration_method(integration_order), + _element(e) + { + unsigned const n_integration_points = + _integration_method.getNumberOfPoints(); + + _ip_data.reserve(n_integration_points); + + _secondary_data.N.resize(n_integration_points); + _secondary_data._sigmaXX.resize(n_integration_points); + _secondary_data._sigmaYY.resize(n_integration_points); + _secondary_data._sigmaXY.resize(n_integration_points); + + auto const shape_matrices = + initShapeMatrices( + e, _integration_method); + + for (unsigned ip = 0; ip < n_integration_points; ip++) + { + _ip_data.emplace_back(*_process_data._material); + auto& ip_data = _ip_data[ip]; + ip_data._detJ = shape_matrices[ip].detJ; + ip_data._b_matrices.resize( + KelvinVectorDimensions::value, + ShapeFunction::NPOINTS * DisplacementDim); + LinearBMatrix::computeBMatrix< + DisplacementDim, ShapeFunction::NPOINTS, + typename ShapeMatricesType::GlobalDimNodalMatrixType, + BMatrixType>(shape_matrices[ip].dNdx, ip_data._b_matrices); + + ip_data._sigma.resize( + KelvinVectorDimensions::value); + ip_data._sigma_prev.resize( + KelvinVectorDimensions::value); + ip_data._eps.resize(KelvinVectorDimensions::value); + ip_data._eps_prev.resize( + KelvinVectorDimensions::value); + ip_data._C.resize(KelvinVectorDimensions::value, + KelvinVectorDimensions::value); + + _secondary_data.N[ip] = shape_matrices[ip].N; + } + } + + void assembleConcrete( + double const t, std::vector const& local_x, + NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices, + GlobalMatrix& /*M*/, GlobalMatrix& /*K*/, GlobalVector& b) override + { + _localRhs.setZero(); + + unsigned const n_integration_points = + _integration_method.getNumberOfPoints(); + + NodalDisplacementVectorType local_displacement(ShapeFunction::NPOINTS * + DisplacementDim); + + SpatialPosition x_position; + x_position.setElementID(_element.getID()); + + for (unsigned ip = 0; ip < n_integration_points; ip++) + { + x_position.setIntegrationPoint(ip); + auto const& wp = _integration_method.getWeightedPoint(ip); + auto const& detJ = _ip_data[ip]._detJ; + + auto const& B = _ip_data[ip]._b_matrices; + auto const& eps_prev = _ip_data[ip]._eps_prev; + auto const& sigma_prev = _ip_data[ip]._sigma_prev; + + auto& eps = _ip_data[ip]._eps; + auto& sigma = _ip_data[ip]._sigma; + auto& C = _ip_data[ip]._C; + auto& material_state_variables = + *_ip_data[ip]._material_state_variables; + + eps.noalias() = + B * + Eigen::Map( + local_x.data(), ShapeFunction::NPOINTS * DisplacementDim); + + _ip_data[ip]._solid_material.computeConstitutiveRelation( + t, x_position, _process_data.dt, eps_prev, eps, sigma_prev, + sigma, C, material_state_variables); + + _localRhs.noalias() -= + B.transpose() * sigma * detJ * wp.getWeight(); + + // TODO: Reuse _ip_data[ip]._sigma + _secondary_data._sigmaXX[ip] = sigma[0]; + _secondary_data._sigmaYY[ip] = sigma[1]; + _secondary_data._sigmaXY[ip] = sigma[3]; + } + + b.add(indices.rows, _localRhs); + } + + void assembleJacobianConcrete( + double const /*t*/, + std::vector const& /*local_x*/, + NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices, + GlobalMatrix& Jac) override + { + _localA.setZero(); + + unsigned const n_integration_points = + _integration_method.getNumberOfPoints(); + + for (unsigned ip = 0; ip < n_integration_points; ip++) + { + auto const& B = _ip_data[ip]._b_matrices; + auto const& wp = _integration_method.getWeightedPoint(ip); + auto& C = _ip_data[ip]._C; + auto const& detJ = _ip_data[ip]._detJ; + + _localA.noalias() += B.transpose() * C * B * detJ * wp.getWeight(); + } + + Jac.add(indices, _localA); + } + + void preTimestepConcrete(std::vector const& /*local_x*/, + double const /*t*/, + double const /*delta_t*/) override + { + unsigned const n_integration_points = + _integration_method.getNumberOfPoints(); + + for (unsigned ip = 0; ip < n_integration_points; ip++) + { + _ip_data[ip].pushBackState(); + } + } + + Eigen::Map getShapeMatrix( + const unsigned integration_point) const override + { + auto const& N = _secondary_data.N[integration_point]; + + // assumes N is stored contiguously in memory + return Eigen::Map(N.data(), N.size()); + } + + std::vector const& getIntPtSigmaXX( + std::vector& /*cache*/) const override + { + return _secondary_data._sigmaXX; + } + + std::vector const& getIntPtSigmaYY( + std::vector& /*cache*/) const override + { + return _secondary_data._sigmaYY; + } + + std::vector const& getIntPtSigmaXY( + std::vector& /*cache*/) const override + { + return _secondary_data._sigmaXY; + } + +private: + SmallDeformationProcessData& _process_data; + + std::vector> _ip_data; + + StiffnessMatrixType _localA; + NodalForceVectorType _localRhs; + + IntegrationMethod _integration_method; + MeshLib::Element const& _element; + SecondaryData _secondary_data; +}; + +template +class LocalAssemblerData final + : public SmallDeformationLocalAssembler +{ +public: + LocalAssemblerData(LocalAssemblerData const&) = delete; + LocalAssemblerData(LocalAssemblerData&&) = delete; + + LocalAssemblerData( + MeshLib::Element const& e, + std::size_t const local_matrix_size, + unsigned const integration_order, + SmallDeformationProcessData& process_data) + : SmallDeformationLocalAssembler( + e, local_matrix_size, integration_order, process_data) + { + } +}; + +} // namespace SmallDeformation +} // namespace ProcessLib + +#endif // PROCESS_LIB_SMALLDEFORMATION_FEM_H_ diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess-fwd.h b/ProcessLib/SmallDeformation/SmallDeformationProcess-fwd.h new file mode 100644 index 00000000000..0cb7eaabdc8 --- /dev/null +++ b/ProcessLib/SmallDeformation/SmallDeformationProcess-fwd.h @@ -0,0 +1,18 @@ +/** + * \copyright + * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef PROCESS_LIB_SMALLDEFORMATIONPROCESS_FWD_H_ +#define PROCESS_LIB_SMALLDEFORMATIONPROCESS_FWD_H_ + +#include "SmallDeformationProcess.h" + +extern template class ProcessLib::SmallDeformation::SmallDeformationProcess<2>; +extern template class ProcessLib::SmallDeformation::SmallDeformationProcess<3>; + +#endif // PROCESS_LIB_SMALLDEFORMATIONPROCESS_FWD_H_ diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp new file mode 100644 index 00000000000..534c17f5306 --- /dev/null +++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp @@ -0,0 +1,22 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include "SmallDeformationProcess-fwd.h" +#include "SmallDeformationProcess.h" + +namespace ProcessLib +{ +namespace SmallDeformation +{ + +template class SmallDeformationProcess<2>; +template class SmallDeformationProcess<3>; + +} // namespace SmallDeformation +} // namespace ProcessLib diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/SmallDeformation/SmallDeformationProcess.h new file mode 100644 index 00000000000..bd1e113d370 --- /dev/null +++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.h @@ -0,0 +1,158 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef PROCESS_LIB_SMALLDEFORMATIONPROCESS_H_ +#define PROCESS_LIB_SMALLDEFORMATIONPROCESS_H_ + +#include + +#include "ProcessLib/Process.h" +#include "ProcessLib/SmallDeformation/CreateLocalAssemblers.h" + +#include "SmallDeformationFEM.h" +#include "SmallDeformationProcessData.h" + +namespace ProcessLib +{ +namespace SmallDeformation +{ +template +class SmallDeformationProcess final : public Process +{ + using Base = Process; + +public: + SmallDeformationProcess( + MeshLib::Mesh& mesh, + std::vector> const& parameters, + std::vector>&& + process_variables, + SmallDeformationProcessData&& process_data, + SecondaryVariableCollection&& secondary_variables, + NumLib::NamedFunctionCaller&& named_function_caller) + : Process(mesh, parameters, std::move(process_variables), + std::move(secondary_variables), + std::move(named_function_caller)), + _process_data(std::move(process_data)) + { + } + + //! \name ODESystem interface + //! @{ + + bool isLinear() const override { return false; } + //! @} + +private: + using LocalAssemblerInterface = SmallDeformationLocalAssemblerInterface; + + void initializeConcreteProcess( + NumLib::LocalToGlobalIndexMap const& dof_table, + MeshLib::Mesh const& mesh, + unsigned const integration_order) override + { + ProcessLib::SmallDeformation::createLocalAssemblers( + mesh.getDimension(), mesh.getElements(), dof_table, + integration_order, _local_assemblers, _process_data); + + // TODO move the two data members somewhere else. + // for extrapolation of secondary variables + std::vector> + all_mesh_subsets_single_component; + all_mesh_subsets_single_component.emplace_back( + new MeshLib::MeshSubsets(_mesh_subset_all_nodes.get())); + _local_to_global_index_map_single_component.reset( + new NumLib::LocalToGlobalIndexMap( + std::move(all_mesh_subsets_single_component), + // by location order is needed for output + NumLib::ComponentOrder::BY_LOCATION)); + + Base::_secondary_variables.addSecondaryVariable( + "sigma_xx", 1, + makeExtrapolator( + getExtrapolator(), _local_assemblers, + &SmallDeformationLocalAssemblerInterface::getIntPtSigmaXX)); + + Base::_secondary_variables.addSecondaryVariable( + "sigma_yy", 1, + makeExtrapolator( + getExtrapolator(), _local_assemblers, + &SmallDeformationLocalAssemblerInterface::getIntPtSigmaYY)); + + Base::_secondary_variables.addSecondaryVariable( + "sigma_xy", 1, + makeExtrapolator( + getExtrapolator(), _local_assemblers, + &SmallDeformationLocalAssemblerInterface::getIntPtSigmaXY)); + } + + void assembleConcreteProcess(const double t, GlobalVector const& x, + GlobalMatrix& M, GlobalMatrix& K, + GlobalVector& b) override + { + DBUG("Assemble SmallDeformationProcess."); + + // Call global assembler for each local assembly item. + GlobalExecutor::executeMemberOnDereferenced( + &SmallDeformationLocalAssemblerInterface::assemble, + _local_assemblers, *_local_to_global_index_map, t, x, M, K, b); + } + + void assembleJacobianConcreteProcess(const double t, GlobalVector const& x, + GlobalVector const& /*xdot*/, + const double /*dxdot_dx*/, + GlobalMatrix const& /*M*/, + const double /*dx_dx*/, + GlobalMatrix const& /*K*/, + GlobalMatrix& Jac) override + { + DBUG("AssembleJacobian SmallDeformationProcess."); + + // Call global assembler for each local assembly item. + GlobalExecutor::executeMemberOnDereferenced( + &SmallDeformationLocalAssemblerInterface::assembleJacobian, + _local_assemblers, *_local_to_global_index_map, t, x, Jac); + } + + void preTimestep(GlobalVector const& x, double const t, + double const dt) override + { + DBUG("PreTimestep SmallDeformationProcess."); + + _process_data.dt = dt; + _process_data.t = t; + + GlobalExecutor::executeMemberOnDereferenced( + &SmallDeformationLocalAssemblerInterface::preTimestep, + _local_assemblers, *_local_to_global_index_map, x, t, dt); + } + + void postTimestep(GlobalVector const& x) override + { + DBUG("PostTimestep SmallDeformationProcess."); + + GlobalExecutor::executeMemberOnDereferenced( + &SmallDeformationLocalAssemblerInterface::postTimestep, + _local_assemblers, *_local_to_global_index_map, x); + } + +private: + SmallDeformationProcessData _process_data; + + std::vector> _local_assemblers; + + std::unique_ptr + _local_to_global_index_map_single_component; +}; + +} // namespace SmallDeformation +} // namespace ProcessLib + +#endif // PROCESS_LIB_SMALLDEFORMATIONPROCESS_H_ diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcessData.h b/ProcessLib/SmallDeformation/SmallDeformationProcessData.h new file mode 100644 index 00000000000..45c11295210 --- /dev/null +++ b/ProcessLib/SmallDeformation/SmallDeformationProcessData.h @@ -0,0 +1,53 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef PROCESSLIB_SMALLDEFORMATION_SMALLDEFORMATIONPROCESSDATA_H_ +#define PROCESSLIB_SMALLDEFORMATION_SMALLDEFORMATIONPROCESSDATA_H_ + +namespace MeshLib +{ +class Element; +} + +namespace ProcessLib +{ +namespace SmallDeformation +{ +template +struct SmallDeformationProcessData +{ + SmallDeformationProcessData( + std::unique_ptr>&& material) + : _material{std::move(material)} + { + } + + SmallDeformationProcessData(SmallDeformationProcessData&& other) + : _material{std::move(other._material)} + { + } + + //! Copies are forbidden. + SmallDeformationProcessData(SmallDeformationProcessData const&) = delete; + + //! Assignments are not needed. + void operator=(SmallDeformationProcessData const&) = delete; + + //! Assignments are not needed. + void operator=(SmallDeformationProcessData&&) = delete; + + std::unique_ptr> _material; + double dt; + double t; +}; + +} // namespace SmallDeformation +} // namespace ProcessLib + +#endif // PROCESSLIB_SMALLDEFORMATION_SMALLDEFORMATIONPROCESSDATA_H_ From 0ba099d3b3163393e380ebbfdeafd2f1243b25b8 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Mon, 4 Jul 2016 16:27:15 +0200 Subject: [PATCH 6/7] [AppL] Instantiate SmallDeformationProcess in prj. --- Applications/ApplicationsLib/ProjectData.cpp | 23 ++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp index 5787f506a2c..c70c0b819d2 100644 --- a/Applications/ApplicationsLib/ProjectData.cpp +++ b/Applications/ApplicationsLib/ProjectData.cpp @@ -32,6 +32,7 @@ #include "ProcessLib/UncoupledProcessesTimeLoop.h" #include "ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h" +#include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h" #include "ProcessLib/TES/CreateTESProcess.h" #include "ProcessLib/HeatConduction/CreateHeatConductionProcess.h" @@ -269,6 +270,28 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config) process = ProcessLib::HeatConduction::createHeatConductionProcess( *_mesh_vec[0], _process_variables, _parameters, process_config); } + else if (type == "SMALL_DEFORMATION") + { + switch (process_config.getConfigParameter("dimension")) + { + case 2: + process = ProcessLib::SmallDeformation:: + createSmallDeformationProcess<2>( + *_mesh_vec[0], _process_variables, _parameters, + process_config); + break; + case 3: + process = ProcessLib::SmallDeformation:: + createSmallDeformationProcess<3>( + *_mesh_vec[0], _process_variables, _parameters, + process_config); + break; + default: + OGS_FATAL( + "SMALL_DEFORMATION process does not support " + "given dimension"); + } + } else { OGS_FATAL("Unknown process type: %s", type.c_str()); From d98f6aaede124134f2089a78a216f4a5959b4c2c Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Tue, 30 Aug 2016 15:59:36 +0200 Subject: [PATCH 7/7] Update Tests/Data. --- Applications/CLI/Tests.cmake | 46 ++++++++++++++++++++++++++++++++++++ Tests/Data | 2 +- 2 files changed, 47 insertions(+), 1 deletion(-) diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake index 5187632ca86..86ad238fc13 100644 --- a/Applications/CLI/Tests.cmake +++ b/Applications/CLI/Tests.cmake @@ -243,6 +243,52 @@ if(NOT OGS_USE_MPI) temperature_analytical.vtu line_60_heat_pcs_0_ts_405_t_31640625.000000.vtu Temperature_Analytical_1year temperature ) + # Mechanics; Small deformations, linear (SDL) + AddTest( + NAME Mechanics_SDL_square_1e0_displacementBC + PATH Mechanics/Linear + EXECUTABLE ogs + EXECUTABLE_ARGS square_1e0.prj + WRAPPER time + TESTER vtkdiff + ABSTOL 1e-16 RELTOL 1e-16 + DIFF_DATA + square_1e0_expected_pcs_0_ts_4_t_1.000000.vtu square_1e0_pcs_0_ts_4_t_1.000000.vtu displacement displacement + ) + AddTest( + NAME Mechanics_SDL_square_1e2_tractionBC + PATH Mechanics/Linear + EXECUTABLE ogs + EXECUTABLE_ARGS square_1e2.prj + WRAPPER time + TESTER vtkdiff + ABSTOL 1e-16 RELTOL 1e-16 + DIFF_DATA + square_1e2_expected_pcs_0_ts_4_t_1.000000.vtu square_1e2_pcs_0_ts_4_t_1.000000.vtu displacement displacement + ) + AddTest( + NAME LARGE_Mechanics_SDL_disc_with_hole + PATH Mechanics/Linear + EXECUTABLE ogs + EXECUTABLE_ARGS disc_with_hole.prj + WRAPPER time + TESTER vtkdiff + ABSTOL 1e-16 RELTOL 1e-16 + DIFF_DATA + disc_with_hole_expected_pcs_0_ts_4_t_1.000000.vtu disc_with_hole_pcs_0_ts_4_t_1.000000.vtu displacement displacement + ) + AddTest( + NAME LARGE_Mechanics_SDL_square_1e5_tractionBC + PATH Mechanics/Linear + EXECUTABLE ogs + EXECUTABLE_ARGS square_1e5.prj + WRAPPER time + TESTER vtkdiff + ABSTOL 1e-16 RELTOL 1e-16 + DIFF_DATA + square_1e5_expected_pcs_0_ts_4_t_1.000000.vtu square_1e5_pcs_0_ts_4_t_1.000000.vtu displacement displacement + ) + else() # MPI groundwater flow tests AddTest( diff --git a/Tests/Data b/Tests/Data index 117219b7e7a..7ac91ea4b42 160000 --- a/Tests/Data +++ b/Tests/Data @@ -1 +1 @@ -Subproject commit 117219b7e7a6ae337257bb50c49730cdc818b26b +Subproject commit 7ac91ea4b4231785326485c42e8ed7720c487e94