diff --git a/mgrid_2/geom.cpp b/mgrid_2/geom.cpp
index a2dbcf2..019f5a5 100644
--- a/mgrid_2/geom.cpp
+++ b/mgrid_2/geom.cpp
@@ -129,6 +129,8 @@ void Mesh::DebugEdgeBased() const
     }
 
     cout << "\n ...............    edges        ...................\n";
+    cout << "_edges.size(): " << _edges.size() << endl;
+    cout << "_nedge: " << _nedge << endl;
     for (int k = 0; k < _nedge; ++k)
     {
         cout << k << " : ";
@@ -147,6 +149,7 @@ void Mesh::DebugEdgeBased() const
         cout << endl;
     }
     cout << "\n ...............    Boundary (edges)    .................\n";
+    cout << "_ebedges.size(): " << _ebedges.size() << endl;
     cout << " _ebedges : " << _ebedges << endl;
 
     return;
@@ -751,6 +754,33 @@ void Mesh::DeriveEdgeFromVertexBased_fast()
 }
 // HG
 
+const std::vector<int> Mesh::BoundaryEdges() const
+{
+    return _ebedges;
+}
+
+
+// Only the outer edges for Robin BC
+const std::vector<int> Mesh::OuterEdges() const
+{
+    vector<int> outerEdges;
+
+    for (int k = 0; k < _ebedges.size()/2; ++k)
+    {
+        if (EdgeSubdomains[k] == 0 || EdgeSubdomains[k + 1] == 0)
+            outerEdges.push_back(_ebedges[k]);
+    }
+
+    cout << "All boundary edges: " << _ebedges.size() << endl;
+    cout << "Outer boundary edges: " << outerEdges.size() << endl;
+
+    cout << _ebedges;
+    cout << endl;
+    cout << outerEdges << endl;
+
+    return outerEdges;
+}
+
 
 #include <utility>             // pair
 
@@ -982,16 +1012,33 @@ Mesh::Mesh(std::string const &fname)
     //cout << " P E R M U T E D !" << endl;
 }
 
+
+// Constructor with subdomain support
 Mesh::Mesh(std::string const &filename, std::string const &subdomain_filename) : Mesh(filename)
 {
     ElementSubdomains = ReadElementSubdomains(subdomain_filename);
+    //EdgeSubdomains = ReadEdgeSubdomains(filename);
 }
 
+// const vector<int> Mesh::ReadEdgeSubdomains(std::string const &filename) const
+// {
+//     vector<int> edgeSubdomains(_nedge);
+
+//     ifstream ifs(filename);
+//     if (!(ifs.is_open() && ifs.good())) {
+//         cerr << "Mesh::ReadEdgeSubdomains: Error cannot open file " << filename << endl;
+//         assert(ifs.is_open());
+//     }
+
+//     return edgeSubdomains;
+// }
+
+
 const vector<int> Mesh::ReadElementSubdomains(string const &dname) const
 {
     ifstream ifs(dname);
     if (!(ifs.is_open() && ifs.good())) {
-        cerr << "ParMesh::ReadElementSubdomain: Error cannot open file " << dname << endl;
+        cerr << "Mesh::ReadElementSubdomain: Error cannot open file " << dname << endl;
         assert(ifs.is_open());
     }
 
@@ -1073,12 +1120,20 @@ void Mesh::ReadVertexBasedMesh(std::string const &fname)
             ifs >> _bedges[k];
             _bedges[k] -= OFFSET;            // Matlab to C indexing
         }
+
+        // Store edge adjacent subdomains
+        EdgeSubdomains.resize(nbedges * 2);
+        for (int k = 0; k < nbedges * 2; ++k)
+        {
+            ifs >> EdgeSubdomains[k];
+        }
     }
     else
     {
         // ToDo: add boundary information to 3D mesh
         cout << std::endl << "NO boundary information available for 3D mesh" << endl;
     }
+
     return;
 }
 
diff --git a/mgrid_2/geom.h b/mgrid_2/geom.h
index 2974c46..ac7787f 100644
--- a/mgrid_2/geom.h
+++ b/mgrid_2/geom.h
@@ -328,6 +328,10 @@ private:
     void Write_ascii_paraview_3D(std::string const &fname, std::vector<double> const &v) const;
     
 public:
+    const std::vector<int> BoundaryEdges() const;
+    const std::vector<int> OuterEdges() const;
+
+
      /**
      * Visualize @p v together with its mesh information via paraview
      *
@@ -546,9 +550,12 @@ public:
      */    
     [[nodiscard]] bool checkObtuseAngles() const;
 
-
+    // Every element belongs to 1 subdomain
     std::vector<int> ElementSubdomains;
 
+    // Every edge has 2 adjacent subdomains
+    std::vector<int> EdgeSubdomains;
+
     /**
      * Reads the global triangle to subdomain mapping.
      *
@@ -558,6 +565,8 @@ public:
     */
     [[nodiscard]] const std::vector<int> ReadElementSubdomains(std::string const &dname) const;
 
+    [[nodiscard]] const std::vector<int> ReadEdgeSubdomains(std::string const &filename) const;
+
 
     /**
      * Calculates the largest inner angle in element @p idx.