diff --git a/MathLib/LinAlg/LinAlg.cpp b/MathLib/LinAlg/LinAlg.cpp
index ceb5fc991b0d7ffaa6c4c3bc75f2c8d1c8b08499..fc6432473e954f9066913522f60f49677c2ee89e 100644
--- a/MathLib/LinAlg/LinAlg.cpp
+++ b/MathLib/LinAlg/LinAlg.cpp
@@ -308,6 +308,8 @@ void finalizeAssembly(EigenMatrix& x)
x.getRawMatrix().makeCompressed();
}
+void finalizeAssembly(EigenVector& x) {}
+
} // namespace LinAlg
} // namespace MathLib
diff --git a/MathLib/LinAlg/LinAlg.h b/MathLib/LinAlg/LinAlg.h
index 16c344376b2c79ce866301af5f9044bb4b893400..f03c4119994e85e1052c1c4cf08bc9bbfb4cb428 100644
--- a/MathLib/LinAlg/LinAlg.h
+++ b/MathLib/LinAlg/LinAlg.h
@@ -98,11 +98,8 @@ double norm(MatrixOrVector const& x, MathLib::VecNormType type)
}
}
-template<typename Matrix>
-void finalizeAssembly(Matrix& /*A*/)
-{
- // By default do nothing.
-}
+template <typename Matrix>
+void finalizeAssembly(Matrix& /*A*/);
// Matrix and Vector
@@ -260,6 +257,7 @@ void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
EigenVector const& v2, EigenVector& v3);
void finalizeAssembly(EigenMatrix& A);
+void finalizeAssembly(EigenVector& A);
} // namespace LinAlg
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
index ad0fc2308430608b225c69ce7cd81e6a20ef3b10..656051d76b55e84fd9c54398f0dcc654fcd5fe16 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
@@ -47,12 +47,39 @@ void LocalLinearLeastSquaresExtrapolator::extrapolate(
{
auto const num_nodal_dof_result =
_dof_table_single_component.dofSizeWithoutGhosts() * num_components;
+
+ std::vector<GlobalIndexType> ghost_indices;
+ { // Create num_components times version of ghost_indices arranged by
+ // location. For example for 3 components and ghost_indices {5,6,10} we
+ // compute {15, 16, 17, 18, 19, 20, 30, 31, 32}.
+ auto const& single_component_ghost_indices =
+ _dof_table_single_component.getGhostIndices();
+ auto const single_component_ghost_indices_size =
+ single_component_ghost_indices.size();
+ ghost_indices.reserve(single_component_ghost_indices_size *
+ num_components);
+ for (unsigned i = 0; i < single_component_ghost_indices_size; ++i)
+ {
+ for (unsigned c = 0; c < num_components; ++c)
+ {
+ ghost_indices.push_back(
+ single_component_ghost_indices[i] * num_components + c);
+ }
+ }
+ }
+
if (!_nodal_values ||
- static_cast<std::size_t>(_nodal_values->size()) != num_nodal_dof_result)
+#ifdef USE_PETSC
+ static_cast<std::size_t>(_nodal_values->getLocalSize() +
+ _nodal_values->getGhostSize())
+#else
+ _nodal_values->size()
+#endif
+ != num_nodal_dof_result)
{
_nodal_values = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
- MathLib::MatrixSpecifications{
- num_nodal_dof_result, num_nodal_dof_result, nullptr, nullptr});
+ {num_nodal_dof_result, num_nodal_dof_result, &ghost_indices,
+ nullptr});
}
_nodal_values->setZero();
@@ -68,6 +95,7 @@ void LocalLinearLeastSquaresExtrapolator::extrapolate(
extrapolateElement(i, num_components, extrapolatables, t,
current_solution, dof_table, *counts);
}
+ MathLib::LinAlg::finalizeAssembly(*_nodal_values);
MathLib::LinAlg::componentwiseDivide(*_nodal_values, *_nodal_values,
*counts);
@@ -104,6 +132,7 @@ void LocalLinearLeastSquaresExtrapolator::calculateResiduals(
calculateResidualElement(i, num_components, extrapolatables, t,
current_solution, dof_table);
}
+ MathLib::LinAlg::finalizeAssembly(*_residuals);
}
void LocalLinearLeastSquaresExtrapolator::extrapolateElement(
@@ -272,6 +301,8 @@ void LocalLinearLeastSquaresExtrapolator::calculateResidualElement(
auto const int_pt_vals_mat =
MathLib::toMatrix(int_pt_vals, num_components, num_int_pts);
+ MathLib::LinAlg::setLocalAccessibleVector(
+ *_nodal_values); // For access in the for-loop.
for (unsigned comp = 0; comp < num_components; ++comp)
{
// filter nodal values of the current element
@@ -279,7 +310,7 @@ void LocalLinearLeastSquaresExtrapolator::calculateResidualElement(
{
// TODO PETSc negative indices?
auto const idx = num_components * global_indices[i] + comp;
- nodal_vals_element[i] = (*_nodal_values)[idx];
+ nodal_vals_element[i] = _nodal_values->get(idx);
}
double const residual = (interpolation_matrix * nodal_vals_element -
diff --git a/ProcessLib/Output/CachedSecondaryVariable.cpp b/ProcessLib/Output/CachedSecondaryVariable.cpp
index c511a3737f6885e9d9d0f16c079a1c9ce7b349ca..f2413cb7c1f0b221d8277b5aa8a6ed716bf8d9a0 100644
--- a/ProcessLib/Output/CachedSecondaryVariable.cpp
+++ b/ProcessLib/Output/CachedSecondaryVariable.cpp
@@ -68,6 +68,9 @@ GlobalVector const& CachedSecondaryVariable::evalFieldNoArgs() const
1, *_extrapolatables, _t, *_current_solution, *_dof_table);
auto const& nodal_values = _extrapolator.getNodalValues();
MathLib::LinAlg::copy(nodal_values, _cached_nodal_values);
+
+ MathLib::LinAlg::setLocalAccessibleVector(
+ _cached_nodal_values); // For access in the getValue()
_needs_recomputation = false;
return nodal_values;
}
diff --git a/ProcessLib/Output/ProcessOutput.cpp b/ProcessLib/Output/ProcessOutput.cpp
index df6fac6affc9d39f92c41d4c9dddc333dc1cd535..2e91b04b03d486d736fe113646ae09e0de079832 100644
--- a/ProcessLib/Output/ProcessOutput.cpp
+++ b/ProcessLib/Output/ProcessOutput.cpp
@@ -10,6 +10,7 @@
#include "ProcessOutput.h"
#include "BaseLib/BuildInfo.h"
+#include "MathLib/LinAlg/LinAlg.h"
#include "MeshLib/IO/VtkIO/VtuInterface.h"
#include "NumLib/DOF/LocalToGlobalIndexMap.h"
@@ -27,7 +28,6 @@ static void addOgsVersion(MeshLib::Mesh& mesh)
BaseLib::BuildInfo::ogs_version.end());
}
-#ifndef USE_PETSC // Not used in PETSc case
static void addSecondaryVariableNodes(
double const t,
GlobalVector const& x,
@@ -55,21 +55,22 @@ static void addSecondaryVariableNodes(
std::unique_ptr<GlobalVector> result_cache;
auto const& nodal_values =
var.fcts.eval_field(t, x, dof_table, result_cache);
- if (nodal_values_mesh.size() !=
- static_cast<std::size_t>(nodal_values.size()))
+#ifdef USE_PETSC
+ std::size_t const global_vector_size =
+ nodal_values.getLocalSize() + nodal_values.getGhostSize();
+#else
+ std::size_t const global_vector_size = nodal_values.size();
+#endif
+ if (nodal_values_mesh.size() != global_vector_size)
{
OGS_FATAL(
"Secondary variable `%s' did not evaluate to the right "
"number of components. Expected: %d, actual: %d.",
- var.name.c_str(), nodal_values_mesh.size(), nodal_values.size());
+ var.name.c_str(), nodal_values_mesh.size(), global_vector_size);
}
// Copy result
- for (GlobalIndexType i = 0; i < nodal_values.size(); ++i)
- {
- assert(!std::isnan(nodal_values[i]));
- nodal_values_mesh[i] = nodal_values[i];
- }
+ nodal_values.copyValues(nodal_values_mesh);
}
static void addSecondaryVariableResiduals(
@@ -103,22 +104,23 @@ static void addSecondaryVariableResiduals(
std::unique_ptr<GlobalVector> result_cache;
auto const& residuals =
var.fcts.eval_residuals(t, x, dof_table, result_cache);
- if (residuals_mesh.size() != static_cast<std::size_t>(residuals.size()))
+#ifdef USE_PETSC
+ std::size_t const global_vector_size =
+ residuals.getLocalSize() + residuals.getGhostSize();
+#else
+ std::size_t const global_vector_size = residuals.size();
+#endif
+ if (residuals_mesh.size() != global_vector_size)
{
OGS_FATAL(
"The residual of secondary variable `%s' did not evaluate to the "
"right number of components. Expected: %d, actual: %d.",
- var.name.c_str(), residuals_mesh.size(), residuals.size());
+ var.name.c_str(), residuals_mesh.size(), global_vector_size);
}
// Copy result
- for (GlobalIndexType i = 0; i < residuals.size(); ++i)
- {
- assert(!std::isnan(residuals[i]));
- residuals_mesh[i] = residuals[i];
- }
+ residuals.copyValues(residuals_mesh);
}
-#endif // USE_PETSC
namespace ProcessLib
{
@@ -212,8 +214,6 @@ void processOutputData(
}
}
-#ifndef USE_PETSC
-
// Secondary variables output
for (auto const& external_variable_name : output_variables)
{
@@ -234,11 +234,6 @@ void processOutputData(
external_variable_name, mesh);
}
}
-#else
- (void)mesh;
- (void)secondary_variables;
- (void)t;
-#endif // USE_PETSC
addIntegrationPointWriter(mesh, integration_point_writer);
}