diff --git a/src/facette.cpp b/src/facette.cpp
index 83c4a1f..c1cf520 100644
--- a/src/facette.cpp
+++ b/src/facette.cpp
@@ -125,3 +125,20 @@ double Fac::potential(std::function<Eigen::Vector3d(Nodes::Node)> getter, int i)
     return 0.5 * dMs * pot;
     }
 
+BasicTri::BasicTri(const Facette::Fac &fac, const bool surface)
+        : nodesInd({fac.ind[0], fac.ind[1], fac.ind[2]}), vRegion(fac.idxPrm), sRegion(surface)
+        { std::sort(nodesInd.begin(), nodesInd.end()); }
+
+bool BasicTri::operator<(const BasicTri &tri) const {
+    return (this->nodesInd[0] < tri.nodesInd[0])
+            || ((this->nodesInd[0] == tri.nodesInd[0]) && (this->nodesInd[1] < tri.nodesInd[1]))
+            || ((this->nodesInd[0] == tri.nodesInd[0]) && (this->nodesInd[1] == tri.nodesInd[1])
+                        && (this->nodesInd[2] < tri.nodesInd[2]));
+}
+
+bool BasicTri::operator==(const BasicTri &tri) const {
+    return (this->nodesInd[0] == tri.nodesInd[0])
+            && (this->nodesInd[1] == tri.nodesInd[1])
+            && (this->nodesInd[2] == tri.nodesInd[2]);
+}
+
diff --git a/src/facette.h b/src/facette.h
index d43e1ee..f4c315f 100644
--- a/src/facette.h
+++ b/src/facette.h
@@ -250,6 +250,53 @@ class Fac : public element<N,NPI>
         }
     };  // end class Fac
 
+/** \class BasicTri
+BasicTri is a simplified version of the class Fac which contains only 
+the nodes, the surface region and the region.
+It is used when checking the validity of the mesh is needed but 
+not all the facet elements are needed.
+*/
+class BasicTri
+    {
+public:
+    /** Constructor using the indices vector and the region. sRegion is false by default. */
+    inline BasicTri(const std::vector<int> &inds, const int idReg, const bool surface = false)
+        : nodesInd({inds[0], inds[1], inds[2]}), vRegion(idReg), sRegion(surface)
+        { std::sort(nodesInd.begin(), nodesInd.end()); }
+
+    /** Constructor used to copy a surface region facette */
+    explicit BasicTri(const Facette::Fac & fac, bool surface = false);
+
+    /** Compares the node indices in order */
+    bool operator<(const BasicTri & tri) const;
+
+    /** Returns true if all the nodes of the same indices are the same */
+    bool operator==(const BasicTri & tri) const;
+
+    /** Returns true if the region is the same */
+    inline bool sameVRegion(const BasicTri & tri) const
+        { return this->vRegion == tri.vRegion; }
+
+    /** getter for sRegion */
+    inline bool getSRegion() const
+        { return sRegion; }
+
+    /** Setter for sRegion */
+    inline void setSRegion(bool sRegion)
+        { this->sRegion = sRegion; }
+
+private:
+    /** The 3 node indices composing the triangle */
+    std::array<int,3> nodesInd;
+
+    /** index of the volume region of the triangle */
+    int vRegion;
+
+    /** If the triangle is in a surface region */
+    bool sRegion;
+
+    }; // class BasicTri
+
     }  // namespace Facette
 
 #endif /* facette_h */
diff --git a/src/mesh.cpp b/src/mesh.cpp
index b6aa8a5..b08d0f1 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -160,6 +160,82 @@ void mesh::sortNodes(Nodes::index long_axis)
                   });
     }
 
+void mesh::checkFacettes() const
+    {
+    std::vector<Facette::BasicTri> allTriCtnr;
+    const int tetraN = 4;
+    // Put all triangles into allTriCtnr after doing the conversion
+    std::for_each(tet.begin(), tet.end(),       // For each tetrahedron
+            [this, &allTriCtnr](const Tetra::Tet &tetrahedron)
+            {
+            for (int i = 0; i < tetraN; i++)    // For each 4 facets
+                {
+                Facette::BasicTri curTri({tetrahedron.ind[i], tetrahedron.ind[(i+1) % tetraN],
+                                                tetrahedron.ind[(i+2) % tetraN]}, 
+                                          tetrahedron.idxPrm);
+                allTriCtnr.push_back(curTri);
+                }
+            });
+    std::for_each(fac.begin(), fac.end(), [this, &allTriCtnr](const Facette::Fac &curFac)
+            {   // For each surface facet
+            Facette::BasicTri curTri(curFac, true);
+            allTriCtnr.push_back(curTri);
+            });
+
+    std::sort(allTriCtnr.begin(), allTriCtnr.end());
+    for (size_t i = 0; i < allTriCtnr.size(); i++)
+        {   // For each triangle
+        if (i < allTriCtnr.size() + 2 &&
+            allTriCtnr[i] == allTriCtnr[i+1] && allTriCtnr[i+1] == allTriCtnr[i+2])
+            {   // If 3 identical triangles are found
+
+            bool sReg1 = allTriCtnr[i].getSRegion();
+            bool sReg2 = allTriCtnr[i+1].getSRegion();
+            bool sReg3 = allTriCtnr[i+2].getSRegion();
+            if ((sReg1 && (sReg2 || sReg3)) || (sReg2 && sReg3)
+                || ((i < allTriCtnr.size() + 3 && allTriCtnr[i] == allTriCtnr[i+3])
+                && (allTriCtnr[i+3].getSRegion() && (sReg1 || sReg2 || sReg3))))
+                {
+                std::cout << "Error: wrong mesh generation. A triangle which belongs to "
+                             "multiple surface regions has been found\n";
+                exit(1);
+                }
+
+            else if (!(sReg1 ^ sReg2 ^ sReg3) || (sReg1 && sReg2 && sReg3) ||
+                (i < allTriCtnr.size() + 3 && allTriCtnr[i+2] == allTriCtnr[i+3]))
+                {
+                std::cout << "Error: wrong mesh generation. Three "
+                             "identical triangles have been found\n";
+                exit(1);
+                }
+
+            else if ((sReg1 && allTriCtnr[i+1].sameVRegion(allTriCtnr[i+2]))
+                    || (sReg2 && allTriCtnr[i].sameVRegion(allTriCtnr[i+2]))
+                    || (sReg3 && allTriCtnr[i].sameVRegion(allTriCtnr[i+1])))
+                {
+                std::cout << "Error: wrong mesh generation. An internal facet "
+                             "belonging to a surface region has been found\n";
+                exit(1);
+                }
+
+            else
+                { i = i + 2; }
+            }
+
+        else if (i < allTriCtnr.size() + 1 && allTriCtnr[i] == allTriCtnr[i+1])
+            {   // If 2 identical triangles are found
+            if (allTriCtnr[i].getSRegion() && allTriCtnr[i+1].getSRegion())
+                {
+                std::cout << "Error: wrong mesh generation. A triangle which belongs to "
+                             "multiple surface regions has been found\n";
+                exit(1);
+                }
+            else
+                { i++; }
+            }
+        }
+    }
+
 double mesh::surface(std::vector<int> &facIndices)
     {
     double S(0);
diff --git a/src/mesh.h b/src/mesh.h
index 914b06a..c0203b0 100644
--- a/src/mesh.h
+++ b/src/mesh.h
@@ -118,6 +118,7 @@ class mesh
                     }
                 });
 
+        checkFacettes();
         for(unsigned int i=0;i<fac.size();i++)
             {
             if (isMagnetic(fac[i]) && !mySets.paramFacette[fac[i].idxPrm].suppress_charges
@@ -430,6 +431,14 @@ class mesh
      * the matrix we will have to solve for. */
     void sortNodes(Nodes::index long_axis /**< [in] */);
 
+    /** Check if all facette verify one of the three following conditions :
+     *      - Is part of 2 tetraedrons
+     *      - Is part of 1 tetraedron and one surface region
+     *      - Is part of 1 tetraedron and no surface region
+     * Throws an error if it finds one which is not in one of those three cases.
+     */
+    void checkFacettes() const;
+
     /** returns the surface defined by the set of facets of indices in facIndices
      * each elementary surface triangle defined by points p0,p1,p2 is computed using
      * norm(cross(p0p1,p0p2))/2, it is always positive  */