AbstractTetrahedralMesh.hpp

/*

Copyright (C) University of Oxford, 2005-2011

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 "ChasteSerializationVersion.hpp"
#include <boost/serialization/vector.hpp>
#include <boost/serialization/base_object.hpp>
#include <boost/serialization/split_member.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 "GenericMeshReader.hpp"


template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
class AbstractConductivityTensors;

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
class AbstractTetrahedralMesh : public AbstractMesh<ELEMENT_DIM, SPACE_DIM>
{
    friend class AbstractConductivityTensors<ELEMENT_DIM, SPACE_DIM>; //A class which needs a global to local element mapping

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();

        // Archive the mesh permutation, so we can just copy the original mesh files whenever possible
        bool permutation_available = (this->rGetNodePermutation().size() != 0);
        archive & permutation_available;

        if (permutation_available)
        {
            const std::vector<unsigned>& rPermutation = this->rGetNodePermutation();
            archive & rPermutation;
        }

        if (!this->IsMeshOnDisk())
        {
            mesh_writer.WriteFilesUsingMesh(*(const_cast<AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>*>(this)));
        }
        else
        {
            unsigned order_of_element = (mMeshIsLinear?1:2);
            unsigned& order_of_boundary_element = order_of_element;

            // Mesh in disc, copy it to the archiving folder
            std::string original_file=this->GetMeshFileBaseName();
            std::auto_ptr<AbstractMeshReader<ELEMENT_DIM, SPACE_DIM> > p_original_mesh_reader
                = GenericMeshReader<ELEMENT_DIM, SPACE_DIM>(original_file, order_of_element, order_of_boundary_element);

            if (p_original_mesh_reader->IsFileFormatBinary())
            {
                // Mesh is in binary format, we can just copy the files across ignoring the mesh reader
                std::stringstream cp_command;
                cp_command << "for filename in `ls " << this->GetMeshFileBaseName() << ".*`; do " <<
                                 "extension=${filename##*.};" << // ##*. deletes basename and . from the filename (i.e. leaves only the extension (http://tldp.org/LDP/Bash-Beginners-Guide/html/sect_10_03.html)
                                 "cp $filename " << ArchiveLocationInfo::GetArchiveDirectory() << ArchiveLocationInfo::GetMeshFilename() <<".$extension;" <<
                              "done";

                if (PetscTools::AmMaster())
                {
                    ABORT_IF_NON0(system, cp_command.str());
                }
            }
            else
            {
                // Mesh in text format, use the mesh writer to "binarise" it
                mesh_writer.WriteFilesUsingMeshReaderAndMesh(*p_original_mesh_reader,
                                                             *(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;

        bool permutation_available=false;
        std::vector<unsigned> permutation;

        if (version > 0)
        {
            archive & permutation_available;

            if (permutation_available)
            {
                archive & permutation;
            }
        }

        // 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());

            if (permutation_available)
            {
                mesh_reader.SetNodePermutation(permutation);
            }

            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;

    void SetElementOwnerships();

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 GetNumLocalElements() const;

    virtual unsigned GetNumBoundaryElements() const;

    unsigned GetNumAllElements() const;

    unsigned GetNumAllBoundaryElements() const;

    virtual unsigned GetNumCableElements() const;

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

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


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

    void ConstructFromMesh(AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>& rOtherMesh);


    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;

    void CheckOutwardNormals();

    virtual void ConstructLinearMesh(unsigned width);



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



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



    void ConstructRegularSlabMesh(double spaceStep, double width, double height=0, double depth=0);



    virtual bool CalculateDesignatedOwnershipOfBoundaryElement( unsigned faceIndex );

    virtual bool CalculateDesignatedOwnershipOfElement( unsigned elementIndex );

    unsigned CalculateMaximumNodeConnectivityPerProcess() const;

    virtual void GetHaloNodeIndices(std::vector<unsigned>& rHaloIndices) const;

     void CalculateNodeExchange( std::vector<std::vector<unsigned> >& rNodesToSendPerProcess,
                                 std::vector<std::vector<unsigned> >& rNodesToReceivePerProcess);


    //                         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)

namespace boost {
namespace serialization {
template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
struct version<AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM> >
{
    CHASTE_VERSION_CONTENT(1);
};
} // namespace serialization
} // namespace boost


//      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_*/