diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 8b3f01f6a63e2c915160de3686935d300a80c3e2..5988610d3f3ec5c1aa4af37852165e814ac69352 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -843,7 +843,7 @@ void signedDistancesPointsTriangle(std::vector<double>&distances,
   if (det = 0.0) return;
 
   const double n2t1 = dot(t1, t1);
-  const double n2t2 = dot(t2, t3);
+  const double n2t2 = dot(t2, t2);
   const double n2t3 = dot(t3, t3);
 
   double u, v, d;
diff --git a/Numeric/cartesian.h b/Numeric/cartesian.h
index 3045c0e3710cd0c8d28fec53c4a726a072fba0bd..23bff5c6b0d0f965e6835a3890c73e5a3c3b8fda 100644
--- a/Numeric/cartesian.h
+++ b/Numeric/cartesian.h
@@ -148,8 +148,11 @@ class cartesianBox {
             _nodalValues[node_index(i+I,j+J,k+K)] = 0.0;
     }
   }
-  void writeMSH (const std::string &filename) const 
+  void writeMSH (const std::string &filename, bool simplex=false) const 
   {
+    
+//     const bool simplex=false;
+    
     FILE *f = fopen (filename.c_str(), "w");
     fprintf(f,"$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
     {
@@ -162,21 +165,76 @@ class cartesianBox {
       fprintf(f,"$EndNodes\n");
     }
     {
-      fprintf(f,"$Elements\n%d\n",_active.size());
+      int coeff=1;
+      if(simplex) coeff=6;
+      fprintf(f,"$Elements\n%d\n",coeff*_active.size());
       std::set<int>::const_iterator it = _active.begin();
       for ( ; it!=_active.end();++it){
-        fprintf(f,"%d 5 3 1 1 1",*it);
-        int i,j,k;
-        element_ijk(*it,i,j,k);
-        fprintf(f," %d",node_index(i,j,k));
-        fprintf(f," %d",node_index(i+1,j,k));
-        fprintf(f," %d",node_index(i+1,j+1,k));
-        fprintf(f," %d",node_index(i,j+1,k));
-        fprintf(f," %d",node_index(i,j,k+1));
-        fprintf(f," %d",node_index(i+1,j,k+1));
-        fprintf(f," %d",node_index(i+1,j+1,k+1));
-        fprintf(f," %d",node_index(i,j+1,k+1));
-        fprintf(f,"\n");
+        if(!simplex){
+          fprintf(f,"%d 5 3 1 1 1",*it);
+          int i,j,k;
+          element_ijk(*it,i,j,k);
+          fprintf(f," %d",node_index(i,j,k));
+          fprintf(f," %d",node_index(i+1,j,k));
+          fprintf(f," %d",node_index(i+1,j+1,k));
+          fprintf(f," %d",node_index(i,j+1,k));
+          fprintf(f," %d",node_index(i,j,k+1));
+          fprintf(f," %d",node_index(i+1,j,k+1));
+          fprintf(f," %d",node_index(i+1,j+1,k+1));
+          fprintf(f," %d",node_index(i,j+1,k+1));
+          fprintf(f,"\n");
+        }else{
+          int i,j,k;
+          element_ijk(*it,i,j,k);
+
+          //Elt1
+          fprintf(f,"%d 4 3 1 1 1",*it);
+          fprintf(f," %d",node_index(i,j+1,k));
+          fprintf(f," %d",node_index(i,j+1,k+1));
+          fprintf(f," %d",node_index(i+1,j,k+1));
+          fprintf(f," %d",node_index(i+1,j+1,k+1));
+          fprintf(f,"\n");
+          
+          //Elt2
+          fprintf(f,"%d 4 3 1 1 1",*it);
+          fprintf(f," %d",node_index(i,j+1,k));
+          fprintf(f," %d",node_index(i+1,j+1,k+1));
+          fprintf(f," %d",node_index(i+1,j,k+1));
+          fprintf(f," %d",node_index(i+1,j+1,k));
+          fprintf(f,"\n");
+          
+          //Elt3
+          fprintf(f,"%d 4 3 1 1 1",*it);
+          fprintf(f," %d",node_index(i,j+1,k));
+          fprintf(f," %d",node_index(i,j,k+1));
+          fprintf(f," %d",node_index(i+1,j,k+1));
+          fprintf(f," %d",node_index(i,j+1,k+1));
+          fprintf(f,"\n");
+          
+          //Elt4
+          fprintf(f,"%d 4 3 1 1 1",*it);
+          fprintf(f," %d",node_index(i,j+1,k));
+          fprintf(f," %d",node_index(i+1,j+1,k));
+          fprintf(f," %d",node_index(i+1,j,k+1));
+          fprintf(f," %d",node_index(i+1,j,k));
+          fprintf(f,"\n");
+          
+          //Elt5
+          fprintf(f,"%d 4 3 1 1 1",*it);
+          fprintf(f," %d",node_index(i,j+1,k));
+          fprintf(f," %d",node_index(i+1,j,k));
+          fprintf(f," %d",node_index(i+1,j,k+1));
+          fprintf(f," %d",node_index(i,j,k));
+          fprintf(f,"\n");
+          
+          //Elt6
+          fprintf(f,"%d 4 3 1 1 1",*it);
+          fprintf(f," %d",node_index(i,j+1,k));
+          fprintf(f," %d",node_index(i,j,k));
+          fprintf(f," %d",node_index(i+1,j,k+1));
+          fprintf(f," %d",node_index(i,j,k+1));
+          fprintf(f,"\n");
+        }
       }    
       fprintf(f,"$EndElements\n");    
     }