diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index de4116f0bab853e0bec6c81347f2e6d1ba9f5985..e10b6b401d4c17ad006196e5379462a574fbf7d1 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -40,16 +40,16 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
assert(perm.rows() == GlobalDim || perm.rows() == 1);
if (perm.size() == 1) // isotropic or 1D problem.
- local_assemble<IsotropicLaplacianAndGravityTermCalculator>(
+ local_assemble<IsotropicCalculator>(
t, local_x, local_M_data, local_K_data, local_b_data, pos, perm);
else
- local_assemble<AnisotropicLaplacianAndGravityTermCalculator>(
+ local_assemble<AnisotropicCalculator>(
t, local_x, local_M_data, local_K_data, local_b_data, pos, perm);
}
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
-template <typename LaplacianAndGravityTermCalculator>
+template <typename Calculator>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
local_assemble(double const t, std::vector<double> const& local_x,
std::vector<double>& local_M_data,
@@ -100,7 +100,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
const double mu = _material_properties.getViscosity(p, _temperature);
// Assemble Laplacian, K, and RHS by the gravitational term
- LaplacianAndGravityTermCalculator::calculate(
+ Calculator::calculateLaplacianAndGravityTerm(
local_K, local_b, sm, perm, integration_factor, mu, rho_g,
_gravitational_axis_id);
}
@@ -109,7 +109,67 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
- IsotropicLaplacianAndGravityTermCalculator::calculate(
+ computeSecondaryVariableConcrete(double const t,
+ std::vector<double> const& local_x)
+{
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+ _material_properties.setMaterialID(pos);
+ const Eigen::MatrixXd& perm =
+ _material_properties.getPermeability(t, pos, _element.getDimension());
+
+ // Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
+ // the assert must be changed to perm.rows() == _element->getDimension()
+ assert(perm.rows() == GlobalDim || perm.rows() == 1);
+
+ if (perm.size() == 1) // isotropic or 1D problem.
+ computeSecondaryVariableLocal<IsotropicCalculator>(t, local_x, pos,
+ perm);
+ else
+ computeSecondaryVariableLocal<AnisotropicCalculator>(t, local_x, pos,
+ perm);
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+template <typename Calculator>
+void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ computeSecondaryVariableLocal(double const /*t*/,
+ std::vector<double> const& local_x,
+ SpatialPosition const& /*pos*/,
+ Eigen::MatrixXd const& perm)
+{
+ auto const local_matrix_size = local_x.size();
+ assert(local_matrix_size == ShapeFunction::NPOINTS);
+
+ const auto local_p_vec =
+ MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& sm = _shape_matrices[ip];
+ double p = 0.;
+ NumLib::shapeFunctionInterpolate(local_x, sm.N, p);
+ // TODO : compute _temperature from the heat transport pcs
+
+ const double rho_g =
+ _material_properties.getLiquidDensity(p, _temperature) *
+ _gravitational_acceleration;
+ // Compute viscosity:
+ const double mu = _material_properties.getViscosity(p, _temperature);
+
+ Calculator::calculateVelocity(_darcy_velocities, local_p_vec, sm, perm,
+ ip, mu, rho_g, _gravitational_axis_id);
+ }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ IsotropicCalculator::calculateLaplacianAndGravityTerm(
Eigen::Map<NodalMatrixType>& local_K,
Eigen::Map<NodalVectorType>& local_b, ShapeMatrices const& sm,
Eigen::MatrixXd const& perm, double const integration_factor,
@@ -133,7 +193,28 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
- AnisotropicLaplacianAndGravityTermCalculator::calculate(
+ IsotropicCalculator::calculateVelocity(
+ std::vector<std::vector<double>>& darcy_velocities,
+ Eigen::Map<const NodalVectorType> const& local_p,
+ ShapeMatrices const& sm, Eigen::MatrixXd const& perm, unsigned const ip,
+ double const mu, double const rho_g, int const gravitational_axis_id)
+{
+ const double K = perm(0, 0) / mu;
+ // Compute the velocity
+ GlobalDimVectorType darcy_velocity = -K * sm.dNdx * local_p;
+ // gravity term
+ if (gravitational_axis_id >= 0)
+ darcy_velocity[gravitational_axis_id] -= K * rho_g;
+ for (unsigned d = 0; d < GlobalDim; ++d)
+ {
+ darcy_velocities[d][ip] = darcy_velocity[d];
+ }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ AnisotropicCalculator::calculateLaplacianAndGravityTerm(
Eigen::Map<NodalMatrixType>& local_K,
Eigen::Map<NodalVectorType>& local_b, ShapeMatrices const& sm,
Eigen::MatrixXd const& perm, double const integration_factor,
@@ -154,70 +235,23 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
- computeSecondaryVariableConcrete(double const t,
- std::vector<double> const& local_x)
+ AnisotropicCalculator::calculateVelocity(
+ std::vector<std::vector<double>>& darcy_velocities,
+ Eigen::Map<const NodalVectorType> const& local_p,
+ ShapeMatrices const& sm, Eigen::MatrixXd const& perm, unsigned const ip,
+ double const mu, double const rho_g, int const gravitational_axis_id)
{
- auto const local_matrix_size = local_x.size();
- assert(local_matrix_size == ShapeFunction::NPOINTS);
-
- const auto local_p_vec =
- MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
-
- unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
-
- SpatialPosition pos;
- pos.setElementID(_element.getID());
- _material_properties.setMaterialID(pos);
- const Eigen::MatrixXd& perm =
- _material_properties.getPermeability(t, pos, _element.getDimension());
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
+ // Compute the velocity
+ GlobalDimVectorType darcy_velocity = -perm * sm.dNdx * local_p / mu;
+ if (gravitational_axis_id >= 0)
{
- auto const& sm = _shape_matrices[ip];
- double p = 0.;
- NumLib::shapeFunctionInterpolate(local_x, sm.N, p);
- // TODO : compute _temperature from the heat transport pcs
-
- const double rho_g =
- _material_properties.getLiquidDensity(p, _temperature) *
- _gravitational_acceleration;
- // Compute viscosity:
- const double mu = _material_properties.getViscosity(p, _temperature);
-
- // Assemble Laplacian, K, and RHS by the gravitational term
- if (perm.size() == 1) // Save time for isotropic permeability.
- {
- // Use scalar number for isotropic permeability
- // to save the computation time.
- const double K = perm(0, 0) / mu;
- // Compute the velocity
- GlobalDimVectorType darcy_velocity = -K * sm.dNdx * local_p_vec;
- // gravity term
- if (_gravitational_axis_id >= 0)
- darcy_velocity[_gravitational_axis_id] -= K * rho_g;
- for (unsigned d = 0; d < GlobalDim; ++d)
- {
- _darcy_velocities[d][ip] = darcy_velocity[d];
- }
- }
- else
- {
- // Compute the velocity
- GlobalDimVectorType darcy_velocity
- = -perm * sm.dNdx * local_p_vec / mu;
- if (_gravitational_axis_id >= 0)
- {
- darcy_velocity.noalias() -=
- rho_g * perm.col(_gravitational_axis_id) / mu;
- }
- for (unsigned d = 0; d < GlobalDim; ++d)
- {
- _darcy_velocities[d][ip] = darcy_velocity[d];
- }
- }
+ darcy_velocity.noalias() -=
+ rho_g * perm.col(gravitational_axis_id) / mu;
+ }
+ for (unsigned d = 0; d < GlobalDim; ++d)
+ {
+ darcy_velocities[d][ip] = darcy_velocity[d];
}
-
}
} // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index ac6264efc881414a8b17ff609a2bf32cdafd6421..d0bcec63de1dfec4dad70d47a4b33a099a757ed3 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -85,8 +85,8 @@ public:
std::vector<double>& local_K_data,
std::vector<double>& local_b_data) override;
- void computeSecondaryVariableConcrete(double const /*t*/,
- std::vector<double> const& local_x) override;
+ void computeSecondaryVariableConcrete(
+ double const /*t*/, std::vector<double> const& local_x) override;
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
@@ -133,33 +133,45 @@ private:
* Calculator of the Laplacian and the gravity term for anisotropic
* permeability tensor
*/
- struct AnisotropicLaplacianAndGravityTermCalculator
+ struct AnisotropicCalculator
{
- static void calculate(Eigen::Map<NodalMatrixType>& local_K,
- Eigen::Map<NodalVectorType>& local_b,
- ShapeMatrices const& sm,
- Eigen::MatrixXd const& perm,
- double const integration_factor, double const mu,
- double const rho_g,
- int const gravitational_axis_id);
+ static void calculateLaplacianAndGravityTerm(
+ Eigen::Map<NodalMatrixType>& local_K,
+ Eigen::Map<NodalVectorType>& local_b, ShapeMatrices const& sm,
+ Eigen::MatrixXd const& perm, double const integration_factor,
+ double const mu, double const rho_g,
+ int const gravitational_axis_id);
+
+ static void calculateVelocity(
+ std::vector<std::vector<double>>& darcy_velocities,
+ Eigen::Map<const NodalVectorType> const& local_p,
+ ShapeMatrices const& sm, Eigen::MatrixXd const& perm,
+ unsigned const ip, double const mu, double const rho_g,
+ int const gravitational_axis_id);
};
/**
* Calculator of the Laplacian and the gravity term for isotropic
* permeability tensor
*/
- struct IsotropicLaplacianAndGravityTermCalculator
+ struct IsotropicCalculator
{
- static void calculate(Eigen::Map<NodalMatrixType>& local_K,
- Eigen::Map<NodalVectorType>& local_b,
- ShapeMatrices const& sm,
- Eigen::MatrixXd const& perm,
- double const integration_factor, double const mu,
- double const rho_g,
- int const gravitational_axis_id);
+ static void calculateLaplacianAndGravityTerm(
+ Eigen::Map<NodalMatrixType>& local_K,
+ Eigen::Map<NodalVectorType>& local_b, ShapeMatrices const& sm,
+ Eigen::MatrixXd const& perm, double const integration_factor,
+ double const mu, double const rho_g,
+ int const gravitational_axis_id);
+
+ static void calculateVelocity(
+ std::vector<std::vector<double>>& darcy_velocities,
+ Eigen::Map<const NodalVectorType> const& local_p,
+ ShapeMatrices const& sm, Eigen::MatrixXd const& perm,
+ unsigned const ip, double const mu, double const rho_g,
+ int const gravitational_axis_id);
};
- template <typename LaplacianAndGravityTermCalculator>
+ template <typename Calculator>
void local_assemble(double const t, std::vector<double> const& local_x,
std::vector<double>& local_M_data,
std::vector<double>& local_K_data,
@@ -167,6 +179,12 @@ private:
SpatialPosition const& pos,
Eigen::MatrixXd const& perm);
+ template <typename Calculator>
+ void computeSecondaryVariableLocal(double const /*t*/,
+ std::vector<double> const& local_x,
+ SpatialPosition const& pos,
+ Eigen::MatrixXd const& perm);
+
const int _gravitational_axis_id;
const double _gravitational_acceleration;
LiquidFlowMaterialProperties& _material_properties;
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 07dbb24c6c80df4251461144e62bfe7e78b6e1cc..795fe907b11880d7c670cce54b7da6937cb31d7f 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -118,9 +118,9 @@ void LiquidFlowProcess::computeSecondaryVariableConcrete(const double t,
{
DBUG("Compute the velocity for LiquidFlowProcess.");
GlobalExecutor::executeMemberOnDereferenced(
- &LiquidFlowLocalAssemblerInterface::computeSecondaryVariable,
- _local_assemblers, *_local_to_global_index_map, t, x);
+ &LiquidFlowLocalAssemblerInterface::computeSecondaryVariable,
+ _local_assemblers, *_local_to_global_index_map, t, x);
}
-
+
} // end of namespace
} // end of namespace
diff --git a/Tests/Data b/Tests/Data
index 8806cf31aa300ee9c631c73674e218e9a60297cc..7a19025a1658a79ca444cee5d12c08a25c1ed6fc 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 8806cf31aa300ee9c631c73674e218e9a60297cc
+Subproject commit 7a19025a1658a79ca444cee5d12c08a25c1ed6fc