diff --git a/m3dc1_scorec/api/m3dc1_scorec.cc b/m3dc1_scorec/api/m3dc1_scorec.cc
index e84aa7fa0..ddaeb9183 100644
--- a/m3dc1_scorec/api/m3dc1_scorec.cc
+++ b/m3dc1_scorec/api/m3dc1_scorec.cc
@@ -869,7 +869,7 @@ int m3dc1_model_getmaxcoord(double* x_max, double* y_max)
 //*******************************************************
 int m3dc1_model_load(char* /* in */ model_file)
 //*******************************************************
-{  
+{ 
   FILE *test_in = fopen (model_file,"r");
   if (!test_in)
   {
@@ -880,19 +880,132 @@ int m3dc1_model_load(char* /* in */ model_file)
   else
     fclose(test_in);
 
-  pGeom g = pumi_geom_load(model_file, "analytic");
-  m3dc1_model::instance()->model = g->getGmi();
-  m3dc1_model::instance()->load_analytic_model(model_file); 
+  pGeom g;
+  std::string model_name = model_file; 
+  if(model_name.substr(model_name.find_last_of(".") + 1) == "dmg")
+  { // modelType = 2 (.dmg model)
+    g = pumi_geom_load(model_file, "mesh");
+    m3dc1_model::instance()->model = g->getGmi();   
+    gmi_register_mesh();
+    m3dc1_model::instance()->model = gmi_load(model_file);
+    m3dc1_model::instance()->modelType = 2;
+  }
+  else
+  { // modelType = 1 (non .dmg models)
+    g = pumi_geom_load(model_file, "analytic");
+    m3dc1_model::instance()->model = g->getGmi();
+    m3dc1_model::instance()->load_analytic_model(model_file); 
+    m3dc1_model::instance()->modelType = 1;
+  }
   pumi_geom_freeze(g);
 
-  m3dc1_model::instance()->caculateBoundingBox();
-  // save the num of geo ent on the oringal plane
+  if (m3dc1_model::instance()->modelType == 1)
+    m3dc1_model::instance()->caculateBoundingBox();
+
+  // save the num of geo ent on the oringal plan
   m3dc1_model::instance()->numEntOrig[0]=m3dc1_model::instance()->model->n[0];
-  //if (m3dc1_model::instance()->model->n[1]==1) assert(m3dc1_model::instance()->numEntOrig[0]==0); // for a smooth loop, there is no geo vtx 
   m3dc1_model::instance()->numEntOrig[1]=m3dc1_model::instance()->model->n[1];
   m3dc1_model::instance()->numEntOrig[2]=m3dc1_model::instance()->model->n[2];
   return M3DC1_SUCCESS;
 }
+//*******************************************************
+void m3dc1_model_getmodeltype(int* modeltype)
+//*******************************************************
+{
+  // 1 for .txt model, 2 for .dmg model
+  *modeltype = m3dc1_model::instance()->modelType;
+}
+
+//*******************************************************
+int m3dc1_modelinfo_load(char* /* in */ modelInfo_file)
+//*******************************************************
+{
+  if (m3dc1_model::instance()->modelType == 1)
+    return M3DC1_SUCCESS;  // only move forward for .dmg models
+
+  assert(m3dc1_model::instance()->modelType == 2);
+  if (std::strcmp(modelInfo_file, "dummyInfo") == 0)
+  {
+    std::cout << "Make sure to provide model_info file since model is .dmg\n";
+    exit(1);  // Need to have model info file
+  }
+
+  std::cout << "Model Info file = " << modelInfo_file << "\n";
+  
+  // Check if the given model info file exists or not:
+  std::ifstream ifs(modelInfo_file); 
+  if (!ifs) 
+  {
+    if (!PCU_Comm_Self()) 
+      std::cout<<"[M3D-C1 ERROR] "<<__func__<<" failed: model file \""<<modelInfo_file<<"\" doesn't exist\n";
+    return M3DC1_FAILURE; 
+  }
+
+  // Load the content of the model info file
+  // First read the bounding box
+  double r_min, z_min, r_max, z_max;
+  ifs >> r_min >> z_min >> r_max >> z_max;
+  m3dc1_model::instance()->boundingBox[0]= r_min;
+  m3dc1_model::instance()->boundingBox[1]= z_min;
+  m3dc1_model::instance()->boundingBox[2]= r_max;
+  m3dc1_model::instance()->boundingBox[3]= z_max;
+
+  // Read total number of edges on the model
+  int numModelEdges = 0;
+  ifs >> numModelEdges;
+  m3dc1_model::instance()->numModelEdges = numModelEdges;
+
+  // Read innermost loop
+  int numEdges = 0;
+  ifs >> numEdges;
+  for (int i = 0; i < numEdges; i++)
+  {
+    int edgeId;
+    ifs >> edgeId;
+    m3dc1_model::instance()->innerLoop.push_back(edgeId);
+  } 
+  
+  // Read outermost loop
+  ifs >> numEdges;
+  for (int i = 0; i < numEdges; i++)
+  {
+    int edgeId;
+    ifs >> edgeId;
+    m3dc1_model::instance()->outerLoop.push_back(edgeId);
+  } 
+}
+
+// ****************************************************
+void m3dc1_model_getnumedges(int* numedges)
+// ****************************************************
+{
+  *numedges = m3dc1_model::instance()->numModelEdges;
+}
+
+// ****************************************************
+void m3dc1_model_getgeometricloop(int* modelEdges, int* numberOfEdges, int* loopType)
+// ****************************************************
+{
+  std::vector <int> geomEdges;
+  if (*loopType == 0)
+    geomEdges = m3dc1_model::instance()->innerLoop;
+  else if (*loopType == 1)
+    geomEdges = m3dc1_model::instance()->outerLoop;
+  else  // might add further types in future
+  {
+    std::cout << "Incorrect Loop Type. Set 0 for inner loop, and 1 for outer loop\n";
+    exit(1);
+  }
+
+  // Print a warning if the requested loop does not exist.
+  if (geomEdges.size() == 0)
+    std::cout << "Warning: The requested loop has zero edges - Loop does not exist in the model\n";
+
+  for (int i = 0; i <  geomEdges.size(); i++)
+    modelEdges[i] = geomEdges[i];
+
+  *numberOfEdges = geomEdges.size();
+}
 
 //*******************************************************
 int m3dc1_model_print()
@@ -1635,6 +1748,36 @@ int m3dc1_ent_getgeomclass (int* /* in */ ent_dim, int* /* in */ ent_id,
   return M3DC1_SUCCESS;
 }
 
+//*******************************************************
+int m3dc1_set_nodes_orientation(int* nodes_array)
+//*******************************************************
+{
+  double node0[3], node1[3], node2[3];
+  m3dc1_node_getcoord(&nodes_array[0], node0);
+  m3dc1_node_getcoord(&nodes_array[1], node1);
+  m3dc1_node_getcoord(&nodes_array[2], node2);
+
+  // Same logic as used in FORTRAN side of the code for consistency
+  double x2 = node1[0] - node0[0];
+  double z2 = node1[1] - node0[1];
+  double x3 = node2[0] - node0[0];
+  double z3 = node2[1] - node0[1];
+
+  double hi = 1.0/sqrt(x2*x2 + z2*z2);
+  double co = x2*hi;
+  double sn = z2*hi;
+  
+  double z3p = -sn*x3 + co*z3;
+
+  if (z3p <= 0.0)
+  {
+    int temp = nodes_array[1];
+    nodes_array[1] = nodes_array[2];
+    nodes_array[2] = temp;
+  }
+  return M3DC1_SUCCESS;
+}
+
 //*******************************************************
 int m3dc1_ent_getadj (int* /* in */ ent_dim, int* /* in */ ent_id, 
                       int* /* in */ adj_dim, int* /* out */ adj_ent, 
@@ -1648,7 +1791,6 @@ int m3dc1_ent_getadj (int* /* in */ ent_dim, int* /* in */ ent_id,
     print_ent_error(__func__, *ent_dim, *ent_id);
 #endif
   assert(e);
-
   if (*adj_dim>*ent_dim) // upward
   {
     apf::Adjacent adjacent;
@@ -1677,6 +1819,11 @@ int m3dc1_ent_getadj (int* /* in */ ent_dim, int* /* in */ ent_id,
     }
     for (int i=0; i<*adj_ent_size; ++i)
       adj_ent[i] = getMdsIndex(mesh, downward[i]);
+ 
+    // Check nodes ordering and fix it if needed
+    if (*ent_dim == 2 && *adj_dim == 0)
+      m3dc1_set_nodes_orientation(adj_ent);
+
     //adjust the order to work with m3dc1
     if (mesh->getDimension()==3 && *ent_dim==3 &&*adj_dim==0 &&adj_ent[0]>adj_ent[3])
     {
@@ -1912,7 +2059,16 @@ int m3dc1_node_getcoord (int* /* in */ node_id, double* /* out */ coord)
   apf::Vector3 xyz;
   m3dc1_mesh::instance()->mesh->getPoint(e, 0, xyz);
   for (int i=0; i<3; ++i)
-    coord[i] = xyz[i]; 
+    coord[i] = xyz[i];
+
+  // Check the plane(If RZ, convert them to XY)
+  double tolerance = 1e-8;
+  if (fabs(coord[1]) <  tolerance)  // Its in RZ plane. change to XY
+  {
+    double temp = coord[1];
+    coord[1] = coord[2];
+    coord[2] = temp;
+  }
   return M3DC1_SUCCESS;
 }
 
@@ -1935,7 +2091,7 @@ int m3dc1_node_getnormvec (int* /* in */ node_id, double* /* out */ xyzt)
   apf::MeshEntity* vt = getMdsEntity(m3dc1_mesh::instance()->mesh, 0, *node_id);
   assert(vt);
   xyzt[2]=0.0;
-  //cout<<"nodnormalvec_ "<<*iNode<<" "<<vt<<endl;
+
   gmi_ent* gent= (gmi_ent*)(m3dc1_mesh::instance()->mesh->toModel(vt));
   int gType = gmi_dim(m3dc1_model::instance()->model,gent);
   if (gType !=  1 && gType !=  0)
@@ -1944,6 +2100,30 @@ int m3dc1_node_getnormvec (int* /* in */ node_id, double* /* out */ xyzt)
     return M3DC1_SUCCESS;
   }
 
+  if (m3dc1_model::instance()->modelType == 2)
+  {
+    apf::MeshTag* norm_curv_tag = m3dc1_mesh::instance()->mesh->findTag("norm_curv");
+    if (norm_curv_tag && m3dc1_mesh::instance()->mesh->hasTag(vt, norm_curv_tag))
+    {
+      double norm_curv[3];
+      m3dc1_mesh::instance()->mesh->getDoubleTag(vt, norm_curv_tag, &norm_curv[0]);
+      xyzt[0] = norm_curv[0]; 
+      xyzt[1] = norm_curv[1];
+      return M3DC1_SUCCESS;
+    }
+    else 
+    {/*
+      std::cout << "Warning: In case of .dmg models, normal vectors are only provided at \n";
+      std::cout << "predefined boundaries while creating mesh, make sure this function is \n";
+      std::cout << "called at the correct model edges at correct boundary loop\n";*/
+
+      // Dummy to return normals so code doesn't crash in case 
+      xyzt[0] = 1.0;
+      xyzt[1] = 0.0;
+      return M3DC1_SUCCESS;
+    }
+  }
+
   apf::MeshTag* norm_curv_tag = m3dc1_mesh::instance()->mesh->findTag("norm_curv");
   if (norm_curv_tag && m3dc1_mesh::instance()->mesh->hasTag(vt, norm_curv_tag))
   {
@@ -1970,7 +2150,6 @@ int m3dc1_node_getnormvec (int* /* in */ node_id, double* /* out */ xyzt)
       int numEdgePlane=0;
       double normalvec[3]={0.,0.,0.};
       xyzt[0]=xyzt[1]=xyzt[2]=0;
-      if (gEdges.size()<2) 
         std::cout<<"["<<PCU_Comm_Self()<<"] "<<__func__<<" ERROR: #adjEdge of gVertex="
 		 <<gEdges.size()<<" (it should be minimum 2) \n";
       assert(gEdges.size()>=2);
@@ -2022,6 +2201,29 @@ int m3dc1_node_getcurv (int* /* in */ node_id, double* /* out */ curv)
   apf::MeshEntity* vt = getMdsEntity(m3dc1_mesh::instance()->mesh, 0, *node_id);
   assert(vt);
 
+  if (m3dc1_model::instance()->modelType == 2)
+  {
+    apf::MeshTag* norm_curv_tag = m3dc1_mesh::instance()->mesh->findTag("norm_curv");
+    if (norm_curv_tag && m3dc1_mesh::instance()->mesh->hasTag(vt, norm_curv_tag))
+    {
+      double norm_curv[3];
+      m3dc1_mesh::instance()->mesh->getDoubleTag(vt, norm_curv_tag, &norm_curv[0]);
+      *curv = norm_curv[2];
+      return M3DC1_SUCCESS;
+    }
+    else 
+    {/*
+      std::cout << "Warning: In case of .dmg models, curvatures are only provided at \n";
+      std::cout << "predefined boundaries while creating mesh, make sure this function is \n";
+      std::cout << "called at the correct model edges at correct boundary loop\n";*/
+
+      // Dummy to return normals so code doesn't crash in case 
+      *curv = 0.0;
+      return M3DC1_SUCCESS;
+    }
+  }
+  
+
   apf::MeshTag* norm_curv_tag = m3dc1_mesh::instance()->mesh->findTag("norm_curv");
   if (norm_curv_tag && m3dc1_mesh::instance()->mesh->hasTag(vt, norm_curv_tag))
   {
@@ -5395,8 +5597,6 @@ void m3dc1_node_getNormVecOnNewVert(apf::MeshEntity* v, double* normalVec)
   
   bool v0_found = false;
   bool v1_found = false;
-  apf::Vector3 param(0,0,0);
-  m->getParam(v,param);
   apf::Adjacent edges;
   m->getAdjacent(v, 1, edges);
   std::vector <apf::MeshEntity*> adjVert;
diff --git a/m3dc1_scorec/api/m3dc1_scorec.h b/m3dc1_scorec/api/m3dc1_scorec.h
index 047c2a763..1f122397b 100644
--- a/m3dc1_scorec/api/m3dc1_scorec.h
+++ b/m3dc1_scorec/api/m3dc1_scorec.h
@@ -60,6 +60,10 @@ int m3dc1_plane_getphi(int* planeid, double* phi);
 
 /** model functions */
 int m3dc1_model_load(char* /* in */ model_file);
+int m3dc1_modelinfo_load(char* /* in */ modelInfo_file);
+void m3dc1_model_getnumedges(int* numedges);
+void m3dc1_model_getgeometricloop(int* modelEdges, int* numberOfEdges, int* loopType);
+void m3dc1_model_getmodeltype(int* modeltype);
 int m3dc1_model_print();
 int m3dc1_model_setnumplane(int*);
 int m3dc1_model_getnumplane(int*);
diff --git a/m3dc1_scorec/api/name_convert.h b/m3dc1_scorec/api/name_convert.h
index f93ab2ad1..67dcce3f9 100644
--- a/m3dc1_scorec/api/name_convert.h
+++ b/m3dc1_scorec/api/name_convert.h
@@ -8,6 +8,10 @@
 #define m3dc1_plane_setphi m3dc1_plane_setphi_
 #define m3dc1_plane_getphi m3dc1_plane_getphi_
 #define m3dc1_model_load m3dc1_model_load_
+#define m3dc1_modelinfo_load m3dc1_modelinfo_load_
+#define m3dc1_model_getnumedges m3dc1_model_getnumedges_
+#define m3dc1_model_getgeometricloop m3dc1_model_getgeometricloop_
+#define m3dc1_model_getmodeltype m3dc1_model_getmodeltype_
 #define m3dc1_model_print m3dc1_model_print_
 #define m3dc1_model_setnumplane m3dc1_model_setnumplane_
 #define m3dc1_model_getnumplane m3dc1_model_getnumplane_
diff --git a/m3dc1_scorec/include/m3dc1_model.h b/m3dc1_scorec/include/m3dc1_model.h
index 8354e9a64..4a5b91a4f 100644
--- a/m3dc1_scorec/include/m3dc1_model.h
+++ b/m3dc1_scorec/include/m3dc1_model.h
@@ -12,6 +12,7 @@
 #include "apf.h"
 #include "apfMesh.h"
 #include "gmi.h"
+#include "gmi_mesh.h"
 #include <mpi.h>
 #include <map>
 #include <vector>
@@ -100,6 +101,10 @@ class m3dc1_model
   int next_plane_partid; // id of corresponding part in the next plane
   bool snapping; // support for snapping
 
+  int modelType = 1;  // = 1 for analytical, = 2 for PUMI (.dmg) model
+  std::vector <int> innerLoop;
+  std::vector <int> outerLoop;
+  int numModelEdges = 0;
   double* phi;
   int numEntOrig[3];
   double boundingBox[4];
diff --git a/m3dc1_scorec/src/m3dc1_matrix.cc b/m3dc1_scorec/src/m3dc1_matrix.cc
index 91ee87e73..662c11e54 100644
--- a/m3dc1_scorec/src/m3dc1_matrix.cc
+++ b/m3dc1_scorec/src/m3dc1_matrix.cc
@@ -126,7 +126,7 @@ int copyField2PetscVec_5(FieldID field_id, Vec petscVec, int scalar_type)
   assert(nodeCounter==num_own_ent);
   ierr=VecAssemblyEnd(petscVec);
   CHKERRQ(ierr);
-  PetscFunctionReturn(PETSC_SUCCESS);
+  PetscFunctionReturn(0);  // PETSC_SUCCESS (0) to indicate success
 //  return 0;
 }
 
diff --git a/m3dc1_scorec/src/m3dc1_model.cc b/m3dc1_scorec/src/m3dc1_model.cc
index ad3ab7e4f..ccc66c8c7 100644
--- a/m3dc1_scorec/src/m3dc1_model.cc
+++ b/m3dc1_scorec/src/m3dc1_model.cc
@@ -237,6 +237,16 @@ void load_model(const char* filename)
   int numL,separatrixLoop, innerWallLoop, outerWallLoop, vacuumLoop;
   fscanf(fp,"%d %d %d %d %d\n", &numL, &separatrixLoop, &innerWallLoop, &outerWallLoop, &vacuumLoop);
 
+  // Identify outermost loop. We will always have inner loop even if numL==1.
+  int outerMostLoop = -1;
+  if (outerWallLoop > 0 || vacuumLoop > 0)
+  {
+    if (vacuumLoop > 0)
+      outerMostLoop = vacuumLoop;
+    else
+      outerMostLoop = outerWallLoop;
+  }
+
   std::vector<int> loop_ids;
   loop_ids.resize(numL);
 
@@ -246,6 +256,7 @@ void load_model(const char* filename)
     int loop; // loop ID
     fscanf(fp,"%d %d\n", &loop, &numE);
     loop_ids[i] = loop;
+    m3dc1_model::instance()->numModelEdges += numE;
     // first read all vtx on the loop
     for( int j=0; j<numE; j++)
     {
@@ -259,7 +270,12 @@ void load_model(const char* filename)
     {
       int edge, beginvtx, endvtx,edgeType;
       fscanf(fp,"%d %d %d %d\n", &edge, &beginvtx, &endvtx, &edgeType);
-            if (edgeContainer.find(edge)!=edgeContainer.end()) // edge already created
+      if (loop == innerWallLoop)
+        m3dc1_model::instance()->innerLoop.push_back(edge);
+      if (loop == outerMostLoop)
+        m3dc1_model::instance()->outerLoop.push_back(edge);
+
+      if (edgeContainer.find(edge)!=edgeContainer.end()) // edge already created
       {
         edges[i]=edge;
         switch (edgeType)
diff --git a/unstructured/input.f90 b/unstructured/input.f90
index e4129ce90..a6fceb790 100644
--- a/unstructured/input.f90
+++ b/unstructured/input.f90
@@ -1200,6 +1200,8 @@ subroutine set_defaults
        "", mesh_grp)
   call add_var_string("mesh_model", mesh_model, 256, "struct.dmg", &
        "", mesh_grp)
+  call add_var_string("model_info", model_info, 256, "dummyInfo", &
+        "", mesh_grp)
   call add_var_int("ipartitioned",ipartitioned,0,&
        "1 = the input mesh is partitioned", mesh_grp)
   call add_var_int("imatassemble", imatassemble, 0, &
diff --git a/unstructured/output.f90 b/unstructured/output.f90
index 50ebd4dcb..ff2c29c19 100644
--- a/unstructured/output.f90
+++ b/unstructured/output.f90
@@ -798,9 +798,9 @@ subroutine output_mesh(time_group_id, nelms, error)
      if(iadapt.eq.0) call boundary_edge(i, is_edge, normal, idim)
 
      bound = 0.
-     if(is_edge(1).ne.0) bound = bound + 1. + (is_edge(1)-1)*2**3
-     if(is_edge(2).ne.0) bound = bound + 2. + (is_edge(2)-1)*2**7
-     if(is_edge(3).ne.0) bound = bound + 4. + (is_edge(3)-1)*2**11
+     if(is_edge(1).ne.0) bound = bound + 1. + boundary_type(is_edge(1))*2**3
+     if(is_edge(2).ne.0) bound = bound + 2. + boundary_type(is_edge(2))*2**7
+     if(is_edge(3).ne.0) bound = bound + 4. + boundary_type(is_edge(3))*2**11
 
      call get_element_data(i, d)
 
diff --git a/unstructured/scorec_mesh.f90 b/unstructured/scorec_mesh.f90
index 6bf808a4d..ef9fa5756 100644
--- a/unstructured/scorec_mesh.f90
+++ b/unstructured/scorec_mesh.f90
@@ -25,6 +25,7 @@ module scorec_mesh_mod
 
   character(len=256) mesh_model
   character(len=256) mesh_filename
+  character(len=256) model_info
   character(len=256) name_buff
   integer :: ipartitioned
   integer :: imatassemble  
@@ -48,8 +49,8 @@ module scorec_mesh_mod
   integer, parameter :: BOUND_FIRSTWALL = 1
   integer, parameter :: BOUND_DOMAIN = 2
 
-  integer, parameter :: max_bounds = 20
-  integer, parameter :: max_zones = 20
+  integer, parameter :: max_bounds = 1000
+  integer, parameter :: max_zones = 100
   integer :: boundary_type(max_bounds)
   integer :: zone_type(max_zones)
 
@@ -80,6 +81,8 @@ subroutine load_mesh
     real, allocatable :: xvals(:)
     integer :: nvals
 #endif
+    integer :: model_type, nedges, i
+    integer, dimension(max_bounds) :: edges
 
     ! load mesh
     call MPI_Comm_size(MPI_COMM_WORLD,maxrank,ier)
@@ -211,11 +214,57 @@ subroutine load_mesh
     ! to get N-part distributed mesh, run split_smb provided as mesh utilities
     write(name_buff,"(A,A)")  mesh_model(1:len_trim(mesh_model)),0
     call m3dc1_model_load(name_buff)
+    write(name_buff,"(A,A)")  model_info(1:len_trim(model_info)),0
+    call m3dc1_modelinfo_load(name_buff)
     write(name_buff,"(A,A)")  mesh_filename(1:len_trim(mesh_filename)),0
     call m3dc1_mesh_load (name_buff)
 #endif
     call update_nodes_owned
 
+    call m3dc1_model_getmodeltype(model_type)
+
+    ! For *.dmg model (model_type == 2),
+    ! find edges associated with first wall and computational domain boundary
+    if(model_type.eq.2) then
+       boundary_type(:) = BOUND_UNKNOWN
+
+       ! Find edges associated with first wall
+       call m3dc1_model_getgeometricloop(edges, nedges, 0)
+       if(myrank.eq.0) &
+            print *, 'nedges for first wall: ', nedges
+       if(nedges.gt.max_bounds) then
+          if(myrank.eq.0) &
+               print *, 'Error: nedges > max_bounds for first wall'
+          call safestop(184)
+       end if
+       do i=1, nedges
+          if(edges(i).gt.max_bounds) then
+             if(myrank.eq.0) &
+                  print *, 'Firstwall edge greater than max_bounds', edges(i)
+             call safestop(185)
+          end if
+          boundary_type(edges(i)) = BOUND_FIRSTWALL
+       end do
+
+       ! Find edges associated with computational domain boundary
+       call m3dc1_model_getgeometricloop(edges, nedges, 1)
+       if(myrank.eq.0) &
+            print *, 'nedges for domain boundary: ', nedges
+       if(nedges.gt.max_bounds) then
+          if(myrank.eq.0) &
+               print *, 'Error: nedges > max_bounds for domain boundary'
+          call safestop(186)
+       end if
+       do i=1, nedges
+          if(edges(i).gt.max_bounds) then
+             if(myrank.eq.0) &
+                  print *, 'Domain edge greater than max_bounds', edges(i)
+             call safestop(187)
+          end if
+          boundary_type(edges(i)) = BOUND_DOMAIN
+       end do
+    end if
+    
     initialized = .true.
   end subroutine load_mesh