TriangulationRExtractor, some debugging, tests

src/db/db/dbPLCTriangulation.cc

@@ -1043,7 +1043,7 @@ Triangulation::search_edges_crossing (Vertex *from, Vertex *to)
         if (e->has_vertex (vv)) {
           return result;
         }
-        if (e->crosses (edge)) {
+        if (e->crosses_including (edge)) {
           result.push_back (e);
           next_edge = e;
           break;

src/pex/pex/pex.pro

@@ -10,12 +10,14 @@ SOURCES = \
   pexForceLink.cc \
   pexRExtractor.cc \
   gsiDeclRExtractor.cc \
-  pexSquareCountingRExtractor.cc
+  pexSquareCountingRExtractor.cc \
+  pexTriangulationRExtractor.cc
 
 HEADERS = \
   pexForceLink.h \
   pexRExtractor.h \
-  pexSquareCountingRExtractor.h
+  pexSquareCountingRExtractor.h \
+  pexTriangulationRExtractor.h
 
 RESOURCES = \
 

src/pex/pex/pexRExtractor.cc

@@ -134,7 +134,12 @@ RNetwork::create_node (RNode::node_type type, unsigned int port_index)
 RElement *
 RNetwork::create_element (double conductivity, RNode *a, RNode *b)
 {
-  auto i = m_elements_by_nodes.find (std::make_pair (a, b));
+  std::pair<RNode *, RNode *> key (a, b);
+  if (size_t (b) < size_t (a)) {
+    std::swap (key.first, key.second);
+  }
+
+  auto i = m_elements_by_nodes.find (key);
   if (i != m_elements_by_nodes.end ()) {
 
     if (conductivity == pex::RElement::short_value () || i->second->conductivity == pex::RElement::short_value ()) {
@@ -149,7 +154,7 @@ RNetwork::create_element (double conductivity, RNode *a, RNode *b)
 
     RElement *element = new RElement (this, conductivity, a, b);
     m_elements.push_back (element);
-    m_elements_by_nodes.insert (std::make_pair (std::make_pair (a, b), element));
+    m_elements_by_nodes.insert (std::make_pair (key, element));
 
     a->m_elements.push_back (element);
     element->m_ia = --a->m_elements.end ();

src/pex/pex/pexRExtractor.h

@@ -142,6 +142,11 @@ class PEX_PUBLIC RNetwork
   : public tl::Object
 {
 public:
+  typedef tl::list<RNode, false> node_list;
+  typedef node_list::const_iterator node_iterator;
+  typedef tl::list<RElement, false> element_list;
+  typedef element_list::const_iterator element_iterator;
+
   RNetwork ();
   ~RNetwork ();
 
@@ -154,9 +159,50 @@ public:
 
   std::string to_string () const;
 
+  node_iterator begin_nodes () const
+  {
+    return m_nodes.begin ();
+  }
+
+  node_iterator end_nodes () const
+  {
+    return m_nodes.end ();
+  }
+
+  size_t num_nodes () const
+  {
+    return m_nodes.size ();
+  }
+
+  size_t num_internal_nodes () const
+  {
+    size_t count = 0;
+    for (auto n = m_nodes.begin (); n != m_nodes.end (); ++n) {
+      if (n->type == pex::RNode::Internal) {
+        ++count;
+      }
+    }
+    return count;
+  }
+
+  element_iterator begin_elements () const
+  {
+    return m_elements.begin ();
+  }
+
+  element_iterator end_elements () const
+  {
+    return m_elements.end ();
+  }
+
+  size_t num_elements () const
+  {
+    return m_elements.size ();
+  }
+
 private:
-  tl::list<RNode, false> m_nodes;
-  tl::list<RElement, false> m_elements;
+  node_list m_nodes;
+  element_list m_elements;
   std::map<std::pair<RNode *, RNode *>, RElement *> m_elements_by_nodes;
   std::map<std::pair<RNode::node_type, unsigned int>, RNode *> m_nodes_by_type;
 

src/pex/pex/pexTriangulationRExtractor.cc

@@ -0,0 +1,221 @@
+
+/*
+
+  KLayout Layout Viewer
+  Copyright (C) 2006-2025 Matthias Koefferlein
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2 of the License, or
+  (at your option) any later version.
+
+  This program 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 General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+
+*/
+
+
+#include "pexTriangulationRExtractor.h"
+#include "dbBoxScanner.h"
+#include "dbPolygonTools.h"
+#include "tlIntervalMap.h"
+
+namespace pex
+{
+
+TriangulationRExtractor::TriangulationRExtractor (double dbu)
+{
+  m_dbu = dbu;
+
+  m_tri_param.min_b = 0.3;
+  m_tri_param.max_area = 0.0;
+}
+
+void
+TriangulationRExtractor::extract (const db::Polygon &polygon, const std::vector<db::Point> &vertex_ports, const std::vector<db::Polygon> &polygon_ports, pex::RNetwork &rnetwork)
+{
+  rnetwork.clear ();
+
+  db::CplxTrans trans = db::CplxTrans (m_dbu) * db::ICplxTrans (db::Trans (db::Point () - polygon.box ().center ()));
+  auto inv_trans = trans.inverted ();
+
+  //  NOTE: currently we treat polygon ports and points where the location is the center of the bounding box
+  std::vector<db::Point> vp = vertex_ports;
+  vp.reserve (vertex_ports.size () + polygon_ports.size ());
+  for (auto pp = polygon_ports.begin (); pp != polygon_ports.end (); ++pp) {
+    vp.push_back (pp->box ().center ());
+  }
+
+
+  db::plc::Graph plc;
+  db::plc::Triangulation tri (&plc);
+
+  tri.triangulate (polygon, vp, m_tri_param, trans);
+
+  //  create a network node for each triangle node
+
+  std::unordered_map<const db::plc::Vertex *, pex::RNode *> vertex2node;
+
+  size_t internal_node_id = 0;
+
+  for (auto p = plc.begin (); p != plc.end (); ++p) {
+
+    for (size_t iv = 0; iv < p->size (); ++iv) {
+
+      const db::plc::Vertex *vertex = p->vertex (iv);
+      if (vertex2node.find (vertex) != vertex2node.end ()) {
+        continue;
+      }
+
+      pex::RNode::node_type type = pex::RNode::Internal;
+      size_t port_index = 0;
+
+      if (vertex->is_precious ()) {
+        size_t idx = vertex->id ();
+        if (idx >= vertex_ports.size ()) {
+          type = pex::RNode::PolygonPort;
+          port_index = size_t (idx) - vertex_ports.size ();
+        } else {
+          type = pex::RNode::VertexPort;
+          port_index = size_t (idx);
+        }
+      } else {
+        port_index = internal_node_id++;
+      }
+
+      pex::RNode *n = rnetwork.create_node (type, port_index);
+      db::DPoint loc = *vertex;
+      n->location = db::DBox (loc, loc);
+
+      vertex2node.insert (std::make_pair (vertex, n));
+
+    }
+
+  }
+
+  //  produce the conductances for each triangle
+
+  for (auto p = plc.begin (); p != plc.end (); ++p) {
+    create_conductances (*p, vertex2node, rnetwork);
+  }
+
+  //  eliminate internal nodes
+
+  eliminate_all (rnetwork);
+}
+
+void
+TriangulationRExtractor::create_conductances (const db::plc::Polygon &tri, const std::unordered_map<const db::plc::Vertex *, pex::RNode *> &vertex2node, RNetwork &rnetwork)
+{
+  tl_assert (tri.size () == 3);
+
+  for (int i = 0; i < 3; ++i) {
+
+    const db::plc::Vertex *pm1 = tri.vertex (i);
+    const db::plc::Vertex *p0 = tri.vertex (i + 1);
+    const db::plc::Vertex *p1 = tri.vertex (i + 2);
+
+    double a = fabs (db::vprod (*pm1 - *p0, *p1 - *p0) * 0.5);
+
+    double lm1 = (*p0 - *pm1).sq_length ();
+    double l0 = (*p1 - *p0).sq_length ();
+    double l1 = (*pm1 - *p1).sq_length ();
+
+    double s = (l0 + l1 - lm1) / (8.0 * a);
+
+    auto i0 = vertex2node.find (p0);
+    auto im1 = vertex2node.find (pm1);
+    rnetwork.create_element (s, i0->second, im1->second);
+
+  }
+}
+
+void
+TriangulationRExtractor::eliminate_all (RNetwork &rnetwork)
+{
+  if (tl::verbosity () >= m_tri_param.base_verbosity + 10) {
+    tl::info << "Staring elimination with " << rnetwork.num_internal_nodes () << " internal nodes and " << rnetwork.num_elements () << " resistors";
+  }
+
+  unsigned int niter = 0;
+  std::vector<pex::RNode *> to_eliminate;
+
+  size_t nmax = 3;
+  while (nmax > 0) {
+
+    bool another_loop = true;
+    while (another_loop) {
+
+      size_t nmax_next = 0;
+      to_eliminate.clear ();
+
+      for (auto n = rnetwork.begin_nodes (); n != rnetwork.end_nodes (); ++n) {
+        if (n->type == pex::RNode::Internal) {
+          size_t nn = n->elements ().size ();
+          if (nn <= nmax) {
+            to_eliminate.push_back (const_cast<pex::RNode *> (n.operator-> ()));
+          } else if (nmax_next == 0 or nn < nmax_next) {
+            nmax_next = nn;
+          }
+        }
+      }
+
+      if (to_eliminate.empty ()) {
+
+        another_loop = false;
+        nmax = nmax_next;
+
+        if (tl::verbosity () >= m_tri_param.base_verbosity + 10) {
+          tl::info << "Nothing left to eliminate with nmax=" << nmax;
+        }
+
+      } else {
+
+        for (auto n = to_eliminate.begin (); n != to_eliminate.end (); ++n) {
+          eliminate_node (*n, rnetwork);
+        }
+
+        niter += 1;
+
+        if (tl::verbosity () >= m_tri_param.base_verbosity + 10) {
+          tl::info << "Nodes left after iteration " << niter << " with nmax=" << nmax << ": " << rnetwork.num_internal_nodes () << " with " << rnetwork.num_elements () << " edges.";
+        }
+
+      }
+
+    }
+
+  }
+}
+
+void
+TriangulationRExtractor::eliminate_node (pex::RNode *node, RNetwork &rnetwork)
+{
+  double s_sum = 0.0;
+  for (auto e = node->elements ().begin (); e != node->elements ().end (); ++e) {
+    s_sum += (*e)->conductivity;
+  }
+
+  if (fabs (s_sum) > 1e-10) {
+    for (auto e = node->elements ().begin (); e != node->elements ().end (); ++e) {
+      auto ee = e;
+      ++ee;
+      for ( ; ee != node->elements ().end (); ++ee) {
+        pex::RNode *n1 = const_cast <pex::RNode *> ((*e)->other (node));
+        pex::RNode *n2 = const_cast <pex::RNode *> ((*ee)->other (node));
+        double c = (*e)->conductivity * (*ee)->conductivity / s_sum;
+        rnetwork.create_element (c, n1, n2);
+      }
+    }
+  }
+
+  rnetwork.remove_node (node);
+}
+
+}

src/pex/pex/pexTriangulationRExtractor.h

@@ -0,0 +1,70 @@
+
+/*
+
+  KLayout Layout Viewer
+  Copyright (C) 2006-2025 Matthias Koefferlein
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2 of the License, or
+  (at your option) any later version.
+
+  This program 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 General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+
+*/
+
+#ifndef HDR_pexTriangulationRExtractor
+#define HDR_pexTriangulationRExtractor
+
+#include "pexCommon.h"
+#include "pexRExtractor.h"
+
+#include "dbPLCTriangulation.h"
+
+namespace pex
+{
+
+// @@@ doc
+class PEX_PUBLIC TriangulationRExtractor
+  : public RExtractor
+{
+public:
+  TriangulationRExtractor (double dbu);
+
+  db::plc::TriangulationParameters &triangulation_parameters ()
+  {
+    return m_tri_param;
+  }
+
+  void set_dbu (double dbu)
+  {
+    m_dbu = dbu;
+  }
+
+  double dbu () const
+  {
+    return m_dbu;
+  }
+
+  virtual void extract (const db::Polygon &polygon, const std::vector<db::Point> &vertex_ports, const std::vector<db::Polygon> &polygon_ports, RNetwork &rnetwork);
+
+private:
+  db::plc::TriangulationParameters m_tri_param;
+  double m_dbu;
+
+  void create_conductances (const db::plc::Polygon &tri, const std::unordered_map<const db::plc::Vertex *, RNode *> &vertex2node, RNetwork &rnetwork);
+  void eliminate_node (pex::RNode *node, RNetwork &rnetwork);
+  void eliminate_all (RNetwork &rnetwork);
+};
+
+}
+
+#endif
+

src/pex/unit_tests/pexRExtractorTests.cc

@@ -56,6 +56,15 @@ TEST(network_basic)
     "R $2 $3 0.2"
   );
 
+  pex::RElement *e23c = rn.create_element (5.0, n3, n2);
+  EXPECT_EQ (e23 == e23c, true);
+
+  EXPECT_EQ (rn.to_string (),
+    "R $1 $2 2\n"
+    "R $1 $3 4\n"
+    "R $2 $3 0.1"
+  );
+
   rn.remove_element (e23);
 
   EXPECT_EQ (rn.to_string (),

src/pex/unit_tests/pexTriangulationRExtractorTests.cc

@@ -0,0 +1,69 @@
+
+/*
+
+  KLayout Layout Viewer
+  Copyright (C) 2006-2025 Matthias Koefferlein
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2 of the License, or
+  (at your option) any later version.
+
+  This program 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 General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+
+*/
+
+
+#include "pexTriangulationRExtractor.h"
+#include "tlUnitTest.h"
+
+namespace
+{
+
+class TestableTriangulationRExtractor
+  : public pex::TriangulationRExtractor
+{
+public:
+  TestableTriangulationRExtractor ()
+    : pex::TriangulationRExtractor (0.001)
+  { }
+};
+
+}
+
+TEST(extraction)
+{
+  db::Point contour[] = {
+    db::Point (0, 0),
+    db::Point (0, 100),
+    db::Point (1000, 100),
+    db::Point (1000, 0)
+  };
+
+  db::Polygon poly;
+  poly.assign_hull (contour + 0, contour + sizeof (contour) / sizeof (contour[0]));
+
+  double dbu = 0.001;
+
+  pex::RNetwork rn;
+  pex::TriangulationRExtractor rex (dbu);
+
+  std::vector<db::Point> vertex_ports;
+  vertex_ports.push_back (db::Point (0, 50));     //  V0
+  vertex_ports.push_back (db::Point (1000, 50));  //  V1
+
+  std::vector<db::Polygon> polygon_ports;
+
+  rex.extract (poly, vertex_ports, polygon_ports, rn);
+
+  EXPECT_EQ (rn.to_string (),
+    "R V1 V0 10.0938"
+  )
+}

src/pex/unit_tests/unit_tests.pro

@@ -8,7 +8,8 @@ include($$PWD/../../lib_ut.pri)
 
 SOURCES = \
   pexRExtractorTests.cc \
-  pexSquareCountingRExtractorTests.cc
+  pexSquareCountingRExtractorTests.cc \
+  pexTriangulationRExtractorTests.cc
 
 INCLUDEPATH += $$TL_INC $$DB_INC $$GSI_INC $$PEX_INC
 DEPENDPATH += $$TL_INC $$DB_INC $$GSI_INC $$PEX_INC