From 5230452d1a6765fd4426dd286efe32732bbcb74e Mon Sep 17 00:00:00 2001
From: xymiao <miaoxingyuan@live.cn>
Date: Mon, 8 May 2017 19:16:45 +0200
Subject: [PATCH] added thermo-mechanical process and test.
---
Applications/ApplicationsLib/ProjectData.cpp | 21 +
ProcessLib/CMakeLists.txt | 2 +
.../ThermoMechanics/CreateLocalAssemblers.h | 100 ++++
.../CreateThermoMechanicsProcess.cpp | 196 ++++++
.../CreateThermoMechanicsProcess.h | 28 +
.../ThermoMechanics/LocalDataInitializer.h | 306 ++++++++++
ProcessLib/ThermoMechanics/Tests.cmake | 16 +
.../ThermoMechanics/ThermoMechanicsFEM.h | 558 ++++++++++++++++++
.../ThermoMechanicsProcess-fwd.h | 15 +
.../ThermoMechanicsProcess.cpp | 21 +
.../ThermoMechanics/ThermoMechanicsProcess.h | 220 +++++++
.../ThermoMechanicsProcessData.h | 80 +++
12 files changed, 1563 insertions(+)
create mode 100644 ProcessLib/ThermoMechanics/CreateLocalAssemblers.h
create mode 100644 ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
create mode 100644 ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.h
create mode 100644 ProcessLib/ThermoMechanics/LocalDataInitializer.h
create mode 100644 ProcessLib/ThermoMechanics/Tests.cmake
create mode 100644 ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
create mode 100644 ProcessLib/ThermoMechanics/ThermoMechanicsProcess-fwd.h
create mode 100644 ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
create mode 100644 ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
create mode 100644 ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index ecfbbc76185..765fd9cebb4 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -46,6 +46,7 @@
#include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h"
#include "ProcessLib/TES/CreateTESProcess.h"
#include "ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.h"
+#include "ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.h"
#include "ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h"
#include "ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowWithPrhoProcess.h"
@@ -434,6 +435,26 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
"given dimension");
}
}
+ else if (type == "THERMOMECHANICS")
+ {
+ switch (process_config.getConfigParameter<int>("dimension"))
+ {
+ case 2:
+ process = ProcessLib::ThermoMechanics::
+ createThermoMechanicsProcess<2>(
+ *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config);
+ break;
+ case 3:
+ process = ProcessLib::ThermoMechanics::
+ createThermoMechanicsProcess<3>(
+ *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config);
+ break;
+ }
+ }
else if (type == "RICHARDS_FLOW")
{
process = ProcessLib::RichardsFlow::createRichardsFlowProcess(
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index daafa478bf5..e64cf8f6782 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -20,6 +20,7 @@ APPEND_SOURCE_FILES(SOURCES Parameter)
APPEND_SOURCE_FILES(SOURCES RichardsFlow)
APPEND_SOURCE_FILES(SOURCES SmallDeformation)
APPEND_SOURCE_FILES(SOURCES TES)
+APPEND_SOURCE_FILES(SOURCES ThermoMechanics)
APPEND_SOURCE_FILES(SOURCES TwoPhaseFlowWithPP)
APPEND_SOURCE_FILES(SOURCES TwoPhaseFlowWithPrho)
APPEND_SOURCE_FILES(SOURCES ThermalTwoPhaseFlowWithPP)
@@ -60,6 +61,7 @@ include(LiquidFlow/Tests.cmake)
include(RichardsFlow/Tests.cmake)
include(SmallDeformation/Tests.cmake)
include(TES/Tests.cmake)
+include(ThermoMechanics/Tests.cmake)
include(TwoPhaseFlowWithPP/Tests.cmake)
include(TwoPhaseFlowWithPrho/Tests.cmake)
include(ThermalTwoPhaseFlowWithPP/Tests.cmake)
diff --git a/ProcessLib/ThermoMechanics/CreateLocalAssemblers.h b/ProcessLib/ThermoMechanics/CreateLocalAssemblers.h
new file mode 100644
index 00000000000..7a4103b84f1
--- /dev/null
+++ b/ProcessLib/ThermoMechanics/CreateLocalAssemblers.h
@@ -0,0 +1,100 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include <vector>
+
+#include <logog/include/logog.hpp>
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+
+#include "LocalDataInitializer.h"
+
+namespace ProcessLib
+{
+namespace ThermoMechanics
+{
+namespace detail
+{
+template <unsigned GlobalDim, int DisplacementDim,
+ template <typename, typename, unsigned, int>
+ class LocalAssemblerImplementation,
+ typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<MeshLib::Element*> const& mesh_elements,
+ std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+ ExtraCtorArgs&&... extra_ctor_args)
+{
+ // Shape matrices initializer
+ using LocalDataInitializer =
+ LocalDataInitializer<LocalAssemblerInterface,
+ LocalAssemblerImplementation, GlobalDim,
+ DisplacementDim, ExtraCtorArgs...>;
+
+ 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,
+ std::forward<ExtraCtorArgs>(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 <int DisplacementDim, template <typename, typename, unsigned, int>
+ class LocalAssemblerImplementation,
+ typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+ const unsigned dimension,
+ std::vector<MeshLib::Element*> const& mesh_elements,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+ ExtraCtorArgs&&... extra_ctor_args)
+{
+ DBUG("Create local assemblers.");
+
+ switch (dimension)
+ {
+ case 2:
+ detail::createLocalAssemblers<2, DisplacementDim,
+ LocalAssemblerImplementation>(
+ dof_table, mesh_elements, local_assemblers,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+ break;
+ case 3:
+ detail::createLocalAssemblers<3, DisplacementDim,
+ LocalAssemblerImplementation>(
+ dof_table, mesh_elements, local_assemblers,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+ break;
+ default:
+ OGS_FATAL(
+ "Meshes with dimension different than two and three are not "
+ "supported.");
+ }
+}
+} // ThermoMechanics
+
+} // ProcessLib
diff --git a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
new file mode 100644
index 00000000000..12320a2d4f9
--- /dev/null
+++ b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
@@ -0,0 +1,196 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "CreateThermoMechanicsProcess.h"
+
+#include <cassert>
+
+#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
+
+#include "ThermoMechanicsProcess.h"
+#include "ThermoMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace ThermoMechanics
+{
+template <int DisplacementDim>
+class ThermoMechanicsProcess;
+
+extern template class ThermoMechanicsProcess<2>;
+extern template class ThermoMechanicsProcess<3>;
+
+template <int DisplacementDim>
+std::unique_ptr<Process> createThermoMechanicsProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{process__type}
+ config.checkConfigParameter("type", "THERMOMECHANICS");
+ DBUG("Create ThermoMechanicsProcess.");
+
+ // Process variable.
+
+ //! \ogs_file_param{prj__processes__process__THERMOMECHANICS__process_variables}
+ auto const pv_config = config.getConfigSubtree("process_variables");
+
+ auto process_variables = findProcessVariables(
+ variables, pv_config,
+ {//! \ogs_file_param_special{prj__processes__process__THERMOMECHANICS__process_variables__temperature}
+ "temperature",
+ //! \ogs_file_param_special{prj__processes__process__THERMOMECHANICS__process_variables__displacement}
+ "displacement"});
+
+ DBUG("Associate displacement with process variable \'%s\'.",
+ process_variables[1].get().getName().c_str());
+
+ if (process_variables[1].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[1].get().getName().c_str(),
+ process_variables[1].get().getNumberOfComponents(),
+ DisplacementDim);
+ }
+
+ DBUG("Associate temperature with process variable \'%s\'.",
+ process_variables[0].get().getName().c_str());
+ if (process_variables[0].get().getNumberOfComponents() != 1)
+ {
+ OGS_FATAL(
+ "Temperature process variable '%s' is not a scalar variable but has"
+ "%d components.",
+ process_variables[0].get().getName().c_str(),
+ process_variables[0].get().getNumberOfComponents(),
+ DisplacementDim);
+ }
+
+ // Constitutive relation.
+ // read type;
+ auto const constitutive_relation_config =
+ //! \ogs_file_param{process__THERMOMECHANICS_constitutive_relation}
+ config.getConfigSubtree("constitutive_relation");
+
+ auto const type =
+ //! \ogs_file_param{prj__processes__process__THERMOMECHANICS__constitutive_relation__type}
+ constitutive_relation_config.peekConfigParameter<std::string>("type");
+
+ std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
+ material = nullptr;
+ if (type == "LinearElasticIsotropic")
+ {
+ material =
+ MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(
+ parameters, constitutive_relation_config);
+ }
+ else
+ {
+ OGS_FATAL(
+ "Cannot construct constitutive relation of given type \'%s\'.",
+ type.c_str());
+ }
+
+ // Solid density
+ auto& solid_density = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__THERMOMECHANICS_solid_density}
+ "solid_density", parameters, 1);
+ DBUG("Use \'%s\' as solid density parameter.", solid_density.name.c_str());
+
+ // Linear thermal expansion coefficient
+ auto& linear_thermal_expansion_coefficient = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__THERMOMECHANICS_linear_thermal_expansion_coefficient}
+ "linear_thermal_expansion_coefficient", parameters, 1);
+ DBUG("Use \'%s\' as linear thermal expansion coefficient.",
+ linear_thermal_expansion_coefficient.name.c_str());
+ // Specific heat capacity
+ auto& specific_heat_capacity = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__THERMOMECHANICS_specific_heat_capacity}
+ "specific_heat_capacity", parameters, 1);
+ DBUG("Use \'%s\' as specific heat capacity parameter.",
+ specific_heat_capacity.name.c_str());
+ // Thermal conductivity
+ auto& thermal_conductivity = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__THERMOMECHANICS_thermal_conductivity}
+ "thermal_conductivity", parameters, 1);
+ DBUG("Use \'%s\' as thermal conductivity parameter.",
+ thermal_conductivity.name.c_str());
+ // Reference temperature
+ const double reference_temperature =
+ //! \ogs_file_param_special{process__THERMOMECHANICS_reference_temperature}
+ config.getConfigParameter<double>("reference_temperature");
+
+ // Specific body force
+ Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
+ {
+ std::vector<double> const b =
+ //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__specific_body_force}
+ config.getConfigParameter<std::vector<double>>(
+ "specific_body_force");
+ if (specific_body_force.size() != DisplacementDim)
+ OGS_FATAL(
+ "The size of the specific body force vector does not match the "
+ "displacement dimension. Vector size is %d, displacement "
+ "dimension is %d",
+ specific_body_force.size(), DisplacementDim);
+
+ std::copy_n(b.data(), b.size(), specific_body_force.data());
+ }
+
+ ThermoMechanicsProcessData<DisplacementDim> process_data{
+ std::move(material),
+ solid_density,
+ linear_thermal_expansion_coefficient,
+ specific_heat_capacity,
+ thermal_conductivity,
+ reference_temperature,
+ specific_body_force};
+
+ SecondaryVariableCollection secondary_variables;
+
+ NumLib::NamedFunctionCaller named_function_caller(
+ {"ThermoMechanics_temperature_displacement"});
+
+ ProcessLib::parseSecondaryVariables(config, secondary_variables,
+ named_function_caller);
+
+ return std::unique_ptr<ThermoMechanicsProcess<DisplacementDim>>{
+ new ThermoMechanicsProcess<DisplacementDim>{
+ mesh, std::move(jacobian_assembler), parameters, integration_order,
+ std::move(process_variables), std::move(process_data),
+ std::move(secondary_variables), std::move(named_function_caller)}};
+}
+
+template std::unique_ptr<Process> createThermoMechanicsProcess<2>(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config);
+
+template std::unique_ptr<Process> createThermoMechanicsProcess<3>(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config);
+
+} // namespace ThermoMechanics
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.h b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.h
new file mode 100644
index 00000000000..5b2849eedd9
--- /dev/null
+++ b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.h
@@ -0,0 +1,28 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace ThermoMechanics
+{
+template <int DisplacementDim>
+std::unique_ptr<Process> createThermoMechanicsProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config);
+
+} // namespace ThermoMechanics
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoMechanics/LocalDataInitializer.h b/ProcessLib/ThermoMechanics/LocalDataInitializer.h
new file mode 100644
index 00000000000..3e766a49756
--- /dev/null
+++ b/ProcessLib/ThermoMechanics/LocalDataInitializer.h
@@ -0,0 +1,306 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <functional>
+#include <memory>
+#include <type_traits>
+#include <typeindex>
+#include <typeinfo>
+#include <unordered_map>
+
+#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
+{
+namespace ThermoMechanics
+{
+/// 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 <typename, typename, unsigned, int> class LocalAssemblerData,
+ unsigned GlobalDim, int DisplacementDim, typename... ConstructorArgs>
+class LocalDataInitializer final
+{
+public:
+ using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
+
+ 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<NumLib::ShapeQuad4>();
+#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<NumLib::ShapeHex8>();
+#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<NumLib::ShapeQuad8>();
+ _builder[std::type_index(typeid(MeshLib::Quad9))] =
+ makeLocalAssemblerBuilder<NumLib::ShapeQuad9>();
+#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<NumLib::ShapeHex20>();
+#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<NumLib::ShapeTri3>();
+#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<NumLib::ShapeTet4>();
+#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<NumLib::ShapeTri6>();
+#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<NumLib::ShapeTet10>();
+#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<NumLib::ShapePrism6>();
+#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<NumLib::ShapePrism15>();
+#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<NumLib::ShapePyra5>();
+#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<NumLib::ShapePyra13>();
+#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,
+ 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,
+ std::forward<ConstructorArgs>(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 or this process requires higher order elements.",
+ type_idx.name());
+ }
+ }
+
+private:
+ using LADataBuilder =
+ std::function<LADataIntfPtr(MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ ConstructorArgs&&...)>;
+
+ template <typename ShapeFunction>
+ using IntegrationMethod = typename NumLib::GaussIntegrationPolicy<
+ typename ShapeFunction::MeshElement>::IntegrationMethod;
+
+ template <typename ShapeFunction>
+ using LAData =
+ LocalAssemblerData<ShapeFunction, IntegrationMethod<ShapeFunction>,
+ GlobalDim, DisplacementDim>;
+
+ /// A helper forwarding to the correct version of makeLocalAssemblerBuilder
+ /// depending whether the global dimension is less than the shape function's
+ /// dimension or not.
+ template <typename ShapeFunction>
+ static LADataBuilder makeLocalAssemblerBuilder()
+ {
+ return makeLocalAssemblerBuilder<ShapeFunction>(
+ static_cast<std::integral_constant<
+ bool, (GlobalDim >= ShapeFunction::DIM)>*>(nullptr));
+ }
+
+ /// Mapping of element types to local assembler constructors.
+ std::unordered_map<std::type_index, LADataBuilder> _builder;
+
+ NumLib::LocalToGlobalIndexMap const& _dof_table;
+
+ // local assembler builder implementations.
+private:
+ /// Generates a function that creates a new LocalAssembler of type
+ /// LAData<ShapeFunction>. Only functions with shape function's dimension
+ /// less or equal to the global dimension are instantiated, e.g. following
+ /// combinations of shape functions and global dimensions: (Line2, 1),
+ /// (Line2, 2), (Line2, 3), (Hex20, 3) but not (Hex20, 2) or (Hex20, 1).
+ template <typename ShapeFunction>
+ static LADataBuilder makeLocalAssemblerBuilder(std::true_type*)
+ {
+ return [](MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ ConstructorArgs&&... args) {
+ return LADataIntfPtr{new LAData<ShapeFunction>{
+ e, local_matrix_size, std::forward<ConstructorArgs>(args)...}};
+ };
+ }
+
+ /// Returns nullptr for shape functions whose dimensions are less than the
+ /// global dimension.
+ template <typename ShapeFunction>
+ static LADataBuilder makeLocalAssemblerBuilder(std::false_type*)
+ {
+ return nullptr;
+ }
+};
+
+} // namespace ThermoMechanics
+} // 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
diff --git a/ProcessLib/ThermoMechanics/Tests.cmake b/ProcessLib/ThermoMechanics/Tests.cmake
new file mode 100644
index 00000000000..4236f2a2cca
--- /dev/null
+++ b/ProcessLib/ThermoMechanics/Tests.cmake
@@ -0,0 +1,16 @@
+AddTest(
+ NAME 3D_ThermoElastic_Stress_Analysis
+ PATH ThermoMechanics
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS cube_1e3.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+ ABSTOL 1e-10 RELTOL 1e-12
+ DIFF_DATA
+ stress_analytical.vtu cube_1e3_tm_pcs_0_ts_80_t_72000.000000.vtu sigma_xx sigma_xx
+ stress_analytical.vtu cube_1e3_tm_pcs_0_ts_80_t_72000.000000.vtu sigma_yy sigma_yy
+ stress_analytical.vtu cube_1e3_tm_pcs_0_ts_80_t_72000.000000.vtu sigma_zz sigma_zz
+ expected_cube_1e3_tm_pcs_0_ts_80_t_72000.000000.vtu cube_1e3_tm_pcs_0_ts_80_t_72000.000000.vtu displacement displacement
+ expected_cube_1e3_tm_pcs_0_ts_80_t_72000.000000.vtu cube_1e3_tm_pcs_0_ts_80_t_72000.000000.vtu temperature temperature
+)
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
new file mode 100644
index 00000000000..48acb1e296f
--- /dev/null
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
@@ -0,0 +1,558 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include <vector>
+
+#include "MaterialLib/SolidModels/KelvinVector.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "NumLib/Function/Interpolation.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 "ThermoMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace ThermoMechanics
+{
+template <typename BMatricesType, typename ShapeMatrixType, int DisplacementDim>
+struct IntegrationPointData final
+{
+ explicit IntegrationPointData(
+ MaterialLib::Solids::MechanicsBase<DisplacementDim>& 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)
+ : sigma(std::move(other.sigma)),
+ sigma_prev(std::move(other.sigma_prev)),
+ eps(std::move(other.eps)),
+ eps_m(std::move(other.eps_m)),
+ eps_m_prev(std::move(other.eps_m_prev)),
+ solid_material(other.solid_material),
+ material_state_variables(std::move(other.material_state_variables)),
+ integration_weight(std::move(other.integration_weight)),
+ {
+ }
+#endif // _MSC_VER
+
+ typename ShapeMatrixType::NodalRowVectorType N;
+ typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;
+ typename BMatricesType::KelvinVectorType sigma, sigma_prev;
+ typename BMatricesType::KelvinVectorType eps;
+ typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;
+
+ MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material;
+ std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
+ DisplacementDim>::MaterialStateVariables>
+ material_state_variables;
+
+ double integration_weight;
+
+ void pushBackState()
+ {
+ eps_m_prev = eps_m;
+ sigma_prev = sigma;
+ material_state_variables->pushBackState();
+ }
+
+ static const int kelvin_vector_size =
+ KelvinVectorDimensions<DisplacementDim>::value;
+ using Invariants = MaterialLib::SolidModels::Invariants<kelvin_vector_size>;
+
+ template <typename DisplacementVectorType>
+ typename BMatricesType::KelvinMatrixType updateConstitutiveRelation(
+ double const t,
+ SpatialPosition const& x_position,
+ double const dt,
+ DisplacementVectorType const& u,
+ double const linear_thermal_strain)
+ {
+ // assume isotropic thermal expansion
+ eps_m.noalias() = eps - linear_thermal_strain * Invariants::identity2;
+ auto&& solution = solid_material.integrateStress(
+ t, x_position, dt, eps_m_prev, eps_m, sigma_prev,
+ *material_state_variables);
+
+ if (!solution)
+ OGS_FATAL("Computation of local constitutive relation failed.");
+
+ KelvinMatrixType<DisplacementDim> C;
+ std::tie(sigma, material_state_variables, C) = std::move(*solution);
+
+ return C;
+ }
+};
+
+/// Used by for extrapolation of the integration point values. It is ordered
+/// (and stored) by integration points.
+template <typename ShapeMatrixType>
+struct SecondaryData
+{
+ std::vector<ShapeMatrixType> N;
+};
+
+struct ThermoMechanicsLocalAssemblerInterface
+ : public ProcessLib::LocalAssemblerInterface,
+ public NumLib::ExtrapolatableElement
+{
+ virtual std::vector<double> const& getIntPtSigmaXX(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaYY(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaZZ(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaXY(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaXZ(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtSigmaYZ(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonXX(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonYY(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonZZ(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonXY(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonXZ(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonYZ(
+ std::vector<double>& cache) const = 0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ int DisplacementDim>
+class ThermoMechanicsLocalAssembler
+ : public ThermoMechanicsLocalAssemblerInterface
+{
+public:
+ using ShapeMatricesType =
+ ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>;
+
+ // Types for displacement.
+ // (Higher order elements = ShapeFunction).
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+ using BMatricesType = BMatrixPolicyType<ShapeFunction, DisplacementDim>;
+
+ using NodalForceVectorType = typename BMatricesType::NodalForceVectorType;
+ using RhsVector = typename ShapeMatricesType::template VectorType<
+ ShapeFunction::NPOINTS + ShapeFunction::NPOINTS * DisplacementDim>;
+ using JacobianMatrix = typename ShapeMatricesType::template MatrixType<
+ ShapeFunction::NPOINTS + ShapeFunction::NPOINTS * DisplacementDim,
+ ShapeFunction::NPOINTS + ShapeFunction::NPOINTS * DisplacementDim>;
+
+ ThermoMechanicsLocalAssembler(ThermoMechanicsLocalAssembler const&) =
+ delete;
+ ThermoMechanicsLocalAssembler(ThermoMechanicsLocalAssembler&&) = delete;
+
+ ThermoMechanicsLocalAssembler(
+ MeshLib::Element const& e,
+ std::size_t const /*local_matrix_size*/,
+ bool const is_axially_symmetric,
+ unsigned const integration_order,
+ ThermoMechanicsProcessData<DisplacementDim>& process_data)
+ : _process_data(process_data),
+ _integration_method(integration_order),
+ _element(e),
+ _is_axially_symmetric(is_axially_symmetric)
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ _ip_data.reserve(n_integration_points);
+ _secondary_data.N.resize(n_integration_points);
+
+ auto const shape_matrices =
+ initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, DisplacementDim>(
+ e, is_axially_symmetric, _integration_method);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ // displacement (subscript u)
+ _ip_data.emplace_back(*_process_data.material);
+ auto& ip_data = _ip_data[ip];
+ _ip_data[ip].integration_weight =
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
+
+ ip_data.sigma.resize(kelvin_vector_size);
+ ip_data.sigma_prev.resize(kelvin_vector_size);
+ ip_data.eps.resize(kelvin_vector_size);
+ ip_data.eps_m.resize(kelvin_vector_size);
+ ip_data.eps_m_prev.resize(kelvin_vector_size);
+
+ ip_data.N = shape_matrices[ip].N;
+ ip_data.dNdx = shape_matrices[ip].dNdx;
+
+ _secondary_data.N[ip] = shape_matrices[ip].N;
+ }
+ }
+
+ void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
+ std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& /*local_rhs_data*/) override
+ {
+ OGS_FATAL(
+ "ThermoMechanicsLocalAssembler: assembly without jacobian is not "
+ "implemented.");
+ }
+
+ void assembleWithJacobian(double const t,
+ std::vector<double> const& local_x,
+ std::vector<double> const& local_xdot,
+ const double /*dxdot_dx*/, const double /*dx_dx*/,
+ std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& local_rhs_data,
+ std::vector<double>& local_Jac_data) override
+ {
+ auto const local_matrix_size = local_x.size();
+ assert(local_matrix_size == temperature_size + displacement_size);
+
+ auto T = Eigen::Map<typename ShapeMatricesType::template VectorType<
+ temperature_size> const>(local_x.data() + temperature_index,
+ temperature_size);
+
+ auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
+ displacement_size> const>(local_x.data() + displacement_index,
+ displacement_size);
+
+ auto T_dot = Eigen::Map<typename ShapeMatricesType::template VectorType<
+ temperature_size> const>(local_xdot.data() + temperature_index,
+ temperature_size);
+
+ auto local_Jac = MathLib::createZeroedMatrix<JacobianMatrix>(
+ local_Jac_data, local_matrix_size, local_matrix_size);
+
+ auto local_rhs = MathLib::createZeroedVector<RhsVector>(
+ local_rhs_data, local_matrix_size);
+
+ typename ShapeMatricesType::template MatrixType<displacement_size,
+ temperature_size>
+ KuT;
+ KuT.setZero(displacement_size, temperature_size);
+
+ typename ShapeMatricesType::NodalMatrixType KTT;
+ KTT.setZero(temperature_size, temperature_size);
+
+ typename ShapeMatricesType::NodalMatrixType DTT;
+ DTT.setZero(temperature_size, temperature_size);
+
+ double const& dt = _process_data.dt;
+
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& w = _ip_data[ip].integration_weight;
+
+ auto const& dNdx = _ip_data[ip].dNdx;
+ auto const& N = _ip_data[ip].N;
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+ _element, N);
+ auto const& B = LinearBMatrix::computeBMatrix<
+ DisplacementDim, ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
+ _is_axially_symmetric);
+ auto const& sigma = _ip_data[ip].sigma;
+ auto& eps = _ip_data[ip].eps;
+
+ double const delta_T =
+ N.dot(T) - _process_data.reference_temperature;
+ // calculate thermally induced strain
+ // assume isotropic thermal expansion
+ auto const alpha =
+ _process_data.linear_thermal_expansion_coefficient(
+ t, x_position)[0];
+ double const linear_thermal_strain = alpha * delta_T;
+
+ //
+ // displacement equation, displacement part
+ //
+ eps.noalias() = B * u;
+ auto C = _ip_data[ip].updateConstitutiveRelation(
+ t, x_position, dt, u, linear_thermal_strain);
+
+ local_Jac
+ .template block<displacement_size, displacement_size>(
+ displacement_index, displacement_index)
+ .noalias() += B.transpose() * C * B * w;
+
+ typename ShapeMatricesType::template MatrixType<DisplacementDim,
+ displacement_size>
+ N_u = ShapeMatricesType::template MatrixType<
+ DisplacementDim,
+ displacement_size>::Zero(DisplacementDim,
+ displacement_size);
+
+ for (int i = 0; i < DisplacementDim; ++i)
+ N_u.template block<1, displacement_size / DisplacementDim>(
+ i, i * displacement_size / DisplacementDim)
+ .noalias() = N;
+
+ using Invariants =
+ MaterialLib::SolidModels::Invariants<kelvin_vector_size>;
+
+ // calculate real density
+ auto const rho_sr = _process_data.solid_density(t, x_position)[0];
+ double const rho_s = rho_sr * (1 - 3 * linear_thermal_strain);
+
+ auto const& b = _process_data.specific_body_force;
+ local_rhs
+ .template block<displacement_size, 1>(displacement_index, 0)
+ .noalias() -=
+ (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
+ //
+ // displacement equation, temperature part
+ //
+ KuT.noalias() +=
+ B.transpose() * C * alpha * Invariants::identity2 * N * w;
+
+ //
+ // temperature equation, temperature part;
+ //
+ auto const lambda =
+ _process_data.thermal_conductivity(t, x_position)[0];
+ KTT.noalias() += dNdx.transpose() * lambda * dNdx * w;
+
+ auto const c =
+ _process_data.specific_heat_capacity(t, x_position)[0];
+ DTT.noalias() += N.transpose() * rho_s * c * N * w;
+ }
+ // temperature equation, temperature part
+ local_Jac
+ .template block<temperature_size, temperature_size>(
+ temperature_index, temperature_index)
+ .noalias() += KTT + DTT / dt;
+
+ // displacement equation, temperature part
+ local_Jac
+ .template block<displacement_size, temperature_size>(
+ displacement_index, temperature_index)
+ .noalias() -= KuT;
+
+ local_rhs.template block<temperature_size, 1>(temperature_index, 0)
+ .noalias() -= KTT * T + DTT * T_dot;
+ }
+
+ void preTimestepConcrete(std::vector<double> 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<const Eigen::RowVectorXd> 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<const Eigen::RowVectorXd>(N.data(), N.size());
+ }
+
+ std::vector<double> const& getIntPtSigmaXX(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtSigma(cache, 0);
+ }
+
+ std::vector<double> const& getIntPtSigmaYY(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtSigma(cache, 1);
+ }
+
+ std::vector<double> const& getIntPtSigmaZZ(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtSigma(cache, 2);
+ }
+
+ std::vector<double> const& getIntPtSigmaXY(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtSigma(cache, 3);
+ }
+
+ std::vector<double> const& getIntPtSigmaYZ(
+ std::vector<double>& cache) const override
+ {
+ assert(DisplacementDim == 3);
+ return getIntPtSigma(cache, 4);
+ }
+
+ std::vector<double> const& getIntPtSigmaXZ(
+ std::vector<double>& cache) const override
+ {
+ assert(DisplacementDim == 3);
+ return getIntPtSigma(cache, 5);
+ }
+
+ std::vector<double> const& getIntPtEpsilonXX(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtEpsilon(cache, 0);
+ }
+
+ std::vector<double> const& getIntPtEpsilonYY(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtEpsilon(cache, 1);
+ }
+
+ std::vector<double> const& getIntPtEpsilonZZ(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtEpsilon(cache, 2);
+ }
+
+ std::vector<double> const& getIntPtEpsilonXY(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtEpsilon(cache, 3);
+ }
+
+ std::vector<double> const& getIntPtEpsilonYZ(
+ std::vector<double>& cache) const override
+ {
+ assert(DisplacementDim == 3);
+ return getIntPtEpsilon(cache, 4);
+ }
+
+ std::vector<double> const& getIntPtEpsilonXZ(
+ std::vector<double>& cache) const override
+ {
+ assert(DisplacementDim == 3);
+ return getIntPtEpsilon(cache, 5);
+ }
+
+private:
+ std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
+ std::size_t const component) const
+ {
+ cache.clear();
+ cache.reserve(_ip_data.size());
+
+ for (auto const& ip_data : _ip_data)
+ {
+ if (component < 3) // xx, yy, zz components
+ cache.push_back(ip_data.sigma[component]);
+ else // mixed xy, yz, xz components
+ cache.push_back(ip_data.sigma[component] / std::sqrt(2));
+ }
+
+ return cache;
+ }
+ std::vector<double> const& getIntPtEpsilon(
+ std::vector<double>& cache, std::size_t const component) const
+ {
+ cache.clear();
+ cache.reserve(_ip_data.size());
+
+ for (auto const& ip_data : _ip_data)
+ {
+ if (component < 3) // xx, yy, zz components
+ cache.push_back(ip_data.eps[component]);
+ else // mixed xy, yz, xz components
+ cache.push_back(ip_data.eps[component] / std::sqrt(2));
+ }
+
+ return cache;
+ }
+
+ ThermoMechanicsProcessData<DisplacementDim>& _process_data;
+
+ std::vector<IntegrationPointData<BMatricesType, ShapeMatricesType,
+ DisplacementDim>> _ip_data;
+
+ IntegrationMethod _integration_method;
+ MeshLib::Element const& _element;
+ bool const _is_axially_symmetric;
+ SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
+
+ static const int temperature_index = 0;
+ static const int temperature_size = ShapeFunction::NPOINTS;
+ static const int displacement_index = ShapeFunction::NPOINTS;
+ static const int displacement_size =
+ ShapeFunction::NPOINTS * DisplacementDim;
+ static const int kelvin_vector_size =
+ KelvinVectorDimensions<DisplacementDim>::value;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim, int DisplacementDim>
+class LocalAssemblerData final
+ : public ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
+ DisplacementDim>
+{
+public:
+ LocalAssemblerData(LocalAssemblerData const&) = delete;
+ LocalAssemblerData(LocalAssemblerData&&) = delete;
+
+ LocalAssemblerData(
+ MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ ThermoMechanicsProcessData<DisplacementDim>& process_data)
+ : ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
+ DisplacementDim>(
+ e, local_matrix_size, is_axially_symmetric, integration_order,
+ process_data)
+ {
+ }
+};
+
+} // namespace ThermoMechanics
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-fwd.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-fwd.h
new file mode 100644
index 00000000000..fc48ba4a62a
--- /dev/null
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess-fwd.h
@@ -0,0 +1,15 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "ThermoMechanicsProcess.h"
+
+extern template class ProcessLib::ThermoMechanics::ThermoMechanicsProcess<2>;
+extern template class ProcessLib::ThermoMechanics::ThermoMechanicsProcess<3>;
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
new file mode 100644
index 00000000000..55fddd02e3d
--- /dev/null
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
@@ -0,0 +1,21 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "ThermoMechanicsProcess.h"
+#include "ThermoMechanicsProcess-fwd.h"
+
+namespace ProcessLib
+{
+namespace ThermoMechanics
+{
+template class ThermoMechanicsProcess<2>;
+template class ThermoMechanicsProcess<3>;
+
+} // namespace ThermoMechanics
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
new file mode 100644
index 00000000000..7e6c4c290da
--- /dev/null
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
@@ -0,0 +1,220 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <cassert>
+
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "ProcessLib/Process.h"
+#include "ProcessLib/ThermoMechanics/CreateLocalAssemblers.h"
+
+#include "ThermoMechanicsFEM.h"
+#include "ThermoMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace ThermoMechanics
+{
+template <int DisplacementDim>
+class ThermoMechanicsProcess final : public Process
+{
+ using Base = Process;
+
+public:
+ ThermoMechanicsProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+ jacobian_assembler,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ std::vector<std::reference_wrapper<ProcessVariable>>&&
+ process_variables,
+ ThermoMechanicsProcessData<DisplacementDim>&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller)
+ : Process(mesh, std::move(jacobian_assembler), parameters,
+ integration_order, 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 = ThermoMechanicsLocalAssemblerInterface;
+
+ void initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order) override
+ {
+ ProcessLib::ThermoMechanics::createLocalAssemblers<DisplacementDim,
+ LocalAssemblerData>(
+ mesh.getDimension(), mesh.getElements(), dof_table,
+ _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
+ _process_data);
+
+ // TODO move the two data members somewhere else.
+ // for extrapolation of secondary variables
+ std::vector<MeshLib::MeshSubsets> all_mesh_subsets_single_component;
+ all_mesh_subsets_single_component.emplace_back(
+ _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,
+ &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaXX));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_yy", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaYY));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_zz", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaZZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_xy", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaXY));
+
+ if (DisplacementDim == 3)
+ {
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_xz", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaXZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_yz", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaYZ));
+ }
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_xx", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonXX));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_yy", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonYY));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_zz", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonZZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_xy", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonXY));
+ if (DisplacementDim == 3)
+ {
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_yz", 1,
+ makeExtrapolator(getExtrapolator(), _local_assemblers,
+ &ThermoMechanicsLocalAssemblerInterface::
+ getIntPtEpsilonYZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_xz", 1,
+ makeExtrapolator(getExtrapolator(), _local_assemblers,
+ &ThermoMechanicsLocalAssemblerInterface::
+ getIntPtEpsilonXZ));
+ }
+ }
+
+ void assembleConcreteProcess(
+ const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b, StaggeredCouplingTerm const& coupling_term) override
+ {
+ DBUG("Assemble ThermoMechanicsProcess.");
+
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assemble,
+ _local_assemblers, *_local_to_global_index_map, t, x, M, K, b,
+ coupling_term);
+ }
+
+ void assembleWithJacobianConcreteProcess(
+ const double t, GlobalVector const& x, GlobalVector const& xdot,
+ const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
+ GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+ StaggeredCouplingTerm const& coupling_term) override
+ {
+ DBUG("AssembleJacobian ThermoMechanicsProcess.");
+
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+ _local_assemblers, *_local_to_global_index_map, t, x, xdot,
+ dxdot_dx, dx_dx, M, K, b, Jac, coupling_term);
+ }
+
+ void preTimestepConcreteProcess(GlobalVector const& x, double const t,
+ double const dt) override
+ {
+ DBUG("PreTimestep ThermoMechanicsProcess.");
+
+ _process_data.dt = dt;
+ _process_data.t = t;
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &ThermoMechanicsLocalAssemblerInterface::preTimestep,
+ _local_assemblers, *_local_to_global_index_map, x, t, dt);
+ }
+
+ void postTimestepConcreteProcess(GlobalVector const& x) override
+ {
+ DBUG("PostTimestep ThermoMechanicsProcess.");
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &ThermoMechanicsLocalAssemblerInterface::postTimestep,
+ _local_assemblers, *_local_to_global_index_map, x);
+ }
+
+private:
+ std::vector<MeshLib::Node*> _base_nodes;
+ std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_base_nodes;
+ ThermoMechanicsProcessData<DisplacementDim> _process_data;
+
+ std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap>
+ _local_to_global_index_map_single_component;
+};
+
+} // namespace ThermoMechanics
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
new file mode 100644
index 00000000000..432db7b1927
--- /dev/null
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
@@ -0,0 +1,80 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+namespace MeshLib
+{
+class Element;
+}
+
+namespace ProcessLib
+{
+namespace ThermoMechanics
+{
+template <int DisplacementDim>
+struct ThermoMechanicsProcessData
+{
+ ThermoMechanicsProcessData(
+ std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>&&
+ material_,
+ Parameter<double> const& solid_density_,
+ Parameter<double> const& linear_thermal_expansion_coefficient_,
+ Parameter<double> const& specific_heat_capacity_,
+ Parameter<double> const& thermal_conductivity_,
+ double const reference_temperature_,
+ Eigen::Matrix<double, DisplacementDim, 1> const& specific_body_force_)
+ : material{std::move(material_)},
+ solid_density(solid_density_),
+ linear_thermal_expansion_coefficient(
+ linear_thermal_expansion_coefficient_),
+ specific_heat_capacity(specific_heat_capacity_),
+ thermal_conductivity(thermal_conductivity_),
+ reference_temperature(reference_temperature_),
+ specific_body_force(specific_body_force_)
+ {
+ }
+
+ ThermoMechanicsProcessData(ThermoMechanicsProcessData&& other)
+ : material{std::move(other.material)},
+ solid_density(other.solid_density),
+ linear_thermal_expansion_coefficient(
+ other.linear_thermal_expansion_coefficient),
+ specific_heat_capacity(other.specific_heat_capacity),
+ thermal_conductivity(other.thermal_conductivity),
+ reference_temperature(other.reference_temperature),
+ specific_body_force(other.specific_body_force),
+ dt(other.dt),
+ t(other.t)
+ {
+ }
+
+ //! Copies are forbidden.
+ ThermoMechanicsProcessData(ThermoMechanicsProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(ThermoMechanicsProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(ThermoMechanicsProcessData&&) = delete;
+
+ std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
+ material;
+ Parameter<double> const& solid_density;
+ Parameter<double> const& linear_thermal_expansion_coefficient;
+ Parameter<double> const& specific_heat_capacity;
+ Parameter<double> const& thermal_conductivity;
+ double const reference_temperature;
+ Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
+ double dt;
+ double t;
+};
+
+} // namespace ThermoMechanics
+} // namespace ProcessLib
--
GitLab