diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
index 41740a1a19e47b1d749bcdf1d8c31b8d0320678c..0ff22ee4240061ac23b90624661d861473051b33 100644
--- a/Applications/CLI/Tests.cmake
+++ b/Applications/CLI/Tests.cmake
@@ -598,7 +598,8 @@ if(NOT OGS_USE_MPI)
TESTER vtkdiff
ABSTOL 1e-8 RELTOL 1e-8
DIFF_DATA
- sat1D.vtu sat_1D_pcs_0_ts_1_t_1.000000.vtu AnalyticSolution pressure
+ sat1D.vtu sat_1D_pcs_0_ts_1_t_1.000000.vtu AnalyticPressure pressure
+ sat1D.vtu sat_1D_pcs_0_ts_1_t_1.000000.vtu AnalyticVec v_x
)
AddTest(
NAME LiquidFlow_PressureBCatCornerOfAnisotropicSquare
@@ -620,7 +621,9 @@ if(NOT OGS_USE_MPI)
TESTER vtkdiff
ABSTOL 1e-8 RELTOL 1e-8
DIFF_DATA
- mesh2D.vtu gravity_driven_pcs_0_ts_1_t_1.000000.vtu AnalyticSolution pressure
+ mesh2D.vtu gravity_driven_pcs_0_ts_1_t_1.000000.vtu AnalyticPressure pressure
+ mesh2D.vtu gravity_driven_pcs_0_ts_1_t_1.000000.vtu AnalyticVx v_x
+ mesh2D.vtu gravity_driven_pcs_0_ts_1_t_1.000000.vtu AnalyticVy v_y
)
AddTest(
NAME LiquidFlow_AxisymTheis
@@ -877,7 +880,8 @@ else()
TESTER vtkdiff
ABSTOL 1e-8 RELTOL 1e-8
DIFF_DATA
- sat1D.vtu sat_1D_pcs_0_ts_1_t_1.000000.vtu AnalyticSolution pressure
+ sat1D.vtu sat_1D_pcs_0_ts_1_t_1.000000.vtu AnalyticPressure pressure
+ sat1D.vtu sat_1D_pcs_0_ts_1_t_1.000000.vtu AnalyticVec v_x
)
AddTest(
NAME LiquidFlow_GravityDriven
@@ -888,7 +892,9 @@ else()
TESTER vtkdiff
ABSTOL 1e-8 RELTOL 1e-8
DIFF_DATA
- mesh2D.vtu gravity_driven_pcs_0_ts_1_t_1.000000.vtu AnalyticSolution pressure
+ mesh2D.vtu gravity_driven_pcs_0_ts_1_t_1.000000.vtu AnalyticPressure pressure
+ mesh2D.vtu gravity_driven_pcs_0_ts_1_t_1.000000.vtu AnalyticVx v_x
+ mesh2D.vtu gravity_driven_pcs_0_ts_1_t_1.000000.vtu AnalyticVy v_y
)
AddTest(
NAME LiquidFlow_PressureBCatCornerOfAnisotropicSquare
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
index 6e102649406bb6902c04d85023d9e856bd7442d5..e0c902dbb204a22ec6cae232039a1256305b26c7 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -105,11 +105,34 @@ public:
local_K.noalias() += sm.dNdx.transpose() * k * sm.dNdx * sm.detJ *
sm.integralMeasure * wp.getWeight();
+ }
+ }
+
+ void computeSecondaryVariableConcrete(
+ const double t,
+ std::vector<double> const& local_x) override
+ {
+ auto const local_matrix_size = local_x.size();
+ // This assertion is valid only if all nodal d.o.f. use the same shape
+ // matrices.
+ assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ const auto local_x_vec =
+ MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
+
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ pos.setIntegrationPoint(ip);
+ auto const& sm = _shape_matrices[ip];
+ auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
// Darcy velocity only computed for output.
- GlobalDimVectorType const darcy_velocity =
- -k * sm.dNdx * Eigen::Map<const NodalVectorType>(
- local_x.data(), ShapeFunction::NPOINTS);
+ GlobalDimVectorType const darcy_velocity = -k * sm.dNdx * local_x_vec;
for (unsigned d = 0; d < GlobalDim; ++d)
{
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index 4bb40f2c8781b0b6442a2435196184d73abd0224..1bc2bebb302e6ef599aebb3e51a0346e8dc8bbce 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -100,5 +100,15 @@ void GroundwaterFlowProcess::assembleWithJacobianConcreteProcess(
dx_dx, M, K, b, Jac);
}
+
+void GroundwaterFlowProcess::computeSecondaryVariableConcrete(const double t,
+ GlobalVector const& x)
+{
+ DBUG("Compute the velocity for GroundwaterFlowProcess.");
+ GlobalExecutor::executeMemberOnDereferenced(
+ &GroundwaterFlowLocalAssemblerInterface::computeSecondaryVariable,
+ _local_assemblers, *_local_to_global_index_map, t, x);
+}
+
} // namespace GroundwaterFlow
} // namespace ProcessLib
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index 6bfa9c65676864bde8d3965dcbaf8ff452eb6ae8..7d2fc9129e573e7279f4f885e31906d82e7f529f 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -88,6 +88,9 @@ public:
}
}
+ void computeSecondaryVariableConcrete(double const t,
+ GlobalVector const& x) override;
+
private:
void initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table,
diff --git a/ProcessLib/HeatConduction/HeatConductionFEM.h b/ProcessLib/HeatConduction/HeatConductionFEM.h
index c9de972da3151b33b560e3198db3383fb993d0bb..e2fad4f8534ed8a3a245519394645aea51fd2ec2 100644
--- a/ProcessLib/HeatConduction/HeatConductionFEM.h
+++ b/ProcessLib/HeatConduction/HeatConductionFEM.h
@@ -116,10 +116,33 @@ public:
local_M.noalias() += sm.N.transpose() * density * heat_capacity *
sm.N * sm.detJ * wp.getWeight() *
sm.integralMeasure;
+ }
+ }
+
+ void computeSecondaryVariableConcrete(
+ const double t,
+ std::vector<double> const& local_x) override
+ {
+ auto const local_matrix_size = local_x.size();
+ // This assertion is valid only if all nodal d.o.f. use the same shape
+ // matrices.
+ assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+ const auto local_x_vec =
+ MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ pos.setIntegrationPoint(ip);
+ auto const& sm = _shape_matrices[ip];
+ auto const k = _process_data.thermal_conductivity(t, pos)[0];
// heat flux only computed for output.
- GlobalDimVectorType const heat_flux =
- -k * sm.dNdx * Eigen::Map<const NodalVectorType>(
- local_x.data(), ShapeFunction::NPOINTS);
+ GlobalDimVectorType const heat_flux = -k * sm.dNdx * local_x_vec;
for (unsigned d = 0; d < GlobalDim; ++d)
{
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index 617b90ff0ca5d6804205c92587cc09e380e1bfbb..4f1801bb034a2a9f5351623049607e05c074fba4 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -96,5 +96,14 @@ void HeatConductionProcess::assembleWithJacobianConcreteProcess(
dx_dx, M, K, b, Jac);
}
+void HeatConductionProcess::computeSecondaryVariableConcrete(const double t,
+ GlobalVector const& x)
+{
+ DBUG("Compute heat flux for HeatConductionProcess.");
+ GlobalExecutor::executeMemberOnDereferenced(
+ &HeatConductionLocalAssemblerInterface::computeSecondaryVariable,
+ _local_assemblers, *_local_to_global_index_map, t, x);
+}
+
} // namespace HeatConduction
} // namespace ProcessLib
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.h b/ProcessLib/HeatConduction/HeatConductionProcess.h
index 7b754a2d36f94812a91398f66a076838f4ba992c..5d180db20656d9b4e3796033dc1ef44002f142d6 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.h
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.h
@@ -38,7 +38,9 @@ public:
//! @{
bool isLinear() const override { return true; }
- //! @}
+
+ void computeSecondaryVariableConcrete(double const t,
+ GlobalVector const& x) override;
private:
void initializeConcreteProcess(
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 34c42024fdb6badc912870991fdd1e4e6221b3b1..3d8720aab9a3afc1edfe4be68e33ecf337104868 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 LaplacianGravityVelocityCalculator>
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(
+ LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
local_K, local_b, sm, perm, integration_factor, mu, rho_g,
_gravitational_axis_id);
}
@@ -109,7 +109,68 @@ 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 LaplacianGravityVelocityCalculator>
+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);
+
+ LaplacianGravityVelocityCalculator
+ ::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 +194,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,
@@ -150,6 +232,29 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
fac * rho_g * sm.dNdx.transpose() * perm.col(gravitational_axis_id);
}
}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ 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)
+{
+ // Compute the velocity
+ GlobalDimVectorType darcy_velocity = -perm * sm.dNdx * local_p / 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];
+ }
+}
+
} // end of namespace
} // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index 6b34deaa09d3a3d943147f8d9bb0672be49daf96..0073c353e923a62e6cddd488fd3e6cc224c00959 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -85,6 +85,9 @@ 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;
+
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
@@ -130,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 LaplacianGravityVelocityCalculator>
void local_assemble(double const t, std::vector<double> const& local_x,
std::vector<double>& local_M_data,
std::vector<double>& local_K_data,
@@ -164,6 +179,12 @@ private:
SpatialPosition const& pos,
Eigen::MatrixXd const& perm);
+ template <typename LaplacianGravityVelocityCalculator>
+ 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 e06ff0db8b1c4ac87d1de67aaf44c293d57b9d39..795fe907b11880d7c670cce54b7da6937cb31d7f 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -113,5 +113,14 @@ void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
dx_dx, M, K, b, Jac);
}
+void LiquidFlowProcess::computeSecondaryVariableConcrete(const double t,
+ GlobalVector const& x)
+{
+ DBUG("Compute the velocity for LiquidFlowProcess.");
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LiquidFlowLocalAssemblerInterface::computeSecondaryVariable,
+ _local_assemblers, *_local_to_global_index_map, t, x);
+}
+
} // end of namespace
} // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index ff1a5aac2052550c6cfd6c6cf592925b79c09ffc..60f0151e3f925b9a32ddacc55da4b2ac23d60421 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -69,6 +69,9 @@ public:
double const gravitational_acceleration,
BaseLib::ConfigTree const& config);
+ void computeSecondaryVariableConcrete(double const t,
+ GlobalVector const& x) override;
+
bool isLinear() const override { return true; }
private:
void initializeConcreteProcess(
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index 5e572b76b0521637cdaedf49b57dfdbfc2f3bd43..c35897c1f4a981ebe1348ede53bdf09486247c69 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -26,6 +26,16 @@ void LocalAssemblerInterface::assembleWithJacobian(
"assembler.");
}
+void LocalAssemblerInterface::computeSecondaryVariable(
+ std::size_t const mesh_item_id,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ double const t, GlobalVector const& x)
+{
+ auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
+ auto const local_x = x.get(indices);
+ computeSecondaryVariableConcrete(t, local_x);
+}
+
void LocalAssemblerInterface::preTimestep(
std::size_t const mesh_item_id,
NumLib::LocalToGlobalIndexMap const& dof_table, GlobalVector const& x,
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index fc0520fd0005a190e9c8a1ec06fbdcdf490f9322..959ab56d9f5ea91fefb62d7660d6797032c2f324 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -45,6 +45,10 @@ public:
std::vector<double>& local_b_data,
std::vector<double>& local_Jac_data);
+ virtual void computeSecondaryVariable(std::size_t const mesh_item_id,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ const double t, GlobalVector const& x);
+
virtual void preTimestep(std::size_t const mesh_item_id,
NumLib::LocalToGlobalIndexMap const& dof_table,
GlobalVector const& x, double const t,
@@ -70,6 +74,9 @@ private:
}
virtual void postTimestepConcrete(std::vector<double> const& /*local_x*/) {}
+
+ virtual void computeSecondaryVariableConcrete
+ (double const /*t*/, std::vector<double> const& /*local_x*/) {}
};
} // namespace ProcessLib
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 27e20af297e6476cb3d73ea6d0237f84d03d32d7..0e906ca47769d21f685cbf24121f2bf4a167544b 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -242,6 +242,11 @@ void Process::postTimestep(GlobalVector const& x)
postTimestepConcreteProcess(x);
}
+void Process::computeSecondaryVariable( const double t, GlobalVector const& x)
+{
+ computeSecondaryVariableConcrete(t, x);
+}
+
void Process::preIteration(const unsigned iter, const GlobalVector &x)
{
// In every new iteration cached values of secondary variables are expired.
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 7ee405f8d4ac0456d12ab1b3e3dd1e42e47bc539..b8ed73c921ca7e7f57b4f687fef8b89cfffd2377 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -59,6 +59,9 @@ public:
void preIteration(const unsigned iter,
GlobalVector const& x) override final;
+ /// compute secondary variables for the coupled equations or for output.
+ void computeSecondaryVariable(const double t, GlobalVector const& x);
+
NumLib::IterationResult postIteration(GlobalVector const& x) override final;
void initialize();
@@ -148,6 +151,9 @@ private:
virtual void preIterationConcreteProcess(const unsigned /*iter*/,
GlobalVector const& /*x*/){}
+ virtual void computeSecondaryVariableConcrete(const double /*t*/,
+ GlobalVector const& /*x*/) {}
+
virtual NumLib::IterationResult postIterationConcreteProcess(
GlobalVector const& /*x*/)
{
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index aa5548d1a6c519be479558af526dd1f677b07da7..8ba27a8827246727f5c1441983b9d7b74ff4de43 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -505,6 +505,8 @@ bool UncoupledProcessesTimeLoop::loop()
nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
x, timestep, t, delta_t, *spd, *_output);
+ pcs.computeSecondaryVariable(t, x);
+
INFO("[time] Solving process #%u took %g s in timestep #%u.",
timestep, time_timestep.elapsed(), pcs_idx);
diff --git a/Tests/Data b/Tests/Data
index 95d086de5e3c6db78c07002cec395809a157ad3b..8254884e43b6ac088e70ce74e23b7c0b72e886e6 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 95d086de5e3c6db78c07002cec395809a157ad3b
+Subproject commit 8254884e43b6ac088e70ce74e23b7c0b72e886e6