Detailed Description

template<typename ShapeFunction, int GlobalDim>
class ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >

Definition at line 27 of file HTFEM.h.

#include <HTFEM.h>

Inheritance diagram for ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >:

Collaboration diagram for ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >:

Public Member Functions
	HTFEM (MeshLib::Element const &element, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, HTProcessData const &process_data, const unsigned dof_per_node)
Eigen::Map< const Eigen::RowVectorXd >	getShapeMatrix (const unsigned integration_point) const override
	Provides the shape matrix at the given integration point.
Eigen::Vector3d	getFlux (MathLib::Point3d const &pnt_local_coords, double const t, std::vector< double > const &local_x) const override
Public Member Functions inherited from ProcessLib::HT::HTLocalAssemblerInterface
	HTLocalAssemblerInterface ()=default
virtual std::vector< double > const &	getIntPtDarcyVelocity (const double, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &) const =0
Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual	~LocalAssemblerInterface ()=default
virtual void	setInitialConditions (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, int const process_id)
virtual void	initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
virtual void	preAssemble (double const, double const, std::vector< double > const &)
virtual void	assemble (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
virtual void	assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
virtual void	assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
virtual void	assembleWithJacobianForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
virtual void	computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
virtual void	preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void	postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
void	postNonLinearSolver (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
virtual Eigen::Vector3d	getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double > > const &) const
	Fits to staggered scheme.
virtual int	getNumberOfVectorElementsForDeformation () const
Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual	~ExtrapolatableElement ()=default

Protected Member Functions
double	getHeatEnergyCoefficient (MaterialPropertyLib::VariableArray const &vars, const double porosity, const double fluid_density, const double specific_heat_capacity_fluid, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
GlobalDimMatrixType	getThermalConductivityDispersivity (MaterialPropertyLib::VariableArray const &vars, const double fluid_density, const double specific_heat_capacity_fluid, const GlobalDimVectorType &velocity, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
std::vector< double > const &	getIntPtDarcyVelocityLocal (const double t, std::vector< double > const &local_x, std::vector< double > &cache) const

Protected Attributes
MeshLib::Element const &	_element
HTProcessData const &	_process_data
NumLib::GenericIntegrationMethod const &	_integration_method
std::vector< IntegrationPointData< GlobalDimNodalMatrixType > >	_ip_data

Static Protected Attributes
static const int	pressure_index = ShapeFunction::NPOINTS
static const int	pressure_size = ShapeFunction::NPOINTS
static const int	temperature_index = 0
static const int	temperature_size = ShapeFunction::NPOINTS

Private Types
using	ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>
using	ShapeMatrices = typename ShapeMatricesType::ShapeMatrices
using	NodalVectorType = typename ShapeMatricesType::NodalVectorType
using	NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType
using	GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType
using	GlobalDimNodalMatrixType
using	GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType

Member Typedef Documentation

◆ GlobalDimMatrixType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType

private

Definition at line 38 of file HTFEM.h.

◆ GlobalDimNodalMatrixType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::GlobalDimNodalMatrixType

private

Initial value:

typename ShapeMatricesType::GlobalDimNodalMatrixType

EigenFixedShapeMatrixPolicy::GlobalDimNodalMatrixType

MatrixType< GlobalDim, ShapeFunction::NPOINTS > GlobalDimNodalMatrixType

Definition ShapeMatrixPolicy.h:75

Definition at line 36 of file HTFEM.h.

◆ GlobalDimVectorType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType

private

Definition at line 35 of file HTFEM.h.

◆ NodalRowVectorType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType

private

Definition at line 33 of file HTFEM.h.

◆ NodalVectorType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::NodalVectorType = typename ShapeMatricesType::NodalVectorType

private

Definition at line 32 of file HTFEM.h.

◆ ShapeMatrices

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::ShapeMatrices = typename ShapeMatricesType::ShapeMatrices

private

Definition at line 30 of file HTFEM.h.

◆ ShapeMatricesType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>

private

Definition at line 29 of file HTFEM.h.

Constructor & Destructor Documentation

◆ HTFEM()

template<typename ShapeFunction, int GlobalDim>

ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::HTFEM	(	MeshLib::Element const &	element,
		std::size_t const	local_matrix_size,
		NumLib::GenericIntegrationMethod const &	integration_method,
		bool const	is_axially_symmetric,
		HTProcessData const &	process_data,
		const unsigned	dof_per_node )

inline

Definition at line 41 of file HTFEM.h.

        : HTLocalAssemblerInterface(),
          _element(element),
          _process_data(process_data),
          _integration_method(integration_method)
    {
        // This assertion is valid only if all nodal d.o.f. use the same shape
        // matrices.
        assert(local_matrix_size == ShapeFunction::NPOINTS * dof_per_node);
        (void)local_matrix_size;
        (void)dof_per_node;
 
        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();
        _ip_data.reserve(n_integration_points);
 
        ParameterLib::SpatialPosition pos;
        pos.setElementID(_element.getID());
 
        double const aperture_size = _process_data.aperture_size(0.0, pos)[0];
 
        auto const shape_matrices =
            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
                                      GlobalDim>(element, is_axially_symmetric,
                                                 _integration_method);
 
        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            _ip_data.emplace_back(
                shape_matrices[ip].dNdx,
                _integration_method.getWeightedPoint(ip).getWeight() *
                    shape_matrices[ip].integralMeasure *
                    shape_matrices[ip].detJ * aperture_size);
        }
    }

References ProcessLib::HT::HTLocalAssemblerInterface::HTLocalAssemblerInterface(), _element, _integration_method, _ip_data, _process_data, NumLib::initShapeMatrices(), and ParameterLib::SpatialPosition::setElementID().

Referenced by ProcessLib::HT::MonolithicHTFEM< ShapeFunction, GlobalDim >::MonolithicHTFEM(), and ProcessLib::HT::StaggeredHTFEM< ShapeFunction, GlobalDim >::StaggeredHTFEM().

Member Function Documentation

◆ getFlux()

template<typename ShapeFunction, int GlobalDim>

Eigen::Vector3d ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::getFlux	(	MathLib::Point3d const &	pnt_local_coords,
		double const	t,
		std::vector< double > const &	local_x ) const

inlineoverridevirtual

Computes the flux in the point pnt_local_coords that is given in local coordinates using the values from local_x.

Implements ProcessLib::HT::HTLocalAssemblerInterface.

Definition at line 94 of file HTFEM.h.

    {
        // Eval shape matrices at given point
        // Note: Axial symmetry is set to false here, because we only need dNdx
        // here, which is not affected by axial symmetry.
        auto const shape_matrices =
            NumLib::computeShapeMatrices<ShapeFunction, ShapeMatricesType,
                                         GlobalDim>(
                _element, false /*is_axially_symmetric*/,
                std::array{pnt_local_coords})[0];
 
        ParameterLib::SpatialPosition pos;
        pos.setElementID(this->_element.getID());
 
        MaterialPropertyLib::VariableArray vars;
 
        // local_x contains the local temperature and pressure values
        double T_int_pt = 0.0;
        double p_int_pt = 0.0;
        NumLib::shapeFunctionInterpolate(local_x, shape_matrices.N, T_int_pt,
                                         p_int_pt);
 
        vars.temperature = T_int_pt;
        vars.liquid_phase_pressure = p_int_pt;
 
        auto const& medium =
            *_process_data.media_map.getMedium(_element.getID());
        auto const& liquid_phase = medium.phase("AqueousLiquid");
 
        // TODO (naumov) Temporary value not used by current material models.
        // Need extension of secondary variables interface.
        double const dt = std::numeric_limits<double>::quiet_NaN();
        // fetch permeability, viscosity, density
        auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
            medium.property(MaterialPropertyLib::PropertyType::permeability)
                .value(vars, pos, t, dt));
 
        auto const mu =
            liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
                .template value<double>(vars, pos, t, dt);
        GlobalDimMatrixType const K_over_mu = K / mu;
 
        auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
            &local_x[local_x.size() / 2], ShapeFunction::NPOINTS);
        GlobalDimVectorType q =
            -K_over_mu * shape_matrices.dNdx * p_nodal_values;
 
        if (this->_process_data.has_gravity)
        {
            auto const rho_w =
                liquid_phase
                    .property(MaterialPropertyLib::PropertyType::density)
                    .template value<double>(vars, pos, t, dt);
            auto const b =
                this->_process_data.projected_specific_body_force_vectors
                    [this->_element.getID()];
            q += K_over_mu * rho_w * b;
        }
 
        Eigen::Vector3d flux;
        flux.head<GlobalDim>() = q;
        return flux;
    }

References _element, _process_data, NumLib::computeShapeMatrices(), MaterialPropertyLib::density, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, ParameterLib::SpatialPosition::setElementID(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::temperature, and MaterialPropertyLib::viscosity.

◆ getHeatEnergyCoefficient()

template<typename ShapeFunction, int GlobalDim>

double ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::getHeatEnergyCoefficient	(	MaterialPropertyLib::VariableArray const &	vars,
		const double	porosity,
		const double	fluid_density,
		const double	specific_heat_capacity_fluid,
		ParameterLib::SpatialPosition const &	pos,
		double const	t,
		double const	dt )

inlineprotected

Definition at line 167 of file HTFEM.h.

    {
        auto const& medium =
            *_process_data.media_map.getMedium(this->_element.getID());
        auto const& solid_phase = medium.phase("Solid");
 
        auto const specific_heat_capacity_solid =
            solid_phase
                .property(
                    MaterialPropertyLib::PropertyType::specific_heat_capacity)
                .template value<double>(vars, pos, t, dt);
 
        auto const solid_density =
            solid_phase.property(MaterialPropertyLib::PropertyType::density)
                .template value<double>(vars, pos, t, dt);
 
        return solid_density * specific_heat_capacity_solid * (1 - porosity) +
               fluid_density * specific_heat_capacity_fluid * porosity;
    }

References _process_data, MaterialPropertyLib::density, MeshLib::Element::getID(), and MaterialPropertyLib::specific_heat_capacity.

Referenced by ProcessLib::HT::MonolithicHTFEM< ShapeFunction, GlobalDim >::assemble(), and ProcessLib::HT::StaggeredHTFEM< ShapeFunction, GlobalDim >::assembleHeatTransportEquation().

◆ getIntPtDarcyVelocityLocal()

template<typename ShapeFunction, int GlobalDim>

std::vector< double > const & ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::getIntPtDarcyVelocityLocal	(	const double	t,
		std::vector< double > const &	local_x,
		std::vector< double > &	cache ) const

inlineprotected

Definition at line 232 of file HTFEM.h.

    {
        std::vector<double> local_p{
            local_x.data() + pressure_index,
            local_x.data() + pressure_index + pressure_size};
        std::vector<double> local_T{
            local_x.data() + temperature_index,
            local_x.data() + temperature_index + temperature_size};
 
        auto const n_integration_points =
            _integration_method.getNumberOfPoints();
 
        cache.clear();
        auto cache_mat = MathLib::createZeroedMatrix<
            Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
            cache, GlobalDim, n_integration_points);
 
        ParameterLib::SpatialPosition pos;
        pos.setElementID(_element.getID());
 
        MaterialPropertyLib::VariableArray vars;
 
        auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
            &local_p[0], ShapeFunction::NPOINTS);
 
        auto const& medium =
            *_process_data.media_map.getMedium(_element.getID());
        auto const& liquid_phase = medium.phase("AqueousLiquid");
 
        auto const& Ns =
            _process_data.shape_matrix_cache
                .NsHigherOrder<typename ShapeFunction::MeshElement>();
 
        for (unsigned ip = 0; ip < n_integration_points; ++ip)
        {
            auto const& ip_data = _ip_data[ip];
            auto const& dNdx = ip_data.dNdx;
            auto const& N = Ns[ip];
 
            double T_int_pt = 0.0;
            double p_int_pt = 0.0;
            NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
            NumLib::shapeFunctionInterpolate(local_T, N, T_int_pt);
 
            vars.temperature = T_int_pt;
            vars.liquid_phase_pressure = p_int_pt;
 
            // TODO (naumov) Temporary value not used by current material
            // models. Need extension of secondary variables interface.
            double const dt = std::numeric_limits<double>::quiet_NaN();
            auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
                medium.property(MaterialPropertyLib::PropertyType::permeability)
                    .value(vars, pos, t, dt));
 
            auto const mu =
                liquid_phase
                    .property(MaterialPropertyLib::PropertyType::viscosity)
                    .template value<double>(vars, pos, t, dt);
            GlobalDimMatrixType const K_over_mu = K / mu;
 
            cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
 
            if (_process_data.has_gravity)
            {
                auto const rho_w =
                    liquid_phase
                        .property(MaterialPropertyLib::PropertyType::density)
                        .template value<double>(vars, pos, t, dt);
                auto const b =
                    _process_data.projected_specific_body_force_vectors
                        [_element.getID()];
                // here it is assumed that the vector b is directed 'downwards'
                cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
            }
        }
 
        return cache;
    }

References _element, _integration_method, _ip_data, _process_data, MathLib::createZeroedMatrix(), MaterialPropertyLib::density, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, pressure_index, pressure_size, ParameterLib::SpatialPosition::setElementID(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::temperature, temperature_index, temperature_size, and MaterialPropertyLib::viscosity.

Referenced by ProcessLib::HT::MonolithicHTFEM< ShapeFunction, GlobalDim >::getIntPtDarcyVelocity(), and ProcessLib::HT::StaggeredHTFEM< ShapeFunction, GlobalDim >::getIntPtDarcyVelocity().

◆ getShapeMatrix()

template<typename ShapeFunction, int GlobalDim>

Eigen::Map< const Eigen::RowVectorXd > ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::getShapeMatrix ( const unsigned integration_point ) const

inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 82 of file HTFEM.h.

    {
        auto const& N = _process_data.shape_matrix_cache.NsHigherOrder<
            typename ShapeFunction::MeshElement>()[integration_point];
 
        // assumes N is stored contiguously in memory
        return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
    }

References _process_data.

◆ getThermalConductivityDispersivity()

template<typename ShapeFunction, int GlobalDim>

GlobalDimMatrixType ProcessLib::HT::HTFEM< ShapeFunction, GlobalDim >::getThermalConductivityDispersivity	(	MaterialPropertyLib::VariableArray const &	vars,
		const double	fluid_density,
		const double	specific_heat_capacity_fluid,
		const GlobalDimVectorType &	velocity,
		ParameterLib::SpatialPosition const &	pos,
		double const	t,
		double const	dt )

inlineprotected

Definition at line 191 of file HTFEM.h.

    {
        auto const& medium =
            *_process_data.media_map.getMedium(_element.getID());
 
        auto thermal_conductivity =
            MaterialPropertyLib::formEigenTensor<GlobalDim>(
                medium
                    .property(
                        MaterialPropertyLib::PropertyType::thermal_conductivity)
                    .value(vars, pos, t, dt));
 
        auto const thermal_dispersivity_transversal =
            medium
                .property(MaterialPropertyLib::PropertyType::
                              thermal_transversal_dispersivity)
                .template value<double>();
 
        auto const thermal_dispersivity_longitudinal =
            medium
                .property(MaterialPropertyLib::PropertyType::
                              thermal_longitudinal_dispersivity)
                .template value<double>();
 
        // Thermal conductivity is moved outside and zero matrix is passed
        // instead due to multiplication with fluid's density times specific
        // heat capacity.
        return thermal_conductivity +
               fluid_density * specific_heat_capacity_fluid *
                   NumLib::computeHydrodynamicDispersion(
                       _process_data.stabilizer, _element.getID(),
                       GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
                       velocity, 0 /* phi */, thermal_dispersivity_transversal,
                       thermal_dispersivity_longitudinal);
    }