Mesh.h

Go to the documentation of this file.

 
#pragma once
 
#include <cstdlib>
#include <memory>
#include <string>
#include <vector>
 
#include "BaseLib/Error.h"
 
#include "MeshEnums.h"
#include "Properties.h"
 
namespace ApplicationUtils
{
    class NodeWiseMeshPartitioner;
}
 
namespace MeshLib
{
    class Node;
    class Element;
 
class Mesh
{
    /* friend functions: */
    friend void removeMeshNodes(Mesh& mesh,
                                const std::vector<std::size_t>& nodes);
 
    friend class ApplicationUtils::NodeWiseMeshPartitioner;
 
public:
    Mesh(std::string name,
         std::vector<Node*>
             nodes,
         std::vector<Element*>
             elements,
         Properties const& properties = Properties());
 
    Mesh(const Mesh &mesh);
 
    virtual ~Mesh();
 
    void addElement(Element* elem);
 
    unsigned getDimension() const { return _mesh_dimension; }
 
    const Node* getNode(std::size_t idx) const { return _nodes[idx]; }
 
    const Element* getElement(std::size_t idx) const { return _elements[idx]; }
 
    double getMinEdgeLength() const { return _edge_length.first; }
 
    double getMaxEdgeLength() const { return _edge_length.second; }
 
    std::size_t getNumberOfElements() const { return _elements.size(); }
 
    std::size_t getNumberOfNodes() const { return _nodes.size(); }
 
    const std::string getName() const { return _name; }
 
    std::vector<Node*> const& getNodes() const { return _nodes; }
 
    std::vector<Element*> const& getElements() const { return _elements; }
 
    void resetElementIDs();
 
    void resetNodeIDs();
 
    void setName(const std::string &name) { this->_name = name; }
 
    std::size_t getID() const {return _id; }
 
    std::size_t getNumberOfBaseNodes() const;
 
    bool hasNonlinearElement() const;
 
    std::vector<Element const*> const& getElementsConnectedToNode(
        std::size_t node_id) const;
    std::vector<Element const*> const& getElementsConnectedToNode(
        Node const& node) const;
 
    Properties& getProperties() { return _properties; }
    Properties const& getProperties() const { return _properties; }
 
    bool isAxiallySymmetric() const { return _is_axially_symmetric; }
    void setAxiallySymmetric(bool is_axial_symmetric) {
        _is_axially_symmetric = is_axial_symmetric;
    }
 
protected:
    void calcEdgeLengthRange();
 
    void setDimension();
 
    void setElementNeighbors();
 
    std::size_t const _id;
    unsigned _mesh_dimension;
    std::pair<double, double> _edge_length = {
        std::numeric_limits<double>::quiet_NaN(),
        std::numeric_limits<double>::quiet_NaN()};
    std::pair<double, double> _node_distance;
    std::string _name;
    std::vector<Node*> _nodes;
    std::vector<Element*> _elements;
    Properties _properties;
 
    std::vector<std::vector<Element const*>> _elements_connected_to_nodes;
 
    bool _is_axially_symmetric = false;
}; /* class */
 
std::vector<std::vector<Node*>> calculateNodesConnectedByElements(
    Mesh const& mesh);
 
inline bool operator==(Mesh const& a, Mesh const& b)
{
    return a.getID() == b.getID();
}
 
inline bool operator!=(Mesh const& a, Mesh const& b)
{
    return !(a == b);
}
 
void scaleMeshPropertyVector(Mesh& mesh,
                             std::string const& property_name,
                             double factor);
 
template <typename T>
void addPropertyToMesh(Mesh& mesh, std::string const& name,
                       MeshItemType item_type, std::size_t number_of_components,
                       std::vector<T> const& values)
{
    if (item_type == MeshItemType::Node)
    {
        if (mesh.getNumberOfNodes() != values.size() / number_of_components)
        {
            OGS_FATAL(
                "Error number of nodes ({:d}) does not match the number of "
                "tuples ({:d}).",
                mesh.getNumberOfNodes(), values.size() / number_of_components);
        }
    }
    if (item_type == MeshItemType::Cell)
    {
        if (mesh.getNumberOfElements() != values.size() / number_of_components)
        {
            OGS_FATAL(
                "Error number of elements ({:d}) does not match the number of "
                "tuples ({:d}).",
                mesh.getNumberOfElements(),
                values.size() / number_of_components);
        }
    }
 
    auto* const property = mesh.getProperties().createNewPropertyVector<T>(
        name, item_type, number_of_components);
    if (!property)
    {
        OGS_FATAL("Error while creating PropertyVector '{:s}'.", name);
    }
    property->reserve(values.size());
    std::copy(values.cbegin(), values.cend(), std::back_inserter(*property));
}
 
template <typename T>
PropertyVector<T>* getOrCreateMeshProperty(Mesh& mesh,
                                           std::string const& property_name,
                                           MeshItemType const item_type,
                                           int const number_of_components)
{
    if (property_name.empty())
    {
        OGS_FATAL(
            "Trying to get or to create a mesh property with empty name.");
    }
 
    auto numberOfMeshItems = [&mesh, &item_type]() -> std::size_t {
        switch (item_type)
        {
            case MeshItemType::Cell:
                return mesh.getNumberOfElements();
            case MeshItemType::Node:
                return mesh.getNumberOfNodes();
            case MeshItemType::IntegrationPoint:
                return 0;  // For the integration point data the size is
                           // variable
            default:
                OGS_FATAL(
                    "getOrCreateMeshProperty cannot handle other "
                    "types than Node, Cell, or IntegrationPoint.");
        }
        return 0;
    };
 
    if (mesh.getProperties().existsPropertyVector<T>(property_name))
    {
        auto result =
            mesh.getProperties().template getPropertyVector<T>(property_name);
        assert(result);
        if (item_type != MeshItemType::IntegrationPoint)
        {
            // Test the size if number of mesh items is known, which is not the
            // case for the integration point data.
            assert(result->size() ==
                   numberOfMeshItems() * number_of_components);
        }
        return result;
    }
 
    auto result = mesh.getProperties().template createNewPropertyVector<T>(
        property_name, item_type, number_of_components);
    assert(result);
    result->resize(numberOfMeshItems() * number_of_components);
    return result;
}
 
PropertyVector<int> const* materialIDs(Mesh const& mesh);
 
std::unique_ptr<Mesh> createMeshFromElementSelection(
    std::string mesh_name, std::vector<Element*> const& elements);
 
bool isBaseNode(Node const& node,
                std::vector<Element const*> const& elements_connected_to_node);
}  // namespace MeshLib