From 58d5c91dd7785c8a148f5adb426984ddb0fe3894 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Wed, 11 Oct 2017 17:41:46 +0200
Subject: [PATCH] [Coupling] Modified the process definitions in Process/HT for
the assembly for both of the monolithic scheme and the staggered scheme.
---
ProcessLib/HT/HTFEM.h | 212 ++-------------------
ProcessLib/HT/HTLocalAssemblerInterface.h | 53 ++++++
ProcessLib/HT/HTProcess.cpp | 22 ++-
ProcessLib/HT/MonolithicHTFEM.h | 218 ++++++++++++++++++++++
ProcessLib/HT/StaggeredHTFEM-impl.h | 198 ++++++++++++++++++++
ProcessLib/HT/StaggeredHTFEM.h | 87 +++++++++
6 files changed, 586 insertions(+), 204 deletions(-)
create mode 100644 ProcessLib/HT/HTLocalAssemblerInterface.h
create mode 100644 ProcessLib/HT/MonolithicHTFEM.h
create mode 100644 ProcessLib/HT/StaggeredHTFEM-impl.h
create mode 100644 ProcessLib/HT/StaggeredHTFEM.h
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 0cfc29612ca..21d082bb794 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -14,66 +14,27 @@
#include "HTMaterialProperties.h"
-#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+
#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
-#include "NumLib/Function/Interpolation.h"
-#include "ProcessLib/LocalAssemblerInterface.h"
#include "ProcessLib/Parameter/Parameter.h"
#include "ProcessLib/Utils/InitShapeMatrices.h"
+#include "HTLocalAssemblerInterface.h"
+
namespace ProcessLib
{
namespace HT
{
-template < typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
-struct IntegrationPointData final
-{
- IntegrationPointData(NodalRowVectorType N_,
- GlobalDimNodalMatrixType dNdx_,
- double const& integration_weight_)
- : N(std::move(N_)),
- dNdx(std::move(dNdx_)),
- integration_weight(integration_weight_)
- {}
-
- NodalRowVectorType const N;
- GlobalDimNodalMatrixType const dNdx;
- double const integration_weight;
-
- EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
-};
-
-const unsigned NUM_NODAL_DOF = 2;
-
-class HTLocalAssemblerInterface
- : public ProcessLib::LocalAssemblerInterface,
- public NumLib::ExtrapolatableElement
-{
-public:
- virtual std::vector<double> const& getIntPtDarcyVelocity(
- const double /*t*/,
- GlobalVector const& /*current_solution*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
- std::vector<double>& /*cache*/) const = 0;
-};
-
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
-class LocalAssemblerData : public HTLocalAssemblerInterface
+class HTFEM : public HTLocalAssemblerInterface
{
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
- using LocalMatrixType = typename ShapeMatricesType::template MatrixType<
- NUM_NODAL_DOF * ShapeFunction::NPOINTS,
- NUM_NODAL_DOF * ShapeFunction::NPOINTS>;
- using LocalVectorType =
- typename ShapeMatricesType::template VectorType<NUM_NODAL_DOF *
- ShapeFunction::NPOINTS>;
-
using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType;
@@ -83,18 +44,19 @@ class LocalAssemblerData : public HTLocalAssemblerInterface
using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
public:
- LocalAssemblerData(MeshLib::Element const& element,
- std::size_t const local_matrix_size,
- bool is_axially_symmetric,
- unsigned const integration_order,
- HTMaterialProperties const& process_data)
+ HTFEM(MeshLib::Element const& element,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ HTMaterialProperties const& process_data,
+ const unsigned dof_per_node)
: _element(element),
_process_data(process_data),
_integration_method(integration_order)
{
// This assertion is valid only if all nodal d.o.f. use the same shape
// matrices.
- assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+ assert(local_matrix_size == ShapeFunction::NPOINTS * dof_per_node);
(void)local_matrix_size;
unsigned const n_integration_points =
@@ -116,156 +78,6 @@ public:
}
}
- 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_b_data) override
- {
- auto const local_matrix_size = local_x.size();
- // This assertion is valid only if all nodal d.o.f. use the same shape
- // matrices.
- assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
-
- auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
- local_M_data, local_matrix_size, local_matrix_size);
- auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
- local_K_data, local_matrix_size, local_matrix_size);
- auto local_b = MathLib::createZeroedVector<LocalVectorType>(
- local_b_data, local_matrix_size);
-
- auto const num_nodes = ShapeFunction::NPOINTS;
-
- auto Ktt = local_K.template block<num_nodes, num_nodes>(0, 0);
- auto Mtt = local_M.template block<num_nodes, num_nodes>(0, 0);
- auto Kpp =
- local_K.template block<num_nodes, num_nodes>(num_nodes, num_nodes);
- auto Mpp =
- local_M.template block<num_nodes, num_nodes>(num_nodes, num_nodes);
- auto Bp = local_b.template block<num_nodes, 1>(num_nodes, 0);
-
- SpatialPosition pos;
- pos.setElementID(_element.getID());
-
- auto p_nodal_values =
- Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
-
- auto const& b = _process_data.specific_body_force;
-
- GlobalDimMatrixType const& I(
- GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
-
- MaterialLib::Fluid::FluidProperty::ArrayType vars;
-
- unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
-
- for (std::size_t ip(0); ip < n_integration_points; ip++)
- {
- pos.setIntegrationPoint(ip);
-
- auto const fluid_reference_density =
- _process_data.fluid_reference_density(t, pos)[0];
-
- auto const density_solid = _process_data.density_solid(t, pos)[0];
- // \todo the argument to getValue() has to be changed for non
- // constant storage model
- auto const specific_storage =
- _process_data.porous_media_properties.getSpecificStorage(t, pos)
- .getValue(0.0);
-
- auto const thermal_conductivity_solid =
- _process_data.thermal_conductivity_solid(t, pos)[0];
- auto const thermal_conductivity_fluid =
- _process_data.thermal_conductivity_fluid(t, pos)[0];
-
- auto const& ip_data = _ip_data[ip];
- auto const& N = ip_data.N;
- auto const& dNdx = ip_data.dNdx;
- auto const& w = ip_data.integration_weight;
-
- double T_int_pt = 0.0;
- double p_int_pt = 0.0;
- // Order matters: First T, then P!
- NumLib::shapeFunctionInterpolate(local_x, N, T_int_pt, p_int_pt);
-
- // \todo the first argument has to be changed for non constant
- // porosity model
- auto const porosity =
- _process_data.porous_media_properties.getPorosity(t, pos)
- .getValue(t, pos, 0.0, T_int_pt);
-
- auto const intrinsic_permeability =
- _process_data.porous_media_properties.getIntrinsicPermeability(
- t, pos).getValue(t, pos, 0.0, T_int_pt);
-
- double const thermal_conductivity =
- thermal_conductivity_solid * (1 - porosity) +
- thermal_conductivity_fluid * porosity;
-
- auto const specific_heat_capacity_solid =
- _process_data.specific_heat_capacity_solid(t, pos)[0];
-
- vars[static_cast<int>(
- MaterialLib::Fluid::PropertyVariableType::T)] = T_int_pt;
- vars[static_cast<int>(
- MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
- auto const specific_heat_capacity_fluid =
- _process_data.fluid_properties->getValue(
- MaterialLib::Fluid::FluidPropertyType::HeatCapacity, vars);
-
- auto const thermal_dispersivity_longitudinal =
- _process_data.thermal_dispersivity_longitudinal(t, pos)[0];
- auto const thermal_dispersivity_transversal =
- _process_data.thermal_dispersivity_transversal(t, pos)[0];
-
- // Use the fluid density model to compute the density
- auto const density = _process_data.fluid_properties->getValue(
- MaterialLib::Fluid::FluidPropertyType::Density, vars);
-
- // Use the viscosity model to compute the viscosity
- auto const viscosity = _process_data.fluid_properties->getValue(
- MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
- GlobalDimMatrixType K_over_mu = intrinsic_permeability / viscosity;
-
- GlobalDimVectorType const velocity =
- _process_data.has_gravity
- ? GlobalDimVectorType(-K_over_mu *
- (dNdx * p_nodal_values - density * b))
- : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
-
-
- double const velocity_magnitude = velocity.norm();
- GlobalDimMatrixType const thermal_dispersivity =
- fluid_reference_density * specific_heat_capacity_fluid *
- (thermal_dispersivity_transversal * velocity_magnitude *
- I +
- (thermal_dispersivity_longitudinal -
- thermal_dispersivity_transversal) /
- velocity_magnitude * velocity * velocity.transpose());
-
- GlobalDimMatrixType const hydrodynamic_thermodispersion =
- thermal_conductivity * I + thermal_dispersivity;
-
- double const heat_capacity =
- density_solid * specific_heat_capacity_solid * (1 - porosity) +
- fluid_reference_density * specific_heat_capacity_fluid * porosity;
-
- // matrix assembly
- Ktt.noalias() +=
- (dNdx.transpose() * hydrodynamic_thermodispersion * dNdx +
- N.transpose() * velocity.transpose() * dNdx *
- fluid_reference_density * specific_heat_capacity_fluid) *
- w;
- Kpp.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
- Mtt.noalias() += w * N.transpose() * heat_capacity * N;
- Mpp.noalias() += w * N.transpose() * specific_storage * N;
- if (_process_data.has_gravity)
- Bp += w * density * dNdx.transpose() * K_over_mu * b;
- /* with Oberbeck-Boussing assumption density difference only exists
- * in buoyancy effects */
- }
- }
-
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
@@ -340,7 +152,7 @@ public:
return cache;
}
-private:
+protected:
MeshLib::Element const& _element;
HTMaterialProperties const& _process_data;
diff --git a/ProcessLib/HT/HTLocalAssemblerInterface.h b/ProcessLib/HT/HTLocalAssemblerInterface.h
new file mode 100644
index 00000000000..82408a4d99f
--- /dev/null
+++ b/ProcessLib/HT/HTLocalAssemblerInterface.h
@@ -0,0 +1,53 @@
+/**
+ * \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
+ *
+ * \file HTLocalAssemblerInterface.h
+ * Created on October 11, 2017, 1:35 PM
+ */
+
+#pragma once
+
+#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+
+namespace ProcessLib
+{
+namespace HT
+{
+template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
+struct IntegrationPointData final
+{
+ IntegrationPointData(NodalRowVectorType N_,
+ GlobalDimNodalMatrixType dNdx_,
+ double const& integration_weight_)
+ : N(std::move(N_)),
+ dNdx(std::move(dNdx_)),
+ integration_weight(integration_weight_)
+ {
+ }
+
+ NodalRowVectorType const N;
+ GlobalDimNodalMatrixType const dNdx;
+ double const integration_weight;
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+class HTLocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
+ public NumLib::ExtrapolatableElement
+{
+public:
+ virtual std::vector<double> const& getIntPtDarcyVelocity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& /*cache*/) const = 0;
+};
+
+} // namespace HT
+} // namespace ProcessLib
\ No newline at end of file
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index bcc58332edc..3b95d1e54b5 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -13,6 +13,9 @@
#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "MonolithicHTFEM.h"
+#include "StaggeredHTFEM.h"
+
namespace ProcessLib
{
namespace HT
@@ -39,10 +42,21 @@ void HTProcess::initializeConcreteProcess(
unsigned const integration_order)
{
ProcessLib::ProcessVariable const& pv = getProcessVariables()[0];
- ProcessLib::createLocalAssemblers<LocalAssemblerData>(
- mesh.getDimension(), mesh.getElements(), dof_table,
- pv.getShapeFunctionOrder(), _local_assemblers,
- mesh.isAxiallySymmetric(), integration_order, *_material_properties);
+
+ if (_is_monolithic_scheme)
+ {
+ ProcessLib::createLocalAssemblers<MonolithicHTFEM>(
+ mesh.getDimension(), mesh.getElements(), dof_table,
+ pv.getShapeFunctionOrder(), _local_assemblers,
+ mesh.isAxiallySymmetric(), integration_order, *_material_properties);
+ }
+ else
+ {
+ ProcessLib::createLocalAssemblers<StaggeredHTFEM>(
+ mesh.getDimension(), mesh.getElements(), dof_table,
+ pv.getShapeFunctionOrder(), _local_assemblers,
+ mesh.isAxiallySymmetric(), integration_order, *_material_properties);
+ }
_secondary_variables.addSecondaryVariable(
"darcy_velocity",
diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h
new file mode 100644
index 00000000000..6efaba1a7e1
--- /dev/null
+++ b/ProcessLib/HT/MonolithicHTFEM.h
@@ -0,0 +1,218 @@
+/**
+ * \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
+ *
+ * \file MonolithicHTFEM.h
+ * Created on October 11, 2017, 2:33 PM
+ */
+
+#pragma once
+
+#include <Eigen/Dense>
+#include <vector>
+
+#include "HTProcessData.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "HTFEM.h"
+
+namespace ProcessLib
+{
+namespace HT
+{
+const unsigned NUM_NODAL_DOF = 2;
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class MonolithicHTFEM
+ : public HTFEM<ShapeFunction, IntegrationMethod, GlobalDim>
+{
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+ using LocalMatrixType = typename ShapeMatricesType::template MatrixType<
+ NUM_NODAL_DOF * ShapeFunction::NPOINTS,
+ NUM_NODAL_DOF * ShapeFunction::NPOINTS>;
+ using LocalVectorType =
+ typename ShapeMatricesType::template VectorType<NUM_NODAL_DOF *
+ ShapeFunction::NPOINTS>;
+
+ using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+ using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType;
+
+ using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+ using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
+
+public:
+ MonolithicHTFEM(MeshLib::Element const& element,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ HTMaterialProperties const& process_data)
+ : HTFEM<ShapeFunction, IntegrationMethod, GlobalDim>(
+ element, local_matrix_size, is_axially_symmetric,
+ integration_order, process_data, NUM_NODAL_DOF)
+ {
+ }
+
+ 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_b_data) override
+ {
+ auto const local_matrix_size = local_x.size();
+ // This assertion is valid only if all nodal d.o.f. use the same shape
+ // matrices.
+ assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+
+ auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
+ local_M_data, local_matrix_size, local_matrix_size);
+ auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
+ local_K_data, local_matrix_size, local_matrix_size);
+ auto local_b = MathLib::createZeroedVector<LocalVectorType>(
+ local_b_data, local_matrix_size);
+
+ auto const num_nodes = ShapeFunction::NPOINTS;
+
+ auto Ktt = local_K.template block<num_nodes, num_nodes>(0, 0);
+ auto Mtt = local_M.template block<num_nodes, num_nodes>(0, 0);
+ auto Kpp =
+ local_K.template block<num_nodes, num_nodes>(num_nodes, num_nodes);
+ auto Mpp =
+ local_M.template block<num_nodes, num_nodes>(num_nodes, num_nodes);
+ auto Bp = local_b.template block<num_nodes, 1>(num_nodes, 0);
+
+ SpatialPosition pos;
+ pos.setElementID(this->_element.getID());
+
+ auto p_nodal_values =
+ Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
+
+ auto const& process_data = this->_process_data;
+
+ auto const& b = process_data.specific_body_force;
+
+ GlobalDimMatrixType const& I(
+ GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
+
+ MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+ unsigned const n_integration_points =
+ this->_integration_method.getNumberOfPoints();
+
+ for (std::size_t ip(0); ip < n_integration_points; ip++)
+ {
+ pos.setIntegrationPoint(ip);
+
+ auto const fluid_reference_density =
+ process_data.fluid_reference_density(t, pos)[0];
+
+ auto const density_solid = process_data.density_solid(t, pos)[0];
+ // \todo the argument to getValue() has to be changed for non
+ // constant storage model
+ auto const specific_storage =
+ process_data.porous_media_properties.getSpecificStorage(t, pos)
+ .getValue(0.0);
+
+ auto const thermal_conductivity_solid =
+ process_data.thermal_conductivity_solid(t, pos)[0];
+ auto const thermal_conductivity_fluid =
+ process_data.thermal_conductivity_fluid(t, pos)[0];
+
+ auto const& ip_data = this->_ip_data[ip];
+ auto const& N = ip_data.N;
+ auto const& dNdx = ip_data.dNdx;
+ auto const& w = ip_data.integration_weight;
+
+ double T_int_pt = 0.0;
+ double p_int_pt = 0.0;
+ // Order matters: First T, then P!
+ NumLib::shapeFunctionInterpolate(local_x, N, T_int_pt, p_int_pt);
+
+ // \todo the first argument has to be changed for non constant
+ // porosity model
+ auto const porosity =
+ process_data.porous_media_properties.getPorosity(t, pos)
+ .getValue(t, pos, 0.0, T_int_pt);
+ auto const intrinsic_permeability =
+ process_data.porous_media_properties.getIntrinsicPermeability(
+ t, pos).getValue(t, pos, 0.0, T_int_pt);
+
+ double const thermal_conductivity =
+ thermal_conductivity_solid * (1 - porosity) +
+ thermal_conductivity_fluid * porosity;
+
+ auto const specific_heat_capacity_solid =
+ process_data.specific_heat_capacity_solid(t, pos)[0];
+
+ vars[static_cast<int>(
+ MaterialLib::Fluid::PropertyVariableType::T)] = T_int_pt;
+ vars[static_cast<int>(
+ MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
+ auto const specific_heat_capacity_fluid =
+ process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::HeatCapacity, vars);
+
+ auto const thermal_dispersivity_longitudinal =
+ process_data.thermal_dispersivity_longitudinal(t, pos)[0];
+ auto const thermal_dispersivity_transversal =
+ process_data.thermal_dispersivity_transversal(t, pos)[0];
+
+ // Use the fluid density model to compute the density
+ auto const density = process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars);
+
+ // Use the viscosity model to compute the viscosity
+ auto const viscosity = process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+ GlobalDimMatrixType K_over_mu = intrinsic_permeability / viscosity;
+
+ GlobalDimVectorType const velocity =
+ process_data.has_gravity
+ ? GlobalDimVectorType(-K_over_mu *
+ (dNdx * p_nodal_values - density * b))
+ : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
+
+ double const velocity_magnitude = velocity.norm();
+ GlobalDimMatrixType const thermal_dispersivity =
+ fluid_reference_density * specific_heat_capacity_fluid *
+ (thermal_dispersivity_transversal * velocity_magnitude * I +
+ (thermal_dispersivity_longitudinal -
+ thermal_dispersivity_transversal) /
+ velocity_magnitude * velocity * velocity.transpose());
+
+ GlobalDimMatrixType const hydrodynamic_thermodispersion =
+ thermal_conductivity * I + thermal_dispersivity;
+
+ double const heat_capacity =
+ density_solid * specific_heat_capacity_solid * (1 - porosity) +
+ fluid_reference_density * specific_heat_capacity_fluid *
+ porosity;
+
+ // matrix assembly
+ Ktt.noalias() +=
+ (dNdx.transpose() * hydrodynamic_thermodispersion * dNdx +
+ N.transpose() * velocity.transpose() * dNdx *
+ fluid_reference_density * specific_heat_capacity_fluid) *
+ w;
+ Kpp.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
+ Mtt.noalias() += w * N.transpose() * heat_capacity * N;
+ Mpp.noalias() += w * N.transpose() * specific_storage * N;
+ if (process_data.has_gravity)
+ Bp += w * density * dNdx.transpose() * K_over_mu * b;
+ /* with Oberbeck-Boussing assumption density difference only exists
+ * in buoyancy effects */
+ }
+ }
+};
+
+} // namespace HT
+} // namespace ProcessLib
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
new file mode 100644
index 00000000000..1d6508049c4
--- /dev/null
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -0,0 +1,198 @@
+/**
+ * \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
+ *
+ * \file StaggeredHTFEM-impl.h
+ * Created on October 13, 2017, 3:52 PM
+ */
+
+#pragma once
+
+#include "StaggeredHTFEM.h"
+
+#include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
+
+namespace ProcessLib
+{
+namespace HT
+{
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
+ assembleWithCoupledTerm(double const t,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data,
+ LocalCoupledSolutions const& coupled_term)
+{
+ if (coupled_term.variable_id == 0)
+ {
+ assembleHydraulicEquation(t, local_M_data, local_K_data, local_b_data,
+ coupled_term);
+ return;
+ }
+
+ assembleHeatTransportEquation(t, local_M_data, local_K_data, local_b_data,
+ coupled_term);
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
+ assembleHydraulicEquation(double const t, std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data,
+ LocalCoupledSolutions const& coupled_term)
+{
+ auto const& local_p = coupled_term.local_coupled_xs[0];
+ auto const& local_T = coupled_term.local_coupled_xs[1];
+ auto const& local_T1 = coupled_term.local_coupled_xs0[1];
+ const double dt = coupled_term.dt;
+
+ auto const local_matrix_size = local_x.size();
+ // This assertion is valid only if all nodal d.o.f. use the same shape
+ // matrices.
+ assert(local_matrix_size == ShapeFunction::NPOINTS);
+
+ auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
+ local_M_data, local_matrix_size, local_matrix_size);
+ auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
+ local_K_data, local_matrix_size, local_matrix_size);
+ auto local_b = MathLib::createZeroedVector<LocalVectorType>(
+ local_b_data, local_matrix_size);
+
+ SpatialPosition pos;
+ pos.setElementID(this->_element.getID());
+
+ auto const& process_data = this->_process_data;
+
+ auto const& b = process_data.specific_body_force;
+
+ GlobalDimMatrixType const& I(
+ GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
+
+ MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+ unsigned const n_integration_points =
+ this->_integration_method.getNumberOfPoints();
+
+ for (std::size_t ip(0); ip < n_integration_points; ip++)
+ {
+ pos.setIntegrationPoint(ip);
+
+ auto const fluid_reference_density =
+ process_data.fluid_reference_density(t, pos)[0];
+
+ auto const density_solid = process_data.density_solid(t, pos)[0];
+ // \todo the argument to getValue() has to be changed for non
+ // constant storage model
+ auto const specific_storage =
+ process_data.porous_media_properties.getSpecificStorage(t, pos)
+ .getValue(0.0);
+
+ auto const thermal_conductivity_solid =
+ process_data.thermal_conductivity_solid(t, pos)[0];
+ auto const thermal_conductivity_fluid =
+ process_data.thermal_conductivity_fluid(t, pos)[0];
+
+ auto const& ip_data = this->_ip_data[ip];
+ auto const& N = ip_data.N;
+ auto const& dNdx = ip_data.dNdx;
+ auto const& w = ip_data.integration_weight;
+
+ double T_int_pt = 0.0;
+ double p_int_pt = 0.0;
+ // Order matters: First T, then P!
+ NumLib::shapeFunctionInterpolate(local_x, N, T_int_pt, p_int_pt);
+
+ // \todo the first argument has to be changed for non constant
+ // porosity model
+ auto const porosity =
+ process_data.porous_media_properties.getPorosity(t, pos).getValue(
+ t, pos, 0.0, T_int_pt);
+ auto const intrinsic_permeability =
+ process_data.porous_media_properties.getIntrinsicPermeability(
+ t, pos).getValue(t, pos, 0.0, T_int_pt);
+
+ double const thermal_conductivity =
+ thermal_conductivity_solid * (1 - porosity) +
+ thermal_conductivity_fluid * porosity;
+
+ auto const specific_heat_capacity_solid =
+ process_data.specific_heat_capacity_solid(t, pos)[0];
+
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] =
+ T_int_pt;
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] =
+ p_int_pt;
+ auto const specific_heat_capacity_fluid =
+ process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::HeatCapacity, vars);
+
+ auto const thermal_dispersivity_longitudinal =
+ process_data.thermal_dispersivity_longitudinal(t, pos)[0];
+ auto const thermal_dispersivity_transversal =
+ process_data.thermal_dispersivity_transversal(t, pos)[0];
+
+ // Use the fluid density model to compute the density
+ auto const density = process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars);
+
+ // Use the viscosity model to compute the viscosity
+ auto const viscosity = process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+ GlobalDimMatrixType K_over_mu = intrinsic_permeability / viscosity;
+
+ GlobalDimVectorType const velocity =
+ process_data.has_gravity
+ ? GlobalDimVectorType(-K_over_mu *
+ (dNdx * p_nodal_values - density * b))
+ : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
+
+ double const velocity_magnitude = velocity.norm();
+ GlobalDimMatrixType const thermal_dispersivity =
+ fluid_reference_density * specific_heat_capacity_fluid *
+ (thermal_dispersivity_transversal * velocity_magnitude * I +
+ (thermal_dispersivity_longitudinal -
+ thermal_dispersivity_transversal) /
+ velocity_magnitude * velocity * velocity.transpose());
+
+ GlobalDimMatrixType const hydrodynamic_thermodispersion =
+ thermal_conductivity * I + thermal_dispersivity;
+
+ double const heat_capacity =
+ density_solid * specific_heat_capacity_solid * (1 - porosity) +
+ fluid_reference_density * specific_heat_capacity_fluid * porosity;
+
+ // matrix assembly
+ Ktt.noalias() +=
+ (dNdx.transpose() * hydrodynamic_thermodispersion * dNdx +
+ N.transpose() * velocity.transpose() * dNdx *
+ fluid_reference_density * specific_heat_capacity_fluid) *
+ w;
+ Kpp.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
+ Mtt.noalias() += w * N.transpose() * heat_capacity * N;
+ Mpp.noalias() += w * N.transpose() * specific_storage * N;
+ if (process_data.has_gravity)
+ Bp += w * density * dNdx.transpose() * K_over_mu * b;
+ /* with Oberbeck-Boussing assumption density difference only exists
+ * in buoyancy effects */
+ }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
+ assembleHeatTransportEquation(double const t,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data,
+ LocalCoupledSolutions const& coupled_term)
+{
+}
+
+} // namespace HT
+} // namespace ProcessLib
diff --git a/ProcessLib/HT/StaggeredHTFEM.h b/ProcessLib/HT/StaggeredHTFEM.h
new file mode 100644
index 00000000000..b6575cc5cb2
--- /dev/null
+++ b/ProcessLib/HT/StaggeredHTFEM.h
@@ -0,0 +1,87 @@
+/**
+ * \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
+ *
+ * \file StaggeredHTFEM.h
+ * Created on October 11, 2017, 1:42 PM
+ */
+
+#pragma once
+
+#include <Eigen/Dense>
+#include <vector>
+
+#include "HTProcessData.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "HTFEM.h"
+
+namespace ProcessLib
+{
+struct LocalCoupledSolutions;
+
+namespace HT
+{
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class StaggeredHTFEM : public HTFEM<ShapeFunction, IntegrationMethod, GlobalDim>
+{
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+ using LocalMatrixType =
+ typename ShapeMatricesType::template MatrixType<ShapeFunction::NPOINTS,
+ ShapeFunction::NPOINTS>;
+ using LocalVectorType =
+ typename ShapeMatricesType::template VectorType<ShapeFunction::NPOINTS>;
+
+ using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+ using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType;
+
+ using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+ using GlobalDimNodalMatrixType =
+ typename ShapeMatricesType::GlobalDimNodalMatrixType;
+ using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
+
+public:
+ StaggeredHTFEM(MeshLib::Element const& element,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ HTMaterialProperties const& process_data)
+ : HTFEM<ShapeFunction, IntegrationMethod, GlobalDim>(
+ element, local_matrix_size, is_axially_symmetric,
+ integration_order, process_data, 1)
+ {
+ }
+
+ void assembleWithCoupledTerm(
+ double const t, std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+ LocalCoupledSolutions const& coupled_term) override;
+
+private:
+ void assembleHydraulicEquation(double const t,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data,
+ LocalCoupledSolutions const& coupled_term);
+
+ void assembleHeatTransportEquation(
+ double const t, std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+ LocalCoupledSolutions const& coupled_term);
+};
+
+} // namespace HT
+} // namespace ProcessLib
+
+#include "StaggeredHTFEM-impl.h"
--
GitLab