From 2f764edf4ca95a1b0d391a4e582b6f6417756817 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Thu, 3 May 2018 17:03:41 +0200
Subject: [PATCH] [PL] HM; Move process implementation into cpp.
---
.../HydroMechanicsProcess-impl.h | 384 ------------------
.../HydroMechanics/HydroMechanicsProcess.cpp | 367 ++++++++++++++++-
2 files changed, 366 insertions(+), 385 deletions(-)
delete mode 100644 ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h
deleted file mode 100644
index 14aec3eaf06..00000000000
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h
+++ /dev/null
@@ -1,384 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2018, 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 "MeshLib/Elements/Utils.h"
-#include "NumLib/DOF/ComputeSparsityPattern.h"
-#include "ProcessLib/HydroMechanics/CreateLocalAssemblers.h"
-#include "ProcessLib/Process.h"
-
-#include "HydroMechanicsFEM.h"
-#include "HydroMechanicsProcessData.h"
-#include "HydroMechanicsProcess.h"
-
-namespace ProcessLib
-{
-namespace HydroMechanics
-{
-template <int DisplacementDim>
-HydroMechanicsProcess<DisplacementDim>::HydroMechanicsProcess(
- 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::vector<std::reference_wrapper<ProcessVariable>>>&&
- process_variables,
- HydroMechanicsProcessData<DisplacementDim>&& process_data,
- SecondaryVariableCollection&& secondary_variables,
- NumLib::NamedFunctionCaller&& named_function_caller,
- bool const use_monolithic_scheme)
- : Process(mesh, std::move(jacobian_assembler), parameters,
- integration_order, std::move(process_variables),
- std::move(secondary_variables), std::move(named_function_caller),
- use_monolithic_scheme),
- _process_data(std::move(process_data))
-{
-}
-
-template <int DisplacementDim>
-bool HydroMechanicsProcess<DisplacementDim>::isLinear() const
-{
- return false;
-}
-
-template <int DisplacementDim>
-MathLib::MatrixSpecifications
-HydroMechanicsProcess<DisplacementDim>::getMatrixSpecifications(
- const int process_id) const
-{
- // For the monolithic scheme or the M process (deformation) in the staggered
- // scheme.
- if (_use_monolithic_scheme || process_id == 1)
- {
- auto const& l = *_local_to_global_index_map;
- return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
- &l.getGhostIndices(), &this->_sparsity_pattern};
- }
-
- // For staggered scheme and H process (pressure).
- auto const& l = *_local_to_global_index_map_with_base_nodes;
- return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
- &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
-}
-
-template <int DisplacementDim>
-void HydroMechanicsProcess<DisplacementDim>::constructDofTable()
-{
- // Create single component dof in every of the mesh's nodes.
- _mesh_subset_all_nodes =
- std::make_unique<MeshLib::MeshSubset>(_mesh, &_mesh.getNodes());
- // Create single component dof in the mesh's base nodes.
- _base_nodes = MeshLib::getBaseNodes(_mesh.getElements());
- _mesh_subset_base_nodes =
- std::make_unique<MeshLib::MeshSubset>(_mesh, &_base_nodes);
-
- // TODO move the two data members somewhere else.
- // for extrapolation of secondary variables of stress or strain
- 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 =
- std::make_unique<NumLib::LocalToGlobalIndexMap>(
- std::move(all_mesh_subsets_single_component),
- // by location order is needed for output
- NumLib::ComponentOrder::BY_LOCATION);
-
- if (_use_monolithic_scheme)
- {
- // For pressure, which is the first
- std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
- all_mesh_subsets.emplace_back(_mesh_subset_base_nodes.get());
-
- // For displacement.
- const int monolithic_process_id = 0;
- std::generate_n(
- std::back_inserter(all_mesh_subsets),
- getProcessVariables(monolithic_process_id)[1]
- .get()
- .getNumberOfComponents(),
- [&]() {
- return MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()};
- });
-
- std::vector<int> const vec_n_components{1, DisplacementDim};
- _local_to_global_index_map =
- std::make_unique<NumLib::LocalToGlobalIndexMap>(
- std::move(all_mesh_subsets), vec_n_components,
- NumLib::ComponentOrder::BY_LOCATION);
- assert(_local_to_global_index_map);
- }
- else
- {
- // For displacement equation.
- const int process_id = 1;
- std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
- std::generate_n(
- std::back_inserter(all_mesh_subsets),
- getProcessVariables(process_id)[0].get().getNumberOfComponents(),
- [&]() {
- return MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()};
- });
-
- std::vector<int> const vec_n_components{DisplacementDim};
- _local_to_global_index_map =
- std::make_unique<NumLib::LocalToGlobalIndexMap>(
- std::move(all_mesh_subsets), vec_n_components,
- NumLib::ComponentOrder::BY_LOCATION);
-
- // For pressure equation.
- // Collect the mesh subsets with base nodes in a vector.
- std::vector<MeshLib::MeshSubsets> all_mesh_subsets_base_nodes;
- all_mesh_subsets_base_nodes.emplace_back(_mesh_subset_base_nodes.get());
- _local_to_global_index_map_with_base_nodes =
- std::make_unique<NumLib::LocalToGlobalIndexMap>(
- std::move(all_mesh_subsets_base_nodes),
- // by location order is needed for output
- NumLib::ComponentOrder::BY_LOCATION);
-
- _sparsity_pattern_with_linear_element = NumLib::computeSparsityPattern(
- *_local_to_global_index_map_with_base_nodes, _mesh);
-
- assert(_local_to_global_index_map);
- assert(_local_to_global_index_map_with_base_nodes);
- }
-}
-
-template <int DisplacementDim>
-void HydroMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
- NumLib::LocalToGlobalIndexMap const& dof_table,
- MeshLib::Mesh const& mesh,
- unsigned const integration_order)
-{
- const int mechanical_process_id = _use_monolithic_scheme ? 0 : 1;
- const int deformation_variable_id = _use_monolithic_scheme ? 1 : 0;
- ProcessLib::HydroMechanics::createLocalAssemblers<
- DisplacementDim, HydroMechanicsLocalAssembler>(
- mesh.getDimension(), mesh.getElements(), dof_table,
- // use displacement process variable to set shape function order
- getProcessVariables(mechanical_process_id)[deformation_variable_id]
- .get()
- .getShapeFunctionOrder(),
- _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
- _process_data);
-
- Base::_secondary_variables.addSecondaryVariable(
- "sigma_xx",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaXX));
-
- Base::_secondary_variables.addSecondaryVariable(
- "sigma_yy",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaYY));
-
- Base::_secondary_variables.addSecondaryVariable(
- "sigma_zz",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaZZ));
-
- Base::_secondary_variables.addSecondaryVariable(
- "sigma_xy",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaXY));
-
- if (DisplacementDim == 3)
- {
- Base::_secondary_variables.addSecondaryVariable(
- "sigma_xz",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaXZ));
-
- Base::_secondary_variables.addSecondaryVariable(
- "sigma_yz",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaYZ));
- }
-
- Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xx",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonXX));
-
- Base::_secondary_variables.addSecondaryVariable(
- "epsilon_yy",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonYY));
-
- Base::_secondary_variables.addSecondaryVariable(
- "epsilon_zz",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonZZ));
-
- Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xy",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonXY));
-
- Base::_secondary_variables.addSecondaryVariable(
- "velocity",
- makeExtrapolator(mesh.getDimension(), getExtrapolator(),
- _local_assemblers,
- &LocalAssemblerInterface::getIntPtDarcyVelocity));
-}
-
-template <int DisplacementDim>
-void HydroMechanicsProcess<DisplacementDim>::initializeBoundaryConditions()
-{
- if (_use_monolithic_scheme)
- {
- const int process_id_of_hydromechancs = 0;
- initializeProcessBoundaryConditionsAndSourceTerms(
- *_local_to_global_index_map, process_id_of_hydromechancs);
- return;
- }
-
- // Staggered scheme:
- // for the equations of pressure
- const int hydraulic_process_id = 0;
- initializeProcessBoundaryConditionsAndSourceTerms(
- *_local_to_global_index_map_with_base_nodes, hydraulic_process_id);
-
- // for the equations of deformation.
- const int mechanical_process_id = 1;
- initializeProcessBoundaryConditionsAndSourceTerms(
- *_local_to_global_index_map, mechanical_process_id);
-}
-
-template <int DisplacementDim>
-void HydroMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
- const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
- GlobalVector& b)
-{
- DBUG("Assemble the equations for HydroMechanics");
-
- // Note: This assembly function is for the Picard nonlinear solver. Since
- // only the Newton-Raphson method is employed to simulate coupled HM
- // processes in this class, this function is actually not used so far.
-
- std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
- dof_table = {std::ref(*_local_to_global_index_map)};
- // Call global assembler for each local assembly item.
- GlobalExecutor::executeMemberDereferenced(
- _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
- dof_table, t, x, M, K, b, _coupled_solutions);
-}
-
-template <int DisplacementDim>
-void HydroMechanicsProcess<DisplacementDim>::
- 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)
-{
- std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
- dof_tables;
- // For the monolithic scheme
- if (_use_monolithic_scheme)
- {
- DBUG(
- "Assemble the Jacobian of HydroMechanics for the monolithic"
- " scheme.");
- dof_tables.emplace_back(*_local_to_global_index_map);
- }
- else
- {
- // For the staggered scheme
- if (_coupled_solutions->process_id == 0)
- {
- DBUG(
- "Assemble the Jacobian equations of liquid fluid process in "
- "HydroMechanics for the staggered scheme.");
- }
- else
- {
- DBUG(
- "Assemble the Jacobian equations of mechanical process in "
- "HydroMechanics for the staggered scheme.");
- }
- dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
- dof_tables.emplace_back(*_local_to_global_index_map);
- }
-
- GlobalExecutor::executeMemberDereferenced(
- _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
- _local_assemblers, dof_tables, t, x, xdot, dxdot_dx, dx_dx, M, K, b,
- Jac, _coupled_solutions);
-}
-
-template <int DisplacementDim>
-void HydroMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
- GlobalVector const& x, double const t, double const dt,
- const int process_id)
-{
- DBUG("PreTimestep HydroMechanicsProcess.");
-
- _process_data.dt = dt;
- _process_data.t = t;
-
- if (hasMechanicalProcess(process_id))
- GlobalExecutor::executeMemberOnDereferenced(
- &LocalAssemblerInterface::preTimestep, _local_assemblers,
- *_local_to_global_index_map, x, t, dt);
-}
-
-template <int DisplacementDim>
-void HydroMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
- GlobalVector const& x, const int process_id)
-{
- DBUG("PostTimestep HydroMechanicsProcess.");
- GlobalExecutor::executeMemberOnDereferenced(
- &LocalAssemblerInterface::postTimestep, _local_assemblers,
- getDOFTable(process_id), x);
-}
-
-template <int DisplacementDim>
-void HydroMechanicsProcess<DisplacementDim>::postNonLinearSolverConcreteProcess(
- GlobalVector const& x, const double t, const int process_id)
-{
- if (!hasMechanicalProcess(process_id))
- {
- return;
- }
-
- DBUG("PostNonLinearSolver HydroMechanicsProcess.");
- // Calculate strain, stress or other internal variables of mechanics.
- GlobalExecutor::executeMemberOnDereferenced(
- &LocalAssemblerInterface::postNonLinearSolver, _local_assemblers,
- getDOFTable(process_id), x, t, _use_monolithic_scheme);
-}
-
-template <int DisplacementDim>
-std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
-HydroMechanicsProcess<DisplacementDim>::getDOFTableForExtrapolatorData() const
-{
- const bool manage_storage = false;
- return std::make_tuple(_local_to_global_index_map_single_component.get(),
- manage_storage);
-}
-
-template <int DisplacementDim>
-NumLib::LocalToGlobalIndexMap const&
-HydroMechanicsProcess<DisplacementDim>::getDOFTable(const int process_id) const
-{
- if (hasMechanicalProcess(process_id))
- {
- return *_local_to_global_index_map;
- }
-
- // For the equation of pressure
- return *_local_to_global_index_map_with_base_nodes;
-}
-
-} // namespace HydroMechanics
-} // namespace ProcessLib
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
index ac2223bae9e..7e7f6611832 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
@@ -8,12 +8,377 @@
*/
#include "HydroMechanicsProcess.h"
-#include "HydroMechanicsProcess-impl.h"
+
+#include <cassert>
+
+#include "MeshLib/Elements/Utils.h"
+#include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "ProcessLib/HydroMechanics/CreateLocalAssemblers.h"
+#include "ProcessLib/Process.h"
+
+#include "HydroMechanicsFEM.h"
+#include "HydroMechanicsProcessData.h"
namespace ProcessLib
{
namespace HydroMechanics
{
+template <int DisplacementDim>
+HydroMechanicsProcess<DisplacementDim>::HydroMechanicsProcess(
+ 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::vector<std::reference_wrapper<ProcessVariable>>>&&
+ process_variables,
+ HydroMechanicsProcessData<DisplacementDim>&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller,
+ bool const use_monolithic_scheme)
+ : Process(mesh, std::move(jacobian_assembler), parameters,
+ integration_order, std::move(process_variables),
+ std::move(secondary_variables), std::move(named_function_caller),
+ use_monolithic_scheme),
+ _process_data(std::move(process_data))
+{
+}
+
+template <int DisplacementDim>
+bool HydroMechanicsProcess<DisplacementDim>::isLinear() const
+{
+ return false;
+}
+
+template <int DisplacementDim>
+MathLib::MatrixSpecifications
+HydroMechanicsProcess<DisplacementDim>::getMatrixSpecifications(
+ const int process_id) const
+{
+ // For the monolithic scheme or the M process (deformation) in the staggered
+ // scheme.
+ if (_use_monolithic_scheme || process_id == 1)
+ {
+ auto const& l = *_local_to_global_index_map;
+ return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+ &l.getGhostIndices(), &this->_sparsity_pattern};
+ }
+
+ // For staggered scheme and H process (pressure).
+ auto const& l = *_local_to_global_index_map_with_base_nodes;
+ return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+ &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
+}
+
+template <int DisplacementDim>
+void HydroMechanicsProcess<DisplacementDim>::constructDofTable()
+{
+ // Create single component dof in every of the mesh's nodes.
+ _mesh_subset_all_nodes =
+ std::make_unique<MeshLib::MeshSubset>(_mesh, &_mesh.getNodes());
+ // Create single component dof in the mesh's base nodes.
+ _base_nodes = MeshLib::getBaseNodes(_mesh.getElements());
+ _mesh_subset_base_nodes =
+ std::make_unique<MeshLib::MeshSubset>(_mesh, &_base_nodes);
+
+ // TODO move the two data members somewhere else.
+ // for extrapolation of secondary variables of stress or strain
+ 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 =
+ std::make_unique<NumLib::LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets_single_component),
+ // by location order is needed for output
+ NumLib::ComponentOrder::BY_LOCATION);
+
+ if (_use_monolithic_scheme)
+ {
+ // For pressure, which is the first
+ std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
+ all_mesh_subsets.emplace_back(_mesh_subset_base_nodes.get());
+
+ // For displacement.
+ const int monolithic_process_id = 0;
+ std::generate_n(
+ std::back_inserter(all_mesh_subsets),
+ getProcessVariables(monolithic_process_id)[1]
+ .get()
+ .getNumberOfComponents(),
+ [&]() {
+ return MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()};
+ });
+
+ std::vector<int> const vec_n_components{1, DisplacementDim};
+ _local_to_global_index_map =
+ std::make_unique<NumLib::LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets), vec_n_components,
+ NumLib::ComponentOrder::BY_LOCATION);
+ assert(_local_to_global_index_map);
+ }
+ else
+ {
+ // For displacement equation.
+ const int process_id = 1;
+ std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
+ std::generate_n(
+ std::back_inserter(all_mesh_subsets),
+ getProcessVariables(process_id)[0].get().getNumberOfComponents(),
+ [&]() {
+ return MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()};
+ });
+
+ std::vector<int> const vec_n_components{DisplacementDim};
+ _local_to_global_index_map =
+ std::make_unique<NumLib::LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets), vec_n_components,
+ NumLib::ComponentOrder::BY_LOCATION);
+
+ // For pressure equation.
+ // Collect the mesh subsets with base nodes in a vector.
+ std::vector<MeshLib::MeshSubsets> all_mesh_subsets_base_nodes;
+ all_mesh_subsets_base_nodes.emplace_back(_mesh_subset_base_nodes.get());
+ _local_to_global_index_map_with_base_nodes =
+ std::make_unique<NumLib::LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets_base_nodes),
+ // by location order is needed for output
+ NumLib::ComponentOrder::BY_LOCATION);
+
+ _sparsity_pattern_with_linear_element = NumLib::computeSparsityPattern(
+ *_local_to_global_index_map_with_base_nodes, _mesh);
+
+ assert(_local_to_global_index_map);
+ assert(_local_to_global_index_map_with_base_nodes);
+ }
+}
+
+template <int DisplacementDim>
+void HydroMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+{
+ const int mechanical_process_id = _use_monolithic_scheme ? 0 : 1;
+ const int deformation_variable_id = _use_monolithic_scheme ? 1 : 0;
+ ProcessLib::HydroMechanics::createLocalAssemblers<
+ DisplacementDim, HydroMechanicsLocalAssembler>(
+ mesh.getDimension(), mesh.getElements(), dof_table,
+ // use displacement process variable to set shape function order
+ getProcessVariables(mechanical_process_id)[deformation_variable_id]
+ .get()
+ .getShapeFunctionOrder(),
+ _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
+ _process_data);
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_xx",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigmaXX));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_yy",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigmaYY));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_zz",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigmaZZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_xy",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigmaXY));
+
+ if (DisplacementDim == 3)
+ {
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_xz",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigmaXZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "sigma_yz",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigmaYZ));
+ }
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_xx",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtEpsilonXX));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_yy",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtEpsilonYY));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_zz",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtEpsilonZZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_xy",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtEpsilonXY));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "velocity",
+ makeExtrapolator(mesh.getDimension(), getExtrapolator(),
+ _local_assemblers,
+ &LocalAssemblerInterface::getIntPtDarcyVelocity));
+}
+
+template <int DisplacementDim>
+void HydroMechanicsProcess<DisplacementDim>::initializeBoundaryConditions()
+{
+ if (_use_monolithic_scheme)
+ {
+ const int process_id_of_hydromechancs = 0;
+ initializeProcessBoundaryConditionsAndSourceTerms(
+ *_local_to_global_index_map, process_id_of_hydromechancs);
+ return;
+ }
+
+ // Staggered scheme:
+ // for the equations of pressure
+ const int hydraulic_process_id = 0;
+ initializeProcessBoundaryConditionsAndSourceTerms(
+ *_local_to_global_index_map_with_base_nodes, hydraulic_process_id);
+
+ // for the equations of deformation.
+ const int mechanical_process_id = 1;
+ initializeProcessBoundaryConditionsAndSourceTerms(
+ *_local_to_global_index_map, mechanical_process_id);
+}
+
+template <int DisplacementDim>
+void HydroMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
+ const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b)
+{
+ DBUG("Assemble the equations for HydroMechanics");
+
+ // Note: This assembly function is for the Picard nonlinear solver. Since
+ // only the Newton-Raphson method is employed to simulate coupled HM
+ // processes in this class, this function is actually not used so far.
+
+ std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+ dof_table = {std::ref(*_local_to_global_index_map)};
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
+ dof_table, t, x, M, K, b, _coupled_solutions);
+}
+
+template <int DisplacementDim>
+void HydroMechanicsProcess<DisplacementDim>::
+ 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)
+{
+ std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+ dof_tables;
+ // For the monolithic scheme
+ if (_use_monolithic_scheme)
+ {
+ DBUG(
+ "Assemble the Jacobian of HydroMechanics for the monolithic"
+ " scheme.");
+ dof_tables.emplace_back(*_local_to_global_index_map);
+ }
+ else
+ {
+ // For the staggered scheme
+ if (_coupled_solutions->process_id == 0)
+ {
+ DBUG(
+ "Assemble the Jacobian equations of liquid fluid process in "
+ "HydroMechanics for the staggered scheme.");
+ }
+ else
+ {
+ DBUG(
+ "Assemble the Jacobian equations of mechanical process in "
+ "HydroMechanics for the staggered scheme.");
+ }
+ dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
+ dof_tables.emplace_back(*_local_to_global_index_map);
+ }
+
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+ _local_assemblers, dof_tables, t, x, xdot, dxdot_dx, dx_dx, M, K, b,
+ Jac, _coupled_solutions);
+}
+
+template <int DisplacementDim>
+void HydroMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
+ GlobalVector const& x, double const t, double const dt,
+ const int process_id)
+{
+ DBUG("PreTimestep HydroMechanicsProcess.");
+
+ _process_data.dt = dt;
+ _process_data.t = t;
+
+ if (hasMechanicalProcess(process_id))
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerInterface::preTimestep, _local_assemblers,
+ *_local_to_global_index_map, x, t, dt);
+}
+
+template <int DisplacementDim>
+void HydroMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
+ GlobalVector const& x, const int process_id)
+{
+ DBUG("PostTimestep HydroMechanicsProcess.");
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerInterface::postTimestep, _local_assemblers,
+ getDOFTable(process_id), x);
+}
+
+template <int DisplacementDim>
+void HydroMechanicsProcess<DisplacementDim>::postNonLinearSolverConcreteProcess(
+ GlobalVector const& x, const double t, const int process_id)
+{
+ if (!hasMechanicalProcess(process_id))
+ {
+ return;
+ }
+
+ DBUG("PostNonLinearSolver HydroMechanicsProcess.");
+ // Calculate strain, stress or other internal variables of mechanics.
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerInterface::postNonLinearSolver, _local_assemblers,
+ getDOFTable(process_id), x, t, _use_monolithic_scheme);
+}
+
+template <int DisplacementDim>
+std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
+HydroMechanicsProcess<DisplacementDim>::getDOFTableForExtrapolatorData() const
+{
+ const bool manage_storage = false;
+ return std::make_tuple(_local_to_global_index_map_single_component.get(),
+ manage_storage);
+}
+
+template <int DisplacementDim>
+NumLib::LocalToGlobalIndexMap const&
+HydroMechanicsProcess<DisplacementDim>::getDOFTable(const int process_id) const
+{
+ if (hasMechanicalProcess(process_id))
+ {
+ return *_local_to_global_index_map;
+ }
+
+ // For the equation of pressure
+ return *_local_to_global_index_map_with_base_nodes;
+}
+
template class HydroMechanicsProcess<2>;
template class HydroMechanicsProcess<3>;
--
GitLab