tessellated normal computation checked and fixed

Unstructured_mesh/main.cpp

@@ -377,7 +377,7 @@ void testMesh2D() {
 
     DBGMSG("2D normals test");
 
-    auto normals = ComputeFaceNormals(mesh);
+    auto normals = ComputeFaceNormals<DEFAULT>(mesh);
     for(auto& edge : mesh.getEdges()){
         DBGVAR(edge.getIndex(),normals.at(edge));
     }
@@ -524,6 +524,13 @@ DBGMSG("tessellated cell volume");
         DBGVAR(face.getIndex(),normals.at(face));
     }
 
+    DBGMSG("3D normals test");
+
+    auto normalsTess = mesh3.computeFaceNormals<TESSELLATED>();
+    for(auto& face : mesh3.getFaces()){
+        DBGVAR(face.getIndex(),normalsTess.at(face));
+    }
+
     DBGMSG("mesh apply test");
     MeshApply<3, 2, 3>::apply(mesh3, [](size_t ori, size_t i){
         DBGVAR(ori,i);

src/NumericStaticArray/GrammSchmidt.h

@@ -9,8 +9,8 @@
  * @brief GrammSchmidt
  * Gramm-Schmidt process without normalization. Can be used to
  * calculate volume in any dimension.
- * @param vectors [in / out] vector containing the vectors and
- * after the process the vectors are changed.
+ * @param vectors [in / out] array containing the vectors to be processed.
+ * After the process the input vectors are changed.
  */
 template <unsigned int NumVec, unsigned int Dimension,typename IndexType, typename Real>
 void grammSchmidt(std::array<Vector<Dimension, Real>, NumVec>& vectors){
@@ -40,4 +40,46 @@ void grammSchmidt(std::array<Vector<Dimension, Real>, NumVec>& vectors){
     }
 }
 
+
+/**
+ * @brief GrammSchmidt
+ * Gramm-Schmidt process with normalization. Can be used to
+ * calculate volume in any dimension or normal vector. The norms of orthogonalized
+ * vectors are returned using norms array. This function is slightly more efficient
+ * than the method without nomalization.
+ * @param vectors [in / out] array containing the vectors to be processed.
+ * After the process the vectors are changed.
+ * @param norms [out] norms of the processed vectors.
+ */
+template <unsigned int NumVec, unsigned int Dimension,typename IndexType, typename Real>
+void grammSchmidt(std::array<Vector<Dimension, Real>, NumVec>& vectors, std::array<Real, NumVec>& norms){
+    /*
+     * Vector of inverse suquare of norm
+     */
+    std::array<Real, NumVec> coef;
+
+    norms.at(0) = vectors.at(0).normEukleid();
+    vectors.at(0) /= norms.at(0);
+
+    for (IndexType i = 1; i < vectors.size(); i++) {
+        /*
+         * Coefitiens of scalar products.
+         * The coefitients are computed in advance for better stability
+         */
+        for (IndexType j = 0; j < i; j++) {
+
+            coef.at(j) = (vectors.at(i)*(vectors.at(j)));
+
+        }
+        for (IndexType j = 0; j < i; j++) {
+
+            vectors.at(i) -= (vectors.at(j) * coef.at(j));
+
+        }
+
+        norms.at(i) = vectors.at(i).normEukleid();
+        vectors.at(i) /= norms.at(i);
+    }
+}
+
 #endif // GRAMMSCHMIDT_H

src/UnstructuredMesh/MeshFunctions/ComputeCenter.h

@@ -125,9 +125,9 @@ struct _ComputeCenters<2, 3, ComputationMethod::TESSELLATED> {
                 IndexType AI = mesh.getEdges().at(sub.index).getVertexAIndex();
                 IndexType BI = mesh.getEdges().at(sub.index).getVertexBIndex();
                 std::array<Vertex<3, Real>, 2> v = {elemCenters.at(i) - mesh.getVertices().at(AI), elemCenters.at(i) - mesh.getVertices().at(BI)};
-
-                grammSchmidt<2, 3, IndexType, Real>(v);
-                Real surf = v.at(0).normEukleid() * 0.5 * v.at(1).normEukleid();
+                std::array<Real, 2> norms;
+                grammSchmidt<2, 3, IndexType, Real>(v, norms);
+                Real surf = norms.at(0) * 0.5 * norms.at(1);
 
                 tempVert += subElemCenters.at(sub.index) * (surf * (2.0 / 3.0));
                 surfTotal += surf;

src/UnstructuredMesh/MeshFunctions/ComputeMeasures.h

@@ -16,7 +16,7 @@ template <unsigned int dim, unsigned int Dimension, ComputationMethod Method = D
 struct _ComputeMeasures{
     template <typename IndexType, typename Real, unsigned int ...Reserve>
     static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, Dimension>>&,MeshElements<Dimension, IndexType, Real, Reserve...>&){
-        static_assert (Dimension > 3,"The measure computation of mesh of dimension higher than 3 is not implemented yet.");
+        static_assert (Dimension <= 3,"The measure computation of mesh of dimension higher than 3 is not implemented yet.");
         throw std::runtime_error("The measure computation of mesh of dimension higher than 3 is not implemented yet.");
     }
 };
@@ -179,9 +179,10 @@ struct _ComputeMeasures<3, 3, TESSELLATED>{
                     Vertex<3,Real>& vertB = mesh.getVertices().at(mesh.getEdges().at(edgeIndex).getVertexBIndex());
 
                     std::array<Vertex<3,Real>, 3> pyramidVec = {vertA - faceCenter, vertB - faceCenter, cellCenter - faceCenter};
-                    grammSchmidt<3, 3, IndexType, Real>(pyramidVec);
+                    std::array<Real, 3> norms;
+                    grammSchmidt<3, 3, IndexType, Real>(pyramidVec, norms);
 
-                    measure += pyramidVec.at(0).normEukleid() * pyramidVec.at(1).normEukleid() * pyramidVec.at(2).normEukleid() * (1.0/6.0);
+                    measure += norms.at(0) * norms.at(1) * norms.at(2) * (1.0/6.0);
 
                 });
                 }
@@ -193,7 +194,7 @@ struct _ComputeMeasures<3, 3, TESSELLATED>{
     }
 };
 
-}
+} // namespace Detail
 
 
 

src/UnstructuredMesh/MeshFunctions/ComputeNormals.h

@@ -4,10 +4,12 @@
 #include "../MeshElements/MeshElement.h"
 #include "../MeshDataContainer/MeshDataContainer.h"
 #include "../../NumericStaticArray/Vector.h"
+#include "../../NumericStaticArray/GrammSchmidt.h"
 #include "MeshApply.h"
+#include "MeshFunctionsDefine.h"
 
 
-template <unsigned int Dimension>
+template <unsigned int Dimension, ComputationMethod Method = DEFAULT>
 struct _ComputeNormals{
     template <typename IndexType, typename Real, unsigned int ...Reserve>
     static void compute(MeshDataContainer<Vector<Dimension, Real>, Dimension-1>&,MeshElements<Dimension, IndexType, Real, Reserve...>&){
@@ -20,8 +22,8 @@ struct _ComputeNormals{
 
 
 
-template <>
-struct _ComputeNormals<2>{
+template <ComputationMethod Method>
+struct _ComputeNormals<2, Method>{
     template <typename IndexType, typename Real, unsigned int ...Reserve>
     static void compute(MeshDataContainer<Vector<2, Real>, 1>& normals,MeshElements<2, IndexType, Real, Reserve...>& mesh){
         for (auto& face : mesh.getEdges()) {
@@ -39,8 +41,8 @@ struct _ComputeNormals<2>{
 
 
 
-template <>
-struct _ComputeNormals<3>{
+template <ComputationMethod Method>
+struct _ComputeNormals<3, Method>{
     template <typename IndexType, typename Real, unsigned int ...Reserve>
     static void compute(MeshDataContainer<Vector<3, Real>, 2>& normals,MeshElements<3, IndexType, Real, Reserve...>& mesh){
         for (auto& face : mesh.getFaces()) {
@@ -99,16 +101,66 @@ struct _ComputeNormals<3>{
 };
 
 
+template <>
+struct _ComputeNormals<3, TESSELLATED>{
+    template <typename IndexType, typename Real, unsigned int ...Reserve>
+    static void compute(MeshDataContainer<Vector<3, Real>, 2>& normals,MeshElements<3, IndexType, Real, Reserve...>& mesh){
+        for (auto& face : mesh.getFaces()) {
+
+            bool vectorSign = true;
+            IndexType cellIndex = face.getCellLeftIndex();
+            if (
+                    cellIndex == INVALID_INDEX(IndexType) ||
+                    (cellIndex & BOUNDARY_INDEX(IndexType)) == BOUNDARY_INDEX(IndexType)
+                ) {
+                vectorSign = false;
+                cellIndex = face.getCellRightIndex();
+            }
+
+            Vertex<3, Real>& cellCenter = mesh.getCells().at(cellIndex).getCenter();
+            Vertex<3, Real>& faceCenter = face.getCenter();
+
+            Vector<3, Real> faceNormal = {};
+
+            Real surfTotal = Real();
+
+            MeshApply<2,1,3>::apply(face.getIndex(), mesh, [&](IndexType , IndexType edgeIndex){
+                Vertex<3,Real>& vertA = mesh.getVertices().at(mesh.getEdges().at(edgeIndex).getVertexAIndex());
+                Vertex<3,Real>& vertB = mesh.getVertices().at(mesh.getEdges().at(edgeIndex).getVertexBIndex());
+
+                std::array<Vertex<3,Real>, 3> pyramidVec = {faceCenter - vertA, faceCenter - vertB, faceCenter - cellCenter};
+                std::array<Real, 3> norms;
+
+                grammSchmidt<3, 3, IndexType, Real>(pyramidVec, norms);
+
+                faceNormal += pyramidVec[2] * 0.5 * norms[0] * norms[1];
+
+                surfTotal += 0.5 * norms[0] * norms[1];
+
+            });
+
+            faceNormal /= surfTotal;
+
+            if (!vectorSign) {
+                faceNormal *= -1;
+            }
+            normals.at(face)[0] = fabs(faceNormal[0]) < 1e-8 ? 0 : faceNormal[0];
+            normals.at(face)[1] = fabs(faceNormal[1]) < 1e-8 ? 0 : faceNormal[1];
+            normals.at(face)[2] = fabs(faceNormal[2]) < 1e-8 ? 0 : faceNormal[2];
+        }
+    }
+};
+
 
 
 
 
-template <unsigned int Dimension,typename IndexType, typename Real, unsigned int ...Reserve>
+template <ComputationMethod Method, unsigned int Dimension,typename IndexType, typename Real, unsigned int ...Reserve>
 MeshDataContainer<Vector<Dimension, Real>, Dimension-1> ComputeFaceNormals(MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
 
     MeshDataContainer<Vector<Dimension, Real>, Dimension-1> normals(mesh);
 
-    _ComputeNormals<Dimension>::compute(normals, mesh);
+    _ComputeNormals<Dimension, Method>::compute(normals, mesh);
 
     return normals;
 }

src/UnstructuredMesh/MeshFunctions/EdgesOrientation.h

@@ -5,6 +5,7 @@
 #include "../../NumericStaticArray/Vector.h"
 #include "../../Debug/Debug.h"
 #include "MeshApply.h"
+#include "ComputeNormals.h"
 #include <valarray>
 #include <set>
 #include <map>
@@ -63,7 +64,7 @@ bool edgeIsLeft(MeshElements<3, IndexType, Real, Reserve...>& mesh, IndexType fa
     typename MeshElements<3, IndexType, Real, Reserve...>::Edge& edge = mesh.getEdges().at(edgeIndex);
     typename MeshElements<3, IndexType, Real, Reserve...>::template ElementType<2>& face = mesh.template getElements<2>().at(faceIndex);
 
-    auto normals = ComputeFaceNormals(mesh);
+    auto normals = ComputeFaceNormals<ComputationMethod::DEFAULT>(mesh);
 
     return edgeIsLeft(mesh, face, edge, normals[face]);
 
@@ -73,7 +74,7 @@ template<typename IndexType, typename Real, unsigned int ...Reserve>
 MeshDataContainer<std::vector<bool>, 2> edgesOrientation(MeshElements<3, IndexType, Real, Reserve...>& mesh) {
 
     MeshDataContainer<std::vector<bool>, 2> orientations(mesh);
-    auto normals = ComputeFaceNormals(mesh);
+    auto normals = ComputeFaceNormals<ComputationMethod::DEFAULT>(mesh);
 
     for (auto& face : mesh.getFaces()) {
         orientations[face].resize(face.getSubelements().getNumberOfSubElements());

src/UnstructuredMesh/UnstructuredMesh.h

@@ -29,8 +29,9 @@ public:
         return ComputeMeasures<Method>(*this);
     }
 
+    template<ComputationMethod Method = ComputationMethod::DEFAULT>
     MeshDataContainer<Vector<Dimension, Real>, Dimension-1> computeFaceNormals() {
-        return ::ComputeFaceNormals(*this);
+        return ::ComputeFaceNormals<Method>(*this);
     }
 
 /*