diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index e8fa33311e569c6fb1bc155984153c7fad4ce883..b2e870bab5e63aa533d29c4f0c7fe4b49f5b3b48 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -16,6 +16,7 @@
#include "ComponentTransportProcessData.h"
#include "MaterialLib/Fluid/FluidProperties/FluidProperties.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
@@ -51,23 +52,11 @@ class ComponentTransportLocalAssemblerInterface
public NumLib::ExtrapolatableElement
{
public:
- virtual std::vector<double> const& getIntPtDarcyVelocityX(
- const double /*t*/,
- GlobalVector const& /*current_solution*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
- std::vector<double>& /*cache*/) const = 0;
-
- virtual std::vector<double> const& getIntPtDarcyVelocityY(
- const double /*t*/,
- GlobalVector const& /*current_solution*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
- std::vector<double>& /*cache*/) const = 0;
-
- virtual std::vector<double> const& getIntPtDarcyVelocityZ(
- const double /*t*/,
- GlobalVector const& /*current_solution*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
- std::vector<double>& /*cache*/) const = 0;
+ virtual std::vector<double> const& getIntPtDarcyVelocity(
+ const double t,
+ GlobalVector const& current_solution,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<double>& cache) const = 0;
};
template <typename ShapeFunction, typename IntegrationMethod,
@@ -100,10 +89,7 @@ public:
ComponentTransportProcessData const& process_data)
: _element(element),
_process_data(process_data),
- _integration_method(integration_order),
- _darcy_velocities(
- GlobalDim,
- std::vector<double>(_integration_method.getNumberOfPoints()))
+ _integration_method(integration_order)
{
// This assertion is valid only if all nodal d.o.f. use the same shape
// matrices.
@@ -267,24 +253,29 @@ public:
}
}
- void computeSecondaryVariableConcrete(
- double const t, std::vector<double> const& local_x) override
+ std::vector<double> const& getIntPtDarcyVelocity(
+ const double t,
+ GlobalVector const& current_solution,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<double>& cache) const override
{
+ auto const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ auto const indices = NumLib::getIndices(_element.getID(), dof_table);
+ assert(!indices.empty());
+ auto const local_x = current_solution.get(indices);
+
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<
+ Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, GlobalDim, n_integration_points);
+
SpatialPosition pos;
pos.setElementID(_element.getID());
- auto const& K =
- _process_data.porous_media_properties.getIntrinsicPermeability(t,
- pos);
MaterialLib::Fluid::FluidProperty::ArrayType vars;
- auto const mu = _process_data.fluid_properties->getValue(
- MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
- GlobalDimMatrixType const K_over_mu = K / mu;
-
- unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
-
auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
&local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
@@ -294,7 +285,14 @@ public:
auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
- GlobalDimVectorType velocity = -K_over_mu * dNdx * p_nodal_values;
+ auto const& K =
+ _process_data.porous_media_properties.getIntrinsicPermeability(
+ t, pos);
+ auto const mu = _process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+ GlobalDimMatrixType const K_over_mu = K / mu;
+
+ cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
if (_process_data.has_gravity)
{
double C_int_pt = 0.0;
@@ -310,14 +308,11 @@ public:
MaterialLib::Fluid::FluidPropertyType::Density, vars);
auto const b = _process_data.specific_body_force;
// here it is assumed that the vector b is directed 'downwards'
- velocity += K_over_mu * rho_w * b;
- }
-
- for (unsigned d = 0; d < GlobalDim; ++d)
- {
- _darcy_velocities[d][ip] = velocity[d];
+ cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
}
}
+
+ return cache;
}
@@ -330,36 +325,6 @@ public:
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
}
- std::vector<double> const& getIntPtDarcyVelocityX(
- const double /*t*/,
- GlobalVector const& /*current_solution*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
- std::vector<double>& /*cache*/) const override
- {
- assert(_darcy_velocities.size() > 0);
- return _darcy_velocities[0];
- }
-
- std::vector<double> const& getIntPtDarcyVelocityY(
- const double /*t*/,
- GlobalVector const& /*current_solution*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
- std::vector<double>& /*cache*/) const override
- {
- assert(_darcy_velocities.size() > 1);
- return _darcy_velocities[1];
- }
-
- std::vector<double> const& getIntPtDarcyVelocityZ(
- const double /*t*/,
- GlobalVector const& /*current_solution*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
- std::vector<double>& /*cache*/) const override
- {
- assert(_darcy_velocities.size() > 2);
- return _darcy_velocities[2];
- }
-
private:
MeshLib::Element const& _element;
ComponentTransportProcessData const& _process_data;
@@ -370,7 +335,6 @@ private:
Eigen::aligned_allocator<
IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>>>
_ip_data;
- std::vector<std::vector<double>> _darcy_velocities;
};
} // namespace ComponentTransport
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index 2ebabe39ebe8fa8b95d86bf10dd2f08864ea0a2d..40826eaa2903d55565ed7a5e0da52b1c822379dd 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -45,27 +45,10 @@ void ComponentTransportProcess::initializeConcreteProcess(
mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
- "darcy_velocity_x",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &ComponentTransportLocalAssemblerInterface::
- getIntPtDarcyVelocityX));
-
- if (mesh.getDimension() > 1)
- {
- _secondary_variables.addSecondaryVariable(
- "darcy_velocity_y",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &ComponentTransportLocalAssemblerInterface::
- getIntPtDarcyVelocityY));
- }
- if (mesh.getDimension() > 2)
- {
- _secondary_variables.addSecondaryVariable(
- "darcy_velocity_z",
- makeExtrapolator(1, getExtrapolator(), _local_assemblers,
- &ComponentTransportLocalAssemblerInterface::
- getIntPtDarcyVelocityZ));
- }
+ "darcy_velocity",
+ makeExtrapolator(
+ mesh.getDimension(), getExtrapolator(), _local_assemblers,
+ &ComponentTransportLocalAssemblerInterface::getIntPtDarcyVelocity));
}
void ComponentTransportProcess::assembleConcreteProcess(
@@ -99,16 +82,6 @@ void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
dx_dx, M, K, b, Jac, coupling_term);
}
-void ComponentTransportProcess::computeSecondaryVariableConcrete(
- double const t, GlobalVector const& x,
- StaggeredCouplingTerm const& coupling_term)
-{
- DBUG("Compute the Darcy velocity for ComponentTransportProcess.");
- GlobalExecutor::executeMemberOnDereferenced(
- &ComponentTransportLocalAssemblerInterface::computeSecondaryVariable,
- _local_assemblers, *_local_to_global_index_map, t, x, coupling_term);
-}
-
} // namespace ComponentTransport
} // namespace ProcessLib
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.h b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
index 71b5dd902a41c0866754ea852cbe3d82bef6ecc3..84437c0f4851d9a6c3289c2c4fb9c663ff66f046 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
@@ -95,10 +95,6 @@ public:
bool isLinear() const override { return false; }
//! @}
- void computeSecondaryVariableConcrete(
- double const t, GlobalVector const& x,
- StaggeredCouplingTerm const& coupling_term) override;
-
private:
void initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table,