*** empty log message ***

Geo/GEdgeCompound.cpp

@@ -79,6 +79,7 @@ void GEdgeCompound::orderEdges()
   else{
     Msg::Error("EdgeCompound %d is wrong (it has %d end points)",tag(),tempv.size());
     Msg::Exit(1);
+
   }
 
   //loop over all segments to order segments and store it in the list _c

Geo/GFaceCompound.cpp

@@ -17,6 +17,10 @@
 #include "gmshLinearSystemFull.h"
 #include "FunctionSpace.h"
 #include "GmshDefines.h"
+#include "GmshPredicates.h"
+#include "meshGFaceOptimize.h"
+#include "MElementCut.h"
+#include "GEntity.h"
 
 class gmshGradientBasedDiffusivity : public gmshFunction<double>
 {
@@ -96,14 +100,160 @@ bool GFaceCompound::trivial() const {
   return false;
 }
 
+bool GFaceCompound::checkOrientation() const
+{
+  double U1, V1, U2,V2, U3,V3;
+  std::list<GFace*>::const_iterator it = _compound.begin();
+  double a_old, a_new;
+  bool oriented = true;
+  int iter = 0;
+  for ( ; it != _compound.end(); ++it){
+    for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+      MTriangle *t = (*it)->triangles[i];
+      SPoint3 v1 = coordinates[t->getVertex(0)];
+      SPoint3 v2 = coordinates[t->getVertex(1)];
+      SPoint3 v3 = coordinates[t->getVertex(2)];
+      double p1[2] = {v1[0],v1[1]};
+      double p2[2] = {v2[0],v2[1]};
+      double p3[2] = {v3[0],v3[1]};
+      a_new = gmsh::orient2d(p1, p2, p3);
+      if (iter == 0) a_old=a_new;
+      if (a_new*a_old < 0.){
+	oriented = false;
+	break;
+      }
+      else{
+	a_old = a_new;
+      }
+      iter ++;
+    }    
+  }  
+
+  if(!oriented ){
+    Msg::Info("*** Could not compute reparametrization with the harmonic map.");
+    Msg::Info("*** Force One-2-One Map.");
+    one2OneMap();
+    checkOrientation();
+  }
+  else{
+    Msg::Info("Harmonic mapping successfully computed ...");
+  }
+
+  return oriented;
+
+}
+
+void GFaceCompound::one2OneMap() const{
+
+  if (!mapv2Tri){
+    std::vector<MTriangle*> allTri;
+    std::list<GFace*>::const_iterator it = _compound.begin();
+    for ( ; it != _compound.end(); ++it){
+      allTri.insert(allTri.end(), (*it)->triangles.begin(), (*it)->triangles.end() );
+    }
+    buildVertexToTriangle(allTri, adjv);
+    //printf("allTri size =%d\n", allTri.size());
+    mapv2Tri=true;
+  }
+
+ for (v2t_cont::iterator it = adjv.begin(); it!= adjv.end(); ++it){
+   MVertex *v = it->first;
+   std::vector<MElement*> vTri = it->second;
+   int nbV = vTri.size();
+   double a_old, a_new;
+   bool badCavity = false;
+
+   std::map<MVertex*,SPoint3>::iterator itV = coordinates.find(v);
+   //printf("*** Node = %d: (%g %g) cavity size = %d \n", v->getNum(), itV->second.x(), itV->second.y(), nbV);
+
+   std::set<MVertex*> cavV;
+
+   for (unsigned int i = 0; i < nbV; ++i){
+     MTriangle *t = (MTriangle*) vTri[i];
+     SPoint3 v1 = coordinates[t->getVertex(0)];
+     SPoint3 v2 = coordinates[t->getVertex(1)];
+     SPoint3 v3 = coordinates[t->getVertex(2)];
+     for (int k=0; k< 3; k++){
+        MVertex *vk = t->getVertex(k);
+        std::set<MVertex*>::const_iterator it = cavV.find(vk);
+        if( it == cavV.end() && v != vk ) {
+	  cavV.insert(vk);
+	  //if (v != vk) cavV.insert(vk);
+	  //if (v->onWhat()->dim() == 1 && v == vk)  cavV.insert(vk);
+	}
+     }
+    
+     double p1[2] = {v1[0],v1[1]};
+     double p2[2] = {v2[0],v2[1]};
+     double p3[2] = {v3[0],v3[1]};
+     //printf("p1=(%g %g) p2=(%g %g) p3=(%g %g)\n",v1[0],v1[1], v2[0],v2[1], v3[0],v3[1] );
+     a_new = gmsh::orient2d(p1, p2, p3);
+     if (i == 0) a_old=a_new;
+     if (a_new*a_old < 0.)   badCavity = true;
+     a_old = a_new;
+   }
+   if(badCavity){
+     Msg::Info("Wrong cavity place vertex (%d onwhat=%d) at center of gravity (sizeCav=%d).",  v->getNum(),  v->onWhat()->dim(), nbV);
+     
+     int nbVert = cavV.size();
+     //printf("nbVert =%d \n", nbVert);
+     double u_cg = 0.0;
+     double v_cg = 0.0;
+     for (std::set<MVertex*>::const_iterator it = cavV.begin() ; it != cavV.end() ; ++it){
+       SPoint3 vsp = coordinates[*it];
+       u_cg+=vsp[0];
+       v_cg+=vsp[1];
+     }
+     u_cg/=nbVert;
+     v_cg/=nbVert;
+     //printf("ucg=%g vcg=%g \n", u_cg, v_cg);
+     SPoint3 p_cg(u_cg,v_cg,0);
+     coordinates[v] = p_cg;
+
+     //  //NEW CAVITY:
+//       printf("//**** NEW CAVITY****** \n " );
+//       for (int i=0; i< nbV; i++){
+// 	MTriangle *t = (MTriangle*) vTri[i];
+// 	SPoint3 v1 = coordinates[t->getVertex(0)];
+// 	SPoint3 v2 = coordinates[t->getVertex(1)];
+// 	SPoint3 v3 = coordinates[t->getVertex(2)];
+
+// 	printf("//////////////////// \n " );
+// 	double p1[2] = {v1[0],v1[1]};
+// 	double p2[2] = {v2[0],v2[1]};
+// 	double p3[2] = {v3[0],v3[1]};
+// 	a_new = gmsh::orient2d(p1, p2, p3);
+// 	MVertex *e1=t->getVertex(0);	MVertex *e2=t->getVertex(1);	MVertex *e3=t->getVertex(2);
+// 	printf("Point(%d)={%g, %g, 0}; \n ",e1->getNum(), v1[0],v1[1] );
+// 	printf("Point(%d)={%g, %g, 0}; \n ",e2->getNum(), v2[0],v2[1] );
+// 	printf("Point(%d)={%g, %g, 0}; \n ",e3->getNum(), v3[0],v3[1] );
+// 	printf("Line(%d)={%d,%d}; \n", e1->getNum()+e2->getNum() , e1->getNum(), e2->getNum());
+// 	printf("Line(%d)={%d,%d}; \n", e2->getNum()+e3->getNum() , e2->getNum(), e3->getNum());
+// 	printf("Line(%d)={%d,%d}; \n", e3->getNum()+e1->getNum() , e3->getNum(), e1->getNum());
+// 	printf("Surface(%d)={%d,%d, %d}; \n", e3->getNum()+e1->getNum() + e2->getNum(), e3->getNum(), e1->getNum()+e2->getNum() , 
+// 	       e2->getNum()+e3->getNum() , e3->getNum()+e1->getNum() );
+// 	printf("//Area=%g \n", a_new);
+//       }
+//       //NEW CAVITY:
+//       exit(1);
+
+   }
+ }
+  
+}
+
 void GFaceCompound::parametrize() const
 {
   if (trivial())return;
 
   if (!oct){
     coordinates.clear(); 
+
     parametrize(ITERU);
     parametrize(ITERV);
+
+    checkOrientation();
+    
     computeNormals();
     buildOct();
   }
@@ -121,13 +271,11 @@ void GFaceCompound::getBoundingEdges()
     std::list<GEdge*> :: const_iterator it = _U0.begin();
     for ( ; it != _U0.end() ; ++it){
       l_edges.push_back(*it);
-      //printf("U0 for edge %d, add face %d \n", (*it)->tag(), this->tag());
       (*it)->addFace(this);
     }
     it = _V0.begin();
     for ( ; it != _V0.end() ; ++it){
       l_edges.push_back(*it);
-      //printf("V0 for edge %d, add face %d \n", (*it)->tag(), this->tag());
       (*it)->addFace(this);
     }
     return;
@@ -165,7 +313,6 @@ void GFaceCompound::getBoundingEdges()
     std::set<GEdge*>::iterator it = _unique.begin();
     GVertex *vB = (*it)->getBeginVertex();
     GVertex *vE = (*it)->getEndVertex();
-    //    printf("boundary add edge=%d (%d,%d)\n", (*it)->tag(),vB->tag(),vE->tag());
     _loop.push_back(*it);
     _unique.erase(it);
     it++;
@@ -232,7 +379,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
   else if (!_V1.size()) _type = UNITCIRCLE;
   else _type = SQUARE;
 
-
+  mapv2Tri = false;
 
 }
 
@@ -412,10 +559,9 @@ void GFaceCompound::parametrize(iterationStep step) const
   Msg::Info("Parametrizing Surface %d coordinate (ITER %d)", tag(), step); 
   
 #ifdef HAVE_GMM
-
   //gmshLinearSystemGmm<double> lsys;
   gmshLinearSystemCSRGmm<double> lsys;
-  lsys.setPrec(1.e-15);
+  lsys.setPrec(1.e-10);
   if(Msg::GetVerbosity() == 99) lsys.setNoisy(2);
 #else
   gmshLinearSystemFull<double> lsys;
@@ -437,25 +583,11 @@ void GFaceCompound::parametrize(iterationStep step) const
       return;
     }
 
-    //    printf("%d v on the boundary\n", ordered.size());
-    
     for (unsigned int i = 0; i < ordered.size(); i++){
       MVertex *v = ordered[i];
       const double theta = 2 * M_PI * coords[i];
       if (step == ITERU) myAssembler.fixVertex(v, 0, 1, cos(theta));
       else if (step == ITERV) myAssembler.fixVertex(v, 0, 1, sin(theta));    
-      else if (step == ITERD) myAssembler.fixVertex(v, 0, 1, 1.0);      
-    }
-    if (step == ITERD && _V0.size()) {
-      success = orderVertices(_V0, ordered1, coords1);
-      if (!success){
-	Msg::Error("Could not order vertices on boundary");
-	return;
-      }
-      for (unsigned int i = 0; i < ordered1.size(); i++){
-	MVertex *v = ordered1[i];
-	myAssembler.fixVertex(v, 0, 1, 0.0);      
-      }    
     }
   }
   else if (_type == SQUARE){
@@ -495,21 +627,8 @@ void GFaceCompound::parametrize(iterationStep step) const
   Msg::Debug("Creating term %d dofs numbered %d fixed",
              myAssembler.sizeOfR(), myAssembler.sizeOfF());
   
-  if (step == ITERD){
-    gmshLaplaceTerm laplace(model(), &ONE, 1);
-    it = _compound.begin();
-    for ( ; it != _compound.end() ; ++it){
-      for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
-	MTriangle *t = (*it)->triangles[i];
-	//	diffusivity.setCurrent(t);
-	laplace.addToMatrix(myAssembler, t);
-      }
-    }    
-  }
-  else{
   //gmshConvexCombinationTerm laplace(model(), &ONE, 1);
   gmshLaplaceTerm laplace(model(), &ONE, 1);
-    //gmshLaplaceTerm laplace(model(), &diffusivity, 1);
   it = _compound.begin();
   for ( ; it != _compound.end() ; ++it){
     for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
@@ -517,11 +636,12 @@ void GFaceCompound::parametrize(iterationStep step) const
       laplace.addToMatrix(myAssembler, t);
     }
   }
-  }
  
   Msg::Debug("Assembly done");
   lsys.systemSolve();
   Msg::Debug("System solved");
+
+
   it = _compound.begin();
   std::set<MVertex *>allNodes;
   for ( ; it != _compound.end() ; ++it){
@@ -536,8 +656,6 @@ void GFaceCompound::parametrize(iterationStep step) const
   for (std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
     MVertex *v = *itv;
     double value = myAssembler.getDofValue(v,0,1);
-    //if (step == 1)
-    //printf("%p node =%g %g %g , value = %g of size %d \n",v,  v->x(), v->y(), v->z(), value, coordinates.size());
     std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
     if (itf == coordinates.end()){
       SPoint3 p(0,0,0);

Geo/GFaceCompound.h

@@ -11,6 +11,7 @@
 #include "GFace.h"
 #include "GEdge.h"
 #include "GEdgeCompound.h"
+#include "meshGFaceOptimize.h"
 
 /*
 A GFaceCompound is a model face that is the compound of model faces.
@@ -52,11 +53,15 @@ class GFaceCompound : public GFace {
   mutable int nbT;
   mutable GFaceCompoundTriangle *_gfct;
   mutable Octree *oct;
+  mutable v2t_cont adjv;
+  mutable bool mapv2Tri;
   mutable std::map<MVertex*,SPoint3> coordinates;
   mutable std::map<MVertex*,SVector3> _normals;
   void buildOct() const ;
   void parametrize() const ;
   void parametrize(iterationStep) const ;
+  bool checkOrientation() const;
+  void one2OneMap() const;
   void computeNormals () const;
   void getBoundingEdges();
   void getTriangle(double u, double v, GFaceCompoundTriangle **lt, 
@@ -64,6 +69,7 @@ class GFaceCompound : public GFace {
   virtual double curvature(MTriangle *t, double u, double v) const;
   void printStuff() const;
   bool trivial() const ;
+
 public:
   typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
   GFaceCompound(GModel *m, int tag, 

Geo/MElementCut.h

@@ -163,6 +163,15 @@ class MPolygon : public MElement {
       _parts.push_back(new MTriangle(v[i * 3], v[i * 3 + 1], v[i * 3 + 2]));
     _initVertices();
   }
+  MPolygon(std::vector<MElement*> vT, int num=0, int part=0)
+    : MElement(num, part), _owner(false), _orig(0)
+  {
+    for(unsigned int i = 0; i < vT.size(); i++){
+      MTriangle *t = (MTriangle*) vT[i];
+      _parts.push_back(t);
+    }
+    _initVertices();
+  }
   MPolygon(std::vector<MTriangle*> vT, int num=0, int part=0)
     : MElement(num, part), _owner(false), _orig(0)
   {

Geo/discreteEdge.cpp

@@ -150,7 +150,7 @@ void discreteEdge::setBoundVertices()
 {
 
   if (boundv.size()==2){
-    printf("Found a regular open Curve \n");
+    //printf("Found a regular open Curve \n");
     std::vector<GVertex*> bound_vertices;
     for (std::map<MVertex*,MLine*>::const_iterator iter = boundv.begin(); iter != boundv.end(); iter++){
       MVertex* vE = (iter)->first;
@@ -447,55 +447,6 @@ GPoint discreteEdge::point(double par) const
     return GPoint(PP.x(),PP.y(),PP.z());
 
   }
-//   else{
-    
-//     //quadratic Lagrange mesh
-//     //-------------------------
-
-//     const SVector3 n1 = _normals[vB];
-//     const SVector3 n2 = _normals[vE];
-
-//     SPoint3 v1(vB->x(),vB->y(),vB->z());  
-//     SPoint3 v2(vE->x(),vE->y(),vE->z());
-  
-//     SVector3 t1, t2, Pij; 
-//     t1 = n2 - n1*dot(n1,n2);
-//     t2 =  n2*(dot(n1,n2)) - n1; 
-//     Pij = v2 - v1;
-    
-//     double a = dot(t1,t2)/dot(t1,t1);
-//     double b = dot(Pij,t1)/dot(t1,t2);
-//     double ap = dot(t1,t2)/dot(t2,t2);
-//     double bp = dot(Pij,t2)/dot(t1,t2);
-    
-//     SPoint3 X3;
-//     if (dot(n1,n2)/(norm(n1)*norm(n2)) > 0.9999){
-//       printf("coucou normals \n");
-//       X3.setPosition(.5*(v1.x()+v2.x()),  .5*(v1.y()+v2.y()), .5*(v1.z()+v2.z()) );
-//     }
-//     else{
-//       SVector3 XX3 = (.5*((ap*bp*t2)-(a*bp*t1)) ) *.25  +  (Pij + .5*((a*b*t1) - (ap*bp*t2)))*.5 +  v1;
-//       X3.setPosition(XX3.x(), XX3.y(), XX3.z());
-//     } 
-    
-//    //   printf("normal x1 num=%d  coord = %g %g %g \n" , vB->getNum(), v1.x(), v1.y(), v1.z());
-// //      printf("normal n1 = %g %g %g \n", n1.x(), n1.y(), n1.z());
-// //      printf("normal x2 num = %d coord = %g %g %g \n", vE->getNum(), v2.x(), v2.y(), v2.z());
-// //      printf("normal n2 = %g %g %g \n", n2.x(), n2.y(), n2.z());
-// //      printf("normal x3 = %g %g %g \n", X3.x(), X3.y(), X3.z());
-// //      printf("normal x3b = %g %g %g \n", X3b.x(), X3b.y(), X3b.z());
-//      //exit(1);
- 
-//     const gmshFunctionSpace& fs = gmshFunctionSpaces::find(MSH_LIN_3);
-//     double f1[3];
-//     SPoint3 p(0, 0, 0);
-//     double U = 2*tLoc -1;
-//     fs.f(U,0,0,f1);
-
-//     p = v1*f1[0] + v2*f1[1] + X3*f1[2];
-//     return GPoint(p.x(),p.y(),p.z());
-  
-//   }
 
 }
 

Numeric/gmshAssembler.h

@@ -9,6 +9,7 @@
 #include <map>
 #include <vector>
 #include "gmshLinearSystem.h"
+#include <stdlib.h>
 
 class MVertex;
 class MElement;
@@ -35,6 +36,7 @@ class gmshAssembler {
  private:
   std::map<gmshDofKey, int> numbering;
   std::map<gmshDofKey, scalar> fixed;
+  std::map<gmshDofKey, scalar> fixedFULL;
   std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, scalar> > > constraints;
   gmshLinearSystem<scalar> *lsys;
  public:
@@ -79,6 +81,10 @@ class gmshAssembler {
   {
     fixed[gmshDofKey(v, iComp, iField)] = val;
   }
+  inline void fixVertexFULL(MVertex*v, int iComp, int iField, scalar val)
+  {
+    fixedFULL[gmshDofKey(v, iComp, iField)] = val;
+  }
   inline scalar getDofValue(MVertex *v, int iComp, int iField) const
   {
     gmshDofKey key(v, iComp, iField);
@@ -117,13 +123,21 @@ class gmshAssembler {
     if (itR != numbering.end()){
       std::map<gmshDofKey, int>::iterator 
         itC = numbering.find(gmshDofKey(vC, iCompC, iFieldC));
-      if (itC != numbering.end()){
+      typename std::map<gmshDofKey, scalar>::iterator 
+	itFF = fixedFULL.find(gmshDofKey(vR, iCompR, iFieldR));
+      if (itC != numbering.end() && itFF ==  fixedFULL.end()){
         lsys->addToMatrix(itR->second, itC->second, val);
       }
+      else if (itFF != fixedFULL.end()){
+	//printf("RHS = %g,  ligne=%d \n",  itFF->second, itR->second);
+	lsys->addToMatrix(itR->second,itR->second, 1. );
+	lsys->addToRightHandSide(itR->second, itFF->second);
+      }
       else {
         typename std::map<gmshDofKey, scalar>::iterator 
           itF = fixed.find(gmshDofKey(vC, iCompC, iFieldC));
         if (itF != fixed.end()){
+	  //printf("RHS = val%g itF=%g \n", val, itF->second);
           lsys->addToRightHandSide(itR->second, -val*itF->second);
         }
         else{

Numeric/gmshConvexCombination.cpp

@@ -11,12 +11,12 @@ void gmshConvexCombinationTerm::elementMatrix(MElement *e, gmshMatrix<double> &m
  
  m.set_all(0.);
  int nbNodes = e->getNumVertices();
- double x = 0.3*(e->getVertex(0)->x()+e->getVertex(1)->x()+e->getVertex(2)->x());
- double y = 0.3*(e->getVertex(0)->y()+e->getVertex(1)->y()+e->getVertex(2)->y());
- double z = 0.3*(e->getVertex(0)->z()+e->getVertex(1)->z()+e->getVertex(2)->z());
- 
- const double _diff = 1.0; 
+ //double x = 0.3*(e->getVertex(0)->x()+e->getVertex(1)->x()+e->getVertex(2)->x());
+ //double y = 0.3*(e->getVertex(0)->y()+e->getVertex(1)->y()+e->getVertex(2)->y());
+ //double z = 0.3*(e->getVertex(0)->z()+e->getVertex(1)->z()+e->getVertex(2)->z());
  //const double _diff = (*_diffusivity)(x, y, z);
+ const double _diff = 1.0; 
+
 
  for (int j = 0; j < nbNodes; j++){
    for (int k = 0; k < nbNodes; k++){

Numeric/gmshLaplace.cpp

@@ -54,25 +54,6 @@ void gmshLaplaceTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
     for (int k = 0; k < j; k++)
       m(k, j) = m(j, k);
   
-  
-  //check for positive scheme
-//   //printf("ELEM\n");
-//   for (int j = 0; j < nbNodes; j++){
-//     double sum = m(j,0)+m(j,1)+m(j,2);
-//     if (sum < 0.0){
-//       printf("Sum LINE gmshLaplace Term NEG <0  %g %g %g\n", m(j,0),m(j,1),m(j,2));
-//       for (int k = 0; k < nbNodes; k++)	m(j,k) = -1.;
-//       m(j,j) = (nbNodes-1);
-  //   }
-//     else{
-//       printf("Sum POS\n");
-//     }
-//     //for (int k = 0; k < nbNodes; k++) {
-//     //printf("m(%d,%d)=%g ", j,k, m(j,k));
-//     //}
-//     //printf("\n");
-//     }
-
 }
 
 void gmshLaplaceTerm2DParametric::elementMatrix(MElement *e, gmshMatrix<double> &m) const

Numeric/gmshLinearSystemCSR.cpp

@@ -291,6 +291,7 @@ void sortColumns (int NbLines,
 template<>
 int gmshLinearSystemCSRGmm<double> :: systemSolve()
 {
+  if(!sorted)
     sortColumns(_b->size(),
 		CSRList_Nbr(a_),
 		(INDEX_TYPE *) ptr_->array,
@@ -299,10 +300,6 @@ int gmshLinearSystemCSRGmm<double> :: systemSolve()
 		(double*) a_->array);
   sorted = true;
 
-  //  for (int i=0;i<_b->size();i++)
-  //    printf("%d ",((INDEX_TYPE *) jptr_->array)[i]);
-  //  printf("\n");
-
   gmm::csr_matrix_ref<double*,INDEX_TYPE *,INDEX_TYPE *, 0>  ref((double*) a_->array, 
 								 (INDEX_TYPE *) ai_->array,
 								 (INDEX_TYPE *) jptr_->array,
@@ -319,11 +316,35 @@ int gmshLinearSystemCSRGmm<double> :: systemSolve()
 }
 #endif
 
+#if defined(HAVE_GMM)
+#include "gmm.h"
+template<>
+int gmshLinearSystemCSRGmm<double> :: checkSystem()
+{
+  if(!sorted)
+    sortColumns(_b->size(),
+		CSRList_Nbr(a_),
+		(INDEX_TYPE *) ptr_->array,
+		(INDEX_TYPE *) jptr_->array, 
+		(INDEX_TYPE *) ai_->array, 
+		(double*) a_->array);
+  sorted = true;
+
+  printf("Coucou check system \n");
+  for (int i=0;i<_b->size();i++)
+    printf("%d ",((INDEX_TYPE *) jptr_->array)[i]);
+  printf("\n");
+
+  return 1;
+}
+#endif
+
 #if defined(HAVE_TAUCSw)
 #include "taucs.h"
 template<>
 int gmshLinearSystemCSRTaucs<double> :: systemSolve()
 {
+  if(!sorted)
     sortColumns(_b->size(),
 		CSRList_Nbr(a_),
 		(INDEX_TYPE *) ptr_->array,

Numeric/gmshLinearSystemCSR.h

@@ -145,6 +145,15 @@ class gmshLinearSystemCSRGmm : public gmshLinearSystemCSR<scalar> {
     Msg::Error("Gmm++ is not available in this version of Gmsh");
   }
 #endif
+  virtual int checkSystem() 
+#if defined(HAVE_GMM)
+    ;
+#else
+  {
+    Msg::Error("Gmm++ is not available in this version of Gmsh");
+  }
+#endif
+
 };
 
 template <class scalar>

Numeric/gmshLinearSystemFull.h

@@ -11,6 +11,8 @@
 #include "GmshMessage.h"
 #include "gmshLinearSystem.h"
 #include "GmshMatrix.h"
+#include <stdlib.h>
+#include <set>
 
 template <class scalar>
 class gmshLinearSystemFull : public gmshLinearSystem<scalar> {
@@ -66,6 +68,42 @@ class gmshLinearSystemFull : public gmshLinearSystem<scalar> {
   virtual int systemSolve() 
   {
     _a->lu_solve(*_b, *_x);
+    // _x->print("********* mySol");
+    return 1;
+  }
+  virtual int checkSystem() 
+  {
+    //a->print("myMatrix");
+    //_b->print("myVector");
+
+    int nbNonConvex=0;
+    for(int i = 0; i < _b->size(); i++){
+
+      double diag = (*_a)(i, i);
+      double offDiag = 0.0;
+      bool convex = true;
+      std::set<int> Ni;
+      for(int j = 0; j < _b->size(); j++){
+	if ( (j !=i) &&  (*_a)(i, j)  > 0.0 ) convex = false;
+	if (  j !=i && (*_a)(i, j) != 0.0){
+	  offDiag += (*_a)(i, j);
+	  Ni.insert(j);
+	}
+      }
+      if (std::abs(offDiag+diag) > 1.e-10 && std::abs(offDiag) > 1.e-10 ) convex= false;
+     
+      if (convex == false){
+	nbNonConvex+=1;
+	printf("*** WARNING NON CONVEX LINE !!!!!\n");
+	printf("*** i=%d : diag=%g offDiag=%g diff=%g convex=%s size Ni=%d\n", i , diag, offDiag,std::abs(offDiag+diag), (convex)?"true":"false", Ni.size());
+	for (std::set<int>::iterator it = Ni.begin(); it != Ni.end(); ++it)   (*_a)(i, *it) = -1.0;
+	(*_a)(i, i) = Ni.size();
+      }
+    }
+
+    printf("nonConvex=%d total=%d ratio=%g\n ",nbNonConvex, _b->size(), nbNonConvex/_b->size());
+    //_a->print("myMatrix AFTER !!! ");
+
     return 1;
   }
 };

benchmarks/stl/PelvisHEMI_CLASS_GEO.geo

@@ -13,9 +13,10 @@ Compound Line(80)={8};
 Compound Line(90)={9};
 
 //Compound Surface(100)={2} Boundary {{}};
-Compound Surface(200)={3} Boundary {{},{},{},{}};
+//Compound Surface(200)={3} Boundary {{},{},{},{}};
+//Compound Surface(200)={3} Boundary {{8},{6,7},{},{}};
 //Compound Surface(300)={4} Boundary {{}};
-//Compound Surface(400)={5} Boundary {{}};
+Compound Surface(400)={5} Boundary {{}};
 
 //Compound Volume(5000)={1000};
 

utils/api_demos/Makefile

@@ -22,9 +22,10 @@ mainGlut: mainGlut.cpp
           -framework ApplicationServices -llapack -lblas -lz -lm
 
 mainPost: mainPost.cpp
-	g++ -g -m64 -o mainPost -Iinclude mainPost.cpp -Llib -lGmsh\
+	g++ -g  -o mainPost -Iinclude mainPost.cpp -Llib -lGmsh\
           -framework ApplicationServices -llapack -lblas -lz -lm
 
+
 levelset: mainLevelset.cpp
 	g++ -g -I${GMSHINCDIR} -DHAVE_GSL -DHAVE_GLEVELSETS -DHAVE_GMM -DHAVE_SPARSKIT \
 	-o levelset mainLevelset.cpp -L${GMSHLIBDIR} -lGmsh -lgsl -lgslcblas -llapack -lblas

utils/api_demos/mainSimple.cpp

@@ -15,8 +15,8 @@ int main(int argc, char **argv)
   GmshInitialize(argc, argv);
   //GmshSetOption("Mesh", "Algorithm", 5);
   GModel *m = new GModel();
-  //m->readGEO("../../tutorial/t5.geo");
-  GmshMergeFile("../../tutorial/t5.geo"); // will also set the bbox
+  m->readGEO("../../tutorial/t5.geo");
+  //GmshMergeFile("../../tutorial/t5.geo"); // will also set the bbox
   m->mesh(3);
   for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
     GRegion *r = *it;