[entitynet] compute/store intersection info upon network build

frackit/entitynetwork/containedentitynetwork.hh

@@ -26,9 +26,6 @@
 #ifndef FRACKIT_CONTAINED_ENTITY_NETWORK_HH
 #define FRACKIT_CONTAINED_ENTITY_NETWORK_HH
 
-#include <vector>
-#include <unordered_map>
-
 #include <TopTools_DataMapOfShapeInteger.hxx>
 #include <TopTools_ListOfShape.hxx>
 #include <TopoDS_Shape.hxx>
@@ -51,8 +48,9 @@ class ContainedEntityNetwork
 {
 
 public:
-    //! The type of underlying map from shape to primary entity index
+    //! Pull up typedefs
     using typename ContainedEntityNetworkInterface::ShapeMap;
+    using typename ContainedEntityNetworkInterface::ConnectivityMap;
 
     /*!
      * \brief Constructor.
@@ -63,10 +61,14 @@ public:
     ContainedEntityNetwork(int entityDim,
                            int domainDim,
                            ShapeMap&& sdFragmentsMap,
-                           ShapeMap&& entityFragmentsMap)
+                           ShapeMap&& entityFragmentsMap,
+                           ShapeMap&& intersectionFragmentMap,
+                           ConnectivityMap&& isConnectivityMap)
     : ContainedEntityNetworkInterface(entityDim, domainDim)
     , subDomainFragmentsMap_(std::move(sdFragmentsMap))
     , entityFragmentsMap_(std::move(entityFragmentsMap))
+    , intersectionFragmentsMap_(std::move(intersectionFragmentMap))
+    , isConnectivityMap_(std::move(isConnectivityMap))
     {}
 
     /*!
@@ -81,9 +83,23 @@ public:
     const ShapeMap& entityFragmentsMap() const override
     { return entityFragmentsMap_; }
 
+    /*!
+     * \brief The map with all intersections, mapping them to unique ids.
+     */
+    virtual const ShapeMap& intersectionFragmentsMap() const override
+    { return intersectionFragmentsMap_; }
+
+    /*!
+     * \brief The map of intersection id to intersecting entity indices.
+     */
+    virtual const ConnectivityMap& intersectionConnectivityMap() const override
+    { return isConnectivityMap_; }
+
 private:
     ShapeMap subDomainFragmentsMap_;
     ShapeMap entityFragmentsMap_;
+    ShapeMap intersectionFragmentsMap_;
+    ConnectivityMap isConnectivityMap_;
 };
 
 } // end namespace Frackit

frackit/entitynetwork/containedentitynetworkinterface.hh

@@ -26,7 +26,7 @@
 #define FRACKIT_CONTAINED_ENTITY_NETWORK_INTERFACE_HH
 
 #include <vector>
-
+#include <unordered_map>
 #include <TopTools_DataMapOfShapeInteger.hxx>
 
 namespace Frackit {
@@ -40,8 +40,10 @@ class ContainedEntityNetworkInterface
 {
 
 public:
-    //! export underlying map type
+    //! export underlying map types
     using ShapeMap = TopTools_DataMapOfShapeInteger;
+    using IndexList = std::vector<std::size_t>;
+    using ConnectivityMap = std::unordered_map<std::size_t, IndexList>;
 
     //! Abstract base classes need virtual destructors
     virtual ~ContainedEntityNetworkInterface() {}
@@ -79,6 +81,16 @@ public:
      */
     virtual const ShapeMap& entityFragmentsMap() const = 0;
 
+    /*!
+     * \brief The map with all intersections, mapping them to unique ids.
+     */
+    virtual const ShapeMap& intersectionFragmentsMap() const = 0;
+
+    /*!
+     * \brief The map of intersection id to intersecting entity indices.
+     */
+    virtual const ConnectivityMap& intersectionConnectivityMap() const = 0;
+
 private:
     int entityDimension_;
     int domainDimension_;

frackit/entitynetwork/entitynetwork.hh

@@ -46,8 +46,9 @@ class EntityNetwork
 {
 
 public:
-    //! The type of underlying map from shape to primary entity index
+    //! pull up types
     using typename EntityNetworkInterface::ShapeMap;
+    using typename EntityNetworkInterface::ConnectivityMap;
 
     /*!
      * \brief Constructor.
@@ -55,11 +56,19 @@ public:
      * \param fragmentMap Map containing the fragments of a network, where each
      *                         fragment is mapped to the index of the primary
      *                         entity of the network (before fragmentation).
+     * \param isFragmentMap Map containing the entity intersection of a network,
+     *                      where each one is mapped to a unique identifier.
+     * \param isConnectivityMap Map intersection identifiers to a list of primary
+     *                          entity indices that produced the intersection.
      */
     EntityNetwork(int entityDim,
-                  ShapeMap&& fragmentMap)
+                  ShapeMap&& fragmentMap,
+                  ShapeMap&& isFragmentMap,
+                  ConnectivityMap&& isConnectivityMap)
     : EntityNetworkInterface(entityDim)
     , entityFragmentsMap_(std::move(fragmentMap))
+    , intersectionFragmentsMap_(std::move(isFragmentMap))
+    , isConnectivityMap_(std::move(isConnectivityMap))
     {}
 
     /*!
@@ -69,8 +78,22 @@ public:
     virtual const ShapeMap& entityFragmentsMap() const override
     { return entityFragmentsMap_; }
 
+    /*!
+     * \brief The map with all intersections, mapping them to unique ids.
+     */
+    virtual const ShapeMap& intersectionFragmentsMap() const override
+    { return intersectionFragmentsMap_; }
+
+    /*!
+     * \brief The map of intersection id to intersecting entity indices.
+     */
+    virtual const ConnectivityMap& intersectionConnectivityMap() const override
+    { return isConnectivityMap_; }
+
 private:
     ShapeMap entityFragmentsMap_;
+    ShapeMap intersectionFragmentsMap_;
+    ConnectivityMap isConnectivityMap_;
 };
 
 } // end namespace Frackit

frackit/entitynetwork/entitynetworkinterface.hh

@@ -24,6 +24,8 @@
 #ifndef FRACKIT_ENTITY_NETWORK_INTERFACE_HH
 #define FRACKIT_ENTITY_NETWORK_INTERFACE_HH
 
+#include <vector>
+#include <unordered_map>
 #include <TopTools_DataMapOfShapeInteger.hxx>
 
 namespace Frackit {
@@ -35,8 +37,10 @@ namespace Frackit {
 class EntityNetworkInterface
 {
 public:
-    //! export underlying map type
+    //! export underlying map types
     using ShapeMap = TopTools_DataMapOfShapeInteger;
+    using IndexList = std::vector<std::size_t>;
+    using ConnectivityMap = std::unordered_map<std::size_t, IndexList>;
 
     //! Abstract base classes need virtual destructors
     virtual ~EntityNetworkInterface() {}
@@ -61,6 +65,16 @@ public:
      */
     virtual const ShapeMap& entityFragmentsMap() const = 0;
 
+    /*!
+     * \brief The map with all intersections, mapping them to unique ids.
+     */
+    virtual const ShapeMap& intersectionFragmentsMap() const = 0;
+
+    /*!
+     * \brief The map of intersection id to intersecting entity indices.
+     */
+    virtual const ConnectivityMap& intersectionConnectivityMap() const = 0;
+
 private:
     int entityDimension_;
 };

frackit/entitynetwork/networkbuilder.hh

@@ -25,6 +25,8 @@
 #define FRACKIT_ENTITY_NETWORK_BUILDER_HH
 
 #include <cmath>
+#include <tuple>
+#include <vector>
 #include <limits>
 #include <stdexcept>
 #include <unordered_map>
@@ -68,6 +70,7 @@ protected:
     using ShapeList = TopTools_ListOfShape;
     using ShapeIdMap = TopTools_DataMapOfShapeInteger;
     using ShapeIdMapIt = TopTools_DataMapIteratorOfDataMapOfShapeInteger;
+    using ConnectivityMap = std::unordered_map<std::size_t, std::vector<std::size_t>>;
 
 public:
 
@@ -369,6 +372,61 @@ protected:
         return {std::move(entityFragments), std::move(domainFragments)};
     }
 
+    /*!
+     * \brief Returns a tuple of maps, where in the first one the entity
+     *        intersections are mapped to corresponding unique identifiers,
+     *        and in the second one, the intersection identigiers map to the
+     *        primary entities that take part in the intersection.
+     */
+    std::tuple<ShapeIdMap, ConnectivityMap>
+    computeIntersectionData_(const ShapeIdMap& fragments)
+    {
+        ShapeIdMap isShapes;
+        ConnectivityMap isConnectivity;
+        std::size_t curIsIdx = 1;
+
+        // TODO: Start second iterator from ++it, somehow. A first implementation
+        //       of that led to some intersection edges not being meshed at the end
+        for (ShapeIdMapIt it(fragments); it.More(); it.Next())
+        {
+            for (ShapeIdMapIt it2(fragments); it2.More(); it2.Next())
+            {
+                const auto eIdx1 = it.Value();
+                const auto eIdx2 = it2.Value();
+                if (eIdx1 == eIdx2)
+                    continue;
+
+                const auto& f1 = TopoDS::Face(it.Key());
+                const auto& f2 = TopoDS::Face(it2.Key());
+                const auto edges = OCCUtilities::getCommonEdges(f1, f2);
+
+                for (const auto& edge : edges)
+                {
+                    int isIdx;
+                    if (!isShapes.Find(OCCUtilities::getShape(edge), isIdx))
+                    {
+                        isIdx = curIsIdx;
+                        isShapes.Bind(edge, isIdx);
+                        curIsIdx++;
+                    }
+
+                    isConnectivity[isIdx].push_back(eIdx1);
+                    isConnectivity[isIdx].push_back(eIdx2);
+                }
+            }
+        }
+
+        // make connectivity maps unique
+        for (auto& pair : isConnectivity)
+        {
+            auto& list = pair.second;
+            std::sort(list.begin(), list.end());
+            list.erase(std::unique(list.begin(), list.end()), list.end());
+        }
+
+        return {std::move(isShapes), std::move(isConnectivity)};
+    }
+
     /*!
      * \brief Returns the tolerance value.
      */
@@ -507,9 +565,14 @@ public:
         // Then fragment the network entities
         auto fragments = this->computeAllFragments_(confined, true)[0];
 
+        // get intersection data
+        auto [isMap, isConnectivity] = this->computeIntersectionData_(fragments);
+
         // create the network
         network_ = std::make_unique<EntityNetwork>(this->entityDimension_,
-                                                   std::move(fragments));
+                                                   std::move(fragments),
+                                                   std::move(isMap),
+                                                   std::move(isConnectivity));
         this->built_ = true;
         return *network_;
     }
@@ -562,14 +625,19 @@ public:
 
         // Then fragment the complete set of shapes
         auto entityAndDomainFragments = this->computeAllFragments_(confined);
-        auto& entities = entityAndDomainFragments[0];
-        auto& subDomains = entityAndDomainFragments[1];
+        auto& entityFragments = entityAndDomainFragments[0];
+        auto& subDomainFragments = entityAndDomainFragments[1];
+
+        // get intersection data
+        auto [isMap, isConnectivity] = this->computeIntersectionData_(entityFragments);
 
         // create the network
         network_ = std::make_unique<EntityNetwork>(this->entityDimension_,
                                                    this->domainDimension_,
-                                                   std::move(subDomains),
-                                                   std::move(entities));
+                                                   std::move(subDomainFragments),
+                                                   std::move(entityFragments),
+                                                   std::move(isMap),
+                                                   std::move(isConnectivity));
         this->built_ = true;
         return *network_;
     }

frackit/io/brepwriter.hh

@@ -65,7 +65,6 @@ public:
         if (network.entityDimension() != 2)
             throw std::runtime_error("BRepWriter only implemented for 2d networks so far");
 
-        findEntityIntersectionShapes_(network);
         makeSubShapeMaps_(network);
         makeCompound_();
     }
@@ -80,7 +79,6 @@ public:
         if (network.domainDimension() != network.entityDimension() + 1)
             throw std::runtime_error("BRepWriter expects entityDim = domainDim - 1");
 
-        findEntityIntersectionShapes_(network);
         makeSubShapeMaps_(network);
         makeCompound_();
     }
@@ -122,45 +120,6 @@ protected:
 
 private:
 
-    /*!
-     * \brief Determins the shapes that describe intersections of entities.
-     */
-    template<class Network>
-    void findEntityIntersectionShapes_(const Network& network)
-    {
-        using Iterator = typename Network::ShapeMap::Iterator;
-        const auto& map = network.entityFragmentsMap();
-
-        // TODO: Start second iterator from ++it, somehow. A first implementation
-        //       of that led to some intersection edges not being meshed at the end
-        for (Iterator it(map); it.More(); it.Next())
-        {
-            for (Iterator it2(map); it2.More(); it2.Next())
-            {
-                if (it.Value() == it2.Value())
-                    continue;
-
-                const auto& f1 = TopoDS::Face(it.Key());
-                const auto& f2 = TopoDS::Face(it2.Key());
-                const auto edges = OCCUtilities::getOverlapEdges(f1, f2);
-
-                std::vector<std::size_t> handled;
-                for (unsigned int eIdx = 0; eIdx < edges.size(); ++eIdx)
-                    for (const auto& isList : entityIntersections_)
-                        if (isList.Contains(edges[eIdx]))
-                            handled.push_back(eIdx);
-
-                if (handled.size() != edges.size())
-                {
-                    entityIntersections_.resize(entityIntersections_.size()+1);
-                    for (unsigned int eIdx = 0; eIdx < edges.size(); ++eIdx)
-                        if (!std::count(handled.begin(), handled.end(), eIdx))
-                            entityIntersections_.back().Append(edges[eIdx]);
-                }
-            }
-        }
-    }
-
     /*!
      * \brief Fill the sub-shape maps with all shapes that describe the network.
      * \note This overload is for contained entity networks.
@@ -174,13 +133,8 @@ private:
         for (Iterator it(network.entityFragmentsMap()); it.More(); it.Next())
             addEntity_(network.entityDimension(), it.Key(), it.Value());
 
-        std::size_t isIdx = 1;
-        for (const auto& isShapeList : entityIntersections_)
-        {
-            for (const auto& isShape : isShapeList)
-                addIntersection_(network.entityDimension()-1, isShape, isIdx);
-            isIdx++;
-        }
+        for (Iterator it(network.intersectionFragmentsMap()); it.More(); it.Next())
+            addIntersection_(network.entityDimension()-1, it.Key(), it.Value());
     }
 
     /*!
@@ -193,13 +147,8 @@ private:
         for (Iterator it(network.entityFragmentsMap()); it.More(); it.Next())
             addEntity_(network.entityDimension(), it.Key(), it.Value());
 
-        std::size_t isIdx = 1;
-        for (const auto& isShapeList : entityIntersections_)
-        {
-            for (const auto& isShape : isShapeList)
-                addIntersection_(network.entityDimension()-1, isShape, isIdx);
-            isIdx++;
-        }
+        for (Iterator it(network.intersectionFragmentsMap()); it.More(); it.Next())
+            addIntersection_(network.entityDimension()-1, it.Key(), it.Value());
     }
 
     /*!
@@ -355,9 +304,6 @@ private:
         for(Standard_Integer i = 1; i <= somap_.Extent(); i++) b.Add(compound_, somap_(i));
     }
 
-    // containers to store the primary entities
-    std::vector<TopTools_ListOfShape> entityIntersections_;
-
     // sub-shape maps
     TopTools_IndexedMapOfShape vmap_, emap_, wmap_, fmap_, shmap_, somap_;