diff --git a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
index ad33b3cca02a12d530a1c174edb9385e6fdf088f..914064a0d54b703726a45c1c2132f0f65547363d 100644
--- a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
+++ b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
@@ -51,6 +51,10 @@ std::unique_ptr<Process> createThermoMechanicsProcess(
//! \ogs_file_param{prj__processes__process__THERMO_MECHANICS__process_variables}
auto const pv_config = config.getConfigSubtree("process_variables");
+ // Process IDs, which are set according to the appearance order of the
+ int heat_conduction_process_id = 0;
+ int mechanics_process_id = 0;
+
ProcessVariable* variable_T;
ProcessVariable* variable_u;
std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
@@ -78,13 +82,11 @@ std::unique_ptr<Process> createThermoMechanicsProcess(
}
variable_T = &process_variables[0][0].get();
variable_u = &process_variables[1][0].get();
+ // process variables. Up to now, the ordering is fixed as:
+ heat_conduction_process_id = 0;
+ mechanics_process_id = 1;
}
- // Process IDs, which are set according to the appearance order of the
- // process variables. Up to now, the ordering is fixed as:
- int heat_conduction_process_id = 0;
- int mechanics_process_id = 1;
-
DBUG("Associate displacement with process variable '%s'.",
variable_u->getName().c_str());
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
index 83b371381ee6136cc00dfca0a814fa28cfb903a1..0e3adb141f65201f48adebe106e0688c153c397b 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
@@ -12,7 +12,9 @@
#include <cassert>
#include "BaseLib/Functional.h"
+#include "NumLib/DOF/ComputeSparsityPattern.h"
#include "NumLib/DOF/DOFTableUtil.h"
+
#include "ProcessLib/SmallDeformation/CreateLocalAssemblers.h"
#include "ThermoMechanicsFEM.h"
@@ -116,18 +118,35 @@ bool ThermoMechanicsProcess<DisplacementDim>::isLinear() const
}
template <int DisplacementDim>
-void ThermoMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
- NumLib::LocalToGlobalIndexMap const& dof_table,
- MeshLib::Mesh const& mesh,
- unsigned const integration_order)
+MathLib::MatrixSpecifications
+ThermoMechanicsProcess<DisplacementDim>::getMatrixSpecifications(
+ const int process_id) const
{
- ProcessLib::SmallDeformation::createLocalAssemblers<
- DisplacementDim, ThermoMechanicsLocalAssembler>(
- mesh.getElements(), dof_table, _local_assemblers,
- mesh.isAxiallySymmetric(), integration_order, _process_data);
+ // For the monolithic scheme or the M process (deformation) in the staggered
+ // scheme.
+ if (_use_monolithic_scheme || process_id == _mechanics_process_id)
+ {
+ auto const& l = *_local_to_global_index_map;
+ return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+ &l.getGhostIndices(), &this->_sparsity_pattern};
+ }
+
+ // For staggered scheme and T process.
+ auto const& l = *_local_to_global_index_map_single_component;
+ return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+ &l.getGhostIndices(), &_sparsity_pattern_with_single_component};
+}
+
+template <int DisplacementDim>
+void ThermoMechanicsProcess<DisplacementDim>::constructDofTable()
+{
+ // Note: the heat conduction process and the mechnical process use the same
+ // order of shape functions.
+
+ // Create single component dof in every of the mesh's nodes.
+ _mesh_subset_all_nodes =
+ std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
- // TODO move the two data members somewhere else.
- // for extrapolation of secondary variables
std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
*_mesh_subset_all_nodes};
_local_to_global_index_map_single_component.reset(
@@ -136,6 +155,64 @@ void ThermoMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
// by location order is needed for output
NumLib::ComponentOrder::BY_LOCATION));
+ // Vector of mesh subsets.
+ std::vector<MeshLib::MeshSubset> all_mesh_subsets;
+
+ // Vector of the number of variable components
+ std::vector<int> vec_var_n_components;
+ if (_use_monolithic_scheme)
+ {
+ // Collect the mesh subsets in a vector for each variables' components.
+ for (ProcessVariable const& pv : _process_variables[0])
+ {
+ std::generate_n(std::back_inserter(all_mesh_subsets),
+ pv.getNumberOfComponents(),
+ [&]() { return *_mesh_subset_all_nodes; });
+ }
+
+ // Create a vector of the number of variable components
+ for (ProcessVariable const& pv : _process_variables[0])
+ {
+ vec_var_n_components.push_back(pv.getNumberOfComponents());
+ }
+ }
+ else // for staggered scheme
+ {
+ // Collect the mesh subsets in a vector for each variable components for
+ // the mechanical process.
+ std::generate_n(std::back_inserter(all_mesh_subsets),
+ _process_variables[_mechanics_process_id][0]
+ .get()
+ .getNumberOfComponents(),
+ [&]() { return *_mesh_subset_all_nodes; });
+
+ // Create a vector of the number of variable components.
+ vec_var_n_components.push_back(
+ _process_variables[_mechanics_process_id][0]
+ .get()
+ .getNumberOfComponents());
+
+ _sparsity_pattern_with_single_component =
+ NumLib::computeSparsityPattern(
+ *_local_to_global_index_map_single_component, _mesh);
+ }
+ _local_to_global_index_map =
+ std::make_unique<NumLib::LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets), vec_var_n_components,
+ NumLib::ComponentOrder::BY_LOCATION);
+}
+
+template <int DisplacementDim>
+void ThermoMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+{
+ ProcessLib::SmallDeformation::createLocalAssemblers<
+ DisplacementDim, ThermoMechanicsLocalAssembler>(
+ mesh.getElements(), dof_table, _local_assemblers,
+ mesh.isAxiallySymmetric(), integration_order, _process_data);
+
_secondary_variables.addSecondaryVariable(
"sigma",
makeExtrapolator(
@@ -203,6 +280,28 @@ void ThermoMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
}
}
+template <int DisplacementDim>
+void ThermoMechanicsProcess<DisplacementDim>::initializeBoundaryConditions()
+{
+ if (_use_monolithic_scheme)
+ {
+ const int process_id_of_thermomechancs = 0;
+ initializeProcessBoundaryConditionsAndSourceTerms(
+ *_local_to_global_index_map, process_id_of_thermomechancs);
+ return;
+ }
+
+ // Staggered scheme:
+ // for the equations of heat conduction
+ initializeProcessBoundaryConditionsAndSourceTerms(
+ *_local_to_global_index_map_single_component,
+ _heat_conduction_process_id);
+
+ // for the equations of deformation.
+ initializeProcessBoundaryConditionsAndSourceTerms(
+ *_local_to_global_index_map, _mechanics_process_id);
+}
+
template <int DisplacementDim>
void ThermoMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
@@ -234,37 +333,80 @@ void ThermoMechanicsProcess<DisplacementDim>::
DBUG("AssembleJacobian ThermoMechanicsProcess.");
std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
- dof_table = {std::ref(*_local_to_global_index_map)};
+ dof_tables;
+ // For the monolithic scheme
+ if (_use_monolithic_scheme)
+ {
+ DBUG(
+ "Assemble the Jacobian of ThermoMechanics for the monolithic"
+ " scheme.");
+ dof_tables.emplace_back(*_local_to_global_index_map);
+ }
+ else
+ {
+ // For the staggered scheme
+ if (_coupled_solutions->process_id == _heat_conduction_process_id)
+ {
+ DBUG(
+ "Assemble the Jacobian equations of heat conduction process in "
+ "ThermoMechanics for the staggered scheme.");
+ }
+ else
+ {
+ DBUG(
+ "Assemble the Jacobian equations of mechanical process in "
+ "ThermoMechanics for the staggered scheme.");
+ }
+
+ // For the flexible appearance order of processes in the coupling.
+ if (_heat_conduction_process_id ==
+ 0) // First: the heat conduction process
+ {
+ dof_tables.emplace_back(
+ *_local_to_global_index_map_single_component);
+ dof_tables.emplace_back(*_local_to_global_index_map);
+ }
+ else // vice versa
+ {
+ dof_tables.emplace_back(*_local_to_global_index_map);
+ dof_tables.emplace_back(
+ *_local_to_global_index_map_single_component);
+ }
+
+ setCoupledSolutionsOfPreviousTimeStep();
+ }
+
const int process_id =
_use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
- // Call global assembler for each local assembly item.
GlobalExecutor::executeSelectedMemberDereferenced(
_global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
- _local_assemblers, pv.getActiveElementIDs(), dof_table, t, x, xdot,
+ _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, x, xdot,
dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
// TODO (naumov): Refactor the copy rhs part. This is copy from HM.
auto copyRhs = [&](int const variable_id, auto& output_vector) {
if (_use_monolithic_scheme)
{
- transformVariableFromGlobalVector(b, variable_id, dof_table[0],
+ transformVariableFromGlobalVector(b, variable_id, dof_tables[0],
output_vector,
std::negate<double>());
}
else
{
transformVariableFromGlobalVector(
- b, 0, dof_table[_coupled_solutions->process_id], output_vector,
+ b, 0, dof_tables[_coupled_solutions->process_id], output_vector,
std::negate<double>());
}
};
- if (_use_monolithic_scheme || _coupled_solutions->process_id == 0)
+ if (_use_monolithic_scheme ||
+ _coupled_solutions->process_id == _heat_conduction_process_id)
{
copyRhs(0, *_heat_flux);
}
- if (_use_monolithic_scheme || _coupled_solutions->process_id == 1)
+ if (_use_monolithic_scheme ||
+ _coupled_solutions->process_id == _mechanics_process_id)
{
copyRhs(1, *_nodal_forces);
}
@@ -282,9 +424,31 @@ void ThermoMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
- GlobalExecutor::executeSelectedMemberOnDereferenced(
- &ThermoMechanicsLocalAssemblerInterface::preTimestep, _local_assemblers,
- pv.getActiveElementIDs(), *_local_to_global_index_map, x, t, dt);
+ assert(process_id < 2);
+
+
+ if (_use_monolithic_scheme || process_id == _mechanics_process_id)
+ {
+ GlobalExecutor::executeSelectedMemberOnDereferenced(
+ &ThermoMechanicsLocalAssemblerInterface::preTimestep,
+ _local_assemblers, pv.getActiveElementIDs(),
+ *_local_to_global_index_map, x, t, dt);
+ return;
+ }
+
+ // For the staggered scheme.
+ if (!_previous_T)
+ {
+ _previous_T = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
+ }
+ else
+ {
+ auto& x0 = *_previous_T;
+ MathLib::LinAlg::copy(x, x0);
+ }
+
+ auto& x0 = *_previous_T;
+ MathLib::LinAlg::setLocalAccessibleVector(x0);
}
template <int DisplacementDim>
@@ -292,6 +456,9 @@ void ThermoMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
GlobalVector const& x, const double /*t*/, const double /*delta_t*/,
int const process_id)
{
+ if (process_id == _heat_conduction_process_id)
+ return;
+
DBUG("PostTimestep ThermoMechanicsProcess.");
ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
@@ -302,6 +469,27 @@ void ThermoMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
*_local_to_global_index_map, x);
}
+template <int DisplacementDim>
+void ThermoMechanicsProcess<
+ DisplacementDim>::setCoupledSolutionsOfPreviousTimeStep()
+{
+ _coupled_solutions->coupled_xs_t0.resize(1);
+ _coupled_solutions->coupled_xs_t0[0] = _previous_T.get();
+}
+
+template <int DisplacementDim>
+NumLib::LocalToGlobalIndexMap const&
+ThermoMechanicsProcess<DisplacementDim>::getDOFTable(const int process_id) const
+{
+ if (_mechanics_process_id == process_id)
+ {
+ return *_local_to_global_index_map;
+ }
+
+ // For the equation of pressure
+ return *_local_to_global_index_map_single_component;
+}
+
template class ThermoMechanicsProcess<2>;
template class ThermoMechanicsProcess<3>;
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
index 19247f6f6a4f31f7cb25ce84a1382669403eb139..17b2663f50af546d1be2408eae3faf0a70faa392 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
@@ -83,12 +83,28 @@ public:
bool isLinear() const override;
//! @}
+ /**
+ * Get the size and the sparse pattern of the global matrix in order to
+ * create the global matrices and vectors for the system equations of this
+ * process.
+ *
+ * @param process_id Process ID. If the monolithic scheme is applied,
+ * process_id = 0.
+ * @return Matrix specifications including size and sparse pattern.
+ */
+ MathLib::MatrixSpecifications getMatrixSpecifications(
+ const int process_id) const override;
+
private:
+ void constructDofTable() override;
+
void initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table,
MeshLib::Mesh const& mesh,
unsigned const integration_order) override;
+ void initializeBoundaryConditions() override;
+
void assembleConcreteProcess(const double t, GlobalVector const& x,
GlobalMatrix& M, GlobalMatrix& K,
GlobalVector& b) override;
@@ -106,6 +122,9 @@ private:
const double delta_t,
int const process_id) override;
+ NumLib::LocalToGlobalIndexMap const& getDOFTable(
+ const int process_id) const override;
+
private:
ThermoMechanicsProcessData<DisplacementDim> _process_data;
@@ -115,6 +134,10 @@ private:
std::unique_ptr<NumLib::LocalToGlobalIndexMap>
_local_to_global_index_map_single_component;
+ /// Sparsity pattern for the heat conduction equation, and it is initialized
+ /// only if the staggered scheme is used.
+ GlobalSparsityPattern _sparsity_pattern_with_single_component;
+
MeshLib::PropertyVector<double>* _nodal_forces = nullptr;
MeshLib::PropertyVector<double>* _heat_flux = nullptr;
@@ -123,6 +146,13 @@ private:
/// ID of heat conduction process.
int const _heat_conduction_process_id;
+
+ /// Temperature of the previous time step for staggered scheme.
+ std::unique_ptr<GlobalVector> _previous_T;
+
+ /// Set the increment solutions of the present time step to the coupled
+ /// term.
+ void setCoupledSolutionsOfPreviousTimeStep();
};
extern template class ThermoMechanicsProcess<2>;