cell_polyhedron.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
#include "libmesh/cell_polyhedron.h"
 
// Local includes
#include "libmesh/face_polygon.h"
#include "libmesh/enum_elem_quality.h"
#include "libmesh/hashword.h"
 
// C++ includes
#include <array>
#include <unordered_map>
#include <unordered_set>
 
 
namespace libMesh
{
 
// ------------------------------------------------------------
// Polyhedron class static member initializations
const int Polyhedron::num_children;
 
// ------------------------------------------------------------
// Polyhedron class member functions
 
 
Polyhedron::Polyhedron (const std::vector<std::shared_ptr<Polygon>> & sides,
                        Elem * p) :
    Cell(/* unused here */ 0, sides.size(), p, nullptr, nullptr),
    _elemlinks_data(sides.size()+2), // neighbors + parent + interior_parent
    _nodelinks_data(0), // We'll have to resize *this* later too!
    _sidelinks_data(sides.size())
  {
    // Set our sides, and while we're at it figure out our node maps
    // and side normal directions and edge lookup table, and count our
    // sides' nodes.  If we have internal nodes too then the subclass
    // will append those afterward.
    unsigned int nn = 0;
    std::unordered_map<Node *, unsigned int> local_node_number;
    std::unordered_set<std::pair<const Node *, const Node *>,
                       libMesh::hash> edges_seen;
    std::unique_ptr<const Elem> edge;
    for (unsigned int s : index_range(sides))
      {
        libmesh_assert(sides[s].get());
        auto & side_tuple = _sidelinks_data[s];
        std::get<0>(side_tuple) = sides[s];
 
        Polygon & side = *sides[s]; // not const, for writeable nodes
        for (auto n : make_range(side.n_nodes()))
          {
            Node * node = side.node_ptr(n);
            if (auto it = local_node_number.find(node);
                it != local_node_number.end())
              {
                std::get<2>(side_tuple).push_back(it->second);
              }
            else
              {
                std::get<2>(side_tuple).push_back(nn);
                local_node_number[node] = nn++;
                _nodelinks_data.push_back(node);
              }
          }
 
        for (unsigned int e : make_range(side.n_edges()))
          {
            side.build_edge_ptr(edge, e);
            auto edge_vertices = std::make_pair(edge->node_ptr(0), edge->node_ptr(1));
            if (edge_vertices.first > edge_vertices.second)
              std::swap(edge_vertices.first, edge_vertices.second);
 
            if (!edges_seen.count(edge_vertices))
              {
                edges_seen.insert(edge_vertices);
                _edge_lookup.emplace_back(s, e);
              }
          }
      }
 
    // Plan room for an extra node so we don't invalidate this->_nodes when we add to
    // nodelinks_data in the derived class
    _nodelinks_data.reserve(_nodelinks_data.size() + 1);
 
    // Do the manual initialization that Elem::Elem and Cell::Cell
    // couldn't, now that we've resized both our vectors.  No need to
    // manually set nullptr, though, since std::vector does that.
    this->_elemlinks = _elemlinks_data.data();
    this->_nodes = _nodelinks_data.data();
    this->_elemlinks[0] = p;
 
    libmesh_assert_equal_to(nn, this->n_nodes());
 
    // Figure out the orientation of our sides, now that we've got our
    // nodes organized enough to find our center.  The algorithm below
    // only works for convex polyhedra, but that's all we're
    // supporting for now.
    Point center;
    for (auto n : make_range(nn))
      center.add (this->point(n));
    center /= static_cast<Real>(nn);
 
    for (unsigned int s : index_range(sides))
      {
        const Polygon & side = *sides[s];
        const Point x_i = side.point(0);
        const Point n_i =
          (side.point(1) - side.point(0)).cross
          (side.point(0) - side.point(side.n_sides()-1)).unit();
 
        bool & inward_normal = std::get<1>(_sidelinks_data[s]);
        inward_normal = (n_i * (center - x_i) > TOLERANCE);
      }
 
    // We're betting a lot on "our polyhedra are all convex", so let's
    // check that if we have time.
#ifdef DEBUG
    for (unsigned int s : index_range(sides))
      {
        const Polygon & side = *sides[s];
        const Point x_i = side.point(0);
        const bool inward_normal = std::get<1>(this->_sidelinks_data[s]);
 
        const Point n_i =
          (side.point(1) - side.point(0)).cross
          (side.point(0) - side.point(side.n_sides()-1)).unit() *
          (inward_normal ? -1 : 1);
 
        for (const Point & node : this->node_ref_range())
          {
            const Point d_n = node - x_i;
            if (d_n * n_i > TOLERANCE * d_n.norm())
              libmesh_not_implemented_msg
                ("Cannot create a non-convex polyhedron");
          }
      }
#endif
 
    // Is this likely to ever be used?  We may do refinement with
    // polyhedra but it's probably not going to have a hierarchy...
    if (p)
      {
        this->subdomain_id() = p->subdomain_id();
        this->processor_id() = p->processor_id();
        _map_type = p->mapping_type();
        _map_data = p->mapping_data();
 
#ifdef LIBMESH_ENABLE_AMR
        this->set_p_level(p->p_level());
#endif
      }
 
    // Make sure the interior parent isn't undefined
    this->set_interior_parent(nullptr);
  }
 
 
 
Point Polyhedron::master_point (const unsigned int i) const
{
  return this->point(i);
}
 
 
 
bool Polyhedron::convex()
{
  for (unsigned int s : make_range(this->n_sides()))
    {
      const Polygon & side = *std::get<0>(this->_sidelinks_data[s]);
      const Point x_i = side.point(0);
      const bool inward_normal = std::get<1>(this->_sidelinks_data[s]);
 
      const Point n_i =
        (side.point(1) - side.point(0)).cross
        (side.point(0) - side.point(side.n_sides()-1)).unit() *
        (inward_normal ? -1 : 1);
 
      for (const Point & node : this->node_ref_range())
        {
          const Point d_n = node - x_i;
          if (d_n * n_i > TOLERANCE * d_n.norm())
            return false;
        }
    }
  return true;
}
 
 
 
bool Polyhedron::on_reference_element(const Point & p,
                                      const Real eps) const
{
  const unsigned int ns = this->n_sides();
 
  // Check that the point is on the same side of all the faces by
  // testing whether:
  //
  // n_i.(p - x_i) <= 0
  //
  // for each i, where:
  //   n_i is the outward normal of face i,
  //   x_i is a point on face i.
 
  for (auto i : make_range(ns))
    {
      const Polygon & face = *std::get<0>(this->_sidelinks_data[i]);
      const bool inward_normal = std::get<1>(this->_sidelinks_data[i]);
 
      const Point x_i = face.point(0);
 
      const Point n_i =
        (face.point(1) - face.point(0)).cross
        (face.point(0) - face.point(face.n_sides()-1)).unit() *
        (inward_normal ? -1 : 1);
 
  // This only works for polyhedra with flat sides.
#ifdef DEBUG
      for (auto j : make_range(face.n_sides()-1))
        {
          const Point x_j = face.point(j+1);
          const Point d_j = x_j - x_i;
          if (std::abs(d_j * n_i) > eps * d_j.norm())
            libmesh_not_implemented_msg
              ("Polyhedra with non-flat sides are not fully supported.");
        }
#endif
 
      if (n_i * (p - x_i) > eps)
        return false;
    }
 
  return true;
}
 
 
 
dof_id_type Polyhedron::key (const unsigned int s) const
{
  libmesh_assert_less (s, this->n_sides());
 
  const Polygon & face = *std::get<0>(this->_sidelinks_data[s]);
 
  return face.key();
}
 
 
 
dof_id_type Polyhedron::low_order_key (const unsigned int s) const
{
  libmesh_assert_less (s, this->n_sides());
 
  const Polygon & face = *std::get<0>(this->_sidelinks_data[s]);
 
  const unsigned int nv = face.n_vertices();
  std::vector<dof_id_type> vertex_ids(nv);
  for (unsigned int v : make_range(nv))
    vertex_ids[v] = face.node_id(v);
 
  return Utility::hashword(vertex_ids);
}
 
 
 
unsigned int Polyhedron::local_side_node(unsigned int side,
                                         unsigned int side_node) const
{
  libmesh_assert_less (side, this->n_sides());
 
  const std::vector<unsigned int> & node_map =
    std::get<2>(this->_sidelinks_data[side]);
  libmesh_assert_less (side_node, node_map.size());
 
  return node_map[side_node];
}
 
 
 
unsigned int Polyhedron::local_edge_node(unsigned int edge,
                                         unsigned int edge_node) const
{
  libmesh_assert_less (edge, this->n_edges());
  libmesh_assert_less (edge, _edge_lookup.size());
 
  auto [side, edge_of_side] = _edge_lookup[edge];
 
  const Polygon & face = *std::get<0>(this->_sidelinks_data[side]);
 
  const std::vector<unsigned int> & node_map =
    std::get<2>(this->_sidelinks_data[side]);
 
  return node_map[face.local_edge_node(edge_of_side, edge_node)];
}
 
 
 
dof_id_type Polyhedron::key () const
{
  std::vector<dof_id_type> node_ids;
  for (const auto & n : this->node_ref_range())
    node_ids.push_back(n.id());
 
  return Utility::hashword(node_ids);
}
 
 
 
std::unique_ptr<Elem> Polyhedron::side_ptr (const unsigned int i)
{
  return const_cast<const Polyhedron *>(this)->side_ptr(i);
}
 
 
 
std::unique_ptr<Elem> Polyhedron::side_ptr (const unsigned int i) const
{
  libmesh_assert_less (i, this->n_sides());
 
  Polygon & face = *std::get<0>(this->_sidelinks_data[i]);
  std::unique_ptr<Elem> face_copy = face.disconnected_clone();
  for (auto n : face.node_index_range())
    face_copy->set_node(n, face.node_ptr(n));
 
  return face_copy;
}
 
 
 
void Polyhedron::side_ptr (std::unique_ptr<Elem> & side,
                           const unsigned int i)
{
  libmesh_assert_less (i, this->n_sides());
 
  // Polyhedra are irregular enough that we're not even going to try
  // and bother optimizing heap access here.
  side = this->side_ptr(i);
}
 
 
 
std::unique_ptr<Elem> Polyhedron::build_side_ptr (const unsigned int i)
{
  auto returnval = this->side_ptr(i);
  returnval->set_interior_parent(this);
  returnval->inherit_data_from(*this);
  return returnval;
}
 
 
 
void Polyhedron::build_side_ptr (std::unique_ptr<Elem> & side,
                                 const unsigned int i)
{
  this->side_ptr(side, i);
  side->set_interior_parent(this);
  side->inherit_data_from(*this);
}
 
 
 
std::unique_ptr<Elem> Polyhedron::build_edge_ptr (const unsigned int i)
{
  auto [s, se] = _edge_lookup[i];
  Polygon & face = *std::get<0>(_sidelinks_data[s]);
  return face.build_edge_ptr(se);
}
 
 
 
void Polyhedron::build_edge_ptr (std::unique_ptr<Elem> & elem,
                                 const unsigned int i)
{
  auto [s, se] = _edge_lookup[i];
  Polygon & face = *std::get<0>(_sidelinks_data[s]);
  face.build_edge_ptr(elem, se);
}
 
 
 
bool Polyhedron::is_child_on_side(const unsigned int /*c*/,
                                  const unsigned int /*s*/) const
{
  libmesh_not_implemented();
  return false;
}
 
 
 
unsigned int Polyhedron::opposite_side(const unsigned int /* side_in */) const
{
  // This is too ambiguous in general.
  libmesh_not_implemented();
  return libMesh::invalid_uint;
}
 
 
 
unsigned int Polyhedron::opposite_node(const unsigned int /* n */,
                                       const unsigned int /* s */) const
{
  // This is too ambiguous in general.
  libmesh_not_implemented();
  return libMesh::invalid_uint;
}
 
 
 
bool Polyhedron::is_flipped() const
{
  if (this->_triangulation.empty())
    return false;
 
  auto & tet = this->_triangulation[0];
 
  const Point v01 = this->point(tet[1]) - this->point(tet[0]);
  const Point v02 = this->point(tet[2]) - this->point(tet[0]);
  const Point v03 = this->point(tet[3]) - this->point(tet[0]);
 
  return (triple_product(v01, v02, v03) < 0);
}
 
 
std::vector<unsigned int>
Polyhedron::edges_adjacent_to_node(const unsigned int n) const
{
  libmesh_assert_less(n, this->n_nodes());
 
  // For mid-edge or mid-face nodes, the subclass had better have
  // overridden this.
  libmesh_assert_less(n, this->n_vertices());
 
  libmesh_assert_equal_to(this->n_edges(), _edge_lookup.size());
 
  std::vector<unsigned int> adjacent_edges;
 
  for (const auto edge : index_range(_edge_lookup))
    {
      const auto [side, side_edge] = _edge_lookup[edge];
      const Polygon & face = *std::get<0>(_sidelinks_data[side]);
      const std::vector<unsigned int> & node_map =
        std::get<2>(_sidelinks_data[side]);
 
      const unsigned int fnv = face.n_vertices();
      libmesh_assert_equal_to(fnv, face.n_edges());
      libmesh_assert_equal_to(fnv, face.n_sides());
      libmesh_assert_less_equal(fnv, node_map.size());
      libmesh_assert_less(side_edge, fnv);
 
      const unsigned int vn = node_map[side_edge];
      const unsigned int vnp = node_map[(side_edge + 1) % fnv];
      libmesh_assert_less(vn, this->n_vertices());
      libmesh_assert_less(vnp, this->n_vertices());
      if (vn == n || vnp == n)
        adjacent_edges.push_back(edge);
    }
 
  return adjacent_edges;
}
 
 
std::pair<Real, Real> Polyhedron::qual_bounds (const ElemQuality q) const
{
  std::pair<Real, Real> bounds;
 
  switch (q)
    {
    case EDGE_LENGTH_RATIO:
      bounds.first  = 1.;
      bounds.second = 4.;
      break;
 
    case MIN_ANGLE:
      bounds.first  = 30.;
      bounds.second = 180.;
      break;
 
    case MAX_ANGLE:
      bounds.first  = 60.;
      bounds.second = 180.;
      break;
 
    // The generic dihedral metrics measure the unoriented angle
    // between adjacent face planes, so their values are in [0, 90].
    // Use the same recommended lower bounds as MIN,MAX_ANGLE.
    case MIN_DIHEDRAL_ANGLE:
      bounds.first  = 30.;
      bounds.second = 90.;
      break;
 
    case MAX_DIHEDRAL_ANGLE:
      bounds.first  = 60.;
      bounds.second = 90.;
      break;
 
    case JACOBIAN:
    case SCALED_JACOBIAN:
      bounds.first  = 0.5;
      bounds.second = 1.;
      break;
 
    default:
      libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
      bounds.first  = -1;
      bounds.second = -1;
    }
 
  return bounds;
}
 
 
 
std::vector<std::shared_ptr<Polygon>>
Polyhedron::side_clones() const
{
  const auto ns = this->n_sides();
 
  libmesh_assert_equal_to(ns, _sidelinks_data.size());
 
  std::vector<std::shared_ptr<Polygon>> cloned_sides(ns);
 
  for (auto i : make_range(ns))
    {
      const Polygon & face = *std::get<0>(this->_sidelinks_data[i]);
 
      Elem * clone = face.disconnected_clone().release();
      Polygon * polygon_clone = cast_ptr<Polygon *>(clone);
      cloned_sides[i] = std::shared_ptr<Polygon>(polygon_clone);
 
      // We can't actually use a *disconnected* clone to reconstruct
      // links between sides, so we'll temporarily give the clone our
      // own nodes; user code that typically replaces the usual
      // nullptr with permanent nodes will then instead place our
      // nodes with permanent nodes.
      for (auto n : make_range(face.n_nodes()))
        cloned_sides[i]->set_node
          (n, const_cast<Node *>(face.node_ptr(n)));
    }
 
  return cloned_sides;
}
 
 
 
bool Polyhedron::side_has_edge_nodes(unsigned int s,
                                     unsigned int min_node,
                                     unsigned int max_node) const
{
  const Polygon & face = *std::get<0>(_sidelinks_data[s]);
  const std::vector<unsigned int> & node_map =
    std::get<2>(this->_sidelinks_data[s]);
 
  for (unsigned int e : make_range(face.n_sides()))
    {
      std::vector<unsigned int> nodes_on_edge =
        face.nodes_on_side(e);
      libmesh_assert_equal_to(nodes_on_edge.size(), 2);
      nodes_on_edge[0] = node_map[nodes_on_edge[0]];
      nodes_on_edge[1] = node_map[nodes_on_edge[1]];
      if ((nodes_on_edge[0] == min_node) &&
          (nodes_on_edge[1] == max_node))
        return true;
      if ((nodes_on_edge[1] == min_node) &&
          (nodes_on_edge[0] == max_node))
        return true;
    }
 
  return false;
}
 
 
 
std::vector<unsigned int>
Polyhedron::sides_on_edge(const unsigned int e) const
{
  std::vector<unsigned int> returnval(2);
  auto [s1, s1e] = _edge_lookup[e];
  returnval[0] = s1;
 
  const Polygon & face1 = *std::get<0>(_sidelinks_data[s1]);
  const std::vector<unsigned int> & node_map =
    std::get<2>(this->_sidelinks_data[s1]);
 
  std::vector<unsigned int> nodes_on_edge =
    face1.nodes_on_side(s1e);
  libmesh_assert_equal_to(nodes_on_edge.size(), 2);
  nodes_on_edge[0] = node_map[nodes_on_edge[0]];
  nodes_on_edge[1] = node_map[nodes_on_edge[1]];
 
  if (nodes_on_edge[0] > nodes_on_edge[1])
    std::swap(nodes_on_edge[0], nodes_on_edge[1]);
 
  for (unsigned int s2 : make_range(this->n_sides()))
    {
      if (s2 == s1)
        continue;
 
      if (this->side_has_edge_nodes(s2, nodes_on_edge[0],
                                    nodes_on_edge[1]))
        {
          returnval[1] = s2;
          return returnval;
        }
    }
 
  libmesh_error();
 
  return returnval;
}
 
 
 
bool Polyhedron::is_edge_on_side(const unsigned int e,
                                 const unsigned int s) const
{
  auto [s1, s1e] = _edge_lookup[e];
 
  // Did we get lucky with our cache?
  if (s1 == s)
    return true;
 
  const Polygon & face1 = *std::get<0>(_sidelinks_data[s1]);
  const std::vector<unsigned int> & node_map =
    std::get<2>(this->_sidelinks_data[s1]);
  std::vector<unsigned int> nodes_on_edge1 =
    face1.nodes_on_side(s1e);
  libmesh_assert_equal_to(nodes_on_edge1.size(), 2);
 
  nodes_on_edge1[0] = node_map[nodes_on_edge1[0]];
  nodes_on_edge1[1] = node_map[nodes_on_edge1[1]];
  if (nodes_on_edge1[0] > nodes_on_edge1[1])
    std::swap(nodes_on_edge1[0], nodes_on_edge1[1]);
 
  return this->side_has_edge_nodes(s,
                                   nodes_on_edge1[0],
                                   nodes_on_edge1[1]);
}
 
 
 
std::array<Point, 4> Polyhedron::master_subelement (unsigned int i) const
{
  libmesh_assert_less(i, this->_triangulation.size());
 
  const auto & tet = this->_triangulation[i];
 
  return { this->master_point(tet[0]),
           this->master_point(tet[1]),
           this->master_point(tet[2]),
           this->master_point(tet[3]) };
}
 
 
 
std::tuple<unsigned int, Real, Real, Real>
Polyhedron::subelement_coordinates (const Point & p, Real tol) const
{
  std::tuple<unsigned int, Real, Real, Real> returnval =
    {libMesh::invalid_uint, -1, -1, -1};
 
  Real best_bad_coord = -1;
 
  for (auto s : make_range(this->n_subelements()))
    {
      const std::array<Point, 4> subtet =
        this->master_subelement(s);
 
      // Find barycentric coordinates in subelem
      const Point v0 = p - subtet[0];
      // const Point v1 = p - subtet[1];
 
      const Point v01 = subtet[1] - subtet[0];
      const Point v02 = subtet[2] - subtet[0];
      const Point v03 = subtet[3] - subtet[0];
 
      // const Point v12 = subtet[2] - subtet[1];
      // const Point v13 = subtet[3] - subtet[1];
 
      // const Real tp0 = triple_product(v1, v13, v12);
      const Real tp1 = triple_product(v0, v02, v03);
      const Real tp2 = triple_product(v0, v03, v01);
      const Real tp3 = triple_product(v0, v01, v02);
 
      const Real six_vol = triple_product(v01, v02, v03);
 
      const Real xi   = tp1 / six_vol;
      const Real eta  = tp2 / six_vol;
      const Real zeta = tp3 / six_vol;
 
      if (xi>=0 && eta>=0 && zeta>=0 && xi+eta+zeta<=1)
        return { s, xi, eta, zeta };
 
      const Real my_best_bad_coord =
        std::min(std::min(std::min(xi, eta), zeta), 1-xi-eta-zeta);
 
      if (my_best_bad_coord > best_bad_coord)
        {
          best_bad_coord = my_best_bad_coord;
          returnval = { s, xi, eta, zeta };
        }
    }
 
  if (best_bad_coord > -tol)
    return returnval;
 
  return {libMesh::invalid_uint, -1, -1, -1};
}
 
 
} // namespace libMesh