diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 56d38818312ab628741143a10597c99d01ccdd05..a683ad2528b6869ee4c8e94644eb07743df67291 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -250,10 +250,9 @@ LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
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);
+ auto const num_intpts = _integration_method.getNumberOfPoints();
veloctiy_cache.clear();
auto veloctiy_cache_vectors = MathLib::createZeroedMatrix<
Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
@@ -262,16 +261,16 @@ LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
SpatialPosition pos;
pos.setElementID(_element.getID());
const int material_id = _material_properties.getMaterialID(pos);
- const Eigen::MatrixXd& perm = _material_properties.getPermeability(
+ const Eigen::MatrixXd& permeability = _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 perm.rows() == _element->getDimension()
- assert(perm.rows() == GlobalDim || perm.rows() == 1);
+ assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
if (!_coupling_term)
{
- computeDarcyVelocity(perm, local_x, veloctiy_cache_vectors);
+ computeDarcyVelocity(permeability, local_x, veloctiy_cache_vectors);
}
else
{
@@ -279,9 +278,10 @@ LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
getCurrentLocalSolutionsOfCoupledProcesses(
_coupling_term->coupled_xs, indices);
if (local_coupled_xs.empty())
- computeDarcyVelocity(perm, local_x, veloctiy_cache_vectors);
+ computeDarcyVelocity(permeability, local_x, veloctiy_cache_vectors);
else
- computeDarcyVelocityWithCoupling(perm, local_x, local_coupled_xs,
+ computeDarcyVelocityWithCoupling(permeability, local_x,
+ local_coupled_xs,
veloctiy_cache_vectors);
}
@@ -291,9 +291,10 @@ LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
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
+ computeDarcyVelocity(
+ Eigen::MatrixXd const& permeability,
+ std::vector<double> const& local_x,
+ MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const
{
if (permeability.size() == 1) // isotropic or 1D problem.
computeDarcyVelocityLocal<IsotropicCalculator>(local_x, permeability,
@@ -307,9 +308,10 @@ template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
template <typename LaplacianGravityVelocityCalculator>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
- computeDarcyVelocityLocal(std::vector<double> const& local_x,
- Eigen::MatrixXd const& permeability,
- LocalMatrixType& darcy_velocity_at_ips) const
+ computeDarcyVelocityLocal(
+ std::vector<double> const& local_x,
+ Eigen::MatrixXd const& permeability,
+ MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const
{
auto const local_matrix_size = local_x.size();
assert(local_matrix_size == ShapeFunction::NPOINTS);
@@ -334,8 +336,8 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
_material_properties.getViscosity(p, _reference_temperature);
LaplacianGravityVelocityCalculator::calculateVelocity(
- ip, darcy_velocity_at_ips, local_p_vec, sm, permeability, mu, rho_g,
- _gravitational_axis_id);
+ ip, local_p_vec, sm, permeability, mu, rho_g,
+ _gravitational_axis_id, darcy_velocity_at_ips);
}
}
@@ -346,7 +348,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
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
+ MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const
{
const auto local_T = coupled_local_solutions.at(std::type_index(
typeid(ProcessLib::HeatConduction::HeatConductionProcess)));
@@ -368,7 +370,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
std::vector<double> const& local_x,
std::vector<double> const& local_T,
Eigen::MatrixXd const& permeability,
- LocalMatrixType& darcy_velocity_at_ips) const
+ MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const
{
auto const local_matrix_size = local_x.size();
assert(local_matrix_size == ShapeFunction::NPOINTS);
@@ -393,8 +395,8 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
const double mu = _material_properties.getViscosity(p, T);
LaplacianGravityVelocityCalculator::calculateVelocity(
- ip, darcy_velocity_at_ips, local_p_vec, sm, permeability, mu, rho_g,
- _gravitational_axis_id);
+ ip, local_p_vec, sm, permeability, mu, rho_g,
+ _gravitational_axis_id, darcy_velocity_at_ips);
}
}
@@ -426,10 +428,10 @@ template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
IsotropicCalculator::calculateVelocity(
- unsigned ip, LocalMatrixType& darcy_velocity_at_ips,
- Eigen::Map<const NodalVectorType> const& local_p,
+ unsigned const ip, Eigen::Map<const NodalVectorType> const& local_p,
ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
- double const mu, double const rho_g, int const gravitational_axis_id)
+ double const mu, double const rho_g, int const gravitational_axis_id,
+ MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips)
{
const double K = permeability(0, 0) / mu;
// Compute the velocity
@@ -465,10 +467,10 @@ template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
AnisotropicCalculator::calculateVelocity(
- unsigned ip, LocalMatrixType& darcy_velocity_at_ips,
- Eigen::Map<const NodalVectorType> const& local_p,
+ unsigned const ip, Eigen::Map<const NodalVectorType> const& local_p,
ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
- double const mu, double const rho_g, int const gravitational_axis_id)
+ double const mu, double const rho_g, int const gravitational_axis_id,
+ MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips)
{
// Compute the velocity
darcy_velocity_at_ips.col(ip).noalias() =
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index c07fbfea53b4dd9637fdf1cb0191ae39bd5c303f..d4cbb1ae839d5ff64c3eed858244f1c5cb599820 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -61,7 +61,10 @@ public:
std::vector<double>& /*cache*/) const = 0;
protected:
- /// Pointer that gets from Process class
+ // TODO: remove _coupling_term or move integration point data from local
+ // assembler class to a new class to make local assembler unique for each
+ //process.
+ /// Pointer that is set from a Process class.
StaggeredCouplingTerm* _coupling_term;
};
@@ -79,7 +82,7 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface
using NodalVectorType = typename LocalAssemblerTraits::LocalVector;
using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
- using LocalMatrixType = Eigen::Map<
+ using MatrixOfVelocityAtIntegrationPoints = Eigen::Map<
Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>;
public:
@@ -153,11 +156,11 @@ private:
double const rho_g, int const gravitational_axis_id);
static void calculateVelocity(
- unsigned ip, LocalMatrixType& darcy_velocity_at_ips,
- Eigen::Map<const NodalVectorType> const& local_p,
+ unsigned const ip, Eigen::Map<const NodalVectorType> const& local_p,
ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
double const mu, double const rho_g,
- int const gravitational_axis_id);
+ int const gravitational_axis_id,
+ MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips);
};
/**
@@ -174,11 +177,11 @@ private:
double const rho_g, int const gravitational_axis_id);
static void calculateVelocity(
- unsigned ip, LocalMatrixType& darcy_velocity_at_ips,
- Eigen::Map<const NodalVectorType> const& local_p,
+ unsigned const ip, Eigen::Map<const NodalVectorType> const& local_p,
ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
double const mu, double const rho_g,
- int const gravitational_axis_id);
+ int const gravitational_axis_id,
+ MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips);
};
template <typename LaplacianGravityVelocityCalculator>
@@ -198,28 +201,29 @@ 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 computeDarcyVelocity(
+ Eigen::MatrixXd const& permeability,
+ std::vector<double> const& local_x,
+ MatrixOfVelocityAtIntegrationPoints& 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;
+ MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;
template <typename LaplacianGravityVelocityCalculator>
void computeDarcyVelocityLocal(
std::vector<double> const& local_x,
Eigen::MatrixXd const& permeability,
- LocalMatrixType& darcy_velocity_at_ips) const;
+ MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;
template <typename LaplacianGravityVelocityCalculator>
void computeDarcyVelocityCoupledWithHeatTransportLocal(
std::vector<double> const& local_x,
std::vector<double> const& local_T,
Eigen::MatrixXd const& permeability,
- LocalMatrixType& darcy_velocity_at_ips) const;
+ MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;
const int _gravitational_axis_id;
const double _gravitational_acceleration;
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index d398f7fc5cbdf18762a377b9c935605a7f41febf..74b635d78eed87825ff90a8c39d878615c3f4538 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -66,7 +66,7 @@ void LiquidFlowProcess::initializeConcreteProcess(
pv.getShapeFunctionOrder(), _local_assemblers,
mesh.isAxiallySymmetric(), integration_order, _gravitational_axis_id,
_gravitational_acceleration, _reference_temperature,
- *_material_properties, this->_coupling_term);
+ *_material_properties, _coupling_term);
_secondary_variables.addSecondaryVariable(
"darcy_velocity",
diff --git a/ProcessLib/RichardsFlow/RichardsFlowFEM.h b/ProcessLib/RichardsFlow/RichardsFlowFEM.h
index 53c3b1c5eb78e3a1e49559326496e49e2af16b8c..6d1ba5fe8b14342127cce3efde0a6526bb6a565f 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowFEM.h
+++ b/ProcessLib/RichardsFlow/RichardsFlowFEM.h
@@ -251,7 +251,6 @@ public:
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(