33template <
typename ShapeFunction,
int GlobalDim>
49 std::size_t
const local_matrix_size,
51 bool const is_axially_symmetric,
53 const unsigned dof_per_node)
61 assert(local_matrix_size == ShapeFunction::NPOINTS * dof_per_node);
62 (void)local_matrix_size;
65 unsigned const n_integration_points =
67 _ip_data.reserve(n_integration_points);
74 auto const shape_matrices =
76 GlobalDim>(element, is_axially_symmetric,
79 for (
unsigned ip = 0; ip < n_integration_points; ip++)
82 shape_matrices[ip].N, shape_matrices[ip].dNdx,
84 shape_matrices[ip].integralMeasure *
85 shape_matrices[ip].detJ * aperture_size);
90 const unsigned integration_point)
const override
92 auto const& N =
_ip_data[integration_point].N;
95 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
102 std::vector<double>
const& local_x)
const override
107 auto const shape_matrices =
111 std::array{pnt_local_coords})[0];
119 double T_int_pt = 0.0;
120 double p_int_pt = 0.0;
129 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
133 double const dt = std::numeric_limits<double>::quiet_NaN();
135 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
137 .value(vars, pos, t, dt));
141 .template value<double>(vars, pos, t, dt);
144 auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
145 &local_x[local_x.size() / 2], ShapeFunction::NPOINTS);
147 -K_over_mu * shape_matrices.dNdx * p_nodal_values;
154 .template value<double>(vars, pos, t, dt);
157 [this->_element.getID()];
158 q += K_over_mu * rho_w * b;
161 Eigen::Vector3d flux;
162 flux.head<GlobalDim>() = q;
173 Eigen::aligned_allocator<
179 const double fluid_density,
const double specific_heat_capacity_fluid,
185 auto const& solid_phase = medium.
phase(
"Solid");
187 auto const specific_heat_capacity_solid =
191 .template value<double>(vars, pos, t, dt);
193 auto const solid_density =
195 .template value<double>(vars, pos, t, dt);
197 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
198 fluid_density * specific_heat_capacity_fluid * porosity;
203 const double fluid_density,
const double specific_heat_capacity_fluid,
211 auto thermal_conductivity =
212 MaterialPropertyLib::formEigenTensor<GlobalDim>(
216 .value(vars, pos, t, dt));
218 auto const thermal_dispersivity_transversal =
221 thermal_transversal_dispersivity)
222 .template value<double>();
224 auto const thermal_dispersivity_longitudinal =
227 thermal_longitudinal_dispersivity)
228 .template value<double>();
233 return thermal_conductivity +
234 fluid_density * specific_heat_capacity_fluid *
237 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
238 velocity, 0 , thermal_dispersivity_transversal,
239 thermal_dispersivity_longitudinal);
243 const double t, std::vector<double>
const& local_x,
244 std::vector<double>& cache)
const
246 std::vector<double> local_p{
249 std::vector<double> local_T{
253 auto const n_integration_points =
258 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
259 cache, GlobalDim, n_integration_points);
266 auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
267 &local_p[0], ShapeFunction::NPOINTS);
271 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
273 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
276 auto const& N = ip_data.N;
277 auto const& dNdx = ip_data.dNdx;
279 pos.setIntegrationPoint(ip);
281 double T_int_pt = 0.0;
282 double p_int_pt = 0.0;
291 double const dt = std::numeric_limits<double>::quiet_NaN();
292 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
294 .value(vars, pos, t, dt));
299 .template value<double>(vars, pos, t, dt);
302 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
309 .template value<double>(vars, pos, t, dt);
314 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
Medium * getMedium(std::size_t element_id)
Phase const & phase(std::size_t index) const
Property const & property(PropertyType const &p) const
virtual std::size_t getID() const final
Returns the ID of the element.
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
unsigned getNumberOfPoints() const
void setElementID(std::size_t element_id)
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
std::vector< double > const & getIntPtDarcyVelocityLocal(const double t, std::vector< double > const &local_x, std::vector< double > &cache) const
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
NumLib::GenericIntegrationMethod const & _integration_method
static const int temperature_index
typename ShapeMatricesType::NodalVectorType NodalVectorType
static const int temperature_size
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)
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)
static const int pressure_size
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
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)
Eigen::Vector3d getFlux(MathLib::Point3d const &pnt_local_coords, double const t, std::vector< double > const &local_x) const override
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
MeshLib::Element const & _element
HTProcessData const & _process_data
typename ShapeMatricesType::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
static const int pressure_index
std::vector< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType > > > _ip_data
typename ShapeMatricesType::NodalRowVectorType NodalRowVectorType
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)
Eigen::MatrixXd computeHydrodynamicDispersion(NumericalStabilization const &stabilizer, std::size_t const element_id, Eigen::MatrixXd const &pore_diffusion_coefficient, Eigen::VectorXd const &velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal)
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
MatrixType< GlobalDim, ShapeFunction::NPOINTS > GlobalDimNodalMatrixType
MatrixType< GlobalDim, GlobalDim > GlobalDimMatrixType
VectorType< GlobalDim > GlobalDimVectorType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
ParameterLib::Parameter< double > const & aperture_size
MaterialPropertyLib::MaterialSpatialDistributionMap media_map
NumLib::NumericalStabilization stabilizer
std::vector< Eigen::VectorXd > const projected_specific_body_force_vectors
Projected specific body force vector: R * R^T * b.