diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 3c7429b4ac12e69d40bc1e894b03bfdf63050ae7..7eae80a42c6e499c29ccd0c04678ebbf24444873 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -250,8 +250,7 @@ public:
: _element(element),
_process_data(process_data),
_integration_method(integration_method),
- _transport_process_variables(transport_process_variables),
- _b(getProjectedBodyForce())
+ _transport_process_variables(transport_process_variables)
{
(void)local_matrix_size;
@@ -468,6 +467,10 @@ public:
auto local_p = Eigen::Map<const NodalVectorType>(
&local_x[pressure_index], pressure_size);
+ auto const& b =
+ _process_data
+ .projected_specific_body_force_vectors[_element.getID()];
+
auto const number_of_components = num_nodal_dof - 1;
for (int component_id = 0; component_id < number_of_components;
++component_id)
@@ -500,8 +503,8 @@ public:
auto local_C = Eigen::Map<const NodalVectorType>(
&local_x[concentration_index], concentration_size);
- assembleBlockMatrices(_b, component_id, t, dt, local_C, local_p,
- KCC, MCC, MCp, MpC, Kpp, Mpp, Bp);
+ assembleBlockMatrices(b, component_id, t, dt, local_C, local_p, KCC,
+ MCC, MCp, MpC, Kpp, Mpp, Bp);
if (_process_data.chemical_solver_interface)
{
@@ -825,6 +828,10 @@ public:
ParameterLib::SpatialPosition pos;
pos.setElementID(_element.getID());
+ auto const& b =
+ _process_data
+ .projected_specific_body_force_vectors[_element.getID()];
+
auto const& medium =
*_process_data.media_map->getMedium(_element.getID());
auto const& phase = medium.phase("AqueousLiquid");
@@ -902,7 +909,7 @@ public:
if (_process_data.has_gravity)
{
local_b.noalias() +=
- w * density * density * dNdx.transpose() * K_over_mu * _b;
+ w * density * density * dNdx.transpose() * K_over_mu * b;
}
// coupling term
@@ -948,6 +955,10 @@ public:
ParameterLib::SpatialPosition pos;
pos.setElementID(_element.getID());
+ auto const& b =
+ _process_data
+ .projected_specific_body_force_vectors[_element.getID()];
+
MaterialPropertyLib::VariableArray vars;
MaterialPropertyLib::VariableArray vars_prev;
@@ -1038,7 +1049,7 @@ public:
GlobalDimVectorType const velocity =
_process_data.has_gravity
? GlobalDimVectorType(-K_over_mu *
- (dNdx * local_p - density * _b))
+ (dNdx * local_p - density * b))
: GlobalDimVectorType(-K_over_mu * dNdx * local_p);
GlobalDimMatrixType const hydrodynamic_dispersion =
@@ -1136,6 +1147,9 @@ public:
ParameterLib::SpatialPosition pos;
pos.setElementID(_element.getID());
+ auto const& b =
+ _process_data
+ .projected_specific_body_force_vectors[_element.getID()];
auto const& medium =
*_process_data.media_map->getMedium(_element.getID());
@@ -1209,7 +1223,7 @@ public:
if (_process_data.has_gravity)
{
local_rhs.noalias() +=
- w * rho * dNdx.transpose() * k / mu * rho * _b;
+ w * rho * dNdx.transpose() * k / mu * rho * b;
}
}
}
@@ -1245,6 +1259,10 @@ public:
ParameterLib::SpatialPosition pos;
pos.setElementID(_element.getID());
+ auto const& b =
+ _process_data
+ .projected_specific_body_force_vectors[_element.getID()];
+
MaterialPropertyLib::VariableArray vars;
MaterialPropertyLib::VariableArray vars_prev;
@@ -1320,7 +1338,7 @@ public:
// Darcy flux
GlobalDimVectorType const q =
_process_data.has_gravity
- ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * _b))
+ ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
: GlobalDimVectorType(-k / mu * dNdx * p);
GlobalDimMatrixType const D =
@@ -1482,6 +1500,10 @@ public:
ParameterLib::SpatialPosition pos;
pos.setElementID(_element.getID());
+ auto const& b =
+ _process_data
+ .projected_specific_body_force_vectors[_element.getID()];
+
MaterialPropertyLib::VariableArray vars;
auto const& medium =
@@ -1524,7 +1546,7 @@ public:
phase[MaterialPropertyLib::PropertyType::density]
.template value<double>(vars, pos, t, dt);
// here it is assumed that the vector b is directed 'downwards'
- cache_mat.col(ip).noalias() += K_over_mu * rho_w * _b;
+ cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
}
}
@@ -1560,6 +1582,9 @@ public:
ParameterLib::SpatialPosition pos;
pos.setElementID(_element.getID());
+ auto const& b =
+ _process_data
+ .projected_specific_body_force_vectors[_element.getID()];
MaterialPropertyLib::VariableArray vars;
@@ -1592,7 +1617,7 @@ public:
.template value<double>(vars, pos, t, dt);
if (_process_data.has_gravity)
{
- q += K_over_mu * rho_w * _b;
+ q += K_over_mu * rho_w * b;
}
Eigen::Vector3d flux(0.0, 0.0, 0.0);
flux.head<GlobalDim>() = rho_w * q;
@@ -1698,6 +1723,10 @@ public:
ParameterLib::SpatialPosition pos;
pos.setElementID(_element.getID());
+ auto const& b =
+ _process_data
+ .projected_specific_body_force_vectors[_element.getID()];
+
MaterialPropertyLib::VariableArray vars;
GlobalDimMatrixType const& I(
@@ -1740,7 +1769,7 @@ public:
// Darcy flux
GlobalDimVectorType const q =
_process_data.has_gravity
- ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * _b))
+ ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
: GlobalDimVectorType(-k / mu * dNdx * p);
auto const alpha_T = medium.template value<double>(
@@ -1791,14 +1820,6 @@ private:
Eigen::aligned_allocator<
IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>>>
_ip_data;
- GlobalDimVectorType const _b;
- GlobalDimVectorType getProjectedBodyForce()
- {
- return _process_data.element_rotation_matrices[_element.getID()] *
- _process_data.element_rotation_matrices[_element.getID()]
- .transpose() *
- _process_data.specific_body_force;
- }
};
} // namespace ComponentTransport
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcessData.h b/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
index a117aa417a1703c36beedf058f4f88f32e7d6ad5..7757358121573ace18dc937ff1ef347bb7dbc732 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
@@ -39,7 +39,6 @@ struct ComponentTransportProcessData
{
std::unique_ptr<MaterialPropertyLib::MaterialSpatialDistributionMap>
media_map;
- Eigen::VectorXd const specific_body_force;
bool const has_gravity;
bool const non_advective_form;
/// This optional tag provides a simple means of considering the temperature
@@ -73,8 +72,8 @@ struct ComponentTransportProcessData
std::unique_ptr<NumLib::NumericalStabilization> stabilizer;
- /// A vector of the rotation matrices for all elements.
- std::vector<Eigen::MatrixXd> const element_rotation_matrices;
+ /// Projected specific body force vector: R * R^T * b.
+ std::vector<Eigen::VectorXd> const projected_specific_body_force_vectors;
int const mesh_space_dimension;
diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
index 4a3a8887c82f43fdcea94a496e292d8ed698c3a5..012f046cd2c9a53f131dc5c2b890e60baecd0409 100644
--- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
@@ -211,7 +211,6 @@ std::unique_ptr<Process> createComponentTransportProcess(
bool const has_gravity = MathLib::toVector(b).norm() > 0;
if (has_gravity)
{
- specific_body_force.resize(b.size());
std::copy_n(b.data(), b.size(), specific_body_force.data());
}
@@ -254,9 +253,18 @@ std::unique_ptr<Process> createComponentTransportProcess(
*aperture_config, "parameter", parameters, 1);
}
+ auto const rotation_matrices = MeshLib::getElementRotationMatrices(
+ mesh_space_dimension, mesh.getDimension(), mesh.getElements());
+ std::vector<Eigen::VectorXd> projected_specific_body_force_vectors;
+ projected_specific_body_force_vectors.reserve(rotation_matrices.size());
+
+ std::transform(rotation_matrices.begin(), rotation_matrices.end(),
+ std::back_inserter(projected_specific_body_force_vectors),
+ [&specific_body_force](const auto& R)
+ { return R * R.transpose() * specific_body_force; });
+
ComponentTransportProcessData process_data{
std::move(media_map),
- specific_body_force,
has_gravity,
non_advective_form,
temperature,
@@ -264,8 +272,7 @@ std::unique_ptr<Process> createComponentTransportProcess(
chemical_solver_interface.get(),
std::move(lookup_table),
std::move(stabilizer),
- MeshLib::getElementRotationMatrices(
- mesh_space_dimension, mesh.getDimension(), mesh.getElements()),
+ projected_specific_body_force_vectors,
mesh_space_dimension,
*aperture_size_parameter};