AbstractTetrahedralMesh.hpp

/*

Copyright (C) University of Oxford, 2005-2010

University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

Chaste is free software: you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation, either version 2.1 of the License, or
(at your option) any later version.

Chaste 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 Lesser General Public
License for more details. The offer of Chaste under the terms of the
License is subject to the License being interpreted in accordance with
English Law and subject to any action against the University of Oxford
being under the jurisdiction of the English Courts.

You should have received a copy of the GNU Lesser General Public License
along with Chaste. If not, see <http://www.gnu.org/licenses/>.

*/

#ifndef ABSTRACTTETRAHEDRALMESH_HPP_
#define ABSTRACTTETRAHEDRALMESH_HPP_

#include "ChasteSerialization.hpp"
#include "ClassIsAbstract.hpp"
#include <boost/serialization/base_object.hpp>

#include <vector>
#include <string>
#include <cassert>

#include "AbstractMesh.hpp"
#include "BoundaryElement.hpp"
#include "Element.hpp"
#include "AbstractMeshReader.hpp"
#include "TrianglesMeshReader.hpp"
#include "TrianglesMeshWriter.hpp"
#include "ArchiveLocationInfo.hpp"

#include <boost/serialization/split_member.hpp>

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
class AbstractTetrahedralMesh : public AbstractMesh<ELEMENT_DIM, SPACE_DIM>
{
protected:
    bool mMeshIsLinear;

private:
    virtual unsigned SolveElementMapping(unsigned index) const = 0;

    virtual unsigned SolveBoundaryElementMapping(unsigned index) const = 0;

    friend class boost::serialization::access;

    template<class Archive>
    void save(Archive & archive, const unsigned int version) const
    {
        archive & boost::serialization::base_object<AbstractMesh<ELEMENT_DIM,SPACE_DIM> >(*this);
        archive & mMeshIsLinear;
        // Create a mesh writer pointing to the correct file and directory.
        TrianglesMeshWriter<ELEMENT_DIM,SPACE_DIM> mesh_writer(ArchiveLocationInfo::GetArchiveRelativePath(),
                                                               ArchiveLocationInfo::GetMeshFilename(),
                                                               false);
        //Binary meshes have similar content to the original Triangle/Tetgen format, but take up less space on disk
        mesh_writer.SetWriteFilesAsBinary();
        mesh_writer.WriteFilesUsingMesh(*(const_cast<AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>*>(this)));
        // Make sure that the files are written before slave processes proceed
        PetscTools::Barrier("AbstractTetrahedralMesh::save");
    }

    template<class Archive>
    void load(Archive & archive, const unsigned int version)
    {
        archive & boost::serialization::base_object<AbstractMesh<ELEMENT_DIM,SPACE_DIM> >(*this);
        archive & mMeshIsLinear;

        // Store the DistributedVectorFactory loaded from the archive
        DistributedVectorFactory* p_factory = this->mpDistributedVectorFactory;
        this->mpDistributedVectorFactory = NULL;
        // Check whether we're migrating, or if we can use the original partition for the mesh
        DistributedVectorFactory* p_our_factory = NULL;
        if (p_factory)
        {
            p_our_factory = p_factory->GetOriginalFactory();
        }
        if (p_our_factory && p_our_factory->GetNumProcs() == p_factory->GetNumProcs())
        {
            // Specify the node distribution
            this->SetDistributedVectorFactory(p_our_factory);
        }
        else
        {
            // Migrating; let the mesh re-partition if it likes
            p_our_factory = NULL;
        }

        if (mMeshIsLinear)
        {
            //I am a linear mesh
            TrianglesMeshReader<ELEMENT_DIM,SPACE_DIM> mesh_reader(ArchiveLocationInfo::GetArchiveDirectory() + ArchiveLocationInfo::GetMeshFilename());
            this->ConstructFromMeshReader(mesh_reader);
        }
        else
        {
            //I am a quadratic mesh and need quadratic information from the reader
            TrianglesMeshReader<ELEMENT_DIM,SPACE_DIM> mesh_reader(ArchiveLocationInfo::GetArchiveDirectory() + ArchiveLocationInfo::GetMeshFilename(), 2, 2);
            this->ConstructFromMeshReader(mesh_reader);
        }

        // Make sure we're using the correct vector factory
        if (p_factory)
        {
            if (!this->mpDistributedVectorFactory)
            {
                // If we're not using a DistributedTetrahedralMesh, ConstructFromMeshReader won't set
                // this->mpDistributedVectorFactory.
                this->mpDistributedVectorFactory = p_factory;
            }
            else
            {
                // We need to update p_factory to match this->mpDistributedVectorFactory, and then use
                // p_factory, since the rest of the code (e.g. AbstractCardiacPde) will be using p_factory.
                p_factory->SetFromFactory(this->mpDistributedVectorFactory);
                if (p_our_factory != this->mpDistributedVectorFactory)
                {
                    // Avoid memory leak
                    delete this->mpDistributedVectorFactory;
                }
                this->mpDistributedVectorFactory = p_factory;
            }
        }
    }
    BOOST_SERIALIZATION_SPLIT_MEMBER()


protected:  // Give access of these variables to subclasses

    std::vector<Element<ELEMENT_DIM, SPACE_DIM> *> mElements;

    std::vector<BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> *> mBoundaryElements;

public:

    //                            Iterators                             //

    typedef typename std::vector<BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> *>::const_iterator BoundaryElementIterator;

    class ElementIterator;

    inline ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true);

    inline ElementIterator GetElementIteratorEnd();

    //                             Methods                              //

    AbstractTetrahedralMesh();

    virtual ~AbstractTetrahedralMesh();


    virtual unsigned GetNumElements() const;

    virtual unsigned GetNumBoundaryElements() const;

    unsigned GetNumAllElements() const;

    unsigned GetNumAllBoundaryElements() const;

    Element<ELEMENT_DIM, SPACE_DIM>* GetElement(unsigned index) const;

    BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* GetBoundaryElement(unsigned index) const;

    virtual void SetElementOwnerships(unsigned lo, unsigned hi);

    virtual void ConstructFromMeshReader(AbstractMeshReader<ELEMENT_DIM, SPACE_DIM>& rMeshReader)=0;

    BoundaryElementIterator GetBoundaryElementIteratorBegin() const;

    BoundaryElementIterator GetBoundaryElementIteratorEnd() const;

    virtual void GetInverseJacobianForElement(unsigned elementIndex, c_matrix<double, SPACE_DIM, ELEMENT_DIM>& rJacobian,
                                              double& rJacobianDeterminant,
                                              c_matrix<double, ELEMENT_DIM, SPACE_DIM>& rInverseJacobian) const;

    virtual void GetWeightedDirectionForBoundaryElement(unsigned elementIndex,
                                                        c_vector<double, SPACE_DIM>& rWeightedDirection,
                                                        double& rJacobianDeterminant) const;



     virtual void ConstructLinearMesh(unsigned width);

    virtual void ConstructRectangularMesh(unsigned width, unsigned height, bool stagger=true);

    virtual void ConstructCuboid(unsigned width, unsigned height, unsigned depth);

    virtual bool CalculateDesignatedOwnershipOfBoundaryElement( unsigned faceIndex );

    virtual bool CalculateDesignatedOwnershipOfElement( unsigned elementIndex );


    //                         Nested classes                           //


    class ElementIterator
    {
    public:
        inline Element<ELEMENT_DIM, SPACE_DIM>& operator*();

        inline Element<ELEMENT_DIM, SPACE_DIM>* operator->();

        inline bool operator!=(const AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator& rOther);

        inline ElementIterator& operator++();

        ElementIterator(AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>& rMesh,
                        typename std::vector<Element<ELEMENT_DIM, SPACE_DIM> *>::iterator elementIter,
                        bool skipDeletedElements=true);

    private:
        AbstractTetrahedralMesh& mrMesh;

        typename std::vector<Element<ELEMENT_DIM, SPACE_DIM> *>::iterator mElementIter;

        bool mSkipDeletedElements;

        inline bool IsAtEnd();

        inline bool IsAllowedElement();
    };

};

TEMPLATED_CLASS_IS_ABSTRACT_2_UNSIGNED(AbstractTetrahedralMesh);

//      ElementIterator class implementation - most methods are inlined     //

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetElementIteratorBegin(
        bool skipDeletedElements)
{
    return ElementIterator(*this, mElements.begin(), skipDeletedElements);
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::GetElementIteratorEnd()
{
    return ElementIterator(*this, mElements.end());
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
Element<ELEMENT_DIM, SPACE_DIM>& AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator::operator*()
{
    assert(!IsAtEnd());
    return **mElementIter;
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
Element<ELEMENT_DIM, SPACE_DIM>* AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator::operator->()
{
    assert(!IsAtEnd());
    return *mElementIter;
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
bool AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator::operator!=(const AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator& rOther)
{
    return mElementIter != rOther.mElementIter;
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator& AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator::operator++()
{
    do
    {
        ++mElementIter;
    }
    while (!IsAtEnd() && !IsAllowedElement());

    return (*this);
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator::ElementIterator(
        AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>& rMesh,
        typename std::vector<Element<ELEMENT_DIM, SPACE_DIM> *>::iterator elementIter,
        bool skipDeletedElements)
    : mrMesh(rMesh),
      mElementIter(elementIter),
      mSkipDeletedElements(skipDeletedElements)
{
    if (mrMesh.mElements.size() == 0)
    {
        // Cope with empty meshes
        mElementIter = mrMesh.mElements.end();
    }
    else
    {
        // Make sure we start at an allowed element
        if (mElementIter == mrMesh.mElements.begin() && !IsAllowedElement())
        {
            ++(*this);
        }
    }
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
bool AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator::IsAtEnd()
{
    return mElementIter == mrMesh.mElements.end();
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
bool AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator::IsAllowedElement()
{
    return !(mSkipDeletedElements && (*this)->IsDeleted());
}


#endif /*ABSTRACTTETRAHEDRALMESH_HPP_*/