diff --git a/NumLib/NumericalStability/AdvectionMatrixAssembler.h b/NumLib/NumericalStability/AdvectionMatrixAssembler.h
index 3d541b5c048d4353da8880728222db07672052ea..ed5c09964d7d2d33bb0008c8dc12b234e887ffe4 100644
--- a/NumLib/NumericalStability/AdvectionMatrixAssembler.h
+++ b/NumLib/NumericalStability/AdvectionMatrixAssembler.h
@@ -14,12 +14,35 @@
#include <memory>
#include <vector>
+#include "NumLib/Fem/ShapeMatrixCache.h"
#include "NumericalStabilization.h"
namespace NumLib
{
namespace detail
{
+template <typename MeshElementType,
+ typename IPData,
+ typename FluxVectorType,
+ typename Derived>
+void assembleAdvectionMatrix(IPData const& ip_data_vector,
+ NumLib::ShapeMatrixCache const& shape_matrix_cache,
+ std::vector<FluxVectorType> const& ip_flux_vector,
+ Eigen::MatrixBase<Derived>& laplacian_matrix)
+{
+ auto const& Ns = shape_matrix_cache.NsHigherOrder<MeshElementType>();
+
+ for (std::size_t ip = 0; ip < ip_flux_vector.size(); ++ip)
+ {
+ auto const& ip_data = ip_data_vector[ip];
+ auto const w = ip_data.integration_weight;
+ auto const& dNdx = ip_data.dNdx;
+ auto const& N = Ns[ip];
+ laplacian_matrix.noalias() +=
+ N.transpose() * ip_flux_vector[ip].transpose() * dNdx * w;
+ }
+}
+
template <typename IPData, typename FluxVectorType, typename Derived>
void assembleAdvectionMatrix(IPData const& ip_data_vector,
std::vector<FluxVectorType> const& ip_flux_vector,
@@ -79,9 +102,13 @@ void applyFullUpwind(IPData const& ip_data_vector,
}
} // namespace detail
-template <typename IPData, typename FluxVectorType, typename Derived>
+template <typename MeshElementType,
+ typename IPData,
+ typename FluxVectorType,
+ typename Derived>
void assembleAdvectionMatrix(NumericalStabilization const& stabilizer,
IPData const& ip_data_vector,
+ NumLib::ShapeMatrixCache const& shape_matrix_cache,
std::vector<FluxVectorType> const& ip_flux_vector,
double const average_velocity,
Eigen::MatrixBase<Derived>& laplacian_matrix)
@@ -100,8 +127,36 @@ void assembleAdvectionMatrix(NumericalStabilization const& stabilizer,
}
}
- detail::assembleAdvectionMatrix(ip_data_vector, ip_flux_vector,
+ detail::assembleAdvectionMatrix<MeshElementType>(
+ ip_data_vector, shape_matrix_cache, ip_flux_vector,
+ laplacian_matrix);
+ },
+ stabilizer);
+}
+
+template <typename IPData, typename FluxVectorType, typename Derived>
+void assembleAdvectionMatrix(NumericalStabilization const& stabilizer,
+ IPData const& ip_data_vector,
+ std::vector<FluxVectorType> const& ip_flux_vector,
+ double const average_velocity,
+ Eigen::MatrixBase<Derived>& laplacian_matrix)
+{
+ std::visit(
+ [&](auto&& stabilizer)
+ {
+ using Stabilizer = std::decay_t<decltype(stabilizer)>;
+ if constexpr (std::is_same_v<Stabilizer, FullUpwind>)
+ {
+ if (average_velocity > stabilizer.getCutoffVelocity())
+ {
+ detail::applyFullUpwind(ip_data_vector, ip_flux_vector,
laplacian_matrix);
+ return;
+ }
+ }
+
+ detail::assembleAdvectionMatrix(
+ ip_data_vector, ip_flux_vector, laplacian_matrix);
},
stabilizer);
}
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 02ff5f8abff59a800bd84de6b475f496eb909a51..d9cedbbc2f5ebc25330f4cfbaf9f6b81f955514f 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -25,6 +25,7 @@
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/Integration/GenericIntegrationMethod.h"
+#include "NumLib/Fem/ShapeMatrixCache.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "NumLib/Function/Interpolation.h"
#include "NumLib/NumericalStability/AdvectionMatrixAssembler.h"
@@ -37,18 +38,16 @@ namespace ProcessLib
{
namespace ComponentTransport
{
-template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
+template <typename GlobalDimNodalMatrixType>
struct IntegrationPointData final
{
- IntegrationPointData(NodalRowVectorType const& N_,
- GlobalDimNodalMatrixType const& dNdx_,
+ IntegrationPointData(GlobalDimNodalMatrixType const& dNdx_,
double const& integration_weight_)
- : N(N_), dNdx(dNdx_), integration_weight(integration_weight_)
+ : dNdx(dNdx_), integration_weight(integration_weight_)
{
}
void pushBackState() { porosity_prev = porosity; }
- NodalRowVectorType const N;
GlobalDimNodalMatrixType const dNdx;
double const integration_weight;
@@ -278,7 +277,7 @@ public:
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
_ip_data.emplace_back(
- shape_matrices[ip].N, shape_matrices[ip].dNdx,
+ shape_matrices[ip].dNdx,
_integration_method.getWeightedPoint(ip).getWeight() *
shape_matrices[ip].integralMeasure *
shape_matrices[ip].detJ * aperture_size);
@@ -325,12 +324,17 @@ public:
ParameterLib::SpatialPosition pos;
pos.setElementID(_element.getID());
+ auto const& Ns =
+ _process_data.shape_matrix_cache
+ .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
+
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
auto& ip_data = _ip_data[ip];
- auto const& N = ip_data.N;
+ auto const& N = Ns[ip];
auto const& chemical_system_id = ip_data.chemical_system_id;
auto const n_component = _transport_process_variables.size();
@@ -369,12 +373,17 @@ public:
ParameterLib::SpatialPosition pos;
pos.setElementID(_element.getID());
+ auto const& Ns =
+ _process_data.shape_matrix_cache
+ .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
+
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
auto& ip_data = _ip_data[ip];
- auto const& N = ip_data.N;
+ auto const& N = Ns[ip];
auto& porosity = ip_data.porosity;
auto const& porosity_prev = ip_data.porosity_prev;
auto const& chemical_system_id = ip_data.chemical_system_id;
@@ -592,13 +601,17 @@ public:
ip_flux_vector.reserve(n_integration_points);
}
+ auto const& Ns =
+ _process_data.shape_matrix_cache
+ .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
for (unsigned ip(0); ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);
auto& ip_data = _ip_data[ip];
- auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
+ auto const& N = Ns[ip];
auto const& w = ip_data.integration_weight;
auto& porosity = ip_data.porosity;
@@ -719,9 +732,11 @@ public:
if (!_process_data.non_advective_form)
{
- NumLib::assembleAdvectionMatrix(
+ NumLib::assembleAdvectionMatrix<
+ typename ShapeFunction::MeshElement>(
_process_data.stabilizer,
_ip_data,
+ _process_data.shape_matrix_cache,
ip_flux_vector,
average_velocity_norm /
static_cast<double>(n_integration_points),
@@ -750,11 +765,15 @@ public:
auto const& component = phase.component(
_transport_process_variables[component_id].get().getName());
+ auto const& Ns =
+ _process_data.shape_matrix_cache
+ .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
for (unsigned ip(0); ip < n_integration_points; ++ip)
{
auto& ip_data = _ip_data[ip];
- auto const& N = ip_data.N;
auto const& w = ip_data.integration_weight;
+ auto const& N = Ns[ip];
auto& porosity = ip_data.porosity;
auto const retardation_factor =
@@ -843,14 +862,18 @@ public:
MaterialPropertyLib::VariableArray vars;
MaterialPropertyLib::VariableArray vars_prev;
+ auto const& Ns =
+ _process_data.shape_matrix_cache
+ .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
for (unsigned ip(0); ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);
auto& ip_data = _ip_data[ip];
- auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
auto const& w = ip_data.integration_weight;
+ auto const& N = Ns[ip];
auto& porosity = ip_data.porosity;
auto const& porosity_prev = ip_data.porosity_prev;
@@ -967,14 +990,18 @@ public:
double average_velocity_norm = 0.0;
ip_flux_vector.reserve(n_integration_points);
+ auto const& Ns =
+ _process_data.shape_matrix_cache
+ .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
for (unsigned ip(0); ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);
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;
+ auto const& N = Ns[ip];
double p_at_xi = 0.;
NumLib::shapeFunctionInterpolate(local_p, N, p_at_xi);
@@ -1043,8 +1070,9 @@ public:
average_velocity_norm += velocity.norm();
}
- NumLib::assembleAdvectionMatrix(
- process_data.stabilizer, this->_ip_data, ip_flux_vector,
+ NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
+ process_data.stabilizer, this->_ip_data,
+ _process_data.shape_matrix_cache, ip_flux_vector,
average_velocity_norm / static_cast<double>(n_integration_points),
local_K);
}
@@ -1109,14 +1137,18 @@ public:
auto const& component = phase.component(
_transport_process_variables[component_id].get().getName());
+ auto const& Ns =
+ _process_data.shape_matrix_cache
+ .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
for (unsigned ip(0); ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);
auto& ip_data = _ip_data[ip];
- auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
auto const& w = ip_data.integration_weight;
+ auto const& N = Ns[ip];
auto& porosity = ip_data.porosity;
auto const& porosity_prev = ip_data.porosity_prev;
@@ -1240,10 +1272,10 @@ public:
if (!_process_data.non_advective_form)
{
- NumLib::assembleAdvectionMatrix(
- _process_data.stabilizer,
- _ip_data,
- ip_flux_vector,
+ NumLib::assembleAdvectionMatrix<
+ typename ShapeFunction::MeshElement>(
+ _process_data.stabilizer, _ip_data,
+ _process_data.shape_matrix_cache, ip_flux_vector,
average_velocity_norm /
static_cast<double>(n_integration_points),
KCC_Laplacian);
@@ -1307,14 +1339,18 @@ public:
MaterialPropertyLib::VariableArray vars;
MaterialPropertyLib::VariableArray vars_prev;
+ auto const& Ns =
+ _process_data.shape_matrix_cache
+ .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
for (unsigned ip(0); ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);
auto& ip_data = _ip_data[ip];
- auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
auto const& w = ip_data.integration_weight;
+ auto const& N = Ns[ip];
auto& phi = ip_data.porosity;
auto const& phi_prev = ip_data.porosity_prev;
@@ -1429,14 +1465,18 @@ public:
auto const& component = phase.component(
_transport_process_variables[component_id].get().getName());
+ auto const& Ns =
+ _process_data.shape_matrix_cache
+ .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
for (unsigned ip(0); ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);
auto& ip_data = _ip_data[ip];
- auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
auto const& w = ip_data.integration_weight;
+ auto const& N = Ns[ip];
auto& phi = ip_data.porosity;
auto const& phi_prev = ip_data.porosity_prev;
@@ -1513,8 +1553,9 @@ public:
average_velocity_norm += q.norm();
}
- NumLib::assembleAdvectionMatrix(
- _process_data.stabilizer, _ip_data, ip_flux_vector,
+ NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
+ _process_data.stabilizer, _ip_data,
+ _process_data.shape_matrix_cache, ip_flux_vector,
average_velocity_norm / static_cast<double>(n_integration_points),
KCC_Laplacian);
@@ -1552,13 +1593,18 @@ public:
auto const& medium =
*_process_data.media_map.getMedium(_element.getID());
auto const component_id = transport_process_id - 1;
+
+ auto const& Ns =
+ _process_data.shape_matrix_cache
+ .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
for (unsigned ip(0); ip < n_integration_points; ++ip)
{
pos.setIntegrationPoint(ip);
auto& ip_data = _ip_data[ip];
- auto const& N = ip_data.N;
auto const w = ip_data.integration_weight;
+ auto const& N = Ns[ip];
auto& porosity = ip_data.porosity;
auto const& porosity_prev = ip_data.porosity_prev;
auto const chemical_system_id = ip_data.chemical_system_id;
@@ -1669,11 +1715,15 @@ public:
*_process_data.media_map.getMedium(_element.getID());
auto const& 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& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
+ auto const& N = Ns[ip];
auto const& porosity = ip_data.porosity;
pos.setIntegrationPoint(ip);
@@ -1715,7 +1765,8 @@ public:
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
- auto const& N = _ip_data[integration_point].N;
+ 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());
@@ -1896,11 +1947,15 @@ public:
auto const& component = phase.component(
_transport_process_variables[component_id].get().getName());
+ 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& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
+ auto const& N = Ns[ip];
auto const& phi = ip_data.porosity;
pos.setIntegrationPoint(ip);
@@ -1969,11 +2024,7 @@ private:
std::vector<std::reference_wrapper<ProcessVariable>> const
_transport_process_variables;
- std::vector<
- IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>,
- Eigen::aligned_allocator<
- IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>>>
- _ip_data;
+ std::vector<IntegrationPointData<GlobalDimNodalMatrixType>> _ip_data;
double getHeatEnergyCoefficient(
MaterialPropertyLib::VariableArray const& vars, const double porosity,
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index 2e81003a219197df83eb54f935964cf7b7ed1243..74cfee620f7af6b2cbb5879de56038548e53bd07 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -23,6 +23,7 @@
#include "MathLib/LinAlg/LinAlg.h"
#include "MeshLib/Utils/getOrCreateMeshProperty.h"
#include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "NumLib/Fem/ShapeMatrixCache.h"
#include "NumLib/NumericalStability/FluxCorrectedTransport.h"
#include "NumLib/NumericalStability/NumericalStabilization.h"
#include "NumLib/ODESolver/NonlinearSystem.h"
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcessData.h b/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
index 2ac5512eed9f6176ad3603b62eb720ce7e711ac8..121c3e7641c99d6d8de72414128a30b927670dc0 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
@@ -17,6 +17,7 @@
#include "LookupTable.h"
#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/Fem/ShapeMatrixCache.h"
#include "NumLib/NumericalStability/NumericalStabilization.h"
#include "ParameterLib/ConstantParameter.h"
#include "ParameterLib/Parameter.h"
@@ -84,6 +85,9 @@ struct ComponentTransportProcessData
bool const isothermal;
+ /// caches for each mesh element type the shape matrix
+ NumLib::ShapeMatrixCache shape_matrix_cache;
+
static const int hydraulic_process_id = 0;
// Thermal process is optional, indicated by -1. If present, it is positive.
const int thermal_process_id = isothermal ? -1 : 1;
diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
index bb69b4fe94eaf9f6296d9c50d96408955fc8e1de..5710f1e4d518c561b22d13102ef3f32c05786742 100644
--- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
@@ -317,7 +317,7 @@ std::unique_ptr<Process> createComponentTransportProcess(
mesh_space_dimension,
*aperture_size_parameter,
isothermal,
- };
+ NumLib::ShapeMatrixCache{integration_order}};
SecondaryVariableCollection secondary_variables;
diff --git a/ProcessLib/HT/CreateHTProcess.cpp b/ProcessLib/HT/CreateHTProcess.cpp
index f8cc32419791b3aec6bdbe3a1be7360ef1cf54a1..8a5a4b6d44d6007740c1491c125959a083b23185 100644
--- a/ProcessLib/HT/CreateHTProcess.cpp
+++ b/ProcessLib/HT/CreateHTProcess.cpp
@@ -218,7 +218,8 @@ std::unique_ptr<Process> createHTProcess(
std::move(stabilizer),
projected_specific_body_force_vectors,
mesh_space_dimension,
- *aperture_size_parameter};
+ *aperture_size_parameter,
+ NumLib::ShapeMatrixCache{integration_order}};
SecondaryVariableCollection secondary_variables;
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index da14725acc227e43f14588f20cda86ea16f12ea9..f0e08207fa6139639579869f136c3952dcfe4c88 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -79,7 +79,7 @@ public:
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
_ip_data.emplace_back(
- shape_matrices[ip].N, shape_matrices[ip].dNdx,
+ shape_matrices[ip].dNdx,
_integration_method.getWeightedPoint(ip).getWeight() *
shape_matrices[ip].integralMeasure *
shape_matrices[ip].detJ * aperture_size);
@@ -89,7 +89,8 @@ public:
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
- auto const& N = _ip_data[integration_point].N;
+ 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());
@@ -168,11 +169,7 @@ protected:
HTProcessData const& _process_data;
NumLib::GenericIntegrationMethod const& _integration_method;
- std::vector<
- IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>,
- Eigen::aligned_allocator<
- IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>>>
- _ip_data;
+ std::vector<IntegrationPointData<GlobalDimNodalMatrixType>> _ip_data;
double getHeatEnergyCoefficient(
MaterialPropertyLib::VariableArray const& vars, const double porosity,
@@ -270,11 +267,15 @@ protected:
*_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& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
+ auto const& N = Ns[ip];
pos.setIntegrationPoint(ip);
diff --git a/ProcessLib/HT/HTLocalAssemblerInterface.h b/ProcessLib/HT/HTLocalAssemblerInterface.h
index c058a23b5c1fa1fbd17f241065aba9654359b314..8fcd126c7deea126c4e659ba6f7699d19b667078 100644
--- a/ProcessLib/HT/HTLocalAssemblerInterface.h
+++ b/ProcessLib/HT/HTLocalAssemblerInterface.h
@@ -21,19 +21,15 @@ namespace ProcessLib
namespace HT
{
-template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
+template <typename GlobalDimNodalMatrixType>
struct IntegrationPointData final
{
- IntegrationPointData(NodalRowVectorType N_,
- GlobalDimNodalMatrixType dNdx_,
+ IntegrationPointData(GlobalDimNodalMatrixType dNdx_,
double const& integration_weight_)
- : N(std::move(N_)),
- dNdx(std::move(dNdx_)),
- integration_weight(integration_weight_)
+ : dNdx(std::move(dNdx_)), integration_weight(integration_weight_)
{
}
- NodalRowVectorType const N;
GlobalDimNodalMatrixType const dNdx;
double const integration_weight;
diff --git a/ProcessLib/HT/HTProcessData.h b/ProcessLib/HT/HTProcessData.h
index f12e2c81c7905b0b9f25bfc7801ca08c69026ba1..c951cf6b98911dd67f1f0a20dd9a351a06ab1efe 100644
--- a/ProcessLib/HT/HTProcessData.h
+++ b/ProcessLib/HT/HTProcessData.h
@@ -14,6 +14,7 @@
#include <utility>
#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+#include "NumLib/Fem/ShapeMatrixCache.h"
#include "NumLib/NumericalStability/NumericalStabilization.h"
#include "ParameterLib/ConstantParameter.h"
#include "ParameterLib/Parameter.h"
@@ -43,6 +44,9 @@ struct HTProcessData final
/// cross section area of 1D element. For 3D element, the value is set to 1.
ParameterLib::Parameter<double> const& aperture_size =
ParameterLib::ConstantParameter<double>("constant_one", 1.0);
+
+ /// caches for each mesh element type the shape matrix
+ NumLib::ShapeMatrixCache shape_matrix_cache;
};
} // namespace HT
diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h
index 902036df8d9c41f16a2f95148a2e2880ab7fdf4b..91ad82498373ee92115067293de0211fc2fda1cb 100644
--- a/ProcessLib/HT/MonolithicHTFEM.h
+++ b/ProcessLib/HT/MonolithicHTFEM.h
@@ -121,13 +121,17 @@ public:
double average_velocity_norm = 0.0;
ip_flux_vector.reserve(n_integration_points);
+ auto const& Ns =
+ process_data.shape_matrix_cache
+ .template NsHigherOrder<typename ShapeFunction::MeshElement>();
+
for (unsigned ip(0); ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);
auto const& ip_data = this->_ip_data[ip];
- auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
+ auto const& N = Ns[ip];
auto const& w = ip_data.integration_weight;
double T_int_pt = 0.0;
@@ -210,8 +214,9 @@ public:
* in buoyancy effects */
}
- NumLib::assembleAdvectionMatrix(
- process_data.stabilizer, this->_ip_data, ip_flux_vector,
+ NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
+ process_data.stabilizer, this->_ip_data,
+ process_data.shape_matrix_cache, ip_flux_vector,
average_velocity_norm / static_cast<double>(n_integration_points),
KTT);
}
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index 65c67a68b1ebc0dbfcf86acbc8aa760b743bfdc0..7f92bddc13beca6190f8a2f1bb7fbe2891b29bdb 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -79,13 +79,17 @@ void StaggeredHTFEM<ShapeFunction, GlobalDim>::assembleHydraulicEquation(
unsigned const n_integration_points =
this->_integration_method.getNumberOfPoints();
+ auto const& Ns =
+ process_data.shape_matrix_cache
+ .template NsHigherOrder<typename ShapeFunction::MeshElement>();
+
for (unsigned ip(0); ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);
auto const& ip_data = this->_ip_data[ip];
- auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
+ auto const& N = Ns[ip];
auto const& w = ip_data.integration_weight;
double p_int_pt = 0.0;
@@ -208,13 +212,17 @@ void StaggeredHTFEM<ShapeFunction, GlobalDim>::assembleHeatTransportEquation(
double average_velocity_norm = 0.0;
ip_flux_vector.reserve(n_integration_points);
+ auto const& Ns =
+ process_data.shape_matrix_cache
+ .template NsHigherOrder<typename ShapeFunction::MeshElement>();
+
for (unsigned ip(0); ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);
auto const& ip_data = this->_ip_data[ip];
- auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
+ auto const& N = Ns[ip];
auto const& w = ip_data.integration_weight;
double p_at_xi = 0.;
@@ -279,8 +287,9 @@ void StaggeredHTFEM<ShapeFunction, GlobalDim>::assembleHeatTransportEquation(
average_velocity_norm += velocity.norm();
}
- NumLib::assembleAdvectionMatrix(
- process_data.stabilizer, this->_ip_data, ip_flux_vector,
+ NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
+ process_data.stabilizer, this->_ip_data,
+ process_data.shape_matrix_cache, ip_flux_vector,
average_velocity_norm / static_cast<double>(n_integration_points),
local_K);
}
diff --git a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp
index c25acd8344f3b6e1a13709d0b20f0927cbe3a999..c494c08b93df3a97cab43d23f6432e7b9a509f38 100644
--- a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp
@@ -152,7 +152,8 @@ std::unique_ptr<Process> createLiquidFlowProcess(
mesh_space_dimension,
std::move(specific_body_force),
has_gravity,
- *aperture_size_parameter};
+ *aperture_size_parameter,
+ NumLib::ShapeMatrixCache{integration_order}};
return std::make_unique<LiquidFlowProcess>(
std::move(name), mesh, std::move(jacobian_assembler), parameters,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowData.h b/ProcessLib/LiquidFlow/LiquidFlowData.h
index b95d1dbad784f9c16a084da2adc273f5ebb6c4ad..b366ccaf2dec593f155fec1db29ebb9fff930e22 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowData.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowData.h
@@ -13,6 +13,7 @@
#include <Eigen/Core>
#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+#include "NumLib/Fem/ShapeMatrixCache.h"
#include "ParameterLib/Parameter.h"
namespace ProcessLib
@@ -34,6 +35,9 @@ struct LiquidFlowData final
/// It stores aperture size, which is the thickness of 2D element or the
/// cross section area of 1D element. For 3D element, the value is set to 1.
ParameterLib::Parameter<double> const& aperture_size;
+
+ /// caches for each mesh element type the shape matrix
+ NumLib::ShapeMatrixCache shape_matrix_cache;
};
} // namespace LiquidFlow
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 4b3e41fffcc5f1655d1d2bb76dd690e035cd2e3f..99bc458dfdc0408c641f56203e39b390f2a9b4c5 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -143,15 +143,16 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
_process_data.element_rotation_matrices[_element.getID()].transpose() *
_process_data.specific_body_force;
- auto const& N = _shape_matrix_cache
- .NsHigherOrder<typename ShapeFunction::MeshElement>();
+ 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& N = Ns[ip];
double p = 0.;
- NumLib::shapeFunctionInterpolate(local_x, N[ip], p);
+ NumLib::shapeFunctionInterpolate(local_x, N, p);
vars.liquid_phase_pressure = p;
// Compute density:
@@ -176,7 +177,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
// Assemble mass matrix, M
local_M.noalias() +=
(porosity * ddensity_dpressure / fluid_density + storage) *
- N[ip].transpose() * N[ip] * ip_data.integration_weight;
+ N.transpose() * N * ip_data.integration_weight;
// Compute viscosity:
auto const viscosity =
@@ -295,14 +296,15 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
_process_data.element_rotation_matrices[_element.getID()].transpose() *
_process_data.specific_body_force;
- auto const& N = _shape_matrix_cache
- .NsHigherOrder<typename ShapeFunction::MeshElement>();
+ 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& N = Ns[ip];
double p = 0.;
- NumLib::shapeFunctionInterpolate(local_x, N[ip], p);
+ NumLib::shapeFunctionInterpolate(local_x, N, p);
vars.liquid_phase_pressure = p;
// Compute density:
@@ -333,8 +335,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::IsotropicCalculator::
calculateLaplacianAndGravityTerm(
Eigen::Map<NodalMatrixType>& local_K,
Eigen::Map<NodalVectorType>& local_b,
- IntegrationPointData<NodalRowVectorType,
- GlobalDimNodalMatrixType> const& ip_data,
+ IntegrationPointData<GlobalDimNodalMatrixType> const& ip_data,
GlobalDimMatrixType const& permeability, double const mu,
double const rho_L, GlobalDimVectorType const& specific_body_force,
bool const has_gravity)
@@ -355,8 +356,7 @@ Eigen::Matrix<double, GlobalDim, 1>
LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::IsotropicCalculator::
calculateVelocity(
Eigen::Map<const NodalVectorType> const& local_p,
- IntegrationPointData<NodalRowVectorType,
- GlobalDimNodalMatrixType> const& ip_data,
+ IntegrationPointData<GlobalDimNodalMatrixType> const& ip_data,
GlobalDimMatrixType const& permeability, double const mu,
double const rho_L, GlobalDimVectorType const& specific_body_force,
bool const has_gravity)
@@ -377,8 +377,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::AnisotropicCalculator::
calculateLaplacianAndGravityTerm(
Eigen::Map<NodalMatrixType>& local_K,
Eigen::Map<NodalVectorType>& local_b,
- IntegrationPointData<NodalRowVectorType,
- GlobalDimNodalMatrixType> const& ip_data,
+ IntegrationPointData<GlobalDimNodalMatrixType> const& ip_data,
GlobalDimMatrixType const& permeability, double const mu,
double const rho_L, GlobalDimVectorType const& specific_body_force,
bool const has_gravity)
@@ -399,8 +398,7 @@ Eigen::Matrix<double, GlobalDim, 1>
LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::AnisotropicCalculator::
calculateVelocity(
Eigen::Map<const NodalVectorType> const& local_p,
- IntegrationPointData<NodalRowVectorType,
- GlobalDimNodalMatrixType> const& ip_data,
+ IntegrationPointData<GlobalDimNodalMatrixType> const& ip_data,
GlobalDimMatrixType const& permeability, double const mu,
double const rho_L, GlobalDimVectorType const& specific_body_force,
bool const has_gravity)
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index 86ae1838af3d58f8a78f66994b12aba1bb890e03..dcdfeb934e40ab9146fadf230c05059aa2a18957 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -34,7 +34,7 @@ namespace ProcessLib
{
namespace LiquidFlow
{
-template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
+template <typename GlobalDimNodalMatrixType>
struct IntegrationPointData final
{
explicit IntegrationPointData(GlobalDimNodalMatrixType const& dNdx_,
@@ -86,12 +86,10 @@ public:
std::size_t const /*local_matrix_size*/,
NumLib::GenericIntegrationMethod const& integration_method,
bool const is_axially_symmetric,
- LiquidFlowData const& process_data,
- NumLib::ShapeMatrixCache const& shape_matrix_cache)
+ LiquidFlowData const& process_data)
: _element(element),
_integration_method(integration_method),
- _process_data(process_data),
- _shape_matrix_cache(shape_matrix_cache)
+ _process_data(process_data)
{
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -137,7 +135,7 @@ public:
const unsigned integration_point) const override
{
auto const& N =
- _shape_matrix_cache
+ _process_data.shape_matrix_cache
.NsHigherOrder<typename ShapeFunction::MeshElement>();
// assumes N is stored contiguously in memory
@@ -155,11 +153,7 @@ private:
MeshLib::Element const& _element;
NumLib::GenericIntegrationMethod const& _integration_method;
- std::vector<
- IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>,
- Eigen::aligned_allocator<
- IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>>>
- _ip_data;
+ std::vector<IntegrationPointData<GlobalDimNodalMatrixType>> _ip_data;
/**
* Calculator of the Laplacian and the gravity term for anisotropic
@@ -170,16 +164,14 @@ private:
static void calculateLaplacianAndGravityTerm(
Eigen::Map<NodalMatrixType>& local_K,
Eigen::Map<NodalVectorType>& local_b,
- IntegrationPointData<NodalRowVectorType,
- GlobalDimNodalMatrixType> const& ip_data,
+ IntegrationPointData<GlobalDimNodalMatrixType> const& ip_data,
GlobalDimMatrixType const& permeability, double const mu,
double const rho_L, GlobalDimVectorType const& specific_body_force,
bool const has_gravity);
static Eigen::Matrix<double, GlobalDim, 1> calculateVelocity(
Eigen::Map<const NodalVectorType> const& local_p,
- IntegrationPointData<NodalRowVectorType,
- GlobalDimNodalMatrixType> const& ip_data,
+ IntegrationPointData<GlobalDimNodalMatrixType> const& ip_data,
GlobalDimMatrixType const& permeability, double const mu,
double const rho_L, GlobalDimVectorType const& specific_body_force,
bool const has_gravity);
@@ -194,16 +186,14 @@ private:
static void calculateLaplacianAndGravityTerm(
Eigen::Map<NodalMatrixType>& local_K,
Eigen::Map<NodalVectorType>& local_b,
- IntegrationPointData<NodalRowVectorType,
- GlobalDimNodalMatrixType> const& ip_data,
+ IntegrationPointData<GlobalDimNodalMatrixType> const& ip_data,
GlobalDimMatrixType const& permeability, double const mu,
double const rho_L, GlobalDimVectorType const& specific_body_force,
bool const has_gravity);
static Eigen::Matrix<double, GlobalDim, 1> calculateVelocity(
Eigen::Map<const NodalVectorType> const& local_p,
- IntegrationPointData<NodalRowVectorType,
- GlobalDimNodalMatrixType> const& ip_data,
+ IntegrationPointData<GlobalDimNodalMatrixType> const& ip_data,
GlobalDimMatrixType const& permeability, double const mu,
double const rho_L, GlobalDimVectorType const& specific_body_force,
bool const has_gravity);
@@ -231,7 +221,6 @@ private:
VelocityCacheType& darcy_velocity_at_ips) const;
const LiquidFlowData& _process_data;
- NumLib::ShapeMatrixCache const& _shape_matrix_cache;
};
} // namespace LiquidFlow
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 3a5d2af3e8b1759314a97a41bc41e49386ac20cd..23f7891cc4dd511ce86612613489196445d0661a 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -41,7 +41,6 @@ LiquidFlowProcess::LiquidFlowProcess(
integration_order, std::move(process_variables),
std::move(secondary_variables)),
_process_data(std::move(process_data)),
- _shape_matrix_cache{integration_order},
_surfaceflux(std::move(surfaceflux)),
_is_linear(is_linear)
{
@@ -60,7 +59,7 @@ void LiquidFlowProcess::initializeConcreteProcess(
ProcessLib::createLocalAssemblers<LiquidFlowLocalAssembler>(
mesh_space_dimension, mesh.getElements(), dof_table, _local_assemblers,
NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
- _process_data, _shape_matrix_cache);
+ _process_data);
_secondary_variables.addSecondaryVariable(
"darcy_velocity",
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index aac257a1e6f3284647df77de2d253459d097875e..8ac5910005eeb4ec6a5c8ea37724280b5ff5cbfa 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -17,7 +17,6 @@
#include "LiquidFlowData.h"
#include "LiquidFlowLocalAssembler.h"
#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "NumLib/Fem/ShapeMatrixCache.h"
#include "ProcessLib/Process.h"
#include "ProcessLib/SurfaceFlux/SurfaceFluxData.h"
@@ -109,8 +108,6 @@ private:
std::vector<std::unique_ptr<LiquidFlowLocalAssemblerInterface>>
_local_assemblers;
- NumLib::ShapeMatrixCache _shape_matrix_cache;
-
std::unique_ptr<ProcessLib::SurfaceFluxData> _surfaceflux;
MeshLib::PropertyVector<double>* _hydraulic_flow = nullptr;
bool _is_linear = false;