adaptiveTriangMesh.h

Go to the documentation of this file.

#ifndef __ADAPTIVETRIANGMESH_H
#define __ADAPTIVETRIANGMESH_H

#include <triangMesh.h>

typedef short LocalIndex;
typedef int GlobalIndex;
const int IndexNotSet = -1;

template <typename DomVecType>
struct ParentInformation {
    int globalIndices[2];
    int ElementIndex;
};

namespace shellFE
{

template< typename DataTypeContainer, typename TriangleType >
class AdaptiveTriangMesh : public TriangMesh<DataTypeContainer, TriangleType>
{
public:
    typedef TriangleType ElementType;
    typedef typename DataTypeContainer::RealType      RealType;
    typedef typename DataTypeContainer::DomVecType    DomVecType;
    typedef typename DataTypeContainer::Point3DType   Point3DType;
    typedef typename DataTypeContainer::Indices3DType Indices3DType;
    typedef typename DataTypeContainer::VectorType    VectorType;
    typedef typename DataTypeContainer::MaskType      MaskType;
    typedef std::vector< Indices3DType >              VertexIndicesType;
    typedef std::vector< Point3DType >                VertexCoordsType;

    class DartIterator
    {
        typedef AdaptiveTriangMesh<DataTypeContainer, TriangleType> MeshType;

    protected:
        const MeshType & _mesh;
        GlobalIndex _triangle;
        LocalIndex  _node;
        LocalIndex  _edge;

    public:
        DartIterator(const MeshType &Mesh, GlobalIndex triangleIndex, LocalIndex localEdgeIndex, LocalIndex localNodeIndex) :
            _mesh(Mesh), _triangle (triangleIndex), _node(localNodeIndex), _edge(localEdgeIndex) {}

        DartIterator(const MeshType &Mesh, GlobalIndex triangleIndex, LocalIndex localEdgeIndex) :
            _mesh(Mesh), _triangle (triangleIndex), _node((localEdgeIndex + 1) % 3), _edge(localEdgeIndex) {}

        void set ( const GlobalIndex triangle, const LocalIndex edge, const LocalIndex node )
        {
            _triangle = triangle;
            _node = node;
            _edge = edge;
        }

        GlobalIndex getGlobalTriangleIndex() const
        {
            return _triangle;
        }

        GlobalIndex getGlobalNodeIndex()     const
        {
            return _mesh.getTriangNodeIdx( _triangle, _node );
        }

        LocalIndex  getLocalNodeIndex()      const
        {
            return _node;
        }

        LocalIndex  getLocalEdgeIndex()      const
        {
            return _edge;
        }

        LocalIndex  getNextNodeLocalIndex()  const
        {
            return 3 - (_node + _edge);
        }

        LocalIndex  getNextEdgeLocalIndex()  const
        {
            return 3 - (_node + _edge);
        }

        GlobalIndex getNextNodeGlobalIndex() const
        {
            return _mesh.getTriangNodeIdx( _triangle, getNextNodeLocalIndex() );
        }

        GlobalIndex getNextTriangleIndex()   const
        {
            return _mesh.getNeighbour( _triangle, _edge );
        }

        template<typename DartIteratorType>
        LocalIndex getCommonNodeLocalIndex(const DartIteratorType & d) const
        {
            if(_triangle != d.getGlobalTriangleIndex())
                throw std::invalid_argument( aol::strprintf("T=%d, but this->tringle = %d. In File %s at line %d.", d.getGlobalTriangleIndex(), _triangle, __FILE__, __LINE__).c_str() );
            if ( _edge == d.getLocalEdgeIndex() )
                throw std::invalid_argument( aol::strprintf("points to same edge=. In File %s at line %d.", __FILE__, __LINE__).c_str() );
            return 3 - ( getLocalEdgeIndex() + d.getLocalEdgeIndex() );
        }

        template<typename DartIteratorType>
        GlobalIndex getCommonNodeGlobalIndex(const DartIteratorType & d) const
        {
            if ( getGlobalNodeIndex() == d.getGlobalNodeIndex() || getGlobalNodeIndex() == d.getNextNodeGlobalIndex() )
                return getGlobalNodeIndex();
            if ( getNextNodeGlobalIndex() == d.getGlobalNodeIndex() || getNextNodeGlobalIndex() == d.getNextNodeGlobalIndex() )
                return getNextNodeGlobalIndex();
            throw std::invalid_argument( aol::strprintf("Dart d and *this have no node in common. In File %s at line %d.", __FILE__, __LINE__).c_str() );
            return -1;
        }

        void print() const
        {
            cerr << "triangle = " << _triangle << ", node = " << _node << ", edge = " << _edge << endl;
        }

        bool canFlipTriangle() const
        {
            return this->getNextTriangleIndex() != IndexNotSet ;
        }

        void flipTriangle()
        {
            if( !canFlipTriangle() )
                throw std::invalid_argument( aol::strprintf("Cannot flip triangle. In File %s at line %d.", __FILE__, __LINE__).c_str() );

            GlobalIndex newTriangleIndex = getNextTriangleIndex();
            GlobalIndex globalNodeIndex = _mesh.getTriangNodeIdx( _triangle, _node );
            GlobalIndex globalOtherNodesIndex = _mesh.getTriangNodeIdx( _triangle, 3 - (_node + _edge) );

            // find out which node of the new triangle is ours
            LocalIndex newNodeIndex = IndexNotSet;
            for (LocalIndex i = 0; i < 3; ++i) {
                if ( _mesh.getTriangNodeIdx( newTriangleIndex, i) == globalNodeIndex ) {
                    newNodeIndex = i;
                    break;
                }
            }

            // find out which node of the new triangle lies on our edge
            LocalIndex newOtherNodesIndex = IndexNotSet;
            for (LocalIndex i = 1; i < 3; ++i) {
                if ( _mesh.getTriangNodeIdx( newTriangleIndex,(newNodeIndex + i) % 3 ) == globalOtherNodesIndex ) {
                    newOtherNodesIndex = (newNodeIndex + i) % 3;
                    break;
                }
            }

            if ( min( newNodeIndex, newOtherNodesIndex ) < 0 )
                throw std::invalid_argument( aol::strprintf( "newTriang should neighbour of _triangle but has no common node .In File %s at line %d.", __FILE__, __LINE__ ).c_str() );

            LocalIndex newEdgeIndex = (-(newNodeIndex + newOtherNodesIndex)) % 3;
            if (newEdgeIndex < 0) newEdgeIndex += 3;

            // all information found, now set:
            _triangle = newTriangleIndex;
            _node = newNodeIndex;
            _edge = newEdgeIndex;
        }

        void flipNode()
        {
            _node = getNextNodeLocalIndex();
        }

        void flipEdge()
        {
            _edge = getNextEdgeLocalIndex();
        }
    };

protected:
    std::vector<bool> _markedForRefinement;
    std::map< int, ParentInformation<DomVecType> > _interpolationMap;

public:
    AdaptiveTriangMesh ( const string& fileName ) :
        TriangMesh<DataTypeContainer, TriangleType >( fileName ),
        _markedForRefinement( this->getNumTriangs(), false )
    {
        this->makeNeighbour();
    }

    const std::map< int, ParentInformation<DomVecType> > & getInterpolationMap ( ) const
    {
        return _interpolationMap;
    }

    void mark(int element)
    {
        _markedForRefinement[ element ] = true;
    }

    void markAll()
    {
        for( unsigned i = 0; i < _markedForRefinement.size(); i++ ) _markedForRefinement[i] = true;
    }

    void unmark( int element )
    {
        _markedForRefinement[ element ] = false;
    }

    void unmarkAll()
    {
        for( unsigned i = 0; i < _markedForRefinement.size(); i++ )
            _markedForRefinement[i] = false;
    }
    bool isMarkedForRefinement(int element) const
    {
        return _markedForRefinement[element];
    }

    int pushBackTriang(const Indices3DType newTriang)
    {
        _markedForRefinement.push_back( false );
        return TriangMesh<DataTypeContainer, TriangleType >::pushBackTriang( newTriang );
    }

    void refineMarkedTriangles()
    {
        const int numElements = this->getNumTriangs();
        for ( int elementIndex = 0; elementIndex < numElements; ++elementIndex  )
            if ( isMarkedForRefinement( elementIndex ) )  refine( elementIndex );
        this->makeOrientationConsistent();
        this->updateAllTriangles();
    }

    void prolongateLinearly(VectorType& function) const
    {
        VectorType oldFunction = function;
        int oldSize = function.size();
        function.resize( this->getNumVertices() );
        for( int i = 0; i < oldSize; ++i ) function[i] = oldFunction[i];
        typename std::map< int, ParentInformation<DomVecType> >::const_iterator iter;
        for( int i = oldSize; i < function.size(); ++i ) {
            iter = _interpolationMap.find( i );
            if( iter == _interpolationMap.end() ) throw std::invalid_argument ( aol::strprintf ( "unknown vertex in File %s at line %d.", __FILE__, __LINE__ ).c_str() );
            function[i] = 0.5 * ( function[iter->second.globalIndices[0]] + function[iter->second.globalIndices[1]] );
        }
    }

    LocalIndex getLongestEdgeIndex(GlobalIndex triangle) const
    {
        // get global node indices of triangle
        const Indices3DType &triangIdx ( this->getTriang( triangle ).getGlobalNodeIdx() );
        // assume that edge 0 is longest edge
        LocalIndex localIdx = 0;
        Point3DType edge = this->getVertex( triangIdx[1] );
        edge -= this->getVertex( triangIdx[2] );
        RealType maxLength = edge.squaredNorm();
        // now check if edge 1 is longer
        edge = this->getVertex( triangIdx[2] );
        edge -= this->getVertex( triangIdx[0] );
        RealType candidate = edge.squaredNorm();
        if( candidate > maxLength ) {
            localIdx = 1;
            maxLength = candidate;
        }
        // now check if edge 2 is longer
        edge = this->getVertex( triangIdx[0] );
        edge -= this->getVertex( triangIdx[1] );
        candidate = edge.squaredNorm();
        if( candidate > maxLength ) localIdx = 2;

        return localIdx;
    }

protected:
    // internal refinement functions
    GlobalIndex refine(GlobalIndex triangleToRefine)
    {
        DartIterator d( *this, triangleToRefine, getLongestEdgeIndex( triangleToRefine ) );
        return refine( d );
    }

    GlobalIndex refine(const DartIterator& d)
    {

        return 0;
    }

    GlobalIndex addEdgeMidpoint(const DartIterator& d)
    {
        DartIterator d1 = d;
        d1.flipNode();

        ParentInformation<DomVecType> parents;
        parents.globalIndices[0] = d.getGlobalNodeIndex();
        parents.globalIndices[1] = d1.getGlobalNodeIndex();
        parents.ElementIndex = d.getGlobalTriangleIndex();

        // get edge midpoint between nodes n1 and n2
        Point3DType newVertex;
        newVertex.setZero();
        newVertex += this->getVertex ( d.getGlobalNodeIndex() );
        newVertex += this->getVertex ( d1.getGlobalNodeIndex() );
        newVertex /= 2.;

        // add new vertex
        int newNodeIdx = this->pushBackVertex( newVertex );
        _interpolationMap.insert( std::pair< int, ParentInformation<DomVecType> >( newNodeIdx, parents ) );
        return newNodeIdx;
    }

};


} // end namespace

#endif