Skip to content
Snippets Groups Projects
Select Git revision
  • 1209cded8b473ee733a6140f5715483e6dfb4985
  • develop default protected
  • feature/variational-hydro
  • origin/stage/bouguettaia
  • feature/gmsh-reader
  • feature/reconstruction
  • save_clemence
  • feature/kinetic-schemes
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • feature/serraille
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • master protected
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

BuiltinFunctionEmbedderUtils.cpp

Blame
  • ConnectivityBuilderBase.cpp 16.92 KiB
    #include <mesh/ConnectivityBuilderBase.hpp>
    
    #include <mesh/Connectivity.hpp>
    #include <mesh/ConnectivityDescriptor.hpp>
    #include <mesh/ItemId.hpp>
    #include <mesh/Mesh.hpp>
    #include <utils/PugsAssert.hpp>
    #include <utils/PugsMacros.hpp>
    
    #include <map>
    #include <unordered_map>
    #include <vector>
    
    template <size_t Dimension>
    void
    ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
    {
      static_assert((Dimension == 2) or (Dimension == 3), "Invalid dimension to compute cell-face connectivities");
      using CellFaceInfo = std::tuple<CellId, unsigned short, bool>;
      using Face         = ConnectivityFace<Dimension>;
    
      const auto& node_number_vector = descriptor.node_number_vector;
      Array<unsigned short> cell_nb_faces(descriptor.cell_to_node_vector.size());
      std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
      for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
        const auto& cell_nodes = descriptor.cell_to_node_vector[j];
    
        if constexpr (Dimension == 2) {
          switch (descriptor.cell_type_vector[j]) {
          case CellType::Triangle: {
            cell_nb_faces[j] = 3;
            // face 0
            Face f0({cell_nodes[1], cell_nodes[2]}, node_number_vector);
            face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
    
            // face 1
            Face f1({cell_nodes[2], cell_nodes[0]}, node_number_vector);
            face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
    
            // face 2
            Face f2({cell_nodes[0], cell_nodes[1]}, node_number_vector);
            face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
            break;
          }
          case CellType::Quadrangle: {
            cell_nb_faces[j] = 4;
            // face 0
            Face f0({cell_nodes[0], cell_nodes[1]}, node_number_vector);
            face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
    
            // face 1
            Face f1({cell_nodes[1], cell_nodes[2]}, node_number_vector);
            face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
    
            // face 2
            Face f2({cell_nodes[2], cell_nodes[3]}, node_number_vector);
            face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
    
            // face 3
            Face f3({cell_nodes[3], cell_nodes[0]}, node_number_vector);
            face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
            break;
          }
          default: {
            std::ostringstream error_msg;
            error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 2";
            throw UnexpectedError(error_msg.str());
          }
          }
        } else if constexpr (Dimension == 3) {
          switch (descriptor.cell_type_vector[j]) {
          case CellType::Hexahedron: {
            // face 0
            Face f0({cell_nodes[3], cell_nodes[2], cell_nodes[1], cell_nodes[0]}, node_number_vector);
            face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
    
            // face 1
            Face f1({cell_nodes[4], cell_nodes[5], cell_nodes[6], cell_nodes[7]}, node_number_vector);
            face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
    
            // face 2
            Face f2({cell_nodes[0], cell_nodes[4], cell_nodes[7], cell_nodes[3]}, node_number_vector);
            face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
    
            // face 3
            Face f3({cell_nodes[1], cell_nodes[2], cell_nodes[6], cell_nodes[5]}, node_number_vector);
            face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
    
            // face 4
            Face f4({cell_nodes[0], cell_nodes[1], cell_nodes[5], cell_nodes[4]}, node_number_vector);
            face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
    
            // face 5
            Face f5({cell_nodes[3], cell_nodes[7], cell_nodes[6], cell_nodes[2]}, node_number_vector);
            face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
    
            cell_nb_faces[j] = 6;
            break;
          }
          case CellType::Tetrahedron: {
            cell_nb_faces[j] = 4;
            // face 0
            Face f0({cell_nodes[1], cell_nodes[2], cell_nodes[3]}, node_number_vector);
            face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
    
            // face 1
            Face f1({cell_nodes[0], cell_nodes[3], cell_nodes[2]}, node_number_vector);
            face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
    
            // face 2
            Face f2({cell_nodes[0], cell_nodes[1], cell_nodes[3]}, node_number_vector);
            face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
    
            // face 3
            Face f3({cell_nodes[0], cell_nodes[2], cell_nodes[1]}, node_number_vector);
            face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
            break;
          }
          case CellType::Pyramid: {
            cell_nb_faces[j] = cell_nodes.size();
            std::vector<unsigned int> base_nodes;
            std::copy_n(cell_nodes.begin(), cell_nodes.size() - 1, std::back_inserter(base_nodes));
    
            // base face
            {
              Face base_face(base_nodes, node_number_vector);
              face_cells_map[base_face].emplace_back(std::make_tuple(j, 0, base_face.reversed()));
            }
            // side faces
            const auto pyramid_vertex = cell_nodes[cell_nodes.size() - 1];
            for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
              Face side_face({base_nodes[(i_node + 1) % base_nodes.size()], base_nodes[i_node], pyramid_vertex},
                             node_number_vector);
              face_cells_map[side_face].emplace_back(std::make_tuple(j, i_node + 1, side_face.reversed()));
            }
            break;
          }
          case CellType::Diamond: {
            cell_nb_faces[j] = 2 * (cell_nodes.size() - 2);
            std::vector<unsigned int> base_nodes;
            std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
    
            {   // top faces
              const auto top_vertex = cell_nodes[cell_nodes.size() - 1];
              for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
                Face top_face({base_nodes[i_node], base_nodes[(i_node + 1) % base_nodes.size()], top_vertex},
                              node_number_vector);
                face_cells_map[top_face].emplace_back(std::make_tuple(j, i_node, top_face.reversed()));
              }
            }
    
            {   // bottom faces
              const auto bottom_vertex = cell_nodes[0];
              for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
                Face bottom_face({base_nodes[(i_node + 1) % base_nodes.size()], base_nodes[i_node], bottom_vertex},
                                 node_number_vector);
                face_cells_map[bottom_face].emplace_back(
                  std::make_tuple(j, i_node + base_nodes.size(), bottom_face.reversed()));
              }
            }
            break;
          }
          default: {
            std::ostringstream error_msg;
            error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
            throw UnexpectedError(error_msg.str());
          }
          }
        }
      }
    
      {
        descriptor.cell_to_face_vector.resize(descriptor.cell_to_node_vector.size());
        for (CellId j = 0; j < descriptor.cell_to_face_vector.size(); ++j) {
          descriptor.cell_to_face_vector[j].resize(cell_nb_faces[j]);
        }
        FaceId l = 0;
        for (const auto& face_cells_vector : face_cells_map) {
          const auto& cells_vector = face_cells_vector.second;
          for (unsigned short lj = 0; lj < cells_vector.size(); ++lj) {
            const auto& [cell_number, cell_local_face, reversed]         = cells_vector[lj];
            descriptor.cell_to_face_vector[cell_number][cell_local_face] = l;
          }
          ++l;
        }
      }
    
      {
        descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_to_node_vector.size());
        for (CellId j = 0; j < descriptor.cell_face_is_reversed_vector.size(); ++j) {
          descriptor.cell_face_is_reversed_vector[j] = Array<bool>(cell_nb_faces[j]);
        }
        for (const auto& face_cells_vector : face_cells_map) {
          const auto& cells_vector = face_cells_vector.second;
          for (unsigned short lj = 0; lj < cells_vector.size(); ++lj) {
            const auto& [cell_number, cell_local_face, reversed]                  = cells_vector[lj];
            descriptor.cell_face_is_reversed_vector[cell_number][cell_local_face] = reversed;
          }
        }
      }
    
      {
        descriptor.face_to_node_vector.resize(face_cells_map.size());
        int l = 0;
        for (const auto& face_info : face_cells_map) {
          const Face& face                  = face_info.first;
          descriptor.face_to_node_vector[l] = face.nodeIdList();
          ++l;
        }
      }
    
      {
        // Face numbers may change if numbers are provided in the file
        descriptor.face_number_vector.resize(face_cells_map.size());
        for (size_t l = 0; l < face_cells_map.size(); ++l) {
          descriptor.face_number_vector[l] = l;
        }
      }
    }
    
    template <size_t Dimension>
    void
    ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor)
    {
      static_assert(Dimension == 3, "Invalid dimension to compute face-edge connectivities");
      using FaceEdgeInfo = std::tuple<FaceId, unsigned short, bool>;
      using Edge         = ConnectivityFace<2>;
    
      const auto& node_number_vector = descriptor.node_number_vector;
      Array<unsigned short> face_nb_edges(descriptor.face_to_node_vector.size());
      std::map<Edge, std::vector<FaceEdgeInfo>> edge_faces_map;
      for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
        const auto& face_nodes = descriptor.face_to_node_vector[l];
    
        face_nb_edges[l] = face_nodes.size();
        for (size_t r = 0; r < face_nodes.size() - 1; ++r) {
          Edge e({face_nodes[r], face_nodes[r + 1]}, node_number_vector);
          edge_faces_map[e].emplace_back(std::make_tuple(l, r, e.reversed()));
        }
        {
          Edge e({face_nodes[face_nodes.size() - 1], face_nodes[0]}, node_number_vector);
          edge_faces_map[e].emplace_back(std::make_tuple(l, face_nodes.size() - 1, e.reversed()));
        }
      }
    
      std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_id_map;
      {
        descriptor.face_to_edge_vector.resize(descriptor.face_to_node_vector.size());
        for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
          descriptor.face_to_edge_vector[l].resize(face_nb_edges[l]);
        }
        EdgeId e = 0;
        for (const auto& edge_faces_vector : edge_faces_map) {
          const auto& faces_vector = edge_faces_vector.second;
          for (unsigned short l = 0; l < faces_vector.size(); ++l) {
            const auto& [face_number, face_local_edge, reversed]         = faces_vector[l];
            descriptor.face_to_edge_vector[face_number][face_local_edge] = e;
          }
          edge_id_map[edge_faces_vector.first] = e;
          ++e;
        }
      }
    
      {
        descriptor.face_edge_is_reversed_vector.resize(descriptor.face_to_node_vector.size());
        for (FaceId j = 0; j < descriptor.face_edge_is_reversed_vector.size(); ++j) {
          descriptor.face_edge_is_reversed_vector[j] = Array<bool>(face_nb_edges[j]);
        }
        for (const auto& edge_faces_vector : edge_faces_map) {
          const auto& faces_vector = edge_faces_vector.second;
          for (unsigned short lj = 0; lj < faces_vector.size(); ++lj) {
            const auto& [face_number, face_local_edge, reversed]                  = faces_vector[lj];
            descriptor.face_edge_is_reversed_vector[face_number][face_local_edge] = reversed;
          }
        }
      }
    
      {
        descriptor.edge_to_node_vector.resize(edge_faces_map.size());
        int e = 0;
        for (const auto& edge_info : edge_faces_map) {
          const Edge& edge                  = edge_info.first;
          descriptor.edge_to_node_vector[e] = edge.nodeIdList();
          ++e;
        }
      }
    
      {
        // Edge numbers may change if numbers are provided in the file
        descriptor.edge_number_vector.resize(edge_faces_map.size());
        for (size_t e = 0; e < edge_faces_map.size(); ++e) {
          descriptor.edge_number_vector[e] = e;
        }
      }
    
      {
        descriptor.cell_to_edge_vector.reserve(descriptor.cell_to_node_vector.size());
        for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
          const auto& cell_nodes = descriptor.cell_to_node_vector[j];
    
          switch (descriptor.cell_type_vector[j]) {
          case CellType::Tetrahedron: {
            constexpr int local_edge[6][2] = {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {2, 3}, {3, 1}};
            std::vector<unsigned int> cell_edge_vector;
            cell_edge_vector.reserve(6);
            for (int i_edge = 0; i_edge < 6; ++i_edge) {
              const auto e = local_edge[i_edge];
              Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector};
              auto i = edge_id_map.find(edge);
              if (i == edge_id_map.end()) {
                throw NormalError("could not find this edge");
              }
              cell_edge_vector.push_back(i->second);
            }
            descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
            break;
          }
          case CellType::Hexahedron: {
            constexpr int local_edge[12][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
                                               {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
            std::vector<unsigned int> cell_edge_vector;
            cell_edge_vector.reserve(12);
            for (int i_edge = 0; i_edge < 12; ++i_edge) {
              const auto e = local_edge[i_edge];
              Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector};
              auto i = edge_id_map.find(edge);
              if (i == edge_id_map.end()) {
                throw NormalError("could not find this edge");
              }
              cell_edge_vector.push_back(i->second);
            }
            descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
            break;
          }
          case CellType::Pyramid: {
            const size_t number_of_edges = 2 * cell_nodes.size();
            std::vector<unsigned int> base_nodes;
            std::copy_n(cell_nodes.begin(), cell_nodes.size() - 1, std::back_inserter(base_nodes));
    
            std::vector<unsigned int> cell_edge_vector;
            cell_edge_vector.reserve(number_of_edges);
            for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
              Edge edge{{base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]}, node_number_vector};
              auto i = edge_id_map.find(edge);
              if (i == edge_id_map.end()) {
                throw NormalError("could not find this edge");
              }
              cell_edge_vector.push_back(i->second);
            }
    
            const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
            for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
              Edge edge{{base_nodes[i_edge], top_vertex}, node_number_vector};
              auto i = edge_id_map.find(edge);
              if (i == edge_id_map.end()) {
                throw NormalError("could not find this edge");
              }
              cell_edge_vector.push_back(i->second);
            }
            descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
            break;
          }
          case CellType::Diamond: {
            const size_t number_of_edges = 3 * cell_nodes.size();
            std::vector<unsigned int> base_nodes;
            std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
    
            std::vector<unsigned int> cell_edge_vector;
            cell_edge_vector.reserve(number_of_edges);
            for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
              Edge edge{{base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]}, node_number_vector};
              auto i = edge_id_map.find(edge);
              if (i == edge_id_map.end()) {
                throw NormalError("could not find this edge");
              }
              cell_edge_vector.push_back(i->second);
            }
    
            const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
            for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
              Edge edge{{base_nodes[i_edge], top_vertex}, node_number_vector};
              auto i = edge_id_map.find(edge);
              if (i == edge_id_map.end()) {
                throw NormalError("could not find this edge");
              }
              cell_edge_vector.push_back(i->second);
            }
    
            const unsigned int bottom_vertex = cell_nodes[0];
            for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
              Edge edge{{base_nodes[i_edge], bottom_vertex}, node_number_vector};
              auto i = edge_id_map.find(edge);
              if (i == edge_id_map.end()) {
                throw NormalError("could not find this edge");
              }
              cell_edge_vector.push_back(i->second);
            }
    
            descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
    
            break;
          }
          default: {
            std::stringstream error_msg;
            error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
            throw NotImplementedError(error_msg.str());
          }
          }
        }
      }
    }
    
    template void ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(ConnectivityDescriptor& descriptor);
    template void ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(ConnectivityDescriptor& descriptor);
    
    template void ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(
      ConnectivityDescriptor& descriptor);