diff --git a/libsrc/geom2d/python_geom2d.cpp b/libsrc/geom2d/python_geom2d.cpp
index e2089d3e..54f28ef1 100644
--- a/libsrc/geom2d/python_geom2d.cpp
+++ b/libsrc/geom2d/python_geom2d.cpp
@@ -63,7 +63,7 @@ DLL_HEADER void ExportGeom2d(py::module &m)
          py::arg("x"), py::arg("y"), py::arg("maxh") = 1e99, py::arg("hpref")=0, py::arg("name")="")
     .def("Append", FunctionPointer([](SplineGeometry2d &self, py::list segment, int leftdomain, int rightdomain,
                                       optional<variant<int, string>> bc, optional<int> copy, double maxh,
-                                      double hpref)
+                                      double hpref, double hprefleft, double hprefright)
 	  {
             auto segtype = py::cast<std::string>(segment[0]);
             
@@ -86,8 +86,8 @@ DLL_HEADER void ExportGeom2d(py::module &m)
             seg->leftdom = leftdomain;
             seg->rightdom = rightdomain;
             seg->hmax = maxh;
-            seg->hpref_left = hpref;
-            seg->hpref_right = hpref;
+            seg->hpref_left = max(hpref, hprefleft);
+            seg->hpref_right = max(hpref,hprefright);
             seg->reffak = 1;
             seg->copyfrom = -1;
             if (copy.has_value())
@@ -110,7 +110,7 @@ DLL_HEADER void ExportGeom2d(py::module &m)
             return self.GetNSplines()-1;
 	  }), py::arg("point_indices"), py::arg("leftdomain") = 1, py::arg("rightdomain") = py::int_(0),
          py::arg("bc")=nullopt, py::arg("copy")=nullopt, py::arg("maxh")=1e99,
-         py::arg("hpref")=0)
+         py::arg("hpref")=0,py::arg("hprefleft")=0,py::arg("hprefright")=0)
 
     
     .def("AppendSegment", FunctionPointer([](SplineGeometry2d &self, py::list point_indices, int leftdomain, int rightdomain)
diff --git a/libsrc/meshing/classifyhpel.hpp b/libsrc/meshing/classifyhpel.hpp
index 9eaea354..259a7cd0 100644
--- a/libsrc/meshing/classifyhpel.hpp
+++ b/libsrc/meshing/classifyhpel.hpp
@@ -664,7 +664,7 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edge
 				INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int dim, const FaceDescriptor & fd)
 
 {
-  HPREF_ELEMENT_TYPE type = HP_NONE; 
+  HPREF_ELEMENT_TYPE type = HP_NONE;
   
   int pnums[3]; 
   int p[3];   
@@ -761,11 +761,10 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edge
 	      int ep1=p[eledges[k][0]-1];  
 	      int ep2=p[eledges[k][1]-1];  
 	     
-	      INDEX_2 i2(el.PNum(ep1),el.PNum(ep2));  
+	      INDEX_2 i2 = INDEX_2::Sort(el.PNum(ep1),el.PNum(ep2));
 	     
 	      if(edges.Used(i2)) 
 		{
-		  
 		  if(edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[ep1-1])) || 
 		     edgepoint_dom.Used(INDEX_2(-1,pnums[ep1-1])) || 
 		     edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[ep2-1])) || 
@@ -783,10 +782,10 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edge
      
 	 
       for(int k=0;k<3;k++) 
-	if(edgepoint.Test(pnums[k])) //edgepoint, but not member of sing_edge on trig -> cp 
+	if(edgepoint.Test(pnums[k]) && (edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[k])) || edgepoint_dom.Used(INDEX_2(-1,pnums[k])))) //edgepoint, but not member of sing_edge on trig -> cp 
 	  {
 	    INDEX_2 i2a=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+1)%3])); 
-	    INDEX_2 i2b=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+2)%3])); 
+            INDEX_2 i2b=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+2)%3]));
 	    
 	    if(!edges.Used(i2a) && !edges.Used(i2b)) 
 	      point_sing[p[k]-1] = 3; 	
@@ -794,7 +793,7 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edge
       
       for(int k=0;k<3;k++) 
 	if(cornerpoint.Test(el.PNum(p[k]))) 
-	  point_sing[p[k]-1] = 3; 
+	  point_sing[p[k]-1] = 3;
       
       *testout << "point_sing = " << point_sing[0] << point_sing[1] << point_sing[2] << endl;
 
diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp
index 8d046c01..3aef7ebd 100644
--- a/libsrc/meshing/hprefinement.cpp
+++ b/libsrc/meshing/hprefinement.cpp
@@ -1560,7 +1560,7 @@ namespace netgen
 bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
 		       NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
 			INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int & levels, int & act_ref)
-{ 
+{
   bool sing = 0; 
   if (mesh.GetDimension() == 3)
     {
@@ -1722,7 +1722,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
 	    
 	    if (seg.singedge_left * levels >= act_ref)
 	      {
-		INDEX_2 i2 (mesh.LineSegment(i)[0], 
+		INDEX_2 i2 = INDEX_2::Sort(mesh.LineSegment(i)[0], 
 			    mesh.LineSegment(i)[1]);
 		edges.Set(i2,1); 
 		edgepoint.Set(i2.I1());
@@ -1738,7 +1738,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
 	    
 	    if (seg.singedge_right * levels >= act_ref)
 	      {
-		INDEX_2 i2 (mesh.LineSegment(i)[1], 
+		INDEX_2 i2 = INDEX_2::Sort(mesh.LineSegment(i)[1], 
 			    mesh.LineSegment(i)[0]);  
 		edges.Set (i2, 1);
 		edgepoint.Set(i2.I1());