Skip to content
Snippets Groups Projects
Commit 58d5c91d authored by wenqing's avatar wenqing
Browse files

[Coupling] Modified the process definitions in Process/HT for the assembly for... 

[Coupling] Modified the process definitions in Process/HT for the assembly for both of the monolithic scheme and the staggered scheme.
parent c2db4b35
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
ProcessLib/HT/HTFEM.h
+ 12
− 200
View file @ 58d5c91d
......@@ -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;
......
ProcessLib/HT/HTLocalAssemblerInterface.h 0 → 100644
+ 53
− 0
View file @ 58d5c91d
/**
* \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
ProcessLib/HT/HTProcess.cpp
+ 18
− 4
View file @ 58d5c91d
......@@ -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",
......
ProcessLib/HT/MonolithicHTFEM.h 0 → 100644
+ 218
− 0
View file @ 58d5c91d
/**
* \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
ProcessLib/HT/StaggeredHTFEM-impl.h 0 → 100644
+ 198
− 0
View file @ 58d5c91d
/**
* \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
ProcessLib/HT/StaggeredHTFEM.h 0 → 100644
+ 87
− 0
View file @ 58d5c91d
/**
* \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"
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment