Commit 6198cdb0 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#531 TransientSource: at the moment roll back to old code. New assemble for...

#531 TransientSource: at the moment roll back to old code. New assemble for matrix: use a separate function to do the operation (so that with overload vectors could also been taken care of). Add static assert to check no stupid operations are required, such as assmbling mass into a vector.
parent f0a2dfc4
......@@ -299,7 +299,7 @@ namespace HappyHeart
local_operator,
std::forward_as_tuple(args...));
Private::InjectInGlobalMatrix531(linear_algebra_tuple, local_felt_space, local_operator);
Private::InjectInGlobalLinearAlgebra531(linear_algebra_tuple, local_felt_space, local_operator);
}
}
......
......@@ -32,20 +32,39 @@ namespace HappyHeart
template<class LinearAlgebraTupleT, std::size_t I, std::size_t TupleSizeT>
struct Recursivity
{
using current_type = typename std::tuple_element<I, LinearAlgebraTupleT>::type;
using current_type = typename std::tuple_element<I, LinearAlgebraTupleT>::type;
static_assert(std::is_same<current_type, GlobalMatrixWithCoefficient>::value
|| std::is_same<current_type, GlobalVectorWithCoefficient>::value
, "Tuple is expected to include only one of those types.");
static void Assembly531(const LinearAlgebraTupleT& linear_algebra_tuple);
template<class LocalVariationalOperatorT>
static void InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_variational_operator);
static void InjectInGlobalLinearAlgebra531(const LinearAlgebraTupleT& linear_algebra_tuple,
const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_variational_operator);
private:
// PREVIOUS COEFFICIENT ASSUME WRONGLY PREVIOUS ELT WAS SAME TYPE (MATRIX OR VECTOR!) \todo #531
template<class LocalVariationalOperatorT>
static void InjectInGlobalLinearAlgebraImpl531(const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_variational_operator,
const GlobalMatrixWithCoefficient& global_matrix,
double previous_coefficient);
template<class LocalVariationalOperatorT>
static void InjectInGlobalLinearAlgebraImpl531(const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_variational_operator,
const GlobalVectorWithCoefficient& global_vector,
double previous_coefficient);
};
......@@ -59,7 +78,7 @@ namespace HappyHeart
static void Assembly531(const LinearAlgebraTupleT& linear_algebra_tuple);
template<class LocalVariationalOperatorT>
static void InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
static void InjectInGlobalLinearAlgebra531(const LinearAlgebraTupleT& linear_algebra_tuple,
const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_variational_operator);
......
......@@ -39,27 +39,61 @@ namespace HappyHeart
template<class LinearAlgebraTupleT, std::size_t I, std::size_t TupleSizeT>
template<class LocalVariationalOperatorT>
void Recursivity<LinearAlgebraTupleT, I, TupleSizeT>
::InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_variational_operator)
::InjectInGlobalLinearAlgebra531(const LinearAlgebraTupleT& linear_algebra_tuple,
const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_variational_operator)
{
static_assert(I < TupleSizeT, "Otherwise endless recursion!");
auto& global_matrix = std::get<I>(linear_algebra_tuple).first;
auto& elementary_data = local_variational_operator.GetNonCstElementaryData();
auto& local_matrix = elementary_data.GetNonCstMatrixResult();
// \todo #531 We might need here to reduce the local matrix!
static_assert(! (
(std::is_same<current_type, GlobalMatrixWithCoefficient>::value
&& LocalVariationalOperatorT::GetMatrixVectorNature() == MatrixVectorNature::vector)
||
(std::is_same<current_type, GlobalVectorWithCoefficient>::value
&& LocalVariationalOperatorT::GetMatrixVectorNature() == MatrixVectorNature::matrix)
)
, "Given linear algebra given for assembling must be compatible with the operator.");
const double coefficient = std::get<I>(linear_algebra_tuple).second;
const double previous_coefficient =
ZeroSpecialCase<LinearAlgebraTupleT, I>::FetchPreviousCoefficient(linear_algebra_tuple);
assert(!Utilities::IsZero(previous_coefficient));
const auto& global_linear_algebra = std::get<I>(linear_algebra_tuple);
InjectInGlobalLinearAlgebraImpl531(local_felt_space,
local_variational_operator,
global_linear_algebra,
previous_coefficient);
Recursivity<LinearAlgebraTupleT, I + 1, TupleSizeT>::InjectInGlobalLinearAlgebra531(linear_algebra_tuple,
local_felt_space,
local_variational_operator);
}
template<class LinearAlgebraTupleT, std::size_t I, std::size_t TupleSizeT>
template<class LocalVariationalOperatorT>
void Recursivity<LinearAlgebraTupleT, I, TupleSizeT>
::InjectInGlobalLinearAlgebraImpl531(const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_variational_operator,
const GlobalMatrixWithCoefficient& global_matrix_with_coefficient,
const double previous_coefficient)
{
auto& global_matrix = global_matrix_with_coefficient.first;
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 = global_matrix_with_coefficient.second;
assert(!Utilities::IsZero(previous_coefficient));
// Trick here: local matrix was multiplied by a coefficient in previous iteration; we here both apply
// the new one and cancel the previous one.
......@@ -106,10 +140,6 @@ namespace HappyHeart
}
Recursivity<LinearAlgebraTupleT, I + 1, TupleSizeT>::InjectInGlobalMatrix531(linear_algebra_tuple,
local_felt_space,
local_variational_operator);
}
......@@ -125,7 +155,7 @@ namespace HappyHeart
template<class LinearAlgebraTupleT, std::size_t TupleSizeT>
template<class LocalVariationalOperatorT>
void Recursivity<LinearAlgebraTupleT, TupleSizeT, TupleSizeT>
::InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
::InjectInGlobalLinearAlgebra531(const LinearAlgebraTupleT& linear_algebra_tuple,
const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_variational_operator)
{
......
......@@ -28,9 +28,9 @@ namespace HappyHeart
template<class LinearAlgebraTupleT, class LocalVariationalOperatorT>
void InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_operator);
void InjectInGlobalLinearAlgebra531(const LinearAlgebraTupleT& linear_algebra_tuple,
const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_operator);
......
......@@ -31,13 +31,13 @@ namespace HappyHeart
template<class LinearAlgebraTupleT, class LocalVariationalOperatorT>
void InjectInGlobalMatrix531(const LinearAlgebraTupleT& linear_algebra_tuple,
const LocalFEltSpace& local_felt_space,
LocalVariationalOperatorT& local_variational_operator)
void InjectInGlobalLinearAlgebra531(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_felt_space, local_variational_operator);
::InjectInGlobalLinearAlgebra531(linear_algebra_tuple, local_felt_space, local_variational_operator);
}
......
......@@ -119,6 +119,9 @@ namespace HappyHeart
double time,
const Domain& domain = Domain()) const;
void Assemble(double coefficient, GlobalVector& global_vector, double time,
const Domain& domain = Domain()) const;
private:
......
......@@ -29,6 +29,16 @@ namespace HappyHeart
}
inline void TransientSource::Assemble(double coefficient,
GlobalVector& global_vector,
double time,
const Domain& domain) const
{
return Parent::AssembleWithVariadicArguments(coefficient, global_vector, domain, time);
}
inline std::tuple<double>
TransientSource
::SupplArgumentsComputeEltArray(const TransientSource::LocalVariationalOperator& local_operator,
......
......@@ -202,22 +202,36 @@ namespace HappyHeart
{
if (newton_iteration == 0)
{
const auto time = transient_parameters.GetTime();
// const auto time = transient_parameters.GetTime();
//
// if (this->template IsOperatorActivated<ForceIndexList::volumic_force>())
// {
// auto&& force_vector = GlobalVectorWithCoefficient(vm.GetNonCstForce(), 1.);
// this->template GetNonCstForceOperator<ForceIndexList::volumic_force>().Assemble(std::make_tuple(std::move(force_vector)),
// time);
// }
//
// if (this->template IsOperatorActivated<ForceIndexList::surfacic_force>())
// {
// auto&& force_vector = GlobalVectorWithCoefficient(vm.GetNonCstForce(), 1.);
// this->template GetNonCstForceOperator<ForceIndexList::surfacic_force>().Assemble(std::make_tuple(std::move(force_vector)),
// time);
//
// }
auto& force_vector = vm.GetNonCstForce();
auto time = transient_parameters.GetTime();
if (this->template IsOperatorActivated<ForceIndexList::volumic_force>())
{
auto&& force_vector = GlobalVectorWithCoefficient(vm.GetNonCstForce(), 1.);
this->template GetNonCstForceOperator<ForceIndexList::volumic_force>().Assemble(std::move(force_vector),
this->template GetNonCstForceOperator<ForceIndexList::volumic_force>().Assemble(1.,
force_vector,
time);
}
if (this->template IsOperatorActivated<ForceIndexList::surfacic_force>())
{
auto&& force_vector = GlobalVectorWithCoefficient(vm.GetNonCstForce(), 1.);
this->template GetNonCstForceOperator<ForceIndexList::surfacic_force>().Assemble(std::move(force_vector),
this->template GetNonCstForceOperator<ForceIndexList::surfacic_force>().Assemble(1.,
force_vector,
time);
}
}
}
else
......
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