diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 5567b0ffd4af4d3b7e72c2e0bbd70ac10831d262..e47bb62d4e5e976c32774b5225e070a1018c909a 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -240,5 +240,49 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
local_b_data, local_Jac_data, local_coupled_solutions);
}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ postNonLinearSolverConcrete(std::vector<double> const& local_x,
+ double const t,
+ bool const use_monolithic_scheme)
+{
+ const int offset = use_monolithic_scheme ? displacement_index : 0;
+
+ auto u =
+ Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
+ displacement_size> const>(local_x.data() + offset,
+ displacement_size);
+ double const& dt = _process_data.dt;
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ int const n_integration_points = _integration_method.getNumberOfPoints();
+ for (int ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& N_u = _ip_data[ip].N_u;
+ auto const& dNdx_u = _ip_data[ip].dNdx_u;
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(_element,
+ N_u);
+ auto const B =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunctionDisplacement::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx_u, N_u, x_coord, _is_axially_symmetric);
+
+ auto& eps = _ip_data[ip].eps;
+ eps.noalias() = B * u;
+
+ _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u);
+ }
+}
+
} // namespace HydroMechanics
} // namespace ProcessLib
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index 4536c4944b379a6121de2c35c928a8e482c2935f..253e4062d24b8b6a61a4d215c52fc21a60d98efc 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -376,6 +376,10 @@ public:
}
}
+ void postNonLinearSolverConcrete(std::vector<double> const& local_x,
+ double const t,
+ bool const use_monolithic_scheme) override;
+
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h
index d9804f7ce7e467155d7fb95c79399fb8c6362ee2..499707da66ca6c3423a50ee59fc9e0a8398b6041 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h
@@ -91,7 +91,7 @@ void HydroMechanicsProcess<DisplacementDim>::constructDofTable()
std::make_unique<NumLib::LocalToGlobalIndexMap>(
std::move(all_mesh_subsets_single_component),
// by location order is needed for output
- NumLib::ComponentOrder::BY_LOCATION);
+ NumLib::ComponentOrder::BY_LOCATION);
if (_use_monolithic_scheme)
{
@@ -135,13 +135,12 @@ void HydroMechanicsProcess<DisplacementDim>::constructDofTable()
// For pressure equation.
// Collect the mesh subsets with base nodes in a vector.
std::vector<MeshLib::MeshSubsets> all_mesh_subsets_base_nodes;
- all_mesh_subsets_base_nodes.emplace_back(
- _mesh_subset_base_nodes.get());
+ all_mesh_subsets_base_nodes.emplace_back(_mesh_subset_base_nodes.get());
_local_to_global_index_map_with_base_nodes =
std::make_unique<NumLib::LocalToGlobalIndexMap>(
std::move(all_mesh_subsets_base_nodes),
// by location order is needed for output
- NumLib::ComponentOrder::BY_LOCATION);
+ NumLib::ComponentOrder::BY_LOCATION);
_sparsity_pattern_with_linear_element = NumLib::computeSparsityPattern(
*_local_to_global_index_map_with_base_nodes, _mesh);
@@ -168,72 +167,62 @@ void HydroMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
Base::_secondary_variables.addSecondaryVariable(
"sigma_xx",
- makeExtrapolator(
- 1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaXX));
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigmaXX));
Base::_secondary_variables.addSecondaryVariable(
"sigma_yy",
- makeExtrapolator(
- 1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaYY));
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigmaYY));
Base::_secondary_variables.addSecondaryVariable(
"sigma_zz",
- makeExtrapolator(
- 1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaZZ));
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigmaZZ));
Base::_secondary_variables.addSecondaryVariable(
"sigma_xy",
- makeExtrapolator(
- 1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaXY));
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigmaXY));
if (DisplacementDim == 3)
{
Base::_secondary_variables.addSecondaryVariable(
"sigma_xz",
- makeExtrapolator(
- 1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaXZ));
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigmaXZ));
Base::_secondary_variables.addSecondaryVariable(
"sigma_yz",
- makeExtrapolator(
- 1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtSigmaYZ));
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtSigmaYZ));
}
Base::_secondary_variables.addSecondaryVariable(
"epsilon_xx",
- makeExtrapolator(
- 1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonXX));
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtEpsilonXX));
Base::_secondary_variables.addSecondaryVariable(
"epsilon_yy",
- makeExtrapolator(
- 1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonYY));
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtEpsilonYY));
Base::_secondary_variables.addSecondaryVariable(
"epsilon_zz",
- makeExtrapolator(
- 1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonZZ));
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtEpsilonZZ));
Base::_secondary_variables.addSecondaryVariable(
"epsilon_xy",
- makeExtrapolator(
- 1, getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtEpsilonXY));
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &LocalAssemblerInterface::getIntPtEpsilonXY));
Base::_secondary_variables.addSecondaryVariable(
"velocity",
- makeExtrapolator(
- mesh.getDimension(), getExtrapolator(), _local_assemblers,
- &LocalAssemblerInterface::getIntPtDarcyVelocity));
+ makeExtrapolator(mesh.getDimension(), getExtrapolator(),
+ _local_assemblers,
+ &LocalAssemblerInterface::getIntPtDarcyVelocity));
}
template <int DisplacementDim>
@@ -302,7 +291,7 @@ void HydroMechanicsProcess<DisplacementDim>::
{
DBUG(
"Assemble the Jacobian equations of liquid fluid process in "
- "HydroMechanics for the staggered scheme.");
+ "HydroMechanics for the staggered scheme.");
}
else
{
@@ -313,8 +302,8 @@ void HydroMechanicsProcess<DisplacementDim>::
GlobalExecutor::executeMemberDereferenced(
_global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
- _local_assemblers, *_local_to_global_index_map, t, x,
- xdot, dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions,
+ _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
+ dx_dx, M, K, b, Jac, _coupled_solutions,
_local_to_global_index_map_with_base_nodes.get());
}
@@ -331,9 +320,9 @@ void HydroMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
// If monolithic scheme is used or the equation of deformation is solved in
// the staggered scheme.
if (_use_monolithic_scheme || process_id == 1)
- GlobalExecutor::executeMemberOnDereferenced(
- &LocalAssemblerInterface::preTimestep, _local_assemblers,
- *_local_to_global_index_map, x, t, dt);
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerInterface::preTimestep, _local_assemblers,
+ *_local_to_global_index_map, x, t, dt);
}
template <int DisplacementDim>
@@ -347,16 +336,32 @@ void HydroMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
}
template <int DisplacementDim>
-NumLib::LocalToGlobalIndexMap* HydroMechanicsProcess<DisplacementDim>::
- getDOFTableForExtrapolatorData(bool& manage_storage) const
+void HydroMechanicsProcess<DisplacementDim>::postNonLinearSolverProcess(
+ GlobalVector const& x, const double t, const int process_id)
+{
+ if (!_use_monolithic_scheme && process_id == 0)
+ {
+ return;
+ }
+
+ DBUG("PostNonLinearSolver HydroMechanicsProcess.");
+ // Calculate strain, stress or other internal variables of mechanics.
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerInterface::postNonLinearSolver, _local_assemblers,
+ getDOFTable(process_id), x, t, _use_monolithic_scheme);
+}
+
+template <int DisplacementDim>
+NumLib::LocalToGlobalIndexMap*
+HydroMechanicsProcess<DisplacementDim>::getDOFTableForExtrapolatorData(
+ bool& manage_storage) const
{
return _local_to_global_index_map_single_component.get();
}
template <int DisplacementDim>
NumLib::LocalToGlobalIndexMap const&
-HydroMechanicsProcess<DisplacementDim>::getDOFTable(
- const int process_id) const
+HydroMechanicsProcess<DisplacementDim>::getDOFTable(const int process_id) const
{
// If monolithic scheme is used or the equation of deformation is solved in
// the staggered scheme.
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
index 1b773f2f90708361374ee8e03340341e0dda78a3..cbe0b1eac11edbd531f9a80d63dd47c72b888e37 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
@@ -77,6 +77,9 @@ private:
void postTimestepConcreteProcess(GlobalVector const& x,
int const process_id) override;
+ void postNonLinearSolverProcess(GlobalVector const& x, const double t,
+ int const process_id) override;
+
NumLib::LocalToGlobalIndexMap const& getDOFTable(
const int process_id) const override;
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index ab89bd751baec57fc6815f19ff1647951ee9ac4f..937cfbd387a74b5516b0ba7608a3d197619ba3a8 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -100,4 +100,15 @@ void LocalAssemblerInterface::postTimestep(
postTimestepConcrete(local_x);
}
+void LocalAssemblerInterface::postNonLinearSolver(
+ std::size_t const mesh_item_id,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ GlobalVector const& x, double const t, bool const use_monolithic_scheme)
+{
+ auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
+ auto const local_x = x.get(indices);
+
+ postNonLinearSolverConcrete(local_x, t, use_monolithic_scheme);
+}
+
} // namespace ProcessLib
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index a99e41e26d1f2a74bdc3e92712de7b2bc7ce2c60..94e0d82dd3ae96c17556c2f35c9f6c82f5890e91 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -78,6 +78,11 @@ public:
NumLib::LocalToGlobalIndexMap const& dof_table,
GlobalVector const& x);
+ virtual void postNonLinearSolver(std::size_t const mesh_item_id,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ GlobalVector const& x, double const t,
+ bool const use_monolithic_scheme);
+
/// Computes the flux in the point \c p_local_coords that is given in local
/// coordinates using the values from \c local_x.
virtual std::vector<double> getFlux(
@@ -94,6 +99,13 @@ private:
}
virtual void postTimestepConcrete(std::vector<double> const& /*local_x*/) {}
+
+ virtual void postNonLinearSolverConcrete(
+ std::vector<double> const& /*local_x*/, double const /*t*/,
+ bool const /*use_monolithic_scheme*/)
+ {
+ }
+
virtual void computeSecondaryVariableConcrete(
double const /*t*/, std::vector<double> const& /*local_x*/)
{
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 644d4575e436c25d257d45be108dd58b808383c1..05f41c12a4a7a1d613772886e4722e58fcb5b0b7 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -129,8 +129,7 @@ void Process::setInitialConditions(const int process_id, double const t,
for (int component_id = 0; component_id < num_comp; ++component_id)
{
auto const& mesh_subsets =
- dof_table_of_process.getMeshSubsets(variable_id,
- component_id);
+ dof_table_of_process.getMeshSubsets(variable_id, component_id);
for (auto const& mesh_subset : mesh_subsets)
{
auto const mesh_id = mesh_subset->getMeshID();
@@ -350,6 +349,13 @@ void Process::postTimestep(GlobalVector const& x, int const process_id)
postTimestepConcreteProcess(x, process_id);
}
+void Process::postNonLinearSolver(GlobalVector const& x, const double t,
+ int const process_id)
+{
+ MathLib::LinAlg::setLocalAccessibleVector(x);
+ postNonLinearSolverProcess(x, t, process_id);
+}
+
void Process::computeSecondaryVariable(const double t, GlobalVector const& x)
{
MathLib::LinAlg::setLocalAccessibleVector(x);
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 731cefeb64b32276ba899e04e618421df4acf666..955f23ca181f43aada7d6975771c310f792198c7 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -59,6 +59,9 @@ public:
/// Postprocessing after a complete timestep.
void postTimestep(GlobalVector const& x, int const process_id);
+ void postNonLinearSolver(GlobalVector const& x, const double t,
+ int const process_id);
+
void preIteration(const unsigned iter, GlobalVector const& x) final;
/// compute secondary variables for the coupled equations or for output.
@@ -179,6 +182,13 @@ private:
int const /*process_id*/)
{
}
+
+ virtual void postNonLinearSolverProcess(GlobalVector const& /*x*/,
+ const double /*t*/,
+ int const /*process_id*/)
+ {
+ }
+
virtual void preIterationConcreteProcess(const unsigned /*iter*/,
GlobalVector const& /*x*/)
{
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index 4df5dbbc22de85f3b4c4b98a36219b582a36f439..52d438ac36c1234bfa38b6f2e3464216d748bcb2 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -422,7 +422,7 @@ std::vector<GlobalVector*> setInitialConditions(
return process_solutions;
}
-bool solveOneTimeStepOneProcess(const unsigned process_index, GlobalVector& x,
+bool solveOneTimeStepOneProcess(int const process_id, GlobalVector& x,
std::size_t const timestep, double const t,
double const delta_t,
SingleProcessData const& process_data,
@@ -448,7 +448,7 @@ bool solveOneTimeStepOneProcess(const unsigned process_index, GlobalVector& x,
auto const post_iteration_callback = [&](unsigned iteration,
GlobalVector const& x) {
- output_control.doOutputNonlinearIteration(process, process_index,
+ output_control.doOutputNonlinearIteration(process, process_id,
process_data.process_output,
timestep, t, x, iteration);
};
@@ -456,6 +456,11 @@ bool solveOneTimeStepOneProcess(const unsigned process_index, GlobalVector& x,
bool nonlinear_solver_succeeded =
nonlinear_solver.solve(x, post_iteration_callback);
+ if (nonlinear_solver_succeeded)
+ {
+ process.postNonLinearSolver(x, t, process_id);
+ }
+
return nonlinear_solver_succeeded;
}