gmsh curvature bifurcation benchmark

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;

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:

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);
   

benchmarks/curvature/Bifurcation/Bifurcation3D.geo

+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}

benchmarks/curvature/Bifurcation/diffuseCurvature.py

+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")

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"