mesh.h

Go to the documentation of this file.

/*!
 *
 * Copyright (C) 2015 Technical University of Liberec.  All rights reserved.
 * 
 * This program is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License version 3 as published by the
 * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
 * 
 * This program 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 General Public License for more details.
 *
 * 
 * @file    mesh.h
 * @brief   
 */
 
#ifndef MAKE_MESH_H
#define MAKE_MESH_H
 
#include <mpi.h>                             // for MPI_Comm, MPI_COMM_WORLD
 
//#include <boost/range.hpp>
#include <memory>                            // for shared_ptr
#include <string>                            // for string
#include <vector>                            // for vector, vector<>::iterator
#include "input/accessors.hh"                // for Record, Array (ptr only)
#include "input/accessors_impl.hh"           // for Record::val
#include "input/storage.hh"                  // for ExcStorageTypeMismatch
#include "input/type_record.hh"              // for Record (ptr only), Recor...
#include "mesh/region.hh"                    // for RegionDB, RegionDB::MapE...
#include "mesh/elements.h"
#include "mesh/bounding_box.hh"              // for BoundingBox
#include "mesh/range_wrapper.hh"
#include "mesh/mesh_data.hh"
#include "tools/bidirectional_map.hh"
#include "tools/general_iterator.hh"
#include "system/index_types.hh"             // for LongIdx
#include "system/exceptions.hh"              // for operator<<, ExcStream, EI
#include "system/file_path.hh"               // for FilePath
#include "system/armor.hh"
 
 
class BIHTree;
class Distribution;
class Partitioning;
class MixedMeshIntersections;
class Neighbour;
class SideIter;
class Boundary;
class Edge;
class BCMesh;
class DuplicateNodes;
template <int spacedim> class ElementAccessor;
template <int spacedim> class NodeAccessor;
 
 
 
#define ELM  0
#define BC  1
#define NODE  2
 
/**
 *  This parameter limits volume of elements from below.
 */
#define MESH_CRITICAL_VOLUME 1.0E-12
 
class BoundarySegment {
public:
    static Input::Type::Record input_type;
};
 
 
/** Auxiliary structure that keeps the separate
 * element maps (bulk and boundary) for reading mesh and elementwise data.
 * The mappings are form computational mesh to the other (field fe data) mesh.
 * Used by @p check_compatible_mesh().
 */
struct EquivalentMeshMap{
    std::vector<LongIdx> bulk;
    std::vector<LongIdx> boundary;
 
    EquivalentMeshMap(){}
    
    EquivalentMeshMap(unsigned int bulk_size, unsigned int boundary_size, LongIdx def_val)
    : bulk(bulk_size, def_val),
      boundary(boundary_size, def_val)
    {}
 
    bool empty()
    { return bulk.empty() && boundary.empty(); }
};
 
 
/// Base class for Mesh and BCMesh.
class MeshBase {
public:
 
    TYPEDEF_ERR_INFO( EI_ElemId, int);
    TYPEDEF_ERR_INFO( EI_ElemIdOther, int);
 
    DECLARE_EXCEPTION(ExcTooMatchingIds,
            << "Mesh: Duplicate dim-join lower dim elements: " << EI_ElemId::val << ", " << EI_ElemIdOther::val << ".\n" );
 
    static const unsigned int undef_idx=-1;
 
    MeshBase();
 
    virtual ~MeshBase();
 
    inline unsigned int n_elements() const
    { return element_vec_.size(); }
 
    inline unsigned int n_edges() const
    { return edges.size(); }
 
    unsigned int n_vb_neighbours() const;
 
    inline Distribution *get_el_ds() const
    { return el_ds; }
 
    inline LongIdx *get_el_4_loc() const
    { return el_4_loc; }
 
    inline LongIdx *get_row_4_el() const
    { return row_4_el; }
 
    const Element &element(unsigned idx) const
    { return element_vec_[idx]; }
 
    /// Return edge with given index.
    Edge edge(uint edge_idx) const;
 
    /// Return neighbour with given index.
    const Neighbour &vb_neighbour(unsigned int nb) const;
 
    /// Return element id (in GMSH file) of element of given position in element vector.
    int find_elem_id(unsigned int pos) const
    { return element_ids_[pos]; }
 
    /**
     * Returns maximal number of sides of one edge, which connects elements of dimension @p dim.
     * @param dim Dimension of elements sharing the edge.
     */
    unsigned int max_edge_sides(unsigned int dim) const { return max_edge_sides_[dim-1]; }
 
    const DuplicateNodes *duplicate_nodes() const
    { return duplicate_nodes_; }
 
    /**
     * Returns nodes_elements vector, if doesn't exist creates its.
     */
    vector<vector<unsigned int> > const & node_elements();
 
    /// Return node id (in GMSH file) of node of given position in node vector.
    inline int find_node_id(unsigned int pos) const
    {
        return (*node_ids_)[pos];
    }
 
    /// Check if given index is in element_vec_
    inline void check_element_size(unsigned int elem_idx) const
    {
        ASSERT(elem_idx < element_vec_.size())(elem_idx)(element_vec_.size()).error("Index of element is out of bound of element vector!");
    }
 
 
    const std::vector<unsigned int> &get_side_nodes(unsigned int dim, unsigned int side) const
    { return side_nodes[dim][side]; }
 
    inline unsigned int n_nodes() const {
        return nodes_->size();
    }
 
    inline const RegionDB &region_db() const {
        return *region_db_;
    }
 
    /// Create and return NodeAccessor to node of given idx
    NodeAccessor<3> node(unsigned int idx) const;
 
    /// Create and return ElementAccessor to element of given idx
    ElementAccessor<3> element_accessor(unsigned int idx) const;
 
 
    // TODO: have also private non-const accessors and ranges
 
    /// Returns range of mesh elements
    Range<ElementAccessor<3>> elements_range() const;
 
    /// Returns range of nodes
    Range<NodeAccessor<3>> node_range() const;
 
    /// Return range of edges.
    Range<Edge> edge_range() const;
 
 
 
    virtual Boundary boundary(uint edge_idx) const = 0;
    virtual BCMesh *bc_mesh() const = 0;
 
 
    /**
     * Returns pointer to partitioning object. Partitioning is created during setup_topology.
     */
    virtual Partitioning *get_part() = 0;
    virtual const LongIdx *get_local_part() = 0;
 
    /*
     * Check if nodes and elements are compatible with \p mesh.
     */
    virtual std::shared_ptr<EquivalentMeshMap> check_compatible_mesh( Mesh & input_mesh ) = 0;
 
    /**
     * Find intersection of element lists given by Mesh::node_elements_ for elements givne by @p nodes_list parameter.
     * The result is placed into vector @p intersection_element_list. If the @p node_list is empty, and empty intersection is
     * returned.
     */
    void intersect_element_lists(vector<unsigned int> const &nodes_list, vector<unsigned int> &intersection_element_list);
 
    /// For element of given elem_id returns index in element_vec_ or (-1) if element doesn't exist.
    inline int elem_index(int elem_id) const {
        return element_ids_.get_position(elem_id);
    }
 
        /// Initialize element_vec_, set size and reset counters of boundary and bulk elements.
    void init_element_vector(unsigned int size);
 
    /// Initialize node_vec_, set size
    void init_node_vector(unsigned int size);
 
    /// Return permutation vector of nodes
    inline const std::vector<unsigned int> &node_permutations() const
    { return node_permutation_; }
 
    /// Return permutation vector of elements
    inline const std::vector<unsigned int> &element_permutations() const
    { return elem_permutation_; }
 
 
    /// For each node the vector contains a list of elements that use this node
    vector<vector<unsigned int> > node_elements_;
 
    // Temporary solution for numbering of nodes on sides.
    // The data are defined in RefElement<dim>::side_nodes,
    // Mesh::side_nodes can be removed as soon as Element
    // is templated by dimension.
    //
    // side_nodes[dim][elm_side_idx][side_node_idx]
    // for every side dimension D = 0 .. 2
    // for every element side 0 .. D+1
    // for every side node 0 .. D
    // index into element node array
    vector< vector< vector<unsigned int> > > side_nodes;
 
protected:
 
    /**
     * Permute nodes of individual elements so that all elements have same edge orientations and aligned sides have same order of their nodes
     * Canonical edge orientation in elements and faces is from nodes of lower local index to higher local index.
     *
     * Algorithm detals:
     * 1. Orient all edges from lowe global node id to higher node id, fictional step. (substantial is orientation of yet non-oriented edges of a node in direction out of the node.
     *    Can be proven (!?) that this prevents edge cycles of the length 3 (faces with cyclic edges).
     * 2. Having all faces non-cyclic there exists a permutation of any element to the reference element.
     *    Pass through the elements. Sort nodes by global ID.
     */
    void canonical_faces();
 
    /**
     * Create element lists for nodes in Mesh::nodes_elements.
     */
    void create_node_element_lists();
 
    /// Adds element to mesh data structures (element_vec_, element_ids_), returns pointer to this element.
    Element * add_element_to_vector(int id, bool is_boundary = false);
 
    /**
     * Remove elements with dimension not equal to @p dim from @p element_list. Index of the first element of dimension @p dim-1,
     * is returned in @p element_idx. If no such element is found the method returns false, if one such element is found the method returns true,
     * if more elements are found we report an user input error.
     */
    bool find_lower_dim_element(vector<unsigned int> &element_list, unsigned int dim, unsigned int &element_idx);
 
    /**
     * Returns true if side @p si has same nodes as in the list @p side_nodes.
     */
    bool same_sides(const SideIter &si, vector<unsigned int> &side_nodes);
 
    /**
     * Vector of elements of the mesh.
     *
     * Store all elements of the mesh in order bulk elements - boundary elements
     */
    vector<Element> element_vec_;
 
    /// Maps element ids to indexes into vector element_vec_
    BidirectionalMap<int> element_ids_;
 
    /// Vector of MH edges, this should not be part of the geometrical mesh
    std::vector<EdgeData> edges;
 
    /// Vector of compatible neighbourings.
    vector<Neighbour> vb_neighbours_;
 
    /// Maximal number of sides per one edge in the actual mesh (set in make_neighbours_and_edges()).
    unsigned int max_edge_sides_[3];
 
    /**
     * Vector of nodes of the mesh.
     */
    shared_ptr<Armor::Array<double>> nodes_;
 
    /// Maps node ids to indexes into vector node_vec_
    shared_ptr<BidirectionalMap<int>> node_ids_;
 
    /// Vector of node permutations of optimized mesh (see class MeshOptimizer)
    std::vector<unsigned int> node_permutation_;
 
    /// Vector of element permutations of optimized mesh (see class MeshOptimizer)
    std::vector<unsigned int> elem_permutation_;
 
 
    ///
    std::array<std::array<uint, 4>, 64> element_nodes_original_;
 
 
    /// Index set assigning to global element index the local index used in parallel vectors.
    LongIdx *row_4_el;
    /// Index set assigning to local element index its global index.
    LongIdx *el_4_loc;
    /// Parallel distribution of elements.
    Distribution *el_ds;
 
 
    DuplicateNodes *duplicate_nodes_;
 
    /**
     * Database of regions (both bulk and boundary) of the mesh. Regions are logical parts of the
     * domain that allows setting of different data and boundary conditions on them.
     */
    std::shared_ptr<RegionDB> region_db_;
 
 
    friend class Edge;
    // friend class Side;
    // friend class RegionSetBase;
    friend class Element;
    // friend class BIHTree;
    // friend class Boundary;
    // friend class BCMesh;
    template <int spacedim> friend class ElementAccessor;
    template <int spacedim> friend class NodeAccessor;
 
};
 
//=============================================================================
// STRUCTURE OF THE MESH
//=============================================================================
 
class Mesh : public MeshBase {
public:
    TYPEDEF_ERR_INFO( EI_ElemLast, int);
    TYPEDEF_ERR_INFO( EI_ElemNew, int);
    TYPEDEF_ERR_INFO( EI_RegLast, std::string);
    TYPEDEF_ERR_INFO( EI_RegNew, std::string);
    TYPEDEF_ERR_INFO( EI_Region, std::string);
    TYPEDEF_ERR_INFO( EI_RegIdx, unsigned int);
    TYPEDEF_ERR_INFO( EI_Dim, unsigned int);
    TYPEDEF_ERR_INFO( EI_DimOther, unsigned int);
    TYPEDEF_ERR_INFO( EI_Quality, double);
 
 
    DECLARE_EXCEPTION(ExcDuplicateBoundary,
            << "Duplicate boundary elements! \n"
            << "Element id: " << EI_ElemLast::val << " on region name: " << EI_RegLast::val << "\n"
            << "Element id: " << EI_ElemNew::val << " on region name: " << EI_RegNew::val << "\n");
    DECLARE_EXCEPTION(ExcElmWrongOrder,
            << "Element IDs in non-increasing order, ID: " << EI_ElemId::val << "\n");
    DECLARE_EXCEPTION(ExcRegionElmDiffDim,
            << "User defined region " << EI_Region::qval << " (id " << EI_RegIdx::val
            << ") by 'From_Elements' cannot have elements of different dimensions.\n"
            << "Thrown due to: dim " << EI_Dim::val << " neq dim " << EI_DimOther::val << " (ele id " << EI_ElemId::val << ").\n"
            << "Split elements by dim, create separate regions and then possibly use Union.\n" );
    DECLARE_EXCEPTION(ExcBadElement,
            << "Extremely bad quality element ID=" << EI_ElemId::val << ",(" << EI_Quality::val << "<4*epsilon).\n");
    DECLARE_EXCEPTION(ExcBdrElemMatchRegular,
            << "Boundary element (id: " << EI_ElemId::val << ") match a regular element (id: " << EI_ElemIdOther::val << ") of lower dimension.\n" );
 
 
    /**
     * \brief Types of search algorithm for finding intersection candidates.
     */
    typedef enum IntersectionSearch {
        BIHsearch  = 1,
        BIHonly = 2,
        BBsearch = 3
    } IntersectionSearch;
    
    /**
     * \brief The definition of input record for selection of variant of file format
     */
    static const Input::Type::Selection & get_input_intersection_variant();
    
    static const Input::Type::Record & get_input_type();
 
 
    /** Labels for coordinate indexes in arma::vec3 representing vectors and points.*/
    enum {x_coord=0, y_coord=1, z_coord=2};
 
    /**
     * Empty constructor.
     *
     * Use only for unit tests!!!
     */
    Mesh();
    /**
     * Constructor from an input record.
     * Do not process input record. That is done in init_from_input.
     */
    Mesh(Input::Record in_record, MPI_Comm com = MPI_COMM_WORLD);
 
    Mesh(Mesh &other);
 
    /// Destructor.
    ~Mesh() override;
 
    inline unsigned int n_boundaries() const {
        return boundary_.size();
    }
 
    Boundary boundary(uint edge_idx) const override;
 
    unsigned int n_corners();
 
    Partitioning *get_part() override;
 
    const LongIdx *get_local_part() override;
 
    Distribution *get_node_ds() const
    { return node_ds_; }
 
    LongIdx *get_node_4_loc() const
    { return node_4_loc_; }
 
    unsigned int n_local_nodes() const
    { return n_local_nodes_; }
 
    /**
     * Returns MPI communicator of the mesh.
     */
    inline MPI_Comm get_comm() const { return comm_; }
 
 
    MixedMeshIntersections &mixed_intersections();
 
    unsigned int n_sides() const;
 
    /**
     * Reads mesh from stream.
     *
     * Method is especially used in unit tests.
     */
    void read_gmsh_from_stream(istream &in);
    /**
     * Reads input record, creates regions, read the mesh, setup topology. creates region sets.
     */
    void init_from_input();
 
    /**
     * Initialize all mesh structures from raw information about nodes and elements (including boundary elements).
     * Namely: create remaining boundary elements and Boundary objects, find edges and compatible neighborings.
     */
    void setup_topology();
    
    /**
     * Returns vector of ID numbers of elements, either bulk or bc elemnts.
     */
    void elements_id_maps( vector<LongIdx> & bulk_elements_id, vector<LongIdx> & boundary_elements_id) const;
 
    /*
     * Check if nodes and elements are compatible with continuous \p input_mesh.
     *
     * Call this method on computational mesh.
     * @param input_mesh data mesh of input fields
     * @return vector that holds mapping between eleemnts of data and computational meshes
     *             for every element in computational mesh hold idx of equivalent element in input mesh.
     *             If element doesn't exist in input mesh value is set to undef_idx.
     *             If meshes are not compatible returns empty vector.
     */
    std::shared_ptr<EquivalentMeshMap> check_compatible_mesh(Mesh & input_mesh) override;
 
 
    /**
     * Reads elements and their affiliation to regions and region sets defined by user in input file
     * Format of input record is defined in method RegionSetBase::get_input_type()
     *
     * @param region_list Array input AbstractRecords which define regions, region sets and elements
     */
    void read_regions_from_input(Input::Array region_list);
 
 
    /// Vector of boundary sides where is prescribed boundary condition.
    /// TODO: apply all boundary conditions in the main assembling cycle over elements and remove this Vector.
    mutable vector<BoundaryData> boundary_;
 
    //flow::VectorId<int> bcd_group_id; // gives a index of group for an id
 
    /**
     * Vector of individual intersections of two elements.
     * This is enough for local mortar.
     */
    std::shared_ptr<MixedMeshIntersections>  intersections;
 
    /**
     * For every element El we have vector of indices into @var intersections array for every intersection in which El is master element.
     * This is necessary for true mortar.
     */
    vector<vector<unsigned int> >  master_elements;
    
 
    int n_insides; // # of internal sides
    int n_exsides; // # of external sides
    mutable int n_sides_; // total number of sides (should be easy to count when we have separated dimensions
 
 
    /**
     * Check usage of regions, set regions to elements defined by user, close RegionDB
     */
    void check_and_finish();
    
    /// Compute bounding boxes of elements contained in mesh.
    std::vector<BoundingBox> get_element_boxes();
 
    /// Getter for BIH. Creates and compute BIH at first call.
    const BIHTree &get_bih_tree();\
 
    /// Add new node of given id and coordinates to mesh
    void add_node(unsigned int node_id, arma::vec3 coords);
 
    /// Add new element of given id to mesh
    void add_element(unsigned int elm_id, unsigned int dim, unsigned int region_id, unsigned int partition_id,
            std::vector<unsigned int> node_ids);
 
    /// Add new node of given id and coordinates to mesh
    void add_physical_name(unsigned int dim, unsigned int id, std::string name);
 
    /// Return FilePath object representing "mesh_file" input key
    inline FilePath mesh_file() {
        return in_record_.val<FilePath>("mesh_file");
    }
 
    /// Getter for input type selection for intersection search algorithm.
    IntersectionSearch get_intersection_search();
 
    /// Maximal distance of observe point from Mesh relative to its size
    double global_snap_radius() const;
 
    /// For node of given node_id returns index in element_vec_ or (-1) if node doesn't exist.
    inline int node_index(int node_id) const
    {
        return node_ids_->get_position(node_id);
    }
 
    /// Implement MeshBase::bc_mesh(), getter of boundary mesh
    BCMesh *bc_mesh() const override {
        return bc_mesh_;
    }
 
 
protected:
 
    /**
     * Part of the constructor whichdoes not depedn on input record.
     * Initializes node-side numbering according to RefElement.
     */
    void init();
 
    /**
     * Allow store boundary element data to temporary structure.
     *
     * We need this structure to preserve correct order of boundary elements.
     */
    struct ElementTmpData {
        /// Constructor
        ElementTmpData(unsigned int e_id, unsigned int dm, RegionIdx reg_idx, unsigned int part_id, std::vector<unsigned int> nodes)
        : elm_id(e_id), dim(dm), region_idx(reg_idx), partition_id(part_id), node_ids(nodes) {}
 
        unsigned int elm_id;
        unsigned int dim;
        RegionIdx region_idx;
        unsigned int partition_id;
        std::vector<unsigned int> node_ids;
    };
 
    /**
     *  This replaces read_neighbours() in order to avoid using NGH preprocessor.
     *
     *  TODO:
     *  - Avoid maps:
     *
     *    5) need not to have temporary array for Edges, only postpone setting pointers in elements and set them
     *       after edges are found; we can temporary save Edge index instead of pointer in Neigbours and elements
     *
     *    6) Try replace Edge * by indexes in Neigbours and elements (anyway we have mesh pointer in elements so it is accessible also from Neigbours)
     *
     */
    void make_neighbours_and_edges();
 
    void element_to_neigh_vb();
 
    void count_element_types();
    void count_side_types();
 
    /**
     * Check the element quality and remove unused nodes.
     */
    void check_mesh_on_read();
 
    /**
     * Possibly modify region id of elements sets by user in "regions" part of input file.
     *
     * TODO: This method needs check in issue 'Review mesh setting'.
     * Changes have been done during generalized region key and may be causing problems
     * during the further development.
     */
    void modify_element_ids(const RegionDB::MapElementIDToRegionID &map);
 
    /// Initialize element
    void init_element(Element *ele, unsigned int elm_id, unsigned int dim, RegionIdx region_idx, unsigned int partition_id,
            std::vector<unsigned int> node_ids);
 
    unsigned int n_bb_neigh, n_vb_neigh;
 
    /// Output of neighboring data into raw output.
    void output_internal_ngh_data();
    
    /**
     * Apply functionality of MeshOptimizer to sort nodes and elements.
     *
     * Use Hilbert curve, need call sort_permuted_nodes_elements method.
     */
    void optimize();
 
    /// Sort elements and nodes by order stored in permutation vectors.
    void sort_permuted_nodes_elements(std::vector<int> new_node_ids, std::vector<int> new_elem_ids);
 
    /**
     * Looks for the same (compatible) elements between the @p source_mesh and @p target_mesh.
     * Auxiliary function for check_compatible_mesh().
     * Uses the nodal mapping @p node_ids.
     * Fills the element mapping @p map.
     * Returns the number of compatible elements.
     */
    unsigned int check_compatible_elements(MeshBase* source_mesh, MeshBase* target_mesh,
                                           const std::vector<unsigned int>& node_ids,
                                           std::vector<LongIdx>& map);
 
    /**
     * Flag for optimization perfomed at the beginning of setup_topology.
     * Default true, can be set to flase by the optimize_mesh key of the input recoed.
     */
    bool optimize_memory_locality;
 
    /**
     * Mesh partitioning. Created in setup_topology.
     */
    std::shared_ptr<Partitioning> part_;
 
    /**
     * BIH Tree for intersection and observe points lookup.
     */
    std::shared_ptr<BIHTree> bih_tree_;
 
 
    /**
     * Accessor to the input record for the mesh.
     */
    Input::Record in_record_;
 
    /**
     * MPI communicator used for partitioning and ...
     */
    MPI_Comm comm_;
 
    friend class Edge;
    friend class Side;
    friend class RegionSetBase;
    friend class Element;
    friend class BIHTree;
    friend class Boundary;
    friend class BCMesh;
    template <int spacedim> friend class ElementAccessor;
    template <int spacedim> friend class NodeAccessor;
 
 
 
private:
 
    /// Fill array node_4_loc_ and create object node_ds_ according to element distribution.
    void distribute_nodes();
 
    /// Index set assigning to local node index its global index.
    LongIdx *node_4_loc_;
    /// Parallel distribution of nodes. Depends on elements distribution.
    Distribution *node_ds_;
    /// Hold number of local nodes (own + ghost), value is equal with size of node_4_loc array.
    unsigned int n_local_nodes_;
    /// Boundary mesh, object is created only if it's necessary
    BCMesh *bc_mesh_;
        
};
 
#endif
//-----------------------------------------------------------------------------
// vim: set cindent: