diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp index 9800d5b3423ad4dfb80989c0b0c118b0c2865587..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();