diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index a7fb7140b891e98f891778bec8cf0eaa66604b6f..56d38818312ab628741143a10597c99d01ccdd05 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -16,6 +16,8 @@
#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/StaggeredCouplingTerm.h"
+
#include "ProcessLib/HeatConduction/HeatConductionProcess.h"
namespace ProcessLib
@@ -241,37 +243,73 @@ 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)
+std::vector<double> const&
+LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ getIntPtDarcyVelocity(const double t,
+ GlobalVector const& current_solution,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<double>& veloctiy_cache) const
{
+ // auto const num_nodes = ShapeFunction_::NPOINTS;
+ auto const num_intpts = _shape_matrices.size();
+ auto const indices = NumLib::getIndices(_element.getID(), dof_table);
+ auto const local_x = current_solution.get(indices);
+ veloctiy_cache.clear();
+ auto veloctiy_cache_vectors = MathLib::createZeroedMatrix<
+ Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+ veloctiy_cache, GlobalDim, num_intpts);
+
SpatialPosition pos;
pos.setElementID(_element.getID());
const int material_id = _material_properties.getMaterialID(pos);
- const Eigen::MatrixXd& permeability = _material_properties.getPermeability(
+ const Eigen::MatrixXd& perm = _material_properties.getPermeability(
material_id, 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 permeability.rows() ==
- // _element->getDimension()
- assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
+ // the assert must be changed to perm.rows() == _element->getDimension()
+ assert(perm.rows() == GlobalDim || perm.rows() == 1);
+ if (!_coupling_term)
+ {
+ computeDarcyVelocity(perm, local_x, veloctiy_cache_vectors);
+ }
+ else
+ {
+ auto const local_coupled_xs =
+ getCurrentLocalSolutionsOfCoupledProcesses(
+ _coupling_term->coupled_xs, indices);
+ if (local_coupled_xs.empty())
+ computeDarcyVelocity(perm, local_x, veloctiy_cache_vectors);
+ else
+ computeDarcyVelocityWithCoupling(perm, local_x, local_coupled_xs,
+ veloctiy_cache_vectors);
+ }
+
+ return veloctiy_cache;
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ computeDarcyVelocity(Eigen::MatrixXd const& permeability,
+ std::vector<double> const& local_x,
+ LocalMatrixType& darcy_velocity_at_ips) const
+{
if (permeability.size() == 1) // isotropic or 1D problem.
- computeSecondaryVariableLocal<IsotropicCalculator>(t, local_x, pos,
- permeability);
+ computeDarcyVelocityLocal<IsotropicCalculator>(local_x, permeability,
+ darcy_velocity_at_ips);
else
- computeSecondaryVariableLocal<AnisotropicCalculator>(t, local_x, pos,
- permeability);
+ computeDarcyVelocityLocal<AnisotropicCalculator>(local_x, permeability,
+ darcy_velocity_at_ips);
}
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
template <typename LaplacianGravityVelocityCalculator>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
- computeSecondaryVariableLocal(double const /*t*/,
- std::vector<double> const& local_x,
- SpatialPosition const& /*pos*/,
- Eigen::MatrixXd const& permeability)
+ computeDarcyVelocityLocal(std::vector<double> const& local_x,
+ Eigen::MatrixXd const& permeability,
+ LocalMatrixType& darcy_velocity_at_ips) const
{
auto const local_matrix_size = local_x.size();
assert(local_matrix_size == ShapeFunction::NPOINTS);
@@ -296,7 +334,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
_material_properties.getViscosity(p, _reference_temperature);
LaplacianGravityVelocityCalculator::calculateVelocity(
- _darcy_velocities, local_p_vec, sm, permeability, ip, mu, rho_g,
+ ip, darcy_velocity_at_ips, local_p_vec, sm, permeability, mu, rho_g,
_gravitational_axis_id);
}
}
@@ -304,42 +342,33 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
- computeSecondaryVariableWithCoupledProcessConcrete(
- double const t, std::vector<double> const& local_x,
+ computeDarcyVelocityWithCoupling(
+ Eigen::MatrixXd const& permeability, std::vector<double> const& local_x,
std::unordered_map<std::type_index, const std::vector<double>> const&
- coupled_local_solutions)
+ coupled_local_solutions,
+ LocalMatrixType& darcy_velocity_at_ips) const
{
- SpatialPosition pos;
- pos.setElementID(_element.getID());
- const int material_id = _material_properties.getMaterialID(pos);
- const Eigen::MatrixXd& perm = _material_properties.getPermeability(
- material_id, t, pos, _element.getDimension());
-
const auto local_T = coupled_local_solutions.at(std::type_index(
typeid(ProcessLib::HeatConduction::HeatConductionProcess)));
- // 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.
- computeSecondaryVariableCoupledWithHeatTransportLocal<
- IsotropicCalculator>(t, local_x, local_T, pos, perm);
+ if (permeability.size() == 1) // isotropic or 1D problem.
+ computeDarcyVelocityCoupledWithHeatTransportLocal<IsotropicCalculator>(
+ local_x, local_T, permeability, darcy_velocity_at_ips);
else
- computeSecondaryVariableCoupledWithHeatTransportLocal<
- AnisotropicCalculator>(t, local_x, local_T, pos, perm);
+ computeDarcyVelocityCoupledWithHeatTransportLocal<
+ AnisotropicCalculator>(local_x, local_T, permeability,
+ darcy_velocity_at_ips);
}
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
template <typename LaplacianGravityVelocityCalculator>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
- computeSecondaryVariableCoupledWithHeatTransportLocal(
- double const /*t*/,
+ computeDarcyVelocityCoupledWithHeatTransportLocal(
std::vector<double> const& local_x,
std::vector<double> const& local_T,
- SpatialPosition const& /*pos*/,
- Eigen::MatrixXd const& permeability)
+ Eigen::MatrixXd const& permeability,
+ LocalMatrixType& darcy_velocity_at_ips) const
{
auto const local_matrix_size = local_x.size();
assert(local_matrix_size == ShapeFunction::NPOINTS);
@@ -364,7 +393,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
const double mu = _material_properties.getViscosity(p, T);
LaplacianGravityVelocityCalculator::calculateVelocity(
- _darcy_velocities, local_p_vec, sm, permeability, ip, mu, rho_g,
+ ip, darcy_velocity_at_ips, local_p_vec, sm, permeability, mu, rho_g,
_gravitational_axis_id);
}
}
@@ -397,22 +426,17 @@ template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
IsotropicCalculator::calculateVelocity(
- std::vector<std::vector<double>>& darcy_velocities,
+ unsigned ip, LocalMatrixType& darcy_velocity_at_ips,
Eigen::Map<const NodalVectorType> const& local_p,
ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
- unsigned const ip, double const mu, double const rho_g,
- int const gravitational_axis_id)
+ double const mu, double const rho_g, int const gravitational_axis_id)
{
const double K = permeability(0, 0) / mu;
// Compute the velocity
- GlobalDimVectorType darcy_velocity = -K * sm.dNdx * local_p;
+ darcy_velocity_at_ips.col(ip).noalias() = -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];
- }
+ darcy_velocity_at_ips.col(ip)[gravitational_axis_id] -= K * rho_g;
}
template <typename ShapeFunction, typename IntegrationMethod,
@@ -441,23 +465,19 @@ template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
AnisotropicCalculator::calculateVelocity(
- std::vector<std::vector<double>>& darcy_velocities,
+ unsigned ip, LocalMatrixType& darcy_velocity_at_ips,
Eigen::Map<const NodalVectorType> const& local_p,
ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
- unsigned const ip, double const mu, double const rho_g,
- int const gravitational_axis_id)
+ double const mu, double const rho_g, int const gravitational_axis_id)
{
// Compute the velocity
- GlobalDimVectorType darcy_velocity = -permeability * sm.dNdx * local_p / mu;
+ darcy_velocity_at_ips.col(ip).noalias() =
+ -permeability * sm.dNdx * local_p / mu;
if (gravitational_axis_id >= 0)
{
- darcy_velocity.noalias() -=
+ darcy_velocity_at_ips.col(ip).noalias() -=
rho_g * permeability.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 d726b76bd1999936036e40aa674d6f774fdfe3a1..c07fbfea53b4dd9637fdf1cb0191ae39bd5c303f 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -31,6 +31,8 @@
namespace ProcessLib
{
+struct StaggeredCouplingTerm;
+
namespace LiquidFlow
{
const unsigned NUM_NODAL_DOF = 1;
@@ -40,11 +42,27 @@ class LiquidFlowLocalAssemblerInterface
public NumLib::ExtrapolatableElement
{
public:
+ LiquidFlowLocalAssemblerInterface(
+ StaggeredCouplingTerm* const coupling_term)
+ : _coupling_term(coupling_term)
+ {
+ }
+
+ void setStaggeredCouplingTerm(std::size_t const /*mesh_item_id*/,
+ StaggeredCouplingTerm* const coupling_term)
+ {
+ _coupling_term = coupling_term;
+ }
+
virtual std::vector<double> const& getIntPtDarcyVelocity(
const double /*t*/,
GlobalVector const& /*current_solution*/,
NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
+
+protected:
+ /// Pointer that gets from Process class
+ StaggeredCouplingTerm* _coupling_term;
};
template <typename ShapeFunction, typename IntegrationMethod,
@@ -61,6 +79,9 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface
using NodalVectorType = typename LocalAssemblerTraits::LocalVector;
using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+ using LocalMatrixType = Eigen::Map<
+ Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>;
+
public:
LiquidFlowLocalAssembler(
MeshLib::Element const& element,
@@ -70,8 +91,10 @@ public:
int const gravitational_axis_id,
double const gravitational_acceleration,
double const reference_temperature,
- LiquidFlowMaterialProperties const& material_propertries)
- : _element(element),
+ LiquidFlowMaterialProperties const& material_propertries,
+ StaggeredCouplingTerm* coupling_term)
+ : LiquidFlowLocalAssemblerInterface(coupling_term),
+ _element(element),
_integration_method(integration_order),
_shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod, GlobalDim>(
@@ -94,14 +117,6 @@ public:
std::vector<double>& local_b_data,
LocalCouplingTerm const& coupled_term) override;
- void computeSecondaryVariableConcrete(
- double const /*t*/, std::vector<double> const& local_x) override;
-
- void computeSecondaryVariableWithCoupledProcessConcrete(
- double const t, std::vector<double> const& local_x,
- std::unordered_map<std::type_index, const std::vector<double>> const&
- coupled_local_solutions) override;
-
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
@@ -115,40 +130,7 @@ public:
const double t,
GlobalVector const& current_solution,
NumLib::LocalToGlobalIndexMap const& dof_table,
- std::vector<double>& cache) const override
- {
- // auto const num_nodes = ShapeFunction_::NPOINTS;
- auto const num_intpts = _shape_matrices.size();
-
- auto const indices = NumLib::getIndices(
- _element.getID(), dof_table);
- assert(!indices.empty());
- auto const local_x = current_solution.get(indices);
- auto const local_x_vec =
- MathLib::toVector<Eigen::Matrix<double, ShapeFunction::NPOINTS, 1>>(
- local_x, ShapeFunction::NPOINTS);
-
- cache.clear();
- auto cache_vec = MathLib::createZeroedMatrix<
- Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, GlobalDim, num_intpts);
-
- SpatialPosition pos;
- pos.setElementID(_element.getID());
-
- // TODO fix
-#if 0
- for (unsigned i = 0; i < num_intpts; ++i) {
- pos.setIntegrationPoint(i);
- auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
- // dimensions: (d x 1) = (d x n) * (n x 1)
- cache_vec.col(i).noalias() =
- -k * _shape_matrices[i].dNdx * local_x_vec;
- }
-#endif
-
- return cache;
- }
+ std::vector<double>& veloctiy_cache) const override;
private:
MeshLib::Element const& _element;
@@ -157,11 +139,6 @@ private:
std::vector<ShapeMatrices, Eigen::aligned_allocator<ShapeMatrices>>
_shape_matrices;
- std::vector<std::vector<double>> _darcy_velocities =
- std::vector<std::vector<double>>(
- GlobalDim,
- std::vector<double>(_integration_method.getNumberOfPoints()));
-
/**
* Calculator of the Laplacian and the gravity term for anisotropic
* permeability tensor
@@ -176,10 +153,10 @@ private:
double const rho_g, int const gravitational_axis_id);
static void calculateVelocity(
- std::vector<std::vector<double>>& darcy_velocities,
+ unsigned ip, LocalMatrixType& darcy_velocity_at_ips,
Eigen::Map<const NodalVectorType> const& local_p,
ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
- unsigned const ip, double const mu, double const rho_g,
+ double const mu, double const rho_g,
int const gravitational_axis_id);
};
@@ -197,10 +174,10 @@ private:
double const rho_g, int const gravitational_axis_id);
static void calculateVelocity(
- std::vector<std::vector<double>>& darcy_velocities,
+ unsigned ip, LocalMatrixType& darcy_velocity_at_ips,
Eigen::Map<const NodalVectorType> const& local_p,
ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
- unsigned const ip, double const mu, double const rho_g,
+ double const mu, double const rho_g,
int const gravitational_axis_id);
};
@@ -221,19 +198,28 @@ private:
std::vector<double>& local_K_data, std::vector<double>& local_b_data,
SpatialPosition const& pos, Eigen::MatrixXd const& permeability);
+ void computeDarcyVelocity(Eigen::MatrixXd const& permeability,
+ std::vector<double> const& local_x,
+ LocalMatrixType& darcy_velocity_at_ips) const;
+
+ void computeDarcyVelocityWithCoupling(
+ Eigen::MatrixXd const& permeability, std::vector<double> const& local_x,
+ std::unordered_map<std::type_index, const std::vector<double>> const&
+ coupled_local_solutions,
+ LocalMatrixType& darcy_velocity_at_ips) const;
+
template <typename LaplacianGravityVelocityCalculator>
- void computeSecondaryVariableLocal(double const /*t*/,
- std::vector<double> const& local_x,
- SpatialPosition const& pos,
- Eigen::MatrixXd const& permeability);
+ void computeDarcyVelocityLocal(
+ std::vector<double> const& local_x,
+ Eigen::MatrixXd const& permeability,
+ LocalMatrixType& darcy_velocity_at_ips) const;
template <typename LaplacianGravityVelocityCalculator>
- void computeSecondaryVariableCoupledWithHeatTransportLocal(
- double const /*t*/,
+ void computeDarcyVelocityCoupledWithHeatTransportLocal(
std::vector<double> const& local_x,
std::vector<double> const& local_T,
- SpatialPosition const& pos,
- Eigen::MatrixXd const& permeability);
+ Eigen::MatrixXd const& permeability,
+ LocalMatrixType& darcy_velocity_at_ips) const;
const int _gravitational_axis_id;
const double _gravitational_acceleration;