hybrid meshes are generated whenever quads are present on the boundary

not yet finished ... mesh smoothing has to be adapted 

Geo/GRegion.cpp

@@ -31,6 +31,17 @@ GRegion::~GRegion()
   deleteMesh();
 }
 
+void GRegion::set(const std::list<GFace*> f) { 
+  for (std::list<GFace*>::iterator it =  l_faces.begin(); it !=  l_faces.end() ; ++it){
+    (*it)->delRegion(this);
+  }
+  l_faces = f; 
+  for (std::list<GFace*>::iterator it =  l_faces.begin(); it !=  l_faces.end() ; ++it){
+    (*it)->addRegion(this);
+  }
+}
+
+
 void GRegion::deleteMesh()
 {
   for(unsigned int i = 0; i < mesh_vertices.size(); i++) delete mesh_vertices[i];

Geo/GRegion.h

@@ -48,7 +48,7 @@ class GRegion : public GEntity {
   // get/set faces that bound the region
   virtual std::list<GFace*> faces() const{ return l_faces; }
   virtual std::list<int> faceOrientations() const{ return l_dirs; }
-  inline void set(const std::list<GFace*> f) { l_faces = f; }
+  void set(const std::list<GFace*> f);
 
   // edges that bound the region
   virtual std::list<GEdge*> edges() const;

Mesh/meshGFaceDelaunayInsertion.cpp

@@ -373,29 +373,6 @@ void connectTris(ITER beg, ITER end)
   }
 }
 
-template <class ITER>
-void connectQuas(ITER beg, ITER end)
-{
-  std::set<edgeXquad> conn;
-  while (beg != end){
-    for (int i = 0; i < 4; i++){
-      edgeXquad fxt(*beg, i);
-      std::set<edgeXquad>::iterator found = conn.find(fxt);
-      if (found == conn.end())
-	conn.insert(fxt);
-      else if (found->t1 != *beg){
-	found->t1->setNeigh(found->i1, *beg);
-	(*beg)->setNeigh(i, found->t1);
-      }
-    }
-    ++beg;
-  }
-}
-
-void connectQuads(std::vector<MQua4*> &l)
-{
-  connectQuas(l.begin(), l.end());
-}
 
 void connectTriangles(std::list<MTri3*> &l)
 {
@@ -583,7 +560,7 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
     // avoid angles that are too obtuse
     double cosv = ((d1*d1+d2*d2-d3*d3)/(2.*d1*d2));
 
-    if ((d1 < LL * .25 || d2 < LL * .25 || cosv < -0.5) && !force) {
+    if ((d1 < LL * .25 || d2 < LL * .25 || cosv < -.9) && !force) {
       onePointIsTooClose = true;
     }
 

Mesh/meshGFaceDelaunayInsertion.h

@@ -81,48 +81,6 @@ class MTri3
   }
 };
 
-/*        2
-  3 +------------+ 2
-    |            |
-    |            |
-  3 |            | 1
-    |            |
-    |            |
-    +------------+
-  0       0        1 
-	  
-   We require that quads are oriented
-   We'd like to walk into the quad mesh and 
-   create sheets
-   
-   We start from a given quad and one edge
-   We give a path as an array of int's
-   If path[i] == 1 --> 
-      
- */
-
-class MQua4
-{
- protected :
-  MQuadrangle *base;
-  MQua4 *neigh[4];
- public :
- MQua4(MQuadrangle *q) : base(q) {
-    neigh[0] = neigh[1] = neigh[2] = neigh[3] = NULL;}
-  inline MQuadrangle *qua() const { return base; }
-  inline void  setNeigh(int iN , MQua4 *n) { neigh[iN] = n; }
-  inline MQua4 *getNeigh(int iN ) const { return neigh[iN]; }
-  inline int getEdge(MVertex *v1, MVertex *v2){
-    for (int i=0;i<4;i++){
-      MEdge e = base->getEdge(i);
-      if (e.getVertex(0) == v1 && e.getVertex(1) == v2)return i;
-      if (e.getVertex(0) == v2 && e.getVertex(1) == v1)return i;      
-    }
-    return -1;
-  }
-};
-
-
 class compareTri3Ptr
 {
  public:
@@ -134,7 +92,6 @@ class compareTri3Ptr
   }
 };
 
-void connectQuads(std::vector<MQua4*> &);
 void connectTriangles(std::list<MTri3*> &);
 void connectTriangles(std::vector<MTri3*> &);
 void connectTriangles(std::set<MTri3*,compareTri3Ptr> &AllTris);
@@ -163,25 +120,5 @@ struct edgeXface
   }
 };
 
-struct edgeXquad
-{
-  MVertex *v[2];
-  MQua4 * t1;
-  int i1;
-  edgeXquad(MQua4 *_t, int iFac) : t1(_t), i1(iFac)
-  {
-    v[0] = t1->qua()->getVertex(iFac);
-    v[1] = t1->qua()->getVertex((iFac+1)%4);
-    std::sort(v, v + 2);
-  }
-  inline bool operator < ( const edgeXquad &other) const
-  {
-    if(v[0] < other.v[0]) return true;
-    if(v[0] > other.v[0]) return false;
-    if(v[1] < other.v[1]) return true;
-    return false;
-  }
-};
-
 
 #endif

Mesh/meshGRegion.cpp

@@ -19,6 +19,7 @@
 #include "MTriangle.h"
 #include "MQuadrangle.h"
 #include "MTetrahedron.h"
+#include "MPyramid.h"
 #include "BDS.h"
 #include "Context.h"
 #include "GFaceCompound.h"
@@ -30,6 +31,64 @@
 #include "ANN/ANN.h"
 #endif
 
+// hybrid mesh recovery structure
+class splitQuadRecovery {
+  std::multimap<GEntity*, std::pair<MVertex*,MFace> >_data;
+ public : 
+  void add (const MFace &f, MVertex *v, GEntity*ge) ;
+  int buildPyramids (GModel *gm);
+};
+
+void splitQuadRecovery::add (const MFace &f, MVertex *v, GEntity*ge) {
+  _data.insert(std::make_pair(ge, std::make_pair(v,f)));
+}
+
+// re-create quads and 
+int splitQuadRecovery::buildPyramids (GModel *gm) {
+  int NBPY = 0;
+  std::set<MFace,Less_Face> allFaces;
+  for (GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){
+    for (unsigned int i=0; i < (*it)->triangles.size();i++){
+      allFaces.insert ((*it)->triangles[i]->getFace(0));
+      delete (*it)->triangles[i];
+    }
+    (*it)->triangles.clear();
+    for ( std::multimap<GEntity*, std::pair<MVertex*,MFace> >::iterator it2 = _data.lower_bound(*it) ; it2 != _data.upper_bound(*it) ; ++it2) {
+      MVertex *v = it2->second.first;
+      v->onWhat()->mesh_vertices.erase
+	(std::find(v->onWhat()->mesh_vertices.begin(),
+		   v->onWhat()->mesh_vertices.end(),
+		   v));
+      const MFace &f = it2->second.second;
+      std::set<MFace,Less_Face> :: iterator itf0 = allFaces.find(MFace(f.getVertex(0),f.getVertex(1),v));
+      std::set<MFace,Less_Face> :: iterator itf1 = allFaces.find(MFace(f.getVertex(1),f.getVertex(2),v));
+      std::set<MFace,Less_Face> :: iterator itf2 = allFaces.find(MFace(f.getVertex(2),f.getVertex(3),v));
+      std::set<MFace,Less_Face> :: iterator itf3 = allFaces.find(MFace(f.getVertex(3),f.getVertex(0),v));
+      if (itf0 != allFaces.end() && itf1 != allFaces.end() && itf2 != allFaces.end() && itf3 != allFaces.end() ){
+	(*it)->quadrangles.push_back(new MQuadrangle(f.getVertex(0),f.getVertex(1),f.getVertex(2),f.getVertex(3)));
+	allFaces.erase(*itf0);
+	allFaces.erase(*itf1);
+	allFaces.erase(*itf2);
+	allFaces.erase(*itf3);
+
+	//	printf("some pyramids should be created %d regions\n",(*it)->numRegions());
+	for (int iReg = 0; iReg < (*it)->numRegions(); iReg++){
+	  //	  printf("one pyramid is created\n");
+	  if (iReg == 1) {Msg::Fatal("cannot build pyramids on non manifold faces ..."); v = new MVertex (v->x(),v->y(),v->z(),(*it)->getRegion(iReg)); }
+	  else v->setEntity ((*it)->getRegion(iReg));
+	  (*it)->getRegion(iReg)->pyramids.push_back(new MPyramid(f.getVertex(0),f.getVertex(1),f.getVertex(2),f.getVertex(3),v));      
+	  (*it)->getRegion(iReg)->mesh_vertices.push_back(v);
+	  NBPY++;
+	}
+      }
+    }
+    for (std::set<MFace,Less_Face> :: iterator itf = allFaces.begin() ; itf != allFaces.end(); ++itf){
+      (*it)->triangles.push_back(new MTriangle(itf->getVertex(0),itf->getVertex(1),itf->getVertex(2)));
+    }        
+  }
+  return NBPY;
+}
+
 void printVoronoi(GRegion *gr,  std::set<SPoint3> &candidates)
 {
   std::vector<MTetrahedron*> elements = gr->tetrahedra;
@@ -254,6 +313,72 @@ void skeletonFromVoronoi(GRegion *gr, std::set<SPoint3> &voronoiPoles)
 #endif
 }
 
+void getBoundingInfoAndSplitQuads (GRegion *gr, 
+				   std::map<MFace,GEntity*,Less_Face> &allBoundingFaces,
+				   std::set<MVertex*> &allBoundingVertices,
+				   splitQuadRecovery &sqr){
+  
+  std::map<MFace,GEntity*,Less_Face> allBoundingFaces_temp;
+
+  // Get all the faces that are on the boundary
+  std::list<GFace*> faces = gr->faces();
+  std::list<GFace*>::iterator it = faces.begin();
+  while (it != faces.end()){
+    GFace *gf = (*it);
+    for(unsigned int i = 0; i < gf->getNumMeshElements(); i++){
+      allBoundingFaces_temp[gf->getMeshElement(i)->getFace(0)] = gf;
+    }
+    ++it;
+  }
+
+  // if some elements pre-exist in the mesh, then use the
+  // internal faces of those
+
+  for (unsigned int i=0;i<gr->getNumMeshElements();i++){
+    MElement *e = gr->getMeshElement(i);
+    for (int j=0;j<e->getNumFaces();j++){
+      std::map<MFace,GEntity*,Less_Face>::iterator it = allBoundingFaces_temp.find(e->getFace(j));
+      if (it == allBoundingFaces_temp.end())allBoundingFaces_temp[e->getFace(j)] = gr;
+      else allBoundingFaces_temp.erase(it);
+    }
+  }
+  
+  std::map<MFace,GEntity*,Less_Face>::iterator itx = allBoundingFaces_temp.begin();
+  for (; itx != allBoundingFaces_temp.end();++itx){
+    const MFace &f = itx->first;
+    // SPLIT that face into 4 triangular faces 
+    if (f.getNumVertices() == 4){
+      MVertex *v0 = f.getVertex(0); 
+      MVertex *v1 = f.getVertex(1); 
+      MVertex *v2 = f.getVertex(2); 
+      MVertex *v3 = f.getVertex(3); 
+
+      MVertex *newv = new MVertex ((v0->x() + v1->x() + v2->x() + v3->x())*0.25,
+				   (v0->y() + v1->y() + v2->y() + v3->y())*0.25,
+				   (v0->z() + v1->z() + v2->z() + v3->z())*0.25,itx->second); 
+      sqr.add(f,newv,itx->second);
+      allBoundingFaces[MFace(v0,v1,newv)] = itx->second;
+      allBoundingFaces[MFace(v1,v2,newv)] = itx->second;
+      allBoundingFaces[MFace(v2,v3,newv)] = itx->second;
+      allBoundingFaces[MFace(v3,v0,newv)] = itx->second;
+      itx->second->mesh_vertices.push_back(newv);
+      //      myGr->pyramids.push_back(new MPyramid(v0,v1,v2,v3,newv));      
+      allBoundingVertices.insert(v0);
+      allBoundingVertices.insert(v1);
+      allBoundingVertices.insert(v2);
+      allBoundingVertices.insert(v3);
+      allBoundingVertices.insert(newv);
+    }
+    else{
+      allBoundingFaces[f] = itx->second;
+      allBoundingVertices.insert(f.getVertex(0));
+      allBoundingVertices.insert(f.getVertex(1));
+      allBoundingVertices.insert(f.getVertex(2));
+    }
+  }  
+}
+
+
 void getAllBoundingVertices(GRegion *gr, std::set<MVertex*> &allBoundingVertices)
 {
   std::list<GFace*> faces = gr->faces();
@@ -273,10 +398,13 @@ void getAllBoundingVertices(GRegion *gr, std::set<MVertex*> &allBoundingVertices
 
 #if defined(HAVE_TETGEN)
 #include "tetgen.h"
-void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numberedV)
+void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numberedV,  splitQuadRecovery &sqr)
 {
   std::set<MVertex*> allBoundingVertices;
-  getAllBoundingVertices(gr, allBoundingVertices);
+  std::map<MFace,GEntity*,Less_Face> allBoundingFaces;
+  getBoundingInfoAndSplitQuads (gr, allBoundingFaces,allBoundingVertices,sqr);
+  
+  //getAllBoundingVertices(gr, allBoundingVertices);
 
   in.mesh_dim = 3;
   in.firstnumber = 1;
@@ -314,25 +442,14 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb
 	I++;
   }
 
-  int nbFace = 0;
-  std::list<GFace*> faces = gr->faces();
-  std::list<GFace*>::iterator it = faces.begin();
-  while(it != faces.end()){
-    GFace *gf = (*it);
-    nbFace += gf->triangles.size();
-    ++it;
-  }
-
-  in.numberoffacets = nbFace;
+  in.numberoffacets = allBoundingFaces.size();
   in.facetlist = new tetgenio::facet[in.numberoffacets];
   in.facetmarkerlist = new int[in.numberoffacets];
   
-  it = faces.begin();
   I = 0;
-  while(it != faces.end()){
-    GFace *gf = (*it);
-    for(unsigned int i = 0; i < gf->triangles.size(); i++){
-      MTriangle *t = gf->triangles[i];
+  std::map<MFace,GEntity*,Less_Face>::iterator it = allBoundingFaces.begin();
+  for (; it != allBoundingFaces.end();++it){
+    const MFace &fac = it->first;
     tetgenio::facet *f = &in.facetlist[I];
     tetgenio::init(f);
     f->numberofholes = 0;
@@ -342,14 +459,12 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb
     tetgenio::init(p);
     p->numberofvertices = 3;
     p->vertexlist = new int[p->numberofvertices];
-      p->vertexlist[0] = t->getVertex(0)->getIndex();
-      p->vertexlist[1] = t->getVertex(1)->getIndex();
-      p->vertexlist[2] = t->getVertex(2)->getIndex();
-      in.facetmarkerlist[I] = gf->tag();
+    p->vertexlist[0] = fac.getVertex(0)->getIndex();
+    p->vertexlist[1] = fac.getVertex(1)->getIndex();
+    p->vertexlist[2] = fac.getVertex(2)->getIndex();
+    in.facetmarkerlist[I] = (it->second->dim() == 3) ? -it->second->tag() : it->second->tag();
     ++I;    
   }  
-    ++it;
-  }
 }
 
 void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
@@ -379,7 +494,10 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
     GFace *gf = (*it);
     for(unsigned int i = 0; i < gf->triangles.size(); i++)
       delete gf->triangles[i];
+    for(unsigned int i = 0; i < gf->quadrangles.size(); i++)
+      delete gf->quadrangles[i];
     gf->triangles.clear();
+    gf->quadrangles.clear();
     gf->deleteVertexArrays();
   }
 
@@ -397,10 +515,12 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
   // nodes for non manifold geometries (single surface connected to
   // volume)
   for(int i = 0; i < out.numberoftrifaces; i++){
+    if (out.trifacemarkerlist[i] > 0) {
       MVertex *v[3];
       v[0] = numberedV[out.trifacelist[i * 3 + 0] - 1];
       v[1] = numberedV[out.trifacelist[i * 3 + 1] - 1];
       v[2] = numberedV[out.trifacelist[i * 3 + 2] - 1];
+      // do not recover prismatic faces !!!
       GFace *gf = gr->model()->getFaceByTag(out.trifacemarkerlist[i]);
       
       double guess[2] = {0, 0};
@@ -476,8 +596,8 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
       }
       MTriangle *t = new MTriangle(v[0], v[1], v[2]);
       gf->triangles.push_back(t);
+    }// do not do anything with prismatic faces
   }
-
   for(int i = 0; i < out.numberoftetrahedra; i++){
     MVertex *v1 = numberedV[out.tetrahedronlist[i * 4 + 0] - 1];
     MVertex *v2 = numberedV[out.tetrahedronlist[i * 4 + 1] - 1];
@@ -494,6 +614,34 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GRegion *gr)
   buildAdditionalPoints3D (gr);
 }
 
+void _relocateVertex(MVertex *ver,
+		     const std::vector<MElement*> &lt)
+{
+  if(ver->onWhat()->dim() != 3) return;
+  double x = 0, y=0, z=0;
+  int N = 0;
+  bool isPyramid = false;
+  for(unsigned int i = 0; i < lt.size(); i++){
+    double XCG=0,YCG=0,ZCG=0;
+    for (int j=0;j<lt[i]->getNumVertices();j++){
+      XCG += lt[i]->getVertex(j)->x();
+      YCG += lt[i]->getVertex(j)->y();
+      ZCG += lt[i]->getVertex(j)->z();
+    }
+    x += XCG;
+    y += YCG;
+    z += ZCG;
+    if (lt[i]->getNumVertices() == 5) isPyramid = true;
+    N += lt[i]->getNumVertices();
+  }
+  if (isPyramid){
+    ver->x() = x / N;
+    ver->y() = y / N;
+    ver->z() = z / N;
+  }
+}
+
+
 
 void MeshDelaunayVolume(std::vector<GRegion*> &regions)
 {
@@ -523,21 +671,26 @@ void MeshDelaunayVolume(std::vector<GRegion*> &regions)
 
   // mesh with tetgen, possibly changing the mesh on boundaries (leave
   // this in block, so in/out are destroyed before we refine the mesh)
+  splitQuadRecovery sqr;
   {
     tetgenio in, out;
     std::vector<MVertex*> numberedV;
     char opts[128];
-    buildTetgenStructure(gr, in, numberedV);
+    buildTetgenStructure(gr, in, numberedV, sqr);
+    // JFR REMOVED THE -q OPTION BECAUSE MESH SIZES AT VERTICES WERE WRONG ...
+    // COULD BE FIXED IF NEEDED
     if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_DEL ||
        CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_HEX ||
        CTX::instance()->mesh.algo3d == ALGO_3D_MMG3D ||
        CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD ||
        CTX::instance()->mesh.algo2d == ALGO_2D_BAMG){
-      sprintf(opts, "-q1.5pY%c",  (Msg::GetVerbosity() < 3) ? 'Q':
+      sprintf(opts, "-pY%c",  (Msg::GetVerbosity() < 3) ? 'Q':
 	      (Msg::GetVerbosity() > 6) ? 'V': '\0');
+      // sprintf(opts, "-q1.5pY%c",  (Msg::GetVerbosity() < 3) ? 'Q':
+      // 	      (Msg::GetVerbosity() > 6) ? 'V': '\0');
     }
     else {
-      sprintf(opts, "-q3.5Ype%c",  (Msg::GetVerbosity() < 3) ? 'Q':
+      sprintf(opts, "-Ype%c",  (Msg::GetVerbosity() < 3) ? 'Q':
       	      (Msg::GetVerbosity() > 6) ? 'V': '\0');
     }
     try{
@@ -610,6 +763,20 @@ void MeshDelaunayVolume(std::vector<GRegion*> &regions)
    if(!Filler::get_nbr_new_vertices() && !LpSmoother::get_nbr_interior_vertices()){
      insertVerticesInRegion(gr);     
    }
+
+ if (sqr.buildPyramids (gr->model())){
+   // relocate vertices if pyramids
+   for(unsigned int k = 0; k < regions.size(); k++){
+     v2t_cont adj;
+     buildVertexToElement(regions[k]->tetrahedra, adj);
+     buildVertexToElement(regions[k]->pyramids, adj);
+     v2t_cont::iterator it = adj.begin();
+     while (it != adj.end()){
+       //_relocateVertex( it->first, it->second);
+       ++it;
+     }
+   }
+ }
 #endif
 }
 
@@ -872,13 +1039,17 @@ void meshGRegion::operator() (GRegion *gr)
 
   std::list<GFace*> faces = gr->faces();
 
+  // REMOVE SANITY CHECK FOR DELAUNAY : PYRAMIDS AVAILABLE
   // sanity check
+
+  if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL){
     for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
       if((*it)->quadrangles.size()){
 	Msg::Error("Cannot tetrahedralize volume with quadrangles on boundary");
 	return;
       }
     }
+  }
 
   // replace discreteFaces by their compounds
   if(gr->geomType() == GEntity::CompoundVolume){

Mesh/meshGRegion.h

@@ -69,4 +69,6 @@ class adaptMeshGRegion {
   void operator () (GRegion *);
 };
 
+
+
 #endif

Mesh/qualityMeasures.cpp

@@ -172,12 +172,42 @@ double qmTet(const double &x1, const double &y1, const double &z1,
       return 2. * sqrt(6.) * rhoin / l;
     }
     break;
+  case QMTET_COND:
+    {
+      /// condition number is defined as (see Knupp & Freitag in IJNME) 
+      double INVW[3][3] = {{1,-1./sqrt(3.),-1./sqrt(6.)},{0,2/sqrt(3.),-1./sqrt(6.)},{0,0,sqrt(1.5)}};
+      double A[3][3] = {{x2-x1,y2-y1,z2-z1},{x3-x1,y3-y1,z3-z1},{x4-x1,y4-y1,z4-z1}};
+      double S[3][3],INVS[3][3];
+      matmat(A,INVW,S);
+      *volume = inv3x3(S,INVS) * 2 / sqrt(2);
+      double normS = norm2 (S);
+      double normINVS = norm2 (INVS);
+      return normS * normINVS;      
+    }
   default:
     Msg::Error("Unknown quality measure");
     return 0.;
   }
 }
 
+/*
+double conditionNumberAndDerivativeOfTet(const double &x1, const double &y1, const double &z1,
+					 const double &x2, const double &y2, const double &z2,
+					 const double &x3, const double &y3, const double &z3,
+					 const double &x4, const double &y4, const double &z4){
+
+  double INVW[3][3] = {{1,-1./sqrt(3.),-1./sqrt(6.)},{0,2/sqrt(3.),-1./sqrt(6.)},{0,0,sqrt(1.5)}};
+      double A[3][3] = {{x2-x1,y2-y1,z2-z1},{x3-x1,y3-y1,z3-z1},{x4-x1,y4-y1,z4-z1}};
+      double S[3][3],INVS[3][3];
+      matmat(A,INVW,S);
+      double sigma = inv3x3(S,INVS);
+      double normS = norm2 (S);
+      double normINVS = norm2 (INVS);
+      conditionNumber = normS * normINVS;      
+  
+}
+*/
+
 double mesh_functional_distorsion(MElement *t, double u, double v)
 {
   // compute uncurved element jacobian d_u x and d_v x

Numeric/Numeric.cpp

@@ -35,6 +35,19 @@ double myacos(double a)
   else
     return acos(a);
 }
+double norm2(double a[3][3]) {
+  double norm2sq = 
+    SQU(a[0][0])+
+    SQU(a[0][1])+
+    SQU(a[0][2])+
+    SQU(a[1][0])+
+    SQU(a[1][1])+
+    SQU(a[1][2])+
+    SQU(a[2][0])+
+    SQU(a[2][1])+
+    SQU(a[2][2]);
+  return sqrt(norm2sq);
+}
 
 void matvec(double mat[3][3], double vec[3], double res[3])
 {

Numeric/Numeric.h

@@ -65,6 +65,8 @@ inline double norme(double a[3])
   }
   return mod;
 }
+double norm2(double a[3][3]);
+
 void normal3points(double x0, double y0, double z0,
                    double x1, double y1, double z1,
                    double x2, double y2, double z2,