Commit 76a6ec29 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#531 GlobalVariationalOperator: injection is also more generic (even if it...

#531 GlobalVariationalOperator: injection is also more generic (even if it does handle only matrices so far). So now we should be able to assemble into several amtrix (this will be tested immediately).
parent a88c16b0
......@@ -408,18 +408,7 @@ namespace HappyHeart
//! Whether the gradient of finite element is required or not.
AllocateGradientFEltPhi DoAllocateGradientFEltPhi() const;
template<class LinearAlgebraTupleT>
void InjectIntoGlobalLinearAlgebra531(const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_operator,
const LinearAlgebraTupleT& global_linear_algebra) const;
void InjectIntoGlobalLinearAlgebra531(const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_operator,
GlobalMatrixWithCoefficient& global_linear_algebra) const; // \todo #531 remove it
private:
......
......@@ -299,9 +299,7 @@ namespace HappyHeart
local_operator,
std::forward_as_tuple(args...));
InjectIntoGlobalLinearAlgebra531(local_felt_space,
local_operator,
std::get<0>(linear_algebra_tuple));
Private::InjectInGlobalMatrix531(linear_algebra_tuple, local_felt_space, local_operator);
}
}
......@@ -511,63 +509,6 @@ namespace HappyHeart
}
template
<
class DerivedT,
class LocalVariationalOperatorT,
class UnknownPolicyT
>
void GlobalVariationalOperator<DerivedT, LocalVariationalOperatorT, UnknownPolicyT>
::InjectIntoGlobalLinearAlgebra531(const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_operator,
GlobalMatrixWithCoefficient& global_linear_algebra) const
{
auto& elementary_data = local_operator.GetNonCstElementaryData();
auto& local_matrix = elementary_data.GetNonCstMatrixResult();
auto& global_matrix = global_linear_algebra.first;
const auto& unknown = local_operator.GetUnknown();
const auto& unknown_and_numbering_subset = elementary_data.GetRefFElt(unknown).GetUnknownAndNumberingSubset();
const auto& local2global = local_felt_space.GetFiniteElt(unknown_and_numbering_subset).template GetLocal2Global<MpiScale::program_wise>();
Private::ReportValuesInGlobalMatrix(local_matrix, global_linear_algebra.second, local2global, local2global, global_matrix);
}
template
<
class DerivedT,
class LocalVariationalOperatorT,
class UnknownPolicyT
>
template<class LinearAlgebraTupleT>
void GlobalVariationalOperator<DerivedT, LocalVariationalOperatorT, UnknownPolicyT>
::InjectIntoGlobalLinearAlgebra531(const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_operator,
const LinearAlgebraTupleT& linear_algebra_tuple) const
{
auto& elementary_data = local_operator.GetNonCstElementaryData();
const auto& unknown = local_operator.GetUnknown();
const auto& unknown_and_numbering_subset = elementary_data.GetRefFElt(unknown).GetUnknownAndNumberingSubset();
const auto& local2global = local_felt_space.GetFiniteElt(unknown_and_numbering_subset).template GetLocal2Global<MpiScale::program_wise>();
auto& local_matrix = elementary_data.GetNonCstMatrixResult();
Private::InjectInGlobalMatrix531(linear_algebra_tuple, local_operator);
// auto& global_matrix = global_linear_algebra.first;
// Private::ReportValuesInGlobalMatrix(local_matrix, global_linear_algebra.second, local2global, local2global, global_matrix);
}
} // namespace HappyHeart
......
......@@ -50,35 +50,60 @@ namespace HappyHeart
auto& elementary_data = local_variational_operator.GetNonCstElementaryData();
auto& local_matrix = elementary_data.GetNonCstMatrixResult();
// \todo #531 We might need here to reduce the local matrix!
const double coefficient = std::get<I>(linear_algebra_tuple).second;
if (!Utilities::AreEqual(coefficient, 1.))
local_matrix *= coefficient;
const auto& unknown = local_variational_operator.GetUnknown();
const auto& unknown_and_numbering_subset = elementary_data.GetRefFElt(unknown).GetUnknownAndNumberingSubset();
const auto& local2global = local_felt_space.GetFiniteElt(unknown_and_numbering_subset).template GetLocal2Global<MpiScale::program_wise>();
const auto& row_numbering_subset = global_matrix.GetRowNumberingSubset();
const auto& col_numbering_subset = global_matrix.GetColNumberingSubset();
const auto& unknown_and_numbering_subset_list =
local_variational_operator.GetUnknownAndNumberingSubsetList();
decltype(auto) row_local_2_global =
local_felt_space.ComputeLocal2Global<MpiScale::program_wise>(unknown_and_numbering_subset_list,
row_numbering_subset);
assert(static_cast<int>(row_local_2_global.size()) == local_matrix.GetM());
if (row_numbering_subset == col_numbering_subset)
{
assert(static_cast<int>(row_local_2_global.size()) == local_matrix.GetN());
global_matrix.SetValues(row_local_2_global,
row_local_2_global,
local_matrix.GetData(),
ADD_VALUES,
__FILE__, __LINE__);
}
else
{
decltype(auto) col_local_2_global =
local_felt_space.ComputeLocal2Global<MpiScale::program_wise>(unknown_and_numbering_subset_list,
col_numbering_subset);
assert(static_cast<int>(col_local_2_global.size()) == local_matrix.GetN());
global_matrix.SetValues(row_local_2_global,
col_local_2_global,
local_matrix.GetData(),
ADD_VALUES,
__FILE__, __LINE__);
}
// #ifndef NDEBUG
// assert(static_cast<int>(row_local_2_global.size()) == local_matrix.GetM());
// assert(static_cast<int>(col_local_2_global.size()) == local_matrix.GetN());
//
// assert(global_matrix.GetRowNumberingSubset() == row_unknown_and_numbering_subset.GetNumberingSubset());
// assert(global_matrix.GetColNumberingSubset() == col_unknown_and_numbering_subset.GetNumberingSubset());
// #endif // NDEBUG
//
//
// global_matrix.SetValues(row_local_2_global,
// column_local_2_global,
// local_matrix.GetData(),
// ADD_VALUES,
// __FILE__, __LINE__);
if (!Utilities::AreEqual(coefficient, 1.)) // \todo #531 Be clever and get rid of it (use value from previous iteration in the first product).
local_matrix *= 1. / coefficient;
Recursivity<LinearAlgebraTupleT, I + 1, TupleSizeT>::InjectInGlobalMatrix531(linear_algebra_tuple);
Recursivity<LinearAlgebraTupleT, I + 1, TupleSizeT>::InjectInGlobalMatrix531(linear_algebra_tuple,
local_felt_space,
local_variational_operator);
}
......@@ -87,6 +112,7 @@ namespace HappyHeart
::Assembly531(const LinearAlgebraTupleT& linear_algebra_tuple)
{
// End recursion.
static_cast<void>(linear_algebra_tuple);
}
......@@ -98,6 +124,9 @@ namespace HappyHeart
LocalVariationalOperatorT& local_variational_operator)
{
// End recursion.
static_cast<void>(linear_algebra_tuple);
static_cast<void>(local_felt_space);
static_cast<void>(local_variational_operator);
}
......
......@@ -30,12 +30,12 @@ namespace HappyHeart
template<class LinearAlgebraTupleT, class LocalVariationalOperatorT>
void InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_variational_operator)
{
Impl::Recursivity<LinearAlgebraTupleT, 0, std::tuple_size<LinearAlgebraTupleT>::value>
::InjectInGlobalMatrix531(linear_algebra_tuple, local_variational_operator);
::InjectInGlobalMatrix531(linear_algebra_tuple, local_felt_space, local_variational_operator);
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment