triangularSurfaceMesh.h

Go to the documentation of this file.

/* This file is part of the Palabos library.
 *
 * Copyright (C) 2011 FlowKit Sarl
 * Avenue de Chailly 23
 * 1012 Lausanne, Switzerland
 * E-mail contact: contact@flowkit.com
 *
 * The most recent release of Palabos can be downloaded at 
 * <http://www.palabos.org/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

/* Main author: Dimitrios Kontaxakis */

#ifndef TRIANGULAR_SURFACE_MESH_H
#define TRIANGULAR_SURFACE_MESH_H

#include "core/globalDefs.h"
#include "core/geometry3D.h"
#include <string>
#include <vector>

namespace plb {


template<typename T>
class TriangularSurfaceMesh {
public:
    struct Edge {
        plint pv;  // Pointing Vertex.
        plint ne;  // Neighboring Edge.
    };
    struct Lid {
        // Comment: if the close-holes algorithm is to be made more
        //   general one day, you will need to (1) assign the proper
        //   value to numAddedVertices, and (2) replace centerVertex
        //   by something else and change all codes which access
        //   centerVertex.
        Lid() : numAddedVertices(1) { }
        plint firstTriangle;
        plint numTriangles;
        std::vector<plint> boundaryVertices;
        plint centerVertex;
        plint numAddedVertices;
    };
public:
    TriangularSurfaceMesh (
            std::vector<Array<T,3> >& vertexList_,
            std::vector<plint>& emanatingEdgeList_,
            std::vector<Edge>& edgeList_,
            plint numVertices_ = -1 );

    plint getNumTriangles() const {
        return numTriangles;
    }
    plint getNumVertices() const {
        return numVertices;
    }
    void replaceVertex(plint iVertex, Array<T,3> const& newPosition);


    void resetVertices(Array<T,3> const& defaultVertex);

    Array<T,3> const& getVertex(plint iTriangle, int localVertex) const;
    Array<T,3> const& getVertex(plint iVertex) const;
    bool isValidVertex(plint iTriangle, int localVertex) const;
    bool isValidVertex(plint iVertex) const;

    void computeBoundingBox (
            Array<T,2>& xRange,
            Array<T,2>& yRange,
            Array<T,2>& zRange ) const;

    void translate(Array<T,3> const& vector);
    void scale(T alpha);
    void rotate(T phi, T theta, T psi);

    void smooth(plint maxiter = 1, T relax = 1.0, bool isMeasureWeighted = false);

    plint getVertexId(plint iTriangle, plint localVertex) const;
    std::vector<plint> getNeighborVertexIds(plint iVertex) const;
    std::vector<plint> getNeighborVertexIds(plint iVertex, plint jVertex) const;
    std::vector<plint> getNeighborTriangleIds(plint iVertex) const;
    std::vector<plint> getAdjacentTriangleIds(plint iTriangle) const;
    std::vector<plint> getAdjacentTriangleIds(plint iVertex, plint jVertex) const;

    Array<T,3> computeTriangleNormal(plint iTriangle, bool isAreaWeighted = false) const;
    Array<T,3> computeTriangleNormal(
            plint iVertex, plint jVertex, plint kVertex, bool isAreaWeighted = false) const;

    Array<T,3> computeEdgeNormal(
            plint iVertex, plint jVertex, bool isAreaWeighted = false) const;
    Array<T,3> computeVertexNormal(plint iVertex, bool isAreaWeighted = false) const;

    Array<T,3> computeContinuousNormal(
            Array<T,3> const& p, plint iTriangle, bool isAreaWeighted = false) const;

    T computeTriangleArea(plint iTriangle) const;
    T computeTriangleArea(plint iVertex, plint jVertex, plint kVertex) const;
    T computeEdgeArea(plint iVertex, plint jVertex) const;
    T computeVertexArea(plint iVertex) const;

    T computeEdgeLength(plint iVertex, plint jVertex) const;

    T computeDihedralAngle(plint iVertex, plint jVertex) const;

    T computeEdgeTileSpan(plint iVertex, plint jVertex) const;

    void writeAsciiSTL(std::string fname) const;
    void writeBinarySTL(std::string fname) const;

    bool isBoundaryVertex(plint iVertex) const;
    bool isInteriorVertex(plint iVertex) const;

    bool isBoundaryEdge(plint iVertex, plint jVertex) const;
    bool isInteriorEdge(plint iVertex, plint jVertex) const;

    int pointOnTriangle (
            Array<T,3> const& point1, Array<T,3> const& point2, int flag, plint iTriangle,
            Array<T,3>& intersection, Array<T,3>& normal, T& distance ) const;

    bool segmentIntersectsTriangle (
            Array<T,3> const& point1, Array<T,3> const& point2, plint iTriangle ) const;

    void distanceToEdgeLine (
            Array<T,3> const& point, plint iTriangle, plint whichEdge,
            T& distance, bool& intersectionIsInside ) const;

    void distanceToTrianglePlane (
            Array<T,3> const& point, plint iTriangle,
            T& distance, bool& intersectionIsInside, bool& pointIsBehind ) const;

    void distanceToTriangle (
            Array<T,3> const& point, plint iTriangle,
            T& distance, bool& pointIsBehind ) const;

    void reverseOrientation();

    std::vector<Lid> closeHoles();
    std::vector<std::vector<plint> > detectHoles();
    void avoidIntegerPositions();
    void avoidIntegerPosition(plint iVertex);
    void inflate(T amount=1.e-3);
public:
    std::vector<Array<T,3> >const& vertices() const {
        PLB_PRECONDITION(vertexList);
        return *vertexList;
    }
    std::vector<plint> const& emanatingEdges() const {
        PLB_PRECONDITION(emanatingEdgeList);
        return *emanatingEdgeList;
    }
    std::vector<Edge> const& edges() const {
        PLB_PRECONDITION(edgeList);
        return *edgeList;
    }
    Array<T,3>& getVertex(plint iTriangle, int localVertex);
    Array<T,3>& getVertex(plint iVertex);
private:
    void assertVertex(Array<T,3> const& vertex) const;
    std::vector<Array<T,3> >& vertices() {
        PLB_PRECONDITION(vertexList);
        return *vertexList;
    }
    std::vector<plint>& emanatingEdges() {
        PLB_PRECONDITION(emanatingEdgeList);
        return *emanatingEdgeList;
    }
    std::vector<Edge>& edges() {
        PLB_PRECONDITION(edgeList);
        return *edgeList;
    }
private:
    Lid closeHole(std::vector<plint> const& hole);
private:
    plint prev(plint iEdge) const;
    plint next(plint iEdge) const;
    plint changeEdgeId(plint iEdge) const;
private:
    std::vector<Array<T,3> >* vertexList;
    std::vector<plint>* emanatingEdgeList;
    std::vector<Edge>* edgeList;
    plint numTriangles, numVertices;
public:
    static const T eps0, eps1;
};

template<typename T>
class LidLessThan {
public:
    LidLessThan(plint mainDirection_, TriangularSurfaceMesh<T> const& mesh_)
        : mainDirection(mainDirection_),
          mesh(mesh_)
    { }
    bool operator()(typename TriangularSurfaceMesh<T>::Lid const& lid1,
                    typename TriangularSurfaceMesh<T>::Lid const& lid2) const
    {
        plint dim1 = mainDirection;
        plint dim2 = (dim1+1)%3;
        plint dim3 = (dim2+1)%3;
        // Compare on bary-center for main coordinate, and, in the unlikely
        //   case that these two floats are equal, on the two subsequent
        //   coordinates.
        Array<T,3> baryCenter1 = computeBaryCenter(mesh, lid1);
        Array<T,3> baryCenter2 = computeBaryCenter(mesh, lid2);
        return
          (baryCenter1[dim1] < baryCenter2[dim1]) || (
                  equiv(baryCenter1[dim1],baryCenter2[dim1]) && (
                      (baryCenter1[dim2] < baryCenter2[dim2]) || (
                          equiv(baryCenter1[dim2],baryCenter2[dim2]) &&
                          (baryCenter1[dim3] < baryCenter2[dim3]) ) ) ) ;
    }
    static bool equiv(T a, T b) {
        return !(a<b || b<a);
    }
private:
    plint mainDirection;
    TriangularSurfaceMesh<T> const& mesh;
};

template<typename T>
T toLatticeUnits (
        TriangularSurfaceMesh<T>& mesh,
        plint resolution, plint referenceDirection );



/* ******* Lid operations ************************************************** */

template<typename T>
Array<T,3> computeBaryCenter (
        TriangularSurfaceMesh<T> const& mesh,
        typename TriangularSurfaceMesh<T>::Lid const& lid );

template<typename T>
void computeBoundingBox (
        TriangularSurfaceMesh<T> const& mesh,
        typename TriangularSurfaceMesh<T>::Lid const& lid,
        Array<T,2>& xLim, Array<T,2>& yLim, Array<T,2>& zLim );

template<typename T>
T computeInnerRadius (
        TriangularSurfaceMesh<T> const& mesh,
        typename TriangularSurfaceMesh<T>::Lid const& lid );

template<typename T>
T computeOuterRadius (
        TriangularSurfaceMesh<T> const& mesh,
        typename TriangularSurfaceMesh<T>::Lid const& lid );

template<typename T>
T computeArea (
        TriangularSurfaceMesh<T> const& mesh,
        typename TriangularSurfaceMesh<T>::Lid const& lid );

template<typename T>
Array<T,3> computeNormal (
        TriangularSurfaceMesh<T> const& mesh,
        typename TriangularSurfaceMesh<T>::Lid const& lid );

template<typename T>
void reCenter (
        TriangularSurfaceMesh<T>& mesh,
        typename TriangularSurfaceMesh<T>::Lid const& lid );

} // namespace plb

#endif  // TRIANGULAR_SURFACE_MESH_H