Commit b5b4d514 authored by GILLES Sebastien's avatar GILLES Sebastien

#1520 Refactor the test about GlobalCoordsQuadPts to make it more readable (it...

#1520 Refactor the test about GlobalCoordsQuadPts to make it more readable (it remains a tad clumsy but this will do). I have also modified the way the ParamAtQuadraturePoint Write() is working, now using the PrintContainer facility (and at the moment expending the precision - maybe that is not pertinent).
parent 697e70bb
......@@ -233,15 +233,13 @@ namespace MoReFEM
<< 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;
out << std::setprecision(16) << 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.precision(16);
Utilities::PrintContainer<>::Do(value_per_quad_pt[i].value, out, " ", "", "");
}
out << std::endl;
......
......@@ -131,12 +131,14 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
{
void FillExpectedGlobalCoords(const QuadratureRule& quadrature_rule,
std::vector<double>& expected_global_coords)
/*!
* \brief Compute analytically the \a GlobalCoords .
*/
std::vector<double> FillExpectedGlobalCoords(const QuadratureRule& quadrature_rule)
{
decltype(auto) quadrature_point_list = quadrature_rule.GetQuadraturePointList();
const auto topology_identifier = quadrature_rule.GetTopologyIdentifier();
std::vector<double> ret;
for (const auto& it_quad_pt : quadrature_point_list)
{
assert(!(!it_quad_pt));
......@@ -165,15 +167,57 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
}
}
expected_global_coords.push_back(expected_value);
ret.push_back(expected_value);
}
// Set to 0 global coords of the dimensions that are not relevant as local coords of the quadrature
// have the same number of components as the dimension of the finite element considered.
const auto Ndimension = it_quad_pt->GetDimension();
for (auto i = Ndimension; i < 3; ++i)
expected_global_coords.push_back(0.);
ret.push_back(0.);
}
return ret;
}
/*!
* \brief Convert the obtained global coords into a vector of double which is filled with same ordering as the expected ones.
*
*/
std::vector<double>
ConvertObtainedGloballCoords(const FilesystemNS::Directory& result_directory_path,
std::string_view quadrature_order,
const ParameterAtQuadraturePoint<ParameterNS::Type::vector>& obtained_global_coords)
{
const std::string output_file = result_directory_path.AddFile(quadrature_order);
obtained_global_coords.Write(output_file);
std::ifstream file_stream;
FilesystemNS::File::Read(file_stream, output_file, __FILE__, __LINE__);
std::string line;
// skip first two lines
for (int i = 0; i < 2; ++i)
getline(file_stream, line);
std::vector<double> ret;
while (getline(file_stream, line))
{
const auto item_list = Utilities::String::Split(line, ";");
BOOST_CHECK_EQUAL(item_list.size(), 4ul);
const auto coords_list = std::string(item_list.back()); // istringstream can't act upon a
// std::string_view hence this copy.
std::istringstream iconv(coords_list);
double value;
while (iconv >> value)
ret.push_back(value);
}
return ret;
}
......@@ -186,71 +230,44 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
.GetDomain(EnumUnderlyingType(DomainIndex::domain), __FILE__, __LINE__);
decltype(auto) mesh = domain.GetMesh();
decltype(auto) ref_geom_elt_list = mesh.BagOfEltType();
std::vector<double> expected_global_coords;
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 topology_identifier = ref_geom_elt.GetTopologyIdentifier();
const auto& quadrature_rule = quadrature_rule_per_topology.GetRule(topology_identifier);
BOOST_CHECK_EQUAL(domain.GetDimensionList().size(), 1ul);
FillExpectedGlobalCoords(quadrature_rule, expected_global_coords);
// For this test, we intend to work on a simple mesh with a single element and we are doing the
// quadrature points computation only on the highest dimension.
decltype(auto) ref_geom_elt_list = mesh.BagOfEltType(domain.GetDimensionList().back());
BOOST_CHECK_EQUAL(ref_geom_elt_list.size(), 1ul);
FilesystemNS::Directory output_directory(result_directory_path,
ref_geom_elt.GetTopologyName(),
__FILE__, __LINE__,
FilesystemNS::behaviour::create_if_necessary);
BOOST_CHECK(ref_geom_elt_list.back() != nullptr);
const std::string output_file = output_directory.AddFile(quadrature_order);
obtained_global_coords.Write(output_file);
const auto& ref_geom_elt = *(ref_geom_elt_list.back());
BOOST_CHECK(domain.DoRefGeomEltMatchCriteria(ref_geom_elt));
// Read from output and compare results.
std::ifstream file_stream;
FilesystemNS::File::Read(file_stream, output_file, __FILE__, __LINE__);
std::string line;
const auto topology_identifier = ref_geom_elt.GetTopologyIdentifier();
const auto& quadrature_rule = quadrature_rule_per_topology.GetRule(topology_identifier);
// skip first two lines
for (int i = 0; i < 2; ++i)
getline(file_stream, line);
const auto expected_global_coords = FillExpectedGlobalCoords(quadrature_rule);
std::vector<double> obtained_values;
const auto obtained_global_coords_as_vector =
ConvertObtainedGloballCoords(result_directory_path,
quadrature_order,
obtained_global_coords);
while (getline(file_stream, line))
{
const auto item_list = Utilities::String::Split(line, ";");
BOOST_CHECK_EQUAL(item_list.size(), 4ul);
const auto end_obtained = obtained_global_coords_as_vector.cend();
const auto coords_list = std::string(item_list.back()); // istringstream can't act upon a
// std::string_view hence this copy.
std::istringstream iconv(coords_list);
double value;
while (iconv >> value)
obtained_values.push_back(value);
}
BOOST_CHECK_EQUAL(obtained_global_coords_as_vector.size(), expected_global_coords.size());
constexpr const auto precision = 1.e-12;
const auto end_obtained = obtained_values.cend();
for (auto it_obtained = obtained_global_coords_as_vector.cbegin(), it_expected = expected_global_coords.cbegin();
it_obtained != end_obtained;
++it_obtained, ++it_expected)
{
const auto obtained_value = *it_obtained;
const auto expected_value = *it_expected;
BOOST_CHECK_EQUAL(obtained_values.size(), expected_global_coords.size());
constexpr const auto ensight_precision = 1.e-6;
for (auto it_obtained = obtained_values.cbegin(), it_expected = expected_global_coords.cbegin();
it_obtained != end_obtained;
++it_obtained, ++it_expected)
{
const auto obtained_value = *it_obtained;
const auto expected_value = *it_expected;
// Takes into account the precision at which we wrote the quadrature point coordinates.
BOOST_CHECK(NumericNS::AreEqual(expected_value, obtained_value, ensight_precision));
}
// Takes into account the precision at which we wrote the quadrature point coordinates.
BOOST_CHECK_CLOSE(expected_value, obtained_value, precision);
}
}
}
......
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