diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index 3faff3c5d87d1942197f533ef5ae493a49b919aa..83689d1e3f318341b3c94db34249a534de290e6b 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -100,10 +100,14 @@ Curvature& Curvature::getInstance()
        std::list<GFace*>::iterator surfIterator;
 
        for(surfIterator = tempcompounds.begin(); surfIterator != tempcompounds.end(); ++surfIterator) {
-         _ptFinalEntityList.push_back(*surfIterator);
+         if ((*surfIterator)->geomType() == GEntity::DiscreteSurface) {
+           _ptFinalEntityList.push_back(*surfIterator);
+         }
        }
      }
-
+     else if (entities[ie]->geomType() == GEntity::DiscreteSurface) {
+         _ptFinalEntityList.push_back(dynamic_cast<GFace*>(entities[ie]));
+     }
    }
 #endif
  }
@@ -985,6 +989,15 @@ void Curvature::edgeNodalValues(MLine* edge, double& c0, double& c1, int isAbs)
 
 //========================================================================================================
 
+double Curvature::getAtVertex(const MVertex *v) const {
+  std::map<int,int>::const_iterator it = _VertexToInt.find(v->getNum());
+  if (it == _VertexToInt.end()) {
+    Msg::Error("curvature has not been computed for vertex %i (%i)", v->getNum(), _VertexToInt.size());
+    return 1;
+  }
+  return _VertexCurve[it->second];
+}
+
 void Curvature::writeToPosFile( const std::string & filename)
 {
   std::ofstream outfile;
diff --git a/Geo/Curvature.h b/Geo/Curvature.h
index 60281189580fc2ffb5193fcc4f335eb98a1f6e62..07a024a412f4d3a3ec2d653a787efd084e1d7677 100644
--- a/Geo/Curvature.h
+++ b/Geo/Curvature.h
@@ -163,6 +163,7 @@ public:
   
   void setGModel(GModel* model);
   void retrieveCompounds();
+  double getAtVertex(const MVertex *v) const;
   //void retrievePhysicalSurfaces(const std::string & face_tag);
 
   /// The following function implements algorithm from:
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 15560e39979da97cb846065c661895158091c34b..a495a7027741e7c98e6533b750857c0cbe09f81d 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -245,7 +245,6 @@ SPoint2 GEdge::reparamOnFace(const GFace *face, double epar,int dir) const
 
 double GEdge::curvature(double par) const
 {
-  printf("in curv edge \n");
   SVector3 d1 = firstDer(par);
   SVector3 d2 = secondDer(par);
   
diff --git a/benchmarks/curvature/Bifurcation/Bifurcation3D.geo b/benchmarks/curvature/Bifurcation/Bifurcation3D.geo
new file mode 100644
index 0000000000000000000000000000000000000000..eae7c2776f8a236a6b44b89e109ea7fbad4be6ca
--- /dev/null
+++ b/benchmarks/curvature/Bifurcation/Bifurcation3D.geo
@@ -0,0 +1,14 @@
+Merge "BifurcationRemeshCurvatureReclassified.msh";
+CreateTopology;
+Line Loop(100) = {3};
+Plane Surface(200) = {100};
+Line Loop(101) = {4};
+Plane Surface(201) = {101};
+Line Loop(102) = {5};
+Plane Surface(202) = {102};
+
+Surface Loop(30) = {21, 200, 201, 202};
+Volume(1) = {30};
+Physical Volume ("Volume") = {1};
+Physical Surface ("Wall") = {21};
+Physical Surface ("Side") = {200, 201, 202}
diff --git a/benchmarks/curvature/Bifurcation/diffuseCurvature.py b/benchmarks/curvature/Bifurcation/diffuseCurvature.py
new file mode 100644
index 0000000000000000000000000000000000000000..6159f4a75b9754062e6a36e775e06143dda20be1
--- /dev/null
+++ b/benchmarks/curvature/Bifurcation/diffuseCurvature.py
@@ -0,0 +1,49 @@
+from dgpy import *
+
+gmodel = GModel()
+gmodel.load("Bifurcation3D.msh")
+
+curvature = Curvature.getInstance()
+curvature.setGModel(gmodel)
+curvature.computeCurvature_Rusinkiewicz()
+
+groups = dgGroupCollection(gmodel, 3, 1)
+
+curvDof = dgDofContainer(groups, 1)
+curvDof.setAll(10)
+
+groups.buildGroupsOfInterfaces()
+
+for iFG in range(groups.getNbFaceGroups()):
+  fGroup = groups.getFaceGroup(iFG)
+  cGroup = fGroup.getGroupOfConnections(0)
+  if (fGroup.getPhysicalTag() == "Wall"):
+    proxy = curvDof.getGroupProxy(cGroup.getGroupOfElements())
+    for iFace in range (fGroup.getNbElements()):
+      face = fGroup.getFace(iFace)
+      elementId = cGroup.getElementId(iFace)
+      element = cGroup.getElement(iFace)
+      closure = cGroup.getClosure(iFace)
+      for iVertex in range (face.getNumVertices()):
+        vertex = face.getVertex(iVertex)
+        curv = curvature.getAtVertex(vertex)
+        proxy.set(closure[iVertex], elementId, curv)
+
+curvDof.exportMsh("curv3D")
+
+
+coefnu = functionConstant(1)
+law = dgConservationLawAdvectionDiffusion.diffusionLaw(coefnu)
+law.addBoundaryCondition('Wall',law.newOutsideValueBoundary("",curvDof.getFunction()))
+law.addBoundaryCondition('Side',law.new0FluxBoundary())
+law.addBoundaryCondition('none',law.new0FluxBoundary())
+
+solution = dgDofContainer(groups, 1)
+sys = linearSystemPETScBlockDouble()
+sys.setParameter("petscOptions", "-pc_type ilu")
+dof = dgDofManager.newDGBlock(groups, 1, sys)
+solver = dgDIRK(law, dof, 2)
+solver.getNewton().setVerb(10)
+solver.getNewton().setAtol(1e-10)
+solver.iterate(solution, 1 , 0)
+solution.exportMsh("curv3D")
diff --git a/gmshpy/gmshGeo.i b/gmshpy/gmshGeo.i
index 95556901efd1aa814c3fd75322e28159b95c062f..e49b8641b710db4e3ad0e8c7b08f7af02b5b34e8 100644
--- a/gmshpy/gmshGeo.i
+++ b/gmshpy/gmshGeo.i
@@ -32,6 +32,7 @@
   #include "SPoint3.h"
   #include "SPoint2.h"
   #include "SBoundingBox3d.h"
+  #include "Curvature.h"
 %}
 
 namespace std {
@@ -72,3 +73,4 @@ namespace std {
 %include "SPoint3.h"
 %include "SPoint2.h"
 %include "SBoundingBox3d.h"
+%include "Curvature.h"