Attention une mise à jour du serveur va être effectuée le lundi 17 mai entre 13h et 13h30. Cette mise à jour va générer une interruption du service de quelques minutes.

Commit 0365ad3d authored by DIAZ Jerome's avatar DIAZ Jerome Committed by GILLES Sebastien

#1520 Added a test on each topology for the coordinates of the quadrature...

#1520 Added a test on each topology for the coordinates of the quadrature points computed using reference elements for the mesh inputs.
parent 16a1820a
Geometry file
Geometry file
node id assign
element id assign
coordinates
2
0.00000e+00 0.00000e+00 0.00000e+00
1.00000e+00 0.00000e+00 0.00000e+00
part 1
MeshLabel_1
bar2
1
1 2
MeshVersionFormatted 2
Dimension
2
Vertices
3
0 0 100
1 0 101
0 1 102
Edges
3
1 2 2
2 3 3
3 1 4
Triangles
1
1 2 3 1
End
......@@ -12,8 +12,8 @@
*/
#ifndef MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_WRITE_AT_QUAD_POINT_HPP_
# define MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_WRITE_AT_QUAD_POINT_HPP_
#ifndef MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_GLOBAL_COORDS_QUAD_POINTS_HPP_
# define MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_GLOBAL_COORDS_QUAD_POINTS_HPP_
# include "Parameters/Parameter.hpp"
# include "Parameters/InitParameterFromInputData/InitParameterFromInputData.hpp"
......@@ -118,7 +118,7 @@ namespace MoReFEM
/*!
* \brief Only defined to respect the generic interface expected by the parent, it is not used as there is
* no extraction of data from a global to a local level needed for this class.
*
*
* \param[in] local_operator Local operator in charge of the elementary computation. A mutable work
* variable is actually set in this call.
* \param[in] local_felt_space List of finite elements being considered; all those related to the
......@@ -146,4 +146,4 @@ namespace MoReFEM
# include "OperatorInstances/ParameterOperator//GlobalCoordsQuadPoints.hxx"
#endif // MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_WRITE_AT_QUAD_POINT_HPP_
#endif // MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_GLOBAL_COORDS_QUAD_POINTS_HPP_
......@@ -12,8 +12,8 @@
*/
#ifndef MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_WRITE_AT_QUAD_POINT_HXX_
# define MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_WRITE_AT_QUAD_POINT_HXX_
#ifndef MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_GLOBAL_COORDS_QUAD_POINTS_HXX_
# define MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_GLOBAL_COORDS_QUAD_POINTS_HXX_
namespace MoReFEM
......@@ -49,4 +49,4 @@ namespace MoReFEM
/// @} // addtogroup OperatorInstancesGroup
#endif // MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_WRITE_AT_QUAD_POINT_HXX_
#endif // MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_GLOBAL_COORDS_QUAD_POINTS_HXX_
......@@ -12,8 +12,8 @@
*/
#ifndef MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_LOCAL_x_WRITE_AT_QUAD_POINT_HPP_
# define MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_LOCAL_x_WRITE_AT_QUAD_POINT_HPP_
#ifndef MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_LOCAL_x_GLOBAL_COORDS_QUAD_POINTS_HPP_
# define MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_LOCAL_x_GLOBAL_COORDS_QUAD_POINTS_HPP_
# include "Utilities/LinearAlgebra/Storage/Local/LocalMatrixStorage.hpp"
# include "Utilities/LinearAlgebra/Storage/Local/LocalVectorStorage.hpp"
......@@ -33,7 +33,7 @@ namespace MoReFEM
/*!
* \brief Implementation of local \a GlobalCoordsQuadPoints operator.
* \brief Implementation of local \a GlobalCoordsQuadPoints operator.
* This local operator is in charge of computing the global coordinates of each quadrature point.
*/
class GlobalCoordsQuadPoints final
......@@ -116,4 +116,4 @@ namespace MoReFEM
# include "OperatorInstances/ParameterOperator/Local/GlobalCoordsQuadPoints.hxx"
#endif // MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_LOCAL_x_WRITE_AT_QUAD_POINT_HPP_
#endif // MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_LOCAL_x_GLOBAL_COORDS_QUAD_POINTS_HPP_
......@@ -12,8 +12,8 @@
*/
#ifndef MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_LOCAL_x_WRITE_AT_QUAD_POINT_HXX_
# define MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_LOCAL_x_WRITE_AT_QUAD_POINT_HXX_
#ifndef MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_LOCAL_x_GLOBAL_COORDS_QUAD_POINTS_HXX_
# define MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_LOCAL_x_GLOBAL_COORDS_QUAD_POINTS_HXX_
namespace MoReFEM
......@@ -33,4 +33,4 @@ namespace MoReFEM
/// @} // addtogroup OperatorInstancesGroup
#endif // MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_LOCAL_x_WRITE_AT_QUAD_POINT_HXX_
#endif // MOREFEM_x_OPERATOR_INSTANCES_x_PARAMETER_OPERATOR_x_LOCAL_x_GLOBAL_COORDS_QUAD_POINTS_HXX_
......@@ -235,7 +235,7 @@ namespace MoReFEM
for (unsigned int i = 0u; i < Nquadrature_point; ++i)
{
out << geom_elt_index << ';' << quadrature_rule_name << '_'
<< i << ';' << value_per_quad_pt[i].last_update_index << ';';
<< i << ';' << value_per_quad_pt[i].last_update_index << ';';
if constexpr (TypeT == ParameterNS::Type::scalar)
out << std::setw(12) << std::scientific << std::setprecision(5) << value_per_quad_pt[i].value;
......@@ -243,11 +243,9 @@ namespace MoReFEM
// To get rid of the braces "{ }" of the xtensor output and have more control over the precision.
if constexpr (TypeT == ParameterNS::Type::vector)
{
auto separator = Utilities::EmptyString();
for (const auto vector_value : value_per_quad_pt[i].value)
{
out << separator << std::setw(12) << std::scientific << std::setprecision(5) << vector_value;
separator = ',';
out << std::setw(12) << std::scientific << std::setprecision(5) << vector_value;
}
}
......
/*!
//
// \file
//
//
// Created by Sebastien Gilles <sebastien.gilles@inria.fr> on the Mon, 15 Jun 2015 11:59:41 +0200
// Copyright (c) Inria. All rights reserved.
//
// \ingroup ParametersGroup
// \addtogroup ParametersGroup
// \{
*/
#ifndef MOREFEM_x_PARAMETERS_x_POLICY_x_AT_QUADRATURE_POINT_x_AT_QUADRATURE_POINT_HXX_
# define MOREFEM_x_PARAMETERS_x_POLICY_x_AT_QUADRATURE_POINT_x_AT_QUADRATURE_POINT_HXX_
namespace MoReFEM
{
namespace ParameterNS
{
namespace Policy
{
template<ParameterNS::Type TypeT>
AtQuadraturePoint<TypeT>
::AtQuadraturePoint(const std::string& name,
const Domain& domain,
const QuadratureRulePerTopology& quadrature_rule_per_topology,
storage_value_type initial_value,
const TimeManager& time_manager)
: domain_(domain),
initial_value_(initial_value),
quadrature_rule_per_topology_(quadrature_rule_per_topology),
time_manager_(time_manager)
{
static_cast<void>(name);
auto& storage = GetNonCstStorage();
storage.max_load_factor(Utilities::DefaultMaxLoadFactor());
assert(domain.template IsConstraintOn<DomainNS::Criterion::mesh>()
&& "It is an assert because this should have been checked earlier in ParameterInstance "
"constructor. At that place it is an exception, as a mere user might provide an invalid Domain.");
const auto& mesh = domain.GetMesh();
storage.reserve(static_cast<std::size_t>(mesh.NgeometricElt()));
const auto& ref_geom_elt_list = mesh.BagOfEltType();
for (const auto& ref_geom_elt_ptr : ref_geom_elt_list)
{
assert(!(!ref_geom_elt_ptr));
const auto& ref_geom_elt = *ref_geom_elt_ptr;
const auto& quadrature_rule = GetQuadratureRule(ref_geom_elt.GetTopologyIdentifier());
const auto Nquadrature_point = static_cast<std::size_t>(quadrature_rule.NquadraturePoint());
decltype(auto) iterator_range = mesh.GetSubsetGeometricEltList(ref_geom_elt);
// Iterate over all geometric elements that share the ref_geom_elt.
for (auto it_geom_elt = iterator_range.first; it_geom_elt != iterator_range.second; ++it_geom_elt)
{
const auto& geom_elt_ptr = *it_geom_elt;
assert(!(!geom_elt_ptr));
const auto geom_elt_index = geom_elt_ptr->GetIndex();
using value_type = Internal::ParameterNS::AtQuadraturePointNS::ValueHolder<storage_value_type>;
value_type default_value(initial_value, time_manager.NtimeModified());
auto check = storage.insert({geom_elt_index, std::vector<value_type>(Nquadrature_point, default_value)});
assert(check.second && "A geometric element should be entered only once.");
static_cast<void>(check);
}
}
}
template<ParameterNS::Type TypeT>
inline typename AtQuadraturePoint<TypeT>::return_type
AtQuadraturePoint<TypeT>::GetValueFromPolicy(const QuadraturePoint& quad_pt,
const GeometricElt& geom_elt) const
{
return FindValue(quad_pt, geom_elt).value;
}
template<ParameterNS::Type TypeT>
inline typename AtQuadraturePoint<TypeT>::return_type
AtQuadraturePoint<TypeT>::GetAnyValueFromPolicy() const
{
decltype(auto) storage = GetStorage();
assert(!storage.empty());
decltype(auto) geom_elt_content = storage.cbegin()->second;
assert(!geom_elt_content.empty());
return geom_elt_content.back().value;
}
template<ParameterNS::Type TypeT>
const typename AtQuadraturePoint<TypeT>::value_holder_type&
AtQuadraturePoint<TypeT>::FindValue(const QuadraturePoint& quad_pt,
const GeometricElt& geom_elt) const
{
auto it = storage_.find(geom_elt.GetIndex());
assert(it != storage_.cend());
const auto& values_in_geom_elt = it->second;
const auto quad_pt_index = static_cast<std::size_t>(quad_pt.GetIndex());
assert(quad_pt_index < values_in_geom_elt.size());
return values_in_geom_elt[quad_pt_index];
}
template<ParameterNS::Type TypeT>
inline typename AtQuadraturePoint<TypeT>::value_holder_type&
AtQuadraturePoint<TypeT>::FindNonCstValue(const QuadraturePoint& quad_pt,
const GeometricElt& geom_elt)
{
return const_cast<typename AtQuadraturePoint<TypeT>::value_holder_type&>(FindValue(quad_pt, geom_elt));
}
template<ParameterNS::Type TypeT>
template<class UpdateFunctorT>
typename AtQuadraturePoint<TypeT>::return_type
AtQuadraturePoint<TypeT>::UpdateAndGetValue(const QuadraturePoint& quad_pt,
const GeometricElt& geom_elt,
const UpdateFunctorT& update_functor)
{
auto& value_holder = FindNonCstValue(quad_pt, geom_elt);
value_holder.last_update_index = GetTimeManager().NtimeModified();
update_functor(value_holder.value);
return value_holder.value;
}
template<ParameterNS::Type TypeT>
template<class UpdateFunctorT>
void AtQuadraturePoint<TypeT>::UpdateValue(const QuadraturePoint& quad_pt,
const GeometricElt& geom_elt,
const UpdateFunctorT& update_functor)
{
auto& value_holder = FindNonCstValue(quad_pt, geom_elt);
value_holder.last_update_index = GetTimeManager().NtimeModified();
update_functor(value_holder.value);
}
template<ParameterNS::Type TypeT>
void AtQuadraturePoint<TypeT>::Copy(const AtQuadraturePoint& parameter_at_quad_point)
{
GetNonCstStorage() = parameter_at_quad_point.GetStorage();
}
template<ParameterNS::Type TypeT>
[[noreturn]] typename AtQuadraturePoint<TypeT>::return_type
AtQuadraturePoint<TypeT>::GetConstantValueFromPolicy() const
{
assert(false && "Should yield IsConstant() = false!");
exit(EXIT_FAILURE);
}
template<ParameterNS::Type TypeT>
bool AtQuadraturePoint<TypeT>::IsConstant() const
{
return false;
}
template<ParameterNS::Type TypeT>
void AtQuadraturePoint<TypeT>::WriteFromPolicy(std::ostream& out) const
{
decltype(auto) domain = this->GetDomain();
const auto& mesh = domain.GetMesh();
const auto& ref_geom_elt_list = mesh.BagOfEltType();
const auto& storage = GetStorage();
# ifndef NDEBUG
const auto end_storage = storage.cend();
# endif // NDEBUG
out << "# Geometric element index; quadrature point; TimeManager::NtimesModified() at last update; value" << std::endl;
for (const auto& ref_geom_elt_ptr : ref_geom_elt_list)
{
assert(!(!ref_geom_elt_ptr));
const auto& ref_geom_elt = *ref_geom_elt_ptr;
if (!domain.DoRefGeomEltMatchCriteria(ref_geom_elt))
continue;
const auto& quadrature_rule = GetQuadratureRule(ref_geom_elt.GetTopologyIdentifier());
decltype(auto) iterator_range = mesh.GetSubsetGeometricEltList(ref_geom_elt);
const auto quadrature_rule_name = quadrature_rule.GetName();
const auto Nquadrature_point = quadrature_rule.NquadraturePoint();
const auto Ngeo_elt = std::distance(iterator_range.first, iterator_range.second);
const auto Nvalues = Ngeo_elt * Nquadrature_point;
out << "# Nvalues" << std::endl;
out << Nvalues << std::endl;
for (auto it_geom_elt = iterator_range.first; it_geom_elt != iterator_range.second; ++it_geom_elt)
{
const auto& geom_elt_ptr = *it_geom_elt;
assert(!(!geom_elt_ptr));
const auto& geom_elt = *geom_elt_ptr;
if (!domain.IsGeometricEltInside(geom_elt))
continue;
const auto geom_elt_index = geom_elt.GetIndex();
auto it_in_storage = storage.find(geom_elt_index);
assert(it_in_storage != end_storage);
const auto& value_per_quad_pt = it_in_storage->second;
for (unsigned int i = 0u; i < Nquadrature_point; ++i)
{
out << geom_elt_index << ';' << quadrature_rule_name << '_'
<< i << ';' << value_per_quad_pt[i].last_update_index << ';';
if constexpr (TypeT == ParameterNS::Type::scalar)
out << std::setw(12) << std::scientific << std::setprecision(5) << value_per_quad_pt[i].value;
// To get rid of the braces "{ }" of the xtensor output and have more control over the precision.
if constexpr (TypeT == ParameterNS::Type::vector)
{
for (const auto vector_value : value_per_quad_pt[i].value)
{
out << std::setw(12) << std::scientific << std::setprecision(5) << vector_value;
}
}
out << std::endl;
}
}
}
}
template<ParameterNS::Type TypeT>
inline const Domain& AtQuadraturePoint<TypeT>::GetDomain() const noexcept
{
return domain_;
}
template<ParameterNS::Type TypeT>
inline const typename AtQuadraturePoint<TypeT>::storage_type&
AtQuadraturePoint<TypeT>::GetStorage() const noexcept
{
return storage_;
}
template<ParameterNS::Type TypeT>
inline typename AtQuadraturePoint<TypeT>::storage_type&
AtQuadraturePoint<TypeT>::GetNonCstStorage() noexcept
{
return storage_;
}
template<ParameterNS::Type TypeT>
inline
const QuadratureRule& AtQuadraturePoint<TypeT>
::GetQuadratureRule(RefGeomEltNS::TopologyNS::Type topology_id) const
{
return GetQuadratureRulePerTopology().GetRule(topology_id);
}
template<ParameterNS::Type TypeT>
inline const TimeManager& AtQuadraturePoint<TypeT>::GetTimeManager() const noexcept
{
return time_manager_;
}
template<ParameterNS::Type TypeT>
inline const QuadratureRulePerTopology& AtQuadraturePoint<TypeT>
::GetQuadratureRulePerTopology() const noexcept
{
return quadrature_rule_per_topology_;
}
template<ParameterNS::Type TypeT>
inline void AtQuadraturePoint<TypeT>::SetConstantValue(double value)
{
static_cast<void>(value);
assert(false && "SetConstantValue() meaningless for current Parameter.");
exit(EXIT_FAILURE);
}
} // namespace Policy
} // namespace ParameterNS
} // namespace MoReFEM
/// @} // addtogroup ParametersGroup
#endif // MOREFEM_x_PARAMETERS_x_POLICY_x_AT_QUADRATURE_POINT_x_AT_QUADRATURE_POINT_HXX_
/*!
//
// \file
//
//
// Created by Sebastien Gilles <sebastien.gilles@inria.fr> on the Mon, 15 Jun 2015 11:59:41 +0200
// Copyright (c) Inria. All rights reserved.
//
// \ingroup ParametersGroup
// \addtogroup ParametersGroup
// \{
*/
#ifndef MOREFEM_x_PARAMETERS_x_POLICY_x_AT_QUADRATURE_POINT_x_AT_QUADRATURE_POINT_HXX_
# define MOREFEM_x_PARAMETERS_x_POLICY_x_AT_QUADRATURE_POINT_x_AT_QUADRATURE_POINT_HXX_
namespace MoReFEM
{
namespace ParameterNS
{
namespace Policy
{
template<ParameterNS::Type TypeT>
AtQuadraturePoint<TypeT>
::AtQuadraturePoint(const std::string& name,
const Domain& domain,
const QuadratureRulePerTopology& quadrature_rule_per_topology,
storage_value_type initial_value,
const TimeManager& time_manager)
: domain_(domain),
initial_value_(initial_value),
quadrature_rule_per_topology_(quadrature_rule_per_topology),
time_manager_(time_manager)
{
static_cast<void>(name);
auto& storage = GetNonCstStorage();
storage.max_load_factor(Utilities::DefaultMaxLoadFactor());
assert(domain.template IsConstraintOn<DomainNS::Criterion::mesh>()
&& "It is an assert because this should have been checked earlier in ParameterInstance "
"constructor. At that place it is an exception, as a mere user might provide an invalid Domain.");
const auto& mesh = domain.GetMesh();
storage.reserve(static_cast<std::size_t>(mesh.NgeometricElt()));
const auto& ref_geom_elt_list = mesh.BagOfEltType();
for (const auto& ref_geom_elt_ptr : ref_geom_elt_list)
{
assert(!(!ref_geom_elt_ptr));
const auto& ref_geom_elt = *ref_geom_elt_ptr;
const auto& quadrature_rule = GetQuadratureRule(ref_geom_elt.GetTopologyIdentifier());
const auto Nquadrature_point = static_cast<std::size_t>(quadrature_rule.NquadraturePoint());
decltype(auto) iterator_range = mesh.GetSubsetGeometricEltList(ref_geom_elt);
// Iterate over all geometric elements that share the ref_geom_elt.
for (auto it_geom_elt = iterator_range.first; it_geom_elt != iterator_range.second; ++it_geom_elt)
{
const auto& geom_elt_ptr = *it_geom_elt;
assert(!(!geom_elt_ptr));
const auto geom_elt_index = geom_elt_ptr->GetIndex();
using value_type = Internal::ParameterNS::AtQuadraturePointNS::ValueHolder<storage_value_type>;
value_type default_value(initial_value, time_manager.NtimeModified());
auto check = storage.insert({geom_elt_index, std::vector<value_type>(Nquadrature_point, default_value)});
assert(check.second && "A geometric element should be entered only once.");
static_cast<void>(check);
}
}
}
template<ParameterNS::Type TypeT>
inline typename AtQuadraturePoint<TypeT>::return_type
AtQuadraturePoint<TypeT>::GetValueFromPolicy(const QuadraturePoint& quad_pt,
const GeometricElt& geom_elt) const
{
return FindValue(quad_pt, geom_elt).value;
}
template<ParameterNS::Type TypeT>
inline typename AtQuadraturePoint<TypeT>::return_type
AtQuadraturePoint<TypeT>::GetAnyValueFromPolicy() const
{
decltype(auto) storage = GetStorage();
assert(!storage.empty());
decltype(auto) geom_elt_content = storage.cbegin()->second;
assert(!geom_elt_content.empty());
return geom_elt_content.back().value;
}
template<ParameterNS::Type TypeT>
const typename AtQuadraturePoint<TypeT>::value_holder_type&
AtQuadraturePoint<TypeT>::FindValue(const QuadraturePoint& quad_pt,
const GeometricElt& geom_elt) const
{
auto it = storage_.find(geom_elt.GetIndex());
assert(it != storage_.cend());
const auto& values_in_geom_elt = it->second;
const auto quad_pt_index = static_cast<std::size_t>(quad_pt.GetIndex());
assert(quad_pt_index < values_in_geom_elt.size());
return values_in_geom_elt[quad_pt_index];
}
template<ParameterNS::Type TypeT>
inline typename AtQuadraturePoint<TypeT>::value_holder_type&
AtQuadraturePoint<TypeT>::FindNonCstValue(const QuadraturePoint& quad_pt,
const GeometricElt& geom_elt)
{
return const_cast<typename AtQuadraturePoint<TypeT>::value_holder_type&>(FindValue(quad_pt, geom_elt));
}
template<ParameterNS::Type TypeT>
template<class UpdateFunctorT>
typename AtQuadraturePoint<TypeT>::return_type
AtQuadraturePoint<TypeT>::UpdateAndGetValue(const QuadraturePoint& quad_pt,
const GeometricElt& geom_elt,
const UpdateFunctorT& update_functor)