diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 60db84517feb9622e7f63620837b465af936ebf9..72eb86bf8050626414cbb5b382942ad8dae1c9aa 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -41,7 +41,7 @@
#include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h"
#include "ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.h"
#include "ProcessLib/TES/CreateTESProcess.h"
-#include "ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.h"
+#include "ProcessLib/HT/CreateHTProcess.h"
namespace detail
{
@@ -352,9 +352,9 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
"dimension");
}
}
- else if (type == "DensityDrivenFlow")
+ else if (type == "HT")
{
- process = ProcessLib::DensityDrivenFlow::createDensityDrivenFlowProcess(
+ process = ProcessLib::HT::createHTProcess(
*_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters, integration_order,
process_config);
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index c7a4810635a3e6b45ce3c7323f388581536c50c8..592ba502580ea8d608bdb92de4a74602bc7cdccd 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -33,7 +33,7 @@ APPEND_SOURCE_FILES(SOURCES HeatConduction)
APPEND_SOURCE_FILES(SOURCES RichardsFlow)
-APPEND_SOURCE_FILES(SOURCES DensityDrivenFlow)
+APPEND_SOURCE_FILES(SOURCES HT)
add_subdirectory(CalculateSurfaceFlux)
APPEND_SOURCE_FILES(SOURCES CalculateSurfaceFlux)
diff --git a/ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.cpp b/ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.cpp
deleted file mode 100644
index 1c3e4ca8f50962eb81ef0afa9e6b28800e0629c5..0000000000000000000000000000000000000000
--- a/ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.cpp
+++ /dev/null
@@ -1,65 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include "CreateDensityDrivenFlowProcess.h"
-
-#include "ProcessLib/Utils/ParseSecondaryVariables.h"
-#include "ProcessLib/Utils/ProcessUtils.h"
-#include "DensityDrivenFlowProcess.h"
-#include "DensityDrivenFlowProcessData.h"
-
-namespace ProcessLib
-{
-namespace DensityDrivenFlow
-{
-std::unique_ptr<Process> createDensityDrivenFlowProcess(
- 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", "DensityDrivenFlow");
-
- DBUG("Create DensityDrivenFlowProcess.");
-
- // Process variable.
- auto process_variables =
- findProcessVariables(variables, config, {
- "temperature_variable" , "pressure_variable"}); // configure two Pcs
-
- // Thermal conductivity parameter.
- auto& thermal_conductivity = findParameter<double>(
- config,
- //! \ogs_file_param_special{process__DensityDrivenFlow__thermal_conductivity}
- "thermal_conductivity",parameters, 1);
-
- DBUG("Use \'%s\' as thermal conductivity parameter.",
- thermal_conductivity.name.c_str());
-
- DensityDrivenFlowProcessData process_data{thermal_conductivity};
-
- SecondaryVariableCollection secondary_variables;
-
- NumLib::NamedFunctionCaller named_function_caller(
- {"DensityDrivenFlow_temperature_pressure"});
-
- ProcessLib::parseSecondaryVariables(config, secondary_variables,
- named_function_caller);
-
- return std::unique_ptr<Process>{new DensityDrivenFlowProcess{
- 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)}};
-}
-
-} // namespace DensityDrivenFlow
-} // namespace ProcessLib
diff --git a/ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h b/ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h
deleted file mode 100644
index 9564c0b3d406a2bce21840cac569a6d1710488e1..0000000000000000000000000000000000000000
--- a/ProcessLib/DensityDrivenFlow/DensityDrivenFlowFEM.h
+++ /dev/null
@@ -1,338 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef PROCESS_LIB_DensityDrivenFlowPROCESS_FEM_H_
-#define PROCESS_LIB_DensityDrivenFlowPROCESS_FEM_H_
-
-#include <iostream>
-#include <Eigen/Dense>
-#include <vector>
-
-
-#include "DensityDrivenFlowMaterialModel.h"
-#include "DensityDrivenFlowProcessData.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/LocalAssemblerInterface.h"
-#include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Parameter/Parameter.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-using namespace Eigen;
-
-namespace ProcessLib
-{
-
-namespace DensityDrivenFlow
-{
-
-const unsigned NUM_NODAL_DOF = 2;
-
-class DensityDrivenFlowLocalAssemblerInterface
- : public ProcessLib::LocalAssemblerInterface
- , public NumLib::ExtrapolatableElement
-{
-public:
- virtual std::vector<double> const& getIntPtDarcyVelocityX(
- std::vector<double>& /*cache*/) const = 0;
-
- virtual std::vector<double> const& getIntPtDarcyVelocityY(
- std::vector<double>& /*cache*/) const = 0;
-
- virtual std::vector<double> const& getIntPtDarcyVelocityZ(
-std::vector<double>& /*cache*/) const = 0;
-};
-
-
-template <typename ShapeFunction,
- typename IntegrationMethod,
- unsigned GlobalDim>
-class LocalAssemblerData
- : public DensityDrivenFlowLocalAssemblerInterface
-{
- using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
- using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
-
- using LocalAssemblerTraits = ProcessLib::LocalAssemblerTraits<
- ShapeMatricesType, ShapeFunction::NPOINTS,
- NUM_NODAL_DOF, /* number of pcs vars */ GlobalDim>;
- // Two Pcs variable T and P
- using NodalMatrixType = typename LocalAssemblerTraits::LocalMatrix;
- using NodalVectorType = typename LocalAssemblerTraits::LocalVector;
- using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
-
-public:
- /// The thermal_conductivity factor is directly integrated into the local
- /// element matrix.
- LocalAssemblerData(MeshLib::Element const& element,
- std::size_t const local_matrix_size,
- bool is_axially_symmetric,
- unsigned const integration_order,
- DensityDrivenFlowProcessData const& process_data)
- : _element(element)
- , _process_data(process_data)
- , _integration_method(integration_order)
- , _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
- IntegrationMethod, GlobalDim>(
- element, is_axially_symmetric, _integration_method))
- ,_darcy_velocities(
- GlobalDim,
- std::vector<double>(_integration_method.getNumberOfPoints()))
- {
- // 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);
- (void)local_matrix_size;
- }
-
- /* commit t? */
- 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<NodalMatrixType>(
- local_M_data, local_matrix_size, local_matrix_size);
- auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
- local_K_data, local_matrix_size, local_matrix_size);
- auto local_b = MathLib::createZeroedVector<NodalVectorType>(local_b_data,
- local_matrix_size);
-
- unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
-
- SpatialPosition pos;
- pos.setElementID(_element.getID());
-
- const double density_water = 1000.0 ;
- const double temperature0 = 20.0 ;
- const double beta = 4.3e-4 ; // thermal expansion coefficient
- const double density_ice = 1000.0 ; // for the mass balance
- const double density_soil = 0.0 ;
- const double specific_heat_capacity_soil = 0;
- const double specific_heat_capacity_ice = 2000.0 ;
- const double specific_heat_capacity_water = 4200.0 ;
- const double thermal_conductivity_ice = 2.0 ;
- const double thermal_conductivity_soil = 3.0 ;
- const double thermal_conductivity_water = 0.65 ;
- const double specific_storage = 0.0 ; // m-1 if 0 then no stoarge term
- const double permeability = 1e-14 ;
- const double viscosity0 = 1e-3 ;
- const double temperature_con = 75 ; // reference temperature for viscosity
- const double g = 9.81 ; // if 0 then do not consider gravity
- double porosity = 0.001 ; // 0.001 is easier for thermal convection converge
- double phi_i = 0.0 ;
- double sigmoid_coeff = 5.0 ;
- double latent_heat = 334000 ;
- double sigmoid_derive = 0.0 ;
- double thermal_conductivity = 0.0 ;
- double heat_capacity = 0.0 ;
- double Real_hydraulic_conductivity = 0.0 ;
- // typedef Matrix<double, 2, 1> Vecterg;
- // Vecterg _gravity_v;
- // _gravity_v[0] = 0.;
- // _gravity_v[1] = 1.;
-
-
- // IntegrationMethod integration_method(_integration_order);
- // unsigned const n_integration_points = integration_method.getNumberOfPoints();
-
- double T_int_pt = 0.0;
- double p_int_pt = 0.0;
-
- for (std::size_t ip(0); ip < n_integration_points; ip++)
- {
- auto const num_nodes = ShapeFunction::NPOINTS; /* difference with n_integration_points?*/
- auto const& sm = _shape_matrices[ip];
- auto const& wp = _integration_method.getWeightedPoint(ip);
- // auto const k = _process_data.thermal_conductivity(_element);
- typedef Matrix<double, num_nodes, num_nodes> MatrixNN;
- MatrixNN _Ktt;
- MatrixNN _Mtt;
- MatrixNN _Kpp;
- MatrixNN _Mpp;
- MatrixNN _Ktp;
- MatrixNN _Mtp;
- MatrixNN _Kpt;
- MatrixNN _Mpt;
- typedef Matrix<double, num_nodes, 1> VecterNN;
- VecterNN _Bpp;
-
- int n = n_integration_points;
- // Order matters: First T, then P!
- NumLib::shapeFunctionInterpolate(local_x, sm.N, T_int_pt, p_int_pt);
-
- // use T_int_pt here ...
-
- // double density_water_T = DensityWater_T(density_water, T_int_pt, temperature0, beta);
- double density_water_T = density_water*(1 - beta*(T_int_pt - temperature0));
- // std::cout << density_water_T << std::endl;
- // double viscosity = Viscosity(viscosity0, T_int_pt, temperature_con, temperature0); // calculate the dynamic viscosity
- double hydraulic_conductivity = permeability/viscosity0 ; // m/s permeability/viscosity
- // phi_i = CalcIceVolFrac(T_int_pt, sigmoid_coeff, porosity);
-
- // Real_hydraulic_conductivity = KozenyKarman(hydraulic_conductivity, porosity, phi_i);
- // Real_hydraulic_conductivity = hydraulic_conductivity ;
- // sigmoid_derive = Calcsigmoidderive(phi_i, sigmoid_coeff, porosity);
-
- /* thermal_conductivity = TotalThermalConductivity(porosity, phi_i, thermal_conductivity_ice,
-thermal_conductivity_soil, thermal_conductivity_water); */
- thermal_conductivity = thermal_conductivity_soil*(1 - porosity)
- + thermal_conductivity_water*porosity; // mixed thermal conductivity
-
- /* heat_capacity = EquaHeatCapacity(phi_i, density_water, density_soil, density_ice,
- specific_heat_capacity_soil, specific_heat_capacity_ice,
- specific_heat_capacity_water, porosity, sigmoid_derive, latent_heat); */
- heat_capacity = density_soil*specific_heat_capacity_soil*(1 - porosity)
- + density_water*specific_heat_capacity_water*porosity; // mixed heat capacity
- auto const detJ_w_NT = (sm.detJ * wp.getWeight() * sm.N.transpose()).eval();
- auto p_nodal_values =
- Eigen::Map<const Eigen::VectorXd>(&local_x[num_nodes], num_nodes);
-
- auto velocity = (-hydraulic_conductivity*
- sm.dNdx*p_nodal_values /*- g*_gravity_v*density_water_T*hydraulic_conductivity*/).eval();
-
- velocity[1] -= density_water_T*g*hydraulic_conductivity ;
- //std::cout << velocity << std::endl;
- auto detJ_w_NT_vT_dNdx =
- (detJ_w_NT * velocity.transpose() * sm.dNdx).eval();
- // matrix assembly
- _Ktt.noalias() += sm.dNdx.transpose() *
- thermal_conductivity * sm.dNdx *
- sm.detJ * wp.getWeight() + detJ_w_NT_vT_dNdx*density_water*specific_heat_capacity_water;
- /*sm.N*density_water*specific_heat_capacity_water*Real_hydraulic_conductivity*
- (sm.dNdx*p_nodal_values).transpose()*sm.dNdx*sm.detJ*wp.getWeight(); */
- _Kpp.noalias() += sm.dNdx.transpose() *
- /* Real_hydraulic_conductivity */ hydraulic_conductivity* sm.dNdx *
- sm.detJ * wp.getWeight();
- _Kpt.noalias() += sm.dNdx.transpose() *
- 0.0 * sm.dNdx *
- sm.detJ * wp.getWeight();
- _Ktp.noalias() += sm.dNdx.transpose() *
- 0.0 * sm.dNdx *
- sm.detJ * wp.getWeight();
- _Mtt.noalias() += sm.N.transpose() *heat_capacity*sm.N *
- sm.detJ * wp.getWeight();
- _Mpp.noalias() += sm.N.transpose() *specific_storage*sm.N *
- sm.detJ * wp.getWeight();
- _Mpt.noalias() += sm.N.transpose() *0.0*sm.N *
- sm.detJ * wp.getWeight();
- _Mtp.noalias() += sm.N.transpose() *(-beta*porosity*0)*sm.N * // beta or -beta or phi*beta or phi*(-beta)
- sm.detJ * wp.getWeight();
- // _Bpp.noalias() += hydraulic_conductivity* sm.detJ * wp.getWeight()*sm.dNdx.transpose().col(_element.getDimension()-1)*g*density_water_T;
- _Bpp.noalias() += hydraulic_conductivity* sm.detJ * wp.getWeight()*sm.dNdx.transpose().col(1)*g*density_water_T;
- /* with Oberbeck-Boussing assumption density difference only exists in buoyancy effects */
-
- local_K.block<num_nodes,num_nodes>(0,0).noalias() += _Ktt;
- local_M.block<num_nodes,num_nodes>(0,0).noalias() += _Mtt;
- local_K.block<num_nodes,num_nodes>(num_nodes,num_nodes).noalias() += _Kpp;
- local_M.block<num_nodes,num_nodes>(num_nodes,num_nodes).noalias() += _Mpp;
- local_K.block<num_nodes,num_nodes>(num_nodes,0).noalias() += _Ktp;
- local_M.block<num_nodes,num_nodes>(num_nodes,0).noalias() += _Mtp;
- local_K.block<num_nodes,num_nodes>(0,num_nodes).noalias() += _Kpt;
- local_M.block<num_nodes,num_nodes>(0,num_nodes).noalias() += _Mpt;
- local_b.block<num_nodes,1>(num_nodes,0).noalias() -= _Bpp ;
- // heat flux only computed for output.
- /* auto const heat_flux = (-thermal_conductivity * sm.dNdx *
- Eigen::Map<const NodalVectorType>(&local_x[0], num_nodes)
- ).eval();
-
- for (unsigned d=0; d<GlobalDim; ++d) {
- _heat_fluxes[d][ip] = heat_flux[d];
- } */
- // velocity computed for output.
- for (unsigned d = 0; d < GlobalDim; ++d){
- _darcy_velocities[d][ip] = velocity[d];
- }
- }
-
-
-// K.add(indices, _localK);
-// M.add(indices, _localM);
-// b.add(indices.rows, _localRhs);
- }
-
- Eigen::Map<const Eigen::RowVectorXd>
- getShapeMatrix(const unsigned integration_point) const override
- {
- auto const& N = _shape_matrices[integration_point].N;
-
- // assumes N is stored contiguously in memory
- return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
- }
-
-// std::vector<double> const&
-// getIntPtHeatFluxX(std::vector<double>& /*cache*/) const override
-// {
-// assert(_heat_fluxes.size() > 0);
-// return _heat_fluxes[0];
-// }
-
-// std::vector<double> const&
-// getIntPtHeatFluxY(std::vector<double>& /*cache*/) const override
-// {
-// assert(_heat_fluxes.size() > 1);
-// return _heat_fluxes[1];
-// }
-
-// std::vector<double> const&
-// getIntPtHeatFluxZ(std::vector<double>& /*cache*/) const override
-// {
-// assert(_heat_fluxes.size() > 2);
-// return _heat_fluxes[2];
- // }
-
- std::vector<double> const&
- getIntPtDarcyVelocityX(std::vector<double>& /*cache*/) const override
- {
- assert(_darcy_velocities.size() > 0);
- return _darcy_velocities[0];
- }
-
- std::vector<double> const&
- getIntPtDarcyVelocityY(std::vector<double>& /*cache*/) const override
- {
- assert(_darcy_velocities.size() > 1);
- return _darcy_velocities[1];
- }
-
- std::vector<double> const&
- getIntPtDarcyVelocityZ(std::vector<double>& /*cache*/) const override
- {
- assert(_darcy_velocities.size() > 2);
- return _darcy_velocities[2];
- }
-
-private:
- MeshLib::Element const& _element;
- DensityDrivenFlowProcessData const& _process_data;
- // DensityDrivenFlowProcessData const& _process_data;
- // NodalVectorType _gravity_v;
-
- IntegrationMethod const _integration_method;
- std::vector<ShapeMatrices> _shape_matrices;
- // std::vector<std::vector<double>> _heat_fluxes
- // = std::vector<std::vector<double>>(
- // GlobalDim, std::vector<double>(ShapeFunction::NPOINTS));
- std::vector<std::vector<double>> _darcy_velocities;
-
-};
-
-
-} // namespace DensityDrivenFlow
-} // namespace ProcessLib
-
-#endif // PROCESS_LIB_DensityDrivenFlow_FEM_H_
diff --git a/ProcessLib/DensityDrivenFlow/DensityDrivenFlowProcess.cpp b/ProcessLib/DensityDrivenFlow/DensityDrivenFlowProcess.cpp
deleted file mode 100644
index 481d676029be5661616bd015e971b050c6529b67..0000000000000000000000000000000000000000
--- a/ProcessLib/DensityDrivenFlow/DensityDrivenFlowProcess.cpp
+++ /dev/null
@@ -1,120 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include "DensityDrivenFlowProcess.h"
-
-#include <cassert>
-
-#include "ProcessLib/Utils/CreateLocalAssemblers.h"
-
-namespace ProcessLib
-{
-namespace DensityDrivenFlow
-{
-DensityDrivenFlowProcess::DensityDrivenFlowProcess(
- 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,
- DensityDrivenFlowProcessData&& 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))
- {
- }
-
-void DensityDrivenFlowProcess::initializeConcreteProcess(
- NumLib::LocalToGlobalIndexMap const& dof_table,
- MeshLib::Mesh const& mesh,
- unsigned const integration_order)
-{
- ProcessLib::createLocalAssemblers<LocalAssemblerData>(
- mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
- mesh.isAxiallySymmetric(), integration_order, _process_data);
-
- /* _secondary_variables.addSecondaryVariable(
- "heat_flux_x", 1,
- makeExtrapolator(
- getExtrapolator(), _local_assemblers,
- &DensityDrivenFlowLocalAssemblerInterface::getIntPtHeatFluxX));
-
- if (mesh.getDimension() > 1) {
- _secondary_variables.addSecondaryVariable(
- "heat_flux_y", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
- &DensityDrivenFlowLocalAssemblerInterface::
- getIntPtHeatFluxY));
- }
- if (mesh.getDimension() > 2) {
- _secondary_variables.addSecondaryVariable(
- "heat_flux_z", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
- &DensityDrivenFlowLocalAssemblerInterface::
- getIntPtHeatFluxZ));
- }*/
-
- _secondary_variables.addSecondaryVariable(
- "darcy_velocity_x", 1,
- makeExtrapolator(
- getExtrapolator(), _local_assemblers,
- &DensityDrivenFlowLocalAssemblerInterface::getIntPtDarcyVelocityX));
-
- if (mesh.getDimension() > 1) {
- _secondary_variables.addSecondaryVariable(
- "darcy_velocity_y", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
- &DensityDrivenFlowLocalAssemblerInterface::
- getIntPtDarcyVelocityY));
- }
- if (mesh.getDimension() > 2) {
- _secondary_variables.addSecondaryVariable(
- "darcy_velocity_z", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
- &DensityDrivenFlowLocalAssemblerInterface::
- getIntPtDarcyVelocityZ));
- }
-
-}
-
-void DensityDrivenFlowProcess::assembleConcreteProcess(const double t,
- GlobalVector const& x,
- GlobalMatrix& M,
- GlobalMatrix& K,
- GlobalVector& b)
-{
- DBUG("Assemble DensityDrivenFlowProcess.");
-
- // 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);
-}
-
-void DensityDrivenFlowProcess::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)
-{
- DBUG("AssembleWithJacobian DensityDrivenFlowProcess.");
-
- // 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);
-}
-
-
-} // namespace DensityDrivenFlow
-} // namespace ProcessLib
-
diff --git a/ProcessLib/DensityDrivenFlow/DensityDrivenFlowProcess.h b/ProcessLib/DensityDrivenFlow/DensityDrivenFlowProcess.h
deleted file mode 100644
index 2e6d176a68d53977a6321f418346b15fe0478b68..0000000000000000000000000000000000000000
--- a/ProcessLib/DensityDrivenFlow/DensityDrivenFlowProcess.h
+++ /dev/null
@@ -1,69 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef PROCESS_LIB_DensityDrivenFlowPROCESS_H_
-#define PROCESS_LIB_DensityDrivenFlowPROCESS_H_
-
-#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
-#include "ProcessLib/Process.h"
-#include "DensityDrivenFlowFEM.h"
-#include "DensityDrivenFlowProcessData.h"
-
-
-namespace ProcessLib
-{
-namespace DensityDrivenFlow
-{
-class DensityDrivenFlowProcess final : public Process
-{
-
-public:
- DensityDrivenFlowProcess(
- 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,
- DensityDrivenFlowProcessData&& process_data,
- SecondaryVariableCollection&& secondary_variables,
- NumLib::NamedFunctionCaller&& named_function_caller);
-
- //! \name ODESystem interface
- //! @{
-
- bool isLinear() const override { return false; }
- //! @}
-
-private:
- void initializeConcreteProcess(
- NumLib::LocalToGlobalIndexMap const& dof_table,
- MeshLib::Mesh const& mesh,
- unsigned const integration_order) override;
-
- void assembleConcreteProcess(const double t, GlobalVector const& x,
- GlobalMatrix& M, GlobalMatrix& K,
- GlobalVector& b) override;
-
- 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) override;
-
- DensityDrivenFlowProcessData _process_data;
-
- std::vector<std::unique_ptr<DensityDrivenFlowLocalAssemblerInterface>>
-_local_assemblers;
-};
-
-} // namespace DensityDrivenFlow
-} // namespace ProcessLib
-
-#endif // PROCESS_LIB_DensityDrivenFlowPROCESS_H_
diff --git a/ProcessLib/DensityDrivenFlow/DensityDrivenFlowProcessData.h b/ProcessLib/DensityDrivenFlow/DensityDrivenFlowProcessData.h
deleted file mode 100644
index 9d78d62a4bc8ceeda5c621e5ac8f76b629c618bc..0000000000000000000000000000000000000000
--- a/ProcessLib/DensityDrivenFlow/DensityDrivenFlowProcessData.h
+++ /dev/null
@@ -1,57 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef PROCESSLIB_DensityDrivenFlow_DensityDrivenFlowPROCESSDATA_H
-#define PROCESSLIB_DensityDrivenFlow_DensityDrivenFlowPROCESSDATA_H
-
-namespace MeshLib
-{
- class Element;
-}
-
-
-namespace ProcessLib
-{
-
-template <typename ReturnType>
-struct Parameter;
-
-namespace DensityDrivenFlow
-{
-
-// copy that and change names first still only thermal conductivity, later on more parameters will be added.
-struct DensityDrivenFlowProcessData
-{
- DensityDrivenFlowProcessData(
- ProcessLib::Parameter<double> const&
- thermal_conductivity_
- )
- : thermal_conductivity(thermal_conductivity_)
- {}
-
- DensityDrivenFlowProcessData(DensityDrivenFlowProcessData&& other)
- : thermal_conductivity(other.thermal_conductivity)
- {}
-
- //! Copies are forbidden.
- DensityDrivenFlowProcessData(DensityDrivenFlowProcessData const&) = delete;
-
- //! Assignments are not needed.
- void operator=(DensityDrivenFlowProcessData const&) = delete;
-
- //! Assignments are not needed.
- void operator=(DensityDrivenFlowProcessData&&) = delete;
-
- Parameter<double> const& thermal_conductivity;
-};
-
-} // namespace DensityDrivenFlow
-} // namespace ProcessLib
-
-#endif // PROCESSLIB_DensityDrivenFlow_DensityDrivenFlowPROCESSDATA_H
diff --git a/ProcessLib/HT/CreateHTProcess.cpp b/ProcessLib/HT/CreateHTProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3ac5094761bbace4c8414aa0d6721ac696654c03
--- /dev/null
+++ b/ProcessLib/HT/CreateHTProcess.cpp
@@ -0,0 +1,176 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreateHTProcess.h"
+
+#include "HTProcess.h"
+#include "HTProcessData.h"
+#include "ProcessLib/Parameter/ConstantParameter.h"
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace ProcessLib
+{
+namespace HT
+{
+std::unique_ptr<Process> createHTProcess(
+ 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", "HT");
+
+ DBUG("Create HTProcess.");
+
+ // Process variable.
+ auto process_variables = findProcessVariables(
+ variables, config,
+ {"temperature", "pressure"}); // configure two Pcs
+
+ // Porosity parameter.
+ auto& porosity = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__HT__porosity}
+ "porosity", parameters, 1);
+ DBUG("Use \'%s\' as porosity parameter.", porosity.name.c_str());
+
+ // Parameter for the intrinsic permeability (only one scalar per element,
+ // i.e., the isotropic case is handled at the moment)
+ auto& intrinsic_permeability = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__HT__intrinsic_permeability}
+ "intrinsic_permeability", parameters, 1);
+ DBUG("Use \'%s\' as intrinsic_permeability parameter.",
+ intrinsic_permeability.name.c_str());
+
+ // Parameter for the specific storage.
+ auto& specific_storage = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__HT__specific_storage}
+ "specific_storage", parameters, 1);
+ DBUG("Use \'%s\' as specific storage parameter.", specific_storage.name.c_str());
+
+ // Parameter for the reference_temperature.
+ auto& reference_temperature_fluid_density_model = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__HT__reference_temperature}
+ "reference_temperature_fluid_density_model", parameters, 1);
+ DBUG(
+ "Use \'%s\' as reference temperature for the fluid density model "
+ "parameter.",
+ reference_temperature_fluid_density_model.name.c_str());
+
+ // Parameter for the viscosity.
+ auto& viscosity = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__HT__viscosity}
+ "viscosity", parameters, 1);
+ DBUG("Use \'%s\' as viscosity parameter.", viscosity.name.c_str());
+
+ // Parameter for the density of the solid.
+ auto& density_solid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__HT__density_solid}
+ "density_solid", parameters, 1);
+ DBUG("Use \'%s\' as density_solid parameter.", density_solid.name.c_str());
+
+ // Parameter for the density of the fluid.
+ auto& density_fluid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__HT__density_fluid}
+ "density_fluid", parameters, 1);
+ DBUG("Use \'%s\' as density_fluid parameter.", density_fluid.name.c_str());
+
+ // Parameter for the specific heat capacity of the solid.
+ auto& specific_heat_capacity_solid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__HT__specific_heat_capacity_solid}
+ "specific_heat_capacity_solid", parameters, 1);
+ DBUG("Use \'%s\' as specific_heat_capacity_solid parameter.",
+ specific_heat_capacity_solid.name.c_str());
+
+ // Parameter for the specific heat capacity of the fluid.
+ auto& specific_heat_capacity_fluid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__HT__specific_heat_capacity_fluid}
+ "specific_heat_capacity_fluid", parameters, 1);
+ DBUG("Use \'%s\' as specific_heat_capacity_fluid parameter.",
+ specific_heat_capacity_fluid.name.c_str());
+
+ // Parameter for the thermal conductivity of the solid (only one scalar per
+ // element, i.e., the isotropic case is handled at the moment)
+ auto& thermal_conductivity_solid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__HT__thermal_conductivity_solid}
+ "thermal_conductivity_solid", parameters, 1);
+ DBUG("Use \'%s\' as thermal_conductivity_solid parameter.",
+ thermal_conductivity_solid.name.c_str());
+
+ // Parameter for the thermal conductivity of the fluid.
+ auto& thermal_conductivity_fluid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__HT__thermal_conductivity_fluid}
+ "thermal_conductivity_fluid", parameters, 1);
+ DBUG("Use \'%s\' as thermal_conductivity_fluid parameter.",
+ thermal_conductivity_fluid.name.c_str());
+
+ // Parameter for the thermal expansion coefficient.
+ auto& thermal_expansion_coefficient_fluid = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{process__HT__thermal_expansion_coefficient_fluid}
+ "thermal_expansion_coefficient_fluid", parameters, 1);
+ DBUG("Use \'%s\' as thermal expansion coefficient of the fluid parameter.",
+ thermal_expansion_coefficient_fluid.name.c_str());
+
+ // Specific body force parameter.
+ Eigen::Vector3d specific_body_force;
+ std::vector<double> const b =
+ //! \ogs_file_param_special{process__HT__specific_body_force}
+ config.getConfigParameter<std::vector<double>>("specific_body_force");
+ assert(b.size() > 0 && b.size() < 4);
+ bool const has_gravity = MathLib::toVector(b).norm() > 0;
+ if (has_gravity)
+ std::copy_n(b.data(), b.size(), specific_body_force.data());
+
+ HTProcessData process_data{
+ porosity,
+ intrinsic_permeability,
+ specific_storage,
+ viscosity,
+ density_solid,
+ density_fluid,
+ specific_heat_capacity_solid,
+ specific_heat_capacity_fluid,
+ thermal_conductivity_solid,
+ thermal_conductivity_fluid,
+ thermal_expansion_coefficient_fluid,
+ reference_temperature_fluid_density_model,
+ specific_body_force,
+ has_gravity};
+
+ SecondaryVariableCollection secondary_variables;
+
+ NumLib::NamedFunctionCaller named_function_caller(
+ {"HT_temperature_pressure"});
+
+ ProcessLib::parseSecondaryVariables(config, secondary_variables,
+ named_function_caller);
+
+ return std::unique_ptr<Process>{new HTProcess{
+ 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)}};
+}
+
+} // namespace HT
+} // namespace ProcessLib
diff --git a/ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.h b/ProcessLib/HT/CreateHTProcess.h
similarity index 67%
rename from ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.h
rename to ProcessLib/HT/CreateHTProcess.h
index 35e1e7c80e906aed348aac6fb2fd35e4f6f1f3f2..e3f6f3cdc264c837a0ef4f0964244893444fdb02 100644
--- a/ProcessLib/DensityDrivenFlow/CreateDensityDrivenFlowProcess.h
+++ b/ProcessLib/HT/CreateHTProcess.h
@@ -7,19 +7,17 @@
*
*/
-#ifndef PROCESS_LIB_CREATE_DensityDrivenFlowPROCESS_H_
-#define PROCESS_LIB_CREATE_DensityDrivenFlowPROCESS_H_
+#ifndef PROCESS_LIB_CREATE_HTPROCESS_H_
+#define PROCESS_LIB_CREATE_HTPROCESS_H_
#include <memory>
#include "ProcessLib/Process.h"
-
namespace ProcessLib
{
-namespace DensityDrivenFlow
+namespace HT
{
-
-std::unique_ptr<Process> createDensityDrivenFlowProcess(
+std::unique_ptr<Process> createHTProcess(
MeshLib::Mesh& mesh,
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<ProcessVariable> const& variables,
@@ -27,7 +25,7 @@ std::unique_ptr<Process> createDensityDrivenFlowProcess(
unsigned const integration_order,
BaseLib::ConfigTree const& config);
-} // namespace DensityDrivenFlow
-} // namespace ProcessLib
+} // namespace HT
+} // namespace ProcessLib
-#endif // PROCESS_LIB_CREATE_DensityDrivenFlowPROCESS_H_
+#endif // PROCESS_LIB_CREATE_HTPROCESS_H_
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
new file mode 100644
index 0000000000000000000000000000000000000000..7878931f6bef843da593ec541822bb0f07dc44b3
--- /dev/null
+++ b/ProcessLib/HT/HTFEM.h
@@ -0,0 +1,244 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_HTPROCESS_FEM_H_
+#define PROCESS_LIB_HTPROCESS_FEM_H_
+
+#include <iostream>
+#include <Eigen/Dense>
+#include <vector>
+
+
+#include "HTProcessData.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/LocalAssemblerInterface.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+namespace ProcessLib
+{
+namespace HT
+{
+const unsigned NUM_NODAL_DOF = 2;
+
+class HTLocalAssemblerInterface
+ : public ProcessLib::LocalAssemblerInterface,
+ public NumLib::ExtrapolatableElement
+{
+public:
+ virtual std::vector<double> const& getIntPtDarcyVelocityX(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtDarcyVelocityY(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+ std::vector<double>& /*cache*/) const = 0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class LocalAssemblerData : public HTLocalAssemblerInterface
+{
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+ using LocalAssemblerTraits = ProcessLib::LocalAssemblerTraits<
+ ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>;
+ using NodalMatrixType = typename LocalAssemblerTraits::LocalMatrix;
+ using NodalVectorType = typename LocalAssemblerTraits::LocalVector;
+ using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+
+public:
+ LocalAssemblerData(MeshLib::Element const& element,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ HTProcessData const& process_data)
+ : _element(element),
+ _process_data(process_data),
+ _integration_method(integration_order),
+ _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ element, is_axially_symmetric, _integration_method)),
+ _darcy_velocities(
+ GlobalDim,
+ std::vector<double>(_integration_method.getNumberOfPoints()))
+ {
+ // 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);
+ (void)local_matrix_size;
+ }
+
+ 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<NodalMatrixType>(
+ local_M_data, local_matrix_size, local_matrix_size);
+ auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
+ local_K_data, local_matrix_size, local_matrix_size);
+ auto local_b = MathLib::createZeroedVector<NodalVectorType>(
+ local_b_data, local_matrix_size);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+
+ auto const num_nodes = ShapeFunction::NPOINTS;
+ auto p_nodal_values =
+ Eigen::Map<const Eigen::VectorXd>(&local_x[num_nodes], num_nodes);
+
+ auto const & b = _process_data.specific_body_force.head(GlobalDim);
+
+ for (std::size_t ip(0); ip < n_integration_points; ip++)
+ {
+ pos.setIntegrationPoint(ip);
+ auto const temperature0 =
+ _process_data.reference_temperature_fluid_density_model(t,
+ pos)[0];
+ auto const beta =
+ _process_data.thermal_expansion_coefficient(t, pos)[0];
+
+ auto const density_fluid = _process_data.density_fluid(t, pos)[0];
+ auto const density_solid = _process_data.density_solid(t, pos)[0];
+ auto const specific_storage =
+ _process_data.specific_storage(t, pos)[0];
+ auto const intrinsic_permeability =
+ _process_data.intrinsic_permeability(t, pos)[0];
+ auto const viscosity0 = _process_data.viscosity(t, pos)[0];
+
+ auto const porosity = _process_data.porosity(t, pos)[0];
+
+ auto const specific_heat_capacity_solid =
+ _process_data.specific_heat_capacity_solid(t, pos)[0];
+ auto const specific_heat_capacity_fluid =
+ _process_data.specific_heat_capacity_fluid(t, pos)[0];
+ double const heat_capacity =
+ density_solid * specific_heat_capacity_solid * (1 - porosity) +
+ density_fluid * specific_heat_capacity_fluid * porosity;
+
+ 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];
+ double const thermal_conductivity =
+ thermal_conductivity_solid * (1 - porosity) +
+ thermal_conductivity_fluid * porosity;
+
+ auto const& sm = _shape_matrices[ip];
+ auto const& wp = _integration_method.getWeightedPoint(ip);
+ 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);
+
+ double T_int_pt = 0.0;
+ double p_int_pt = 0.0;
+ // Order matters: First T, then P!
+ NumLib::shapeFunctionInterpolate(local_x, sm.N, T_int_pt, p_int_pt);
+
+ double const delta_T(T_int_pt - temperature0);
+ // TODO include this via material lib
+ double density_water_T = density_fluid * (1 - beta * delta_T);
+ // TODO include viscosity computations via material lib
+ // double const viscosity = viscosity0 * std::exp(- delta_T/75.0);
+ double const perm_visc =
+ intrinsic_permeability / viscosity0;
+
+ Eigen::Matrix<double, -1, 1, 0, -1, 1> const velocity =
+ -perm_visc * (sm.dNdx * p_nodal_values - density_water_T * b);
+
+ auto const integral_term =
+ sm.integralMeasure * sm.detJ * wp.getWeight();
+ // matrix assembly
+ Ktt.noalias() +=
+ integral_term *
+ (sm.dNdx.transpose() * thermal_conductivity * sm.dNdx +
+ sm.N.transpose() * velocity.transpose() * sm.dNdx *
+ density_fluid * specific_heat_capacity_fluid);
+ Kpp.noalias() +=
+ integral_term * sm.dNdx.transpose() * perm_visc * sm.dNdx;
+ Mtt.noalias() +=
+ integral_term * sm.N.transpose() * heat_capacity * sm.N;
+ Mpp.noalias() +=
+ integral_term * sm.N.transpose() * specific_storage * sm.N;
+ Bp += perm_visc * integral_term * sm.dNdx.transpose() * b *
+ density_water_T;
+ /* with Oberbeck-Boussing assumption density difference only exists
+ * in buoyancy effects */
+
+ // velocity computed for output.
+ for (unsigned d = 0; d < GlobalDim; ++d)
+ {
+ _darcy_velocities[d][ip] = velocity[d];
+ }
+ }
+ }
+
+ Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+ const unsigned integration_point) const override
+ {
+ auto const& N = _shape_matrices[integration_point].N;
+
+ // assumes N is stored contiguously in memory
+ return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+ }
+
+ std::vector<double> const& getIntPtDarcyVelocityX(
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(_darcy_velocities.size() > 0);
+ return _darcy_velocities[0];
+ }
+
+ std::vector<double> const& getIntPtDarcyVelocityY(
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(_darcy_velocities.size() > 1);
+ return _darcy_velocities[1];
+ }
+
+ std::vector<double> const& getIntPtDarcyVelocityZ(
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(_darcy_velocities.size() > 2);
+ return _darcy_velocities[2];
+ }
+
+private:
+ MeshLib::Element const& _element;
+ HTProcessData const& _process_data;
+
+ IntegrationMethod const _integration_method;
+ std::vector<ShapeMatrices> _shape_matrices;
+ std::vector<std::vector<double>> _darcy_velocities;
+};
+
+} // namespace HT
+} // namespace ProcessLib
+
+#endif // PROCESS_LIB_HT_FEM_H_
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..52ace60069225e4f005a348ad190828a23a8acb5
--- /dev/null
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -0,0 +1,101 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "HTProcess.h"
+
+#include <cassert>
+
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+
+namespace ProcessLib
+{
+namespace HT
+{
+HTProcess::HTProcess(
+ 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,
+ HTProcessData&& 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))
+{
+}
+
+void HTProcess::initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ 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, _process_data);
+
+ _secondary_variables.addSecondaryVariable(
+ "darcy_velocity_x", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &HTLocalAssemblerInterface::getIntPtDarcyVelocityX));
+
+ if (mesh.getDimension() > 1)
+ {
+ _secondary_variables.addSecondaryVariable(
+ "darcy_velocity_y", 1,
+ makeExtrapolator(getExtrapolator(), _local_assemblers,
+ &HTLocalAssemblerInterface::
+ getIntPtDarcyVelocityY));
+ }
+ if (mesh.getDimension() > 2)
+ {
+ _secondary_variables.addSecondaryVariable(
+ "darcy_velocity_z", 1,
+ makeExtrapolator(getExtrapolator(), _local_assemblers,
+ &HTLocalAssemblerInterface::
+ getIntPtDarcyVelocityZ));
+ }
+}
+
+void HTProcess::assembleConcreteProcess(const double t,
+ GlobalVector const& x,
+ GlobalMatrix& M,
+ GlobalMatrix& K,
+ GlobalVector& b)
+{
+ DBUG("Assemble HTProcess.");
+
+ // 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);
+}
+
+void HTProcess::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)
+{
+ DBUG("AssembleWithJacobian HTProcess.");
+
+ // 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);
+}
+
+} // namespace HT
+} // namespace ProcessLib
+
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..a2eaa4e50d9c7c5bc7ab07b7fcd78bef2f7e2f7d
--- /dev/null
+++ b/ProcessLib/HT/HTProcess.h
@@ -0,0 +1,67 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_HTPROCESS_H_
+#define PROCESS_LIB_HTPROCESS_H_
+
+#include "HTFEM.h"
+#include "HTProcessData.h"
+#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace HT
+{
+class HTProcess final : public Process
+{
+public:
+ HTProcess(
+ 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,
+ HTProcessData&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller);
+
+ //! \name ODESystem interface
+ //! @{
+
+ bool isLinear() const override { return false; }
+ //! @}
+
+private:
+ void initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order) override;
+
+ void assembleConcreteProcess(const double t, GlobalVector const& x,
+ GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b) override;
+
+ 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) override;
+
+ HTProcessData _process_data;
+
+ std::vector<std::unique_ptr<HTLocalAssemblerInterface>>
+ _local_assemblers;
+};
+
+} // namespace HT
+} // namespace ProcessLib
+
+#endif // PROCESS_LIB_HTPROCESS_H_
diff --git a/ProcessLib/HT/HTProcessData.h b/ProcessLib/HT/HTProcessData.h
new file mode 100644
index 0000000000000000000000000000000000000000..0270fd520c658bc4f8e89841381da0ab9bb67552
--- /dev/null
+++ b/ProcessLib/HT/HTProcessData.h
@@ -0,0 +1,108 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_HT_HTPROCESSDATA_H
+#define PROCESSLIB_HT_HTPROCESSDATA_H
+
+namespace MeshLib
+{
+class Element;
+}
+
+namespace ProcessLib
+{
+template <typename ReturnType>
+struct Parameter;
+
+namespace HT
+{
+struct HTProcessData
+{
+ HTProcessData(
+ ProcessLib::Parameter<double> const& porosity_,
+ ProcessLib::Parameter<double> const& intrinsic_permeability_,
+ ProcessLib::Parameter<double> const& specific_storage_,
+ ProcessLib::Parameter<double> const& viscosity_,
+ ProcessLib::Parameter<double> const& density_solid_,
+ ProcessLib::Parameter<double> const& density_fluid_,
+ ProcessLib::Parameter<double> const& specific_heat_capacity_solid_,
+ ProcessLib::Parameter<double> const& specific_heat_capacity_fluid_,
+ ProcessLib::Parameter<double> const& thermal_conductivity_solid_,
+ ProcessLib::Parameter<double> const& thermal_conductivity_fluid_,
+ ProcessLib::Parameter<double> const& thermal_expansion_coefficient_,
+ ProcessLib::Parameter<double> const&
+ reference_temperature_fluid_density_model_,
+ Eigen::Vector3d const& specific_body_force_,
+ bool const has_gravity_)
+ : porosity(porosity_),
+ intrinsic_permeability(intrinsic_permeability_),
+ specific_storage(specific_storage_),
+ viscosity(viscosity_),
+ density_solid(density_solid_),
+ density_fluid(density_fluid_),
+ specific_heat_capacity_solid(specific_heat_capacity_solid_),
+ specific_heat_capacity_fluid(specific_heat_capacity_fluid_),
+ thermal_conductivity_solid(thermal_conductivity_solid_),
+ thermal_conductivity_fluid(thermal_conductivity_fluid_),
+ thermal_expansion_coefficient(thermal_expansion_coefficient_),
+ reference_temperature_fluid_density_model(
+ reference_temperature_fluid_density_model_),
+ specific_body_force(specific_body_force_),
+ has_gravity(has_gravity_)
+ {
+ }
+
+ HTProcessData(HTProcessData&& other)
+ : porosity(other.porosity),
+ intrinsic_permeability(other.intrinsic_permeability),
+ specific_storage(other.specific_storage),
+ viscosity(other.viscosity),
+ density_solid(other.density_solid),
+ density_fluid(other.density_fluid),
+ specific_heat_capacity_solid(other.specific_heat_capacity_solid),
+ specific_heat_capacity_fluid(other.specific_heat_capacity_fluid),
+ thermal_conductivity_solid(other.thermal_conductivity_solid),
+ thermal_conductivity_fluid(other.thermal_conductivity_fluid),
+ thermal_expansion_coefficient(other.thermal_expansion_coefficient),
+ reference_temperature_fluid_density_model(
+ other.reference_temperature_fluid_density_model),
+ specific_body_force(other.specific_body_force),
+ has_gravity(other.has_gravity)
+ {
+ }
+
+ //! Copies are forbidden.
+ HTProcessData(HTProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(HTProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(HTProcessData&&) = delete;
+
+ Parameter<double> const& porosity;
+ Parameter<double> const& intrinsic_permeability;
+ Parameter<double> const& specific_storage;
+ Parameter<double> const& viscosity;
+ Parameter<double> const& density_solid;
+ Parameter<double> const& density_fluid;
+ Parameter<double> const& specific_heat_capacity_solid;
+ Parameter<double> const& specific_heat_capacity_fluid;
+ Parameter<double> const& thermal_conductivity_solid;
+ Parameter<double> const& thermal_conductivity_fluid;
+ Parameter<double> const& thermal_expansion_coefficient;
+ Parameter<double> const& reference_temperature_fluid_density_model;
+ Eigen::Vector3d const specific_body_force;
+ bool const has_gravity;
+};
+
+} // namespace HT
+} // namespace ProcessLib
+
+#endif // PROCESSLIB_HT_HTPROCESSDATA_H