diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 7cec15a3911e5e74310ee16d36128ff3e5bcb803..be43cff10566024edc68fc8937becb067f4d2a34 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -28,6 +28,20 @@ namespace ProcessLib
{
namespace HT
{
+template < typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
+struct IntegrationPointData final
+{
+ IntegrationPointData(NodalRowVectorType const& N_,
+ GlobalDimNodalMatrixType const& dNdx_,
+ double const& integration_weight_)
+ : N(N_), dNdx(dNdx_), integration_weight(integration_weight_)
+ {}
+
+ NodalRowVectorType const N;
+ GlobalDimNodalMatrixType const dNdx;
+ double const integration_weight;
+};
+
const unsigned NUM_NODAL_DOF = 2;
class HTLocalAssemblerInterface
@@ -52,10 +66,19 @@ class LocalAssemblerData : public HTLocalAssemblerInterface
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
- using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+ 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 NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType;
+
using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+ using GlobalDimNodalMatrixType =
+ typename ShapeMatricesType::GlobalDimNodalMatrixType;
using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
public:
@@ -67,9 +90,6 @@ public:
: _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()))
@@ -78,6 +98,24 @@ public:
// matrices.
assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
(void)local_matrix_size;
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ _ip_data.reserve(n_integration_points);
+
+ auto const shape_matrices =
+ initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ element, is_axially_symmetric, _integration_method);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ _ip_data.emplace_back(
+ shape_matrices[ip].N, shape_matrices[ip].dNdx,
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ shape_matrices[ip].integralMeasure *
+ shape_matrices[ip].detJ);
+ }
}
void assemble(double const t, std::vector<double> const& local_x,
@@ -90,11 +128,11 @@ public:
// matrices.
assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
- auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
+ auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
local_M_data, local_matrix_size, local_matrix_size);
- auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
+ auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
local_K_data, local_matrix_size, local_matrix_size);
- auto local_b = MathLib::createZeroedVector<NodalVectorType>(
+ auto local_b = MathLib::createZeroedVector<LocalVectorType>(
local_b_data, local_matrix_size);
unsigned const n_integration_points =
@@ -131,11 +169,15 @@ public:
auto const thermal_conductivity_fluid =
_process_data.thermal_conductivity_fluid(t, pos)[0];
- auto const& sm = _shape_matrices[ip];
+ 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, sm.N, T_int_pt, p_int_pt);
+ NumLib::shapeFunctionInterpolate(local_x, N, T_int_pt, p_int_pt);
// TODO the first argument has to be changed for non constant
// porosity model
@@ -156,7 +198,6 @@ public:
auto const thermal_dispersivity_transversal =
_process_data.thermal_dispersivity_transversal(t, pos)[0];
- 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,
@@ -180,8 +221,7 @@ public:
intrinsic_permeability / viscosity;
GlobalDimVectorType const velocity =
- -perm_over_visc *
- (sm.dNdx * p_nodal_values - density_water_T * b);
+ -perm_over_visc * (dNdx * p_nodal_values - density_water_T * b);
double const velocity_magnitude = velocity.norm();
GlobalDimMatrixType const& I(
@@ -201,22 +241,16 @@ public:
density_solid * specific_heat_capacity_solid * (1 - porosity) +
fluid_reference_density * specific_heat_capacity_fluid * porosity;
- auto const integral_term =
- sm.integralMeasure * sm.detJ * wp.getWeight();
// matrix assembly
Ktt.noalias() +=
- integral_term *
- (sm.dNdx.transpose() * hydrodynamic_thermodispersion * sm.dNdx +
- sm.N.transpose() * velocity.transpose() * sm.dNdx *
- fluid_reference_density * specific_heat_capacity_fluid);
- Kpp.noalias() +=
- integral_term * sm.dNdx.transpose() * perm_over_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 += integral_term * density_water_T * sm.dNdx.transpose() *
- perm_over_visc * b;
+ (dNdx.transpose() * hydrodynamic_thermodispersion * dNdx +
+ N.transpose() * velocity.transpose() * dNdx *
+ fluid_reference_density * specific_heat_capacity_fluid) *
+ w;
+ Kpp.noalias() += w * dNdx.transpose() * perm_over_visc * dNdx;
+ Mtt.noalias() += w * N.transpose() * heat_capacity * N;
+ Mpp.noalias() += w * N.transpose() * specific_storage * N;
+ Bp += w * density_water_T * dNdx.transpose() * perm_over_visc * b;
/* with Oberbeck-Boussing assumption density difference only exists
* in buoyancy effects */
@@ -231,7 +265,7 @@ public:
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
- auto const& N = _shape_matrices[integration_point].N;
+ auto const& N = _ip_data[integration_point].N;
// assumes N is stored contiguously in memory
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
@@ -263,7 +297,9 @@ private:
HTProcessData const& _process_data;
IntegrationMethod const _integration_method;
- std::vector<ShapeMatrices> _shape_matrices;
+ std::vector<
+ IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>>
+ _ip_data;
std::vector<std::vector<double>> _darcy_velocities;
};