Skip to content
Snippets Groups Projects
Commit f3150c5b authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Improve MeshNodeBoundary and MeshFlatNodeBoundary 

- improve error messages
- fixes #8 (parallel/sequential reproducibility and straight boundaries)
parent 82cb1358
No related branches found
No related tags found
1 merge request!96Add random number engine encapsulation.
Show whitespace changes
Inline Side-by-side
src/mesh/MeshNodeBoundary.hpp
+ 140
− 108
View file @ f3150c5b
...@@ -16,6 +16,7 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic ...@@ -16,6 +16,7 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic
{ {
protected: protected:
Array<const NodeId> m_node_list; Array<const NodeId> m_node_list;
std::string m_boundary_name;
public: public:
MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default; MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
...@@ -29,6 +30,7 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic ...@@ -29,6 +30,7 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic
template <typename MeshType> template <typename MeshType>
MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list) MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
: m_boundary_name(ref_face_list.refId().tagName())
{ {
static_assert(Dimension == MeshType::Dimension); static_assert(Dimension == MeshType::Dimension);
const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
...@@ -38,7 +40,10 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic ...@@ -38,7 +40,10 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic
face_list.size(), PUGS_LAMBDA(int l) { face_list.size(), PUGS_LAMBDA(int l) {
const auto& face_cells = face_to_cell_matrix[face_list[l]]; const auto& face_cells = face_to_cell_matrix[face_list[l]];
if (face_cells.size() > 1) { if (face_cells.size() > 1) {
throw NormalError("internal faces cannot be used to define mesh boundaries"); std::ostringstream ost;
ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
<< ": inner faces cannot be used to define mesh boundaries";
throw NormalError(ost.str());
} }
}); });
...@@ -66,7 +71,9 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic ...@@ -66,7 +71,9 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic
} }
template <typename MeshType> template <typename MeshType>
MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list) : m_node_list(ref_node_list.list()) MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list)
: m_node_list(ref_node_list.list()), m_boundary_name(ref_node_list.refId().tagName())
{ {
static_assert(Dimension == MeshType::Dimension); static_assert(Dimension == MeshType::Dimension);
} }
...@@ -92,10 +99,27 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension> // clazy:exclu ...@@ -92,10 +99,27 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension> // clazy:exclu
PUGS_INLINE Rd _getNormal(const MeshType& mesh); PUGS_INLINE Rd _getNormal(const MeshType& mesh);
template <typename MeshType> template <typename MeshType>
PUGS_INLINE void _checkBoundaryIsFlat(TinyVector<2, double> normal, PUGS_INLINE void
TinyVector<2, double> xmin, _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
TinyVector<2, double> xmax, const TinyVector<Dimension, double>& origin,
const MeshType& mesh) const; const double length,
const MeshType& mesh) const
{
static_assert(MeshType::Dimension == Dimension);
const NodeValue<const Rd>& xr = mesh.xr();
parallel_for(
this->m_node_list.size(), PUGS_LAMBDA(const size_t r) {
const Rd& x = xr[this->m_node_list[r]];
if (dot(x - origin, normal) > 1E-13 * length) {
std::ostringstream ost;
ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
<< ": boundary is not flat!";
throw NormalError(ost.str());
}
});
}
template <typename MeshType> template <typename MeshType>
PUGS_INLINE Rd _getOutgoingNormal(const MeshType& mesh); PUGS_INLINE Rd _getOutgoingNormal(const MeshType& mesh);
...@@ -130,31 +154,6 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension> // clazy:exclu ...@@ -130,31 +154,6 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension> // clazy:exclu
virtual ~MeshFlatNodeBoundary() = default; virtual ~MeshFlatNodeBoundary() = default;
}; };
template <>
template <typename MeshType>
void
MeshFlatNodeBoundary<2>::_checkBoundaryIsFlat(TinyVector<2, double> normal,
TinyVector<2, double> xmin,
TinyVector<2, double> xmax,
const MeshType& mesh) const
{
static_assert(MeshType::Dimension == 2);
using R2 = TinyVector<2, double>;
const R2 origin = 0.5 * (xmin + xmax);
const double length = l2Norm(xmax - xmin);
const NodeValue<const R2>& xr = mesh.xr();
parallel_for(
m_node_list.size(), PUGS_LAMBDA(size_t r) {
const R2& x = xr[m_node_list[r]];
if (dot(x - origin, normal) > 1E-13 * length) {
throw NormalError("this FlatBoundary is not flat!");
}
});
}
template <> template <>
template <typename MeshType> template <typename MeshType>
PUGS_INLINE TinyVector<1, double> PUGS_INLINE TinyVector<1, double>
...@@ -173,7 +172,10 @@ MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh) ...@@ -173,7 +172,10 @@ MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh)
}(); }();
if (number_of_bc_nodes != 1) { if (number_of_bc_nodes != 1) {
throw NormalError("Node boundaries in 1D require to have exactly 1 node"); std::ostringstream ost;
ost << "invalid boundary " << rang::fgB::yellow << m_boundary_name << rang::style::reset
<< ": node boundaries in 1D require to have exactly 1 node";
throw NormalError(ost.str());
} }
return R{1}; return R{1};
...@@ -190,46 +192,50 @@ MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh) ...@@ -190,46 +192,50 @@ MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh)
const NodeValue<const R2>& xr = mesh.xr(); const NodeValue<const R2>& xr = mesh.xr();
R2 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max()); R2 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
R2 xmax(-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()); R2 xmax(-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max());
for (size_t r = 0; r < m_node_list.size(); ++r) { auto update_xmin = [](const R2& x, R2& xmin) {
const R2& x = xr[m_node_list[r]];
if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) { if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
xmin = x; xmin = x;
} }
};
auto update_xmax = [](const R2& x, R2& xmax) {
if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) { if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
xmax = x; xmax = x;
} }
};
for (size_t r = 0; r < m_node_list.size(); ++r) {
const R2& x = xr[m_node_list[r]];
update_xmin(x, xmin);
update_xmax(x, xmax);
} }
if (parallel::size() > 1) {
Array<R2> xmin_array = parallel::allGather(xmin); Array<R2> xmin_array = parallel::allGather(xmin);
Array<R2> xmax_array = parallel::allGather(xmax); Array<R2> xmax_array = parallel::allGather(xmax);
for (size_t i = 0; i < xmin_array.size(); ++i) { for (size_t i = 0; i < xmin_array.size(); ++i) {
const R2& x = xmin_array[i]; update_xmin(xmin_array[i], xmin);
if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
xmin = x;
}
} }
for (size_t i = 0; i < xmax_array.size(); ++i) { for (size_t i = 0; i < xmax_array.size(); ++i) {
const R2& x = xmax_array[i]; update_xmax(xmax_array[i], xmax);
if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
xmax = x;
} }
} }
if (xmin == xmax) { if (xmin == xmax) {
std::stringstream os; std::ostringstream ost;
os << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal"; ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
throw NormalError(os.str()); << ": unable to compute normal";
throw NormalError(ost.str());
} }
R2 dx = xmax - xmin; R2 dx = xmax - xmin;
dx *= 1. / l2Norm(dx); dx *= 1. / l2Norm(dx);
R2 normal(-dx[1], dx[0]); const R2 normal(-dx[1], dx[0]);
this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh); this->_checkBoundaryIsFlat(normal, 0.5 * (xmin + xmax), l2Norm(xmax - xmin), mesh);
return normal; return normal;
} }
...@@ -242,6 +248,54 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh) ...@@ -242,6 +248,54 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
static_assert(MeshType::Dimension == 3); static_assert(MeshType::Dimension == 3);
using R3 = TinyVector<3, double>; using R3 = TinyVector<3, double>;
auto update_xmin = [](const R3& x, R3& xmin) {
// XMIN: X.xmin X.ymax X.zmax
if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] > xmin[1])) or
((x[0] == xmin[0]) and (x[1] == xmin[1]) and (x[2] > xmin[2]))) {
xmin = x;
}
};
auto update_xmax = [](const R3& x, R3& xmax) {
// XMAX: X.xmax X.ymin X.zmin
if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] < xmax[1])) or
((x[0] == xmax[0]) and (x[1] == xmax[1]) and (x[2] < xmax[2]))) {
xmax = x;
}
};
auto update_ymin = [](const R3& x, R3& ymin) {
// YMIN: X.ymin X.zmax X.xmin
if ((x[1] < ymin[1]) or ((x[1] == ymin[1]) and (x[2] > ymin[2])) or
((x[1] == ymin[1]) and (x[2] == ymin[2]) and (x[0] < ymin[0]))) {
ymin = x;
}
};
auto update_ymax = [](const R3& x, R3& ymax) {
// YMAX: X.ymax X.zmin X.xmax
if ((x[1] > ymax[1]) or ((x[1] == ymax[1]) and (x[2] < ymax[2])) or
((x[1] == ymax[1]) and (x[2] == ymax[1]) and (x[0] > ymax[0]))) {
ymax = x;
}
};
auto update_zmin = [](const R3& x, R3& zmin) {
// ZMIN: X.zmin X.xmin X.ymin
if ((x[2] < zmin[2]) or ((x[2] == zmin[2]) and (x[2] < zmin[2])) or
((x[1] == zmin[1]) and (x[2] == zmin[1]) and (x[0] < zmin[0]))) {
zmin = x;
}
};
auto update_zmax = [](const R3& x, R3& zmax) {
// ZMAX: X.zmax X.xmax X.ymax
if ((x[2] > zmax[2]) or ((x[2] == zmax[2]) and (x[2] > zmax[2])) or
((x[1] == zmax[1]) and (x[2] == zmax[1]) and (x[0] > zmax[0]))) {
zmax = x;
}
};
R3 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()); R3 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
R3 ymin = xmin; R3 ymin = xmin;
R3 zmin = xmin; R3 zmin = xmin;
...@@ -254,66 +308,39 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh) ...@@ -254,66 +308,39 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
for (size_t r = 0; r < m_node_list.size(); ++r) { for (size_t r = 0; r < m_node_list.size(); ++r) {
const R3& x = xr[m_node_list[r]]; const R3& x = xr[m_node_list[r]];
if (x[0] < xmin[0]) { update_xmin(x, xmin);
xmin = x; update_xmax(x, xmax);
} update_ymin(x, ymin);
if (x[1] < ymin[1]) { update_ymax(x, ymax);
ymin = x; update_zmin(x, zmin);
} update_zmax(x, zmax);
if (x[2] < zmin[2]) { }
zmin = x;
} if (parallel::size() > 0) {
if (x[0] > xmax[0]) { Array<const R3> xmin_array = parallel::allGather(xmin);
xmax = x; Array<const R3> xmax_array = parallel::allGather(xmax);
} Array<const R3> ymin_array = parallel::allGather(ymin);
if (x[1] > ymax[1]) { Array<const R3> ymax_array = parallel::allGather(ymax);
ymax = x; Array<const R3> zmin_array = parallel::allGather(zmin);
} Array<const R3> zmax_array = parallel::allGather(zmax);
if (x[2] > zmax[2]) {
zmax = x;
}
}
Array<R3> xmin_array = parallel::allGather(xmin);
Array<R3> xmax_array = parallel::allGather(xmax);
Array<R3> ymin_array = parallel::allGather(ymin);
Array<R3> ymax_array = parallel::allGather(ymax);
Array<R3> zmin_array = parallel::allGather(zmin);
Array<R3> zmax_array = parallel::allGather(zmax);
for (size_t i = 0; i < xmin_array.size(); ++i) { for (size_t i = 0; i < xmin_array.size(); ++i) {
const R3& x = xmin_array[i]; update_xmin(xmin_array[i], xmin);
if (x[0] < xmin[0]) {
xmin = x;
}
} }
for (size_t i = 0; i < ymin_array.size(); ++i) { for (size_t i = 0; i < ymin_array.size(); ++i) {
const R3& x = ymin_array[i]; update_ymin(ymin_array[i], ymin);
if (x[1] < ymin[1]) {
ymin = x;
}
} }
for (size_t i = 0; i < zmin_array.size(); ++i) { for (size_t i = 0; i < zmin_array.size(); ++i) {
const R3& x = zmin_array[i]; update_zmin(zmin_array[i], zmin);
if (x[2] < zmin[2]) {
zmin = x;
}
} }
for (size_t i = 0; i < xmax_array.size(); ++i) { for (size_t i = 0; i < xmax_array.size(); ++i) {
const R3& x = xmax_array[i]; update_xmax(xmax_array[i], xmax);
if (x[0] > xmax[0]) {
xmax = x;
}
} }
for (size_t i = 0; i < ymax_array.size(); ++i) { for (size_t i = 0; i < ymax_array.size(); ++i) {
const R3& x = ymax_array[i]; update_ymax(ymax_array[i], ymax);
if (x[1] > ymax[1]) {
ymax = x;
}
} }
for (size_t i = 0; i < zmax_array.size(); ++i) { for (size_t i = 0; i < zmax_array.size(); ++i) {
const R3& x = zmax_array[i]; update_zmax(zmax_array[i], zmax);
if (x[2] > zmax[2]) {
zmax = x;
} }
} }
...@@ -344,12 +371,17 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh) ...@@ -344,12 +371,17 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
} }
if (normal_l2 == 0) { if (normal_l2 == 0) {
throw NormalError("cannot to compute normal!"); std::ostringstream ost;
ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
<< ": unable to compute normal";
throw NormalError(ost.str());
} }
normal *= 1. / sqrt(normal_l2); const double length = sqrt(normal_l2);
normal *= 1. / length;
// this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh); this->_checkBoundaryIsFlat(normal, 1. / 6. * (xmin + xmax + ymin + ymax + zmin + zmax), length, mesh);
return normal; return normal;
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment