diff --git a/mgrid_2/geom.cpp b/mgrid_2/geom.cpp
index 08bf123..7ed1ef0 100644
--- a/mgrid_2/geom.cpp
+++ b/mgrid_2/geom.cpp
@@ -77,7 +77,7 @@ void Mesh::Init_Solution_mult(std::vector<double> &v,
 
     for (int e = 0; e < Nelems(); ++e)      // loop over all elements
     {  
-        int sd = ElementSubdomains[e];      // get subdomain of element e
+        int sd = _elementSubdomains[e];      // get subdomain of element e
         if (sd == target_sd)                // if is target subdomain then
         {
             int base = e * _nvert_e;        // get starting index of element in coordinate vector
@@ -506,7 +506,7 @@ double Mesh::AverageVectorFunction_perSubdomain(std::vector<double> &v, int targ
     int subdomain_element_counter = 0;
     for (int e = 0; e < Nelems(); ++e)      // loop over all elements
     {  
-        int sd = ElementSubdomains[e];      // get subdomain of element e
+        int sd = _elementSubdomains[e];      // get subdomain of element e
         if (sd == target_sd)                // if is target subdomain then
         {
             subdomain_element_counter++;
@@ -534,7 +534,7 @@ double Mesh::CheckTemp_mult(std::vector<double> &v, int target_sd, double goal_t
     int subdomain_element_counter = 0;
     for (int e = 0; e < Nelems(); ++e)      // loop over all elements
     {  
-        int sd = ElementSubdomains[e];      // get subdomain of element e
+        int sd = _elementSubdomains[e];      // get subdomain of element e
         if (sd == target_sd)                // if is target subdomain then
         {
             subdomain_element_counter++;
@@ -822,31 +822,56 @@ const std::vector<int> Mesh::BoundaryEdgeNodes() const
     return _bedges;
 }
 
+const std::vector<int> Mesh::OuterEdges() const
+{
+    return _outerEdges;
+}
+
+const std::vector<int> Mesh::OuterEdgesNodes() const
+{
+    return _outerEdgesNodes;
+}
+
+const std::vector<int> Mesh::OuterEdgesSubdomains() const
+{
+    return _outerEdgesSubdomains;
+}
+
+const std::vector<int> Mesh::ElementSubdomains() const
+{
+    return _elementSubdomains;
+}
+
+const std::vector<int> Mesh::EdgeSubdomains() const
+{
+    return _edgeSubdomains;
+}
+
 
 // Only the outer edges for Robin BC
 void Mesh::InitializeOuterEdges()
 {
-    assert(EdgeSubdomains.size() == 2*_ebedges.size());
+    assert(_edgeSubdomains.size() == 2*_ebedges.size());
     for (size_t k = 0; k < _ebedges.size(); ++k) // loop over all boundary edges
     {
-        if (EdgeSubdomains[2*k] == -1)
+        if (_edgeSubdomains[2*k] == -1)
         {
-            OuterEdges.push_back(_ebedges[k]);
-            OuterEdgesSubdomains.push_back(EdgeSubdomains[2*k + 1]);
-            OuterEdgesNodes.push_back(_bedges[2*k]);
-            OuterEdgesNodes.push_back(_bedges[2*k + 1]);
+            _outerEdges.push_back(_ebedges[k]);
+            _outerEdgesSubdomains.push_back(_edgeSubdomains[2*k + 1]);
+            _outerEdgesNodes.push_back(_bedges[2*k]);
+            _outerEdgesNodes.push_back(_bedges[2*k + 1]);
         }
-        if (EdgeSubdomains[2*k + 1] == -1)
+        if (_edgeSubdomains[2*k + 1] == -1)
         {            
-            OuterEdges.push_back(_ebedges[k]);
-            OuterEdgesSubdomains.push_back(EdgeSubdomains[2*k]);
-            OuterEdgesNodes.push_back(_bedges[2*k]);
-            OuterEdgesNodes.push_back(_bedges[2*k + 1]);
+            _outerEdges.push_back(_ebedges[k]);
+            _outerEdgesSubdomains.push_back(_edgeSubdomains[2*k]);
+            _outerEdgesNodes.push_back(_bedges[2*k]);
+            _outerEdgesNodes.push_back(_bedges[2*k + 1]);
         }
     }
 
     cout << "All boundary edges: " << _ebedges.size() << endl;
-    cout << "Outer boundary edges: " << OuterEdges.size() << endl;
+    cout << "Outer boundary edges: " << _outerEdges.size() << endl;
 }
 
 
@@ -1088,8 +1113,9 @@ Mesh::Mesh(std::string const &fname)
 // Constructor with subdomain support
 Mesh::Mesh(std::string const &filename, std::string const &subdomain_filename) : Mesh(filename)
 {
-    ElementSubdomains = ReadElementSubdomains(subdomain_filename);
+    _elementSubdomains = ReadElementSubdomains(subdomain_filename);
     InitializeOuterEdges();
+    
 }
 
 
@@ -1181,11 +1207,11 @@ void Mesh::ReadVertexBasedMesh(std::string const &fname)
         }
 
         // Store edge adjacent subdomains
-        EdgeSubdomains.resize(nbedges * 2);
+        _edgeSubdomains.resize(nbedges * 2);
         for (int k = 0; k < nbedges * 2; ++k)
         {
-            ifs >> EdgeSubdomains[k];
-            EdgeSubdomains[k] -= OFFSET;
+            ifs >> _edgeSubdomains[k];
+            _edgeSubdomains[k] -= OFFSET;
         }
     }
     else
diff --git a/mgrid_2/geom.h b/mgrid_2/geom.h
index 84717c5..70698b3 100644
--- a/mgrid_2/geom.h
+++ b/mgrid_2/geom.h
@@ -332,14 +332,27 @@ private:
     void Write_ascii_paraview_3D(std::string const &fname, std::vector<double> const &v) const;
 
     void InitializeOuterEdges();
+
+    std::vector<int> _outerEdges;
+    std::vector<int> _outerEdgesSubdomains;
+    std::vector<int> _outerEdgesNodes;
+
+    // Every element belongs to 1 subdomain
+    std::vector<int> _elementSubdomains;
+
+    // Every edge has 2 adjacent subdomains
+    std::vector<int> _edgeSubdomains;
+
     
 public:
     const std::vector<int> BoundaryEdges() const;
     const std::vector<int> BoundaryEdgeNodes() const;
+    const std::vector<int> OuterEdges() const;
+    const std::vector<int> OuterEdgesSubdomains() const;
+    const std::vector<int> OuterEdgesNodes() const;
+    const std::vector<int> ElementSubdomains() const;
+    const std::vector<int> EdgeSubdomains() const;
 
-    std::vector<int> OuterEdges;
-    std::vector<int> OuterEdgesSubdomains;
-    std::vector<int> OuterEdgesNodes;
 
 
 
@@ -561,11 +574,7 @@ 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.
diff --git a/mgrid_2/getmatrix.cpp b/mgrid_2/getmatrix.cpp
index 394adf1..9d5378a 100644
--- a/mgrid_2/getmatrix.cpp
+++ b/mgrid_2/getmatrix.cpp
@@ -404,7 +404,7 @@ void FEM_Matrix::CalculateLaplace_mult(vector<double> &f)
     auto const &ia   = _mesh.GetConnectivity();
     auto const &xc   = _mesh.GetCoords();
 
-    const vector<int> sd_vec = _mesh.ElementSubdomains;
+    const vector<int> sd_vec = _mesh.ElementSubdomains();
 
     #pragma omp parallel for private(ske,fe)
     for (int i = 0; i < nelem; ++i) {
@@ -498,7 +498,7 @@ void FEM_Matrix::AddMass_mult(vector<double> &f, const double scale_factor)
     auto const &ia   = _mesh.GetConnectivity();
     auto const &xc   = _mesh.GetCoords();
 
-    const vector<int> sd_vec = _mesh.ElementSubdomains;
+    const vector<int> sd_vec = _mesh.ElementSubdomains();
 
     #pragma omp parallel for private(ske,fe)
     for (int i = 0; i < nelem; ++i) {
@@ -620,9 +620,9 @@ void FEM_Matrix::ApplyDirichletBC(std::vector<double> const &u, std::vector<doub
 
 void FEM_Matrix::ApplyRobinBC_mult(std::vector<double> &f, const double u_out)
 {
-    auto const RobinEdges = _mesh.OuterEdges;
-    auto const RobinEdgesSubdomains = _mesh.OuterEdgesSubdomains;
-    auto const RobinEdgeNodes = _mesh.OuterEdgesNodes;
+    auto const RobinEdges = _mesh.OuterEdges();
+    auto const RobinEdgesSubdomains = _mesh.OuterEdgesSubdomains();
+    auto const RobinEdgeNodes = _mesh.OuterEdgesNodes();
     assert(RobinEdgeNodes.size() == 2 * RobinEdges.size());
 
     vector<double> const Coordinates = _mesh.GetCoords();
@@ -632,7 +632,7 @@ void FEM_Matrix::ApplyRobinBC_mult(std::vector<double> &f, const double u_out)
         int const subdomain = RobinEdgesSubdomains[i];
         double const alpha = Heat_transfer_coefficient(subdomain);
 
-        int const EdgeNumber = RobinEdges[i];
+        //int const EdgeNumber = RobinEdges[i];
 
         int const EdgeNode1 = RobinEdgeNodes[2*i];
         int const EdgeNode2 = RobinEdgeNodes[2*i + 1];