swh:1:snp:6088ab52ef49920e01e3f334cdf4d5d6c8a822b9
Raw File
Tip revision: 805d917b35b70ff779a3963230bcdd3eb289525d authored by Lars Bilke on 25 March 2021, 07:30:19 UTC
[cmake] Replace file(REAL_PATH ..) with get_filename_component.
Tip revision: 805d917
GmshReader.cpp
/**
 *
 * \copyright
 * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 */

#include "GmshReader.h"

#include "MeshLib/Elements/Element.h"
#include "MeshLib/Elements/Line.h"
#include "MeshLib/Elements/Tri.h"
#include "MeshLib/Elements/Quad.h"
#include "MeshLib/Elements/Tet.h"
#include "MeshLib/Elements/Hex.h"
#include "MeshLib/Elements/Prism.h"
#include "MeshLib/Elements/Pyramid.h"

#include "MeshLib/Node.h"
#include "MeshLib/Mesh.h"

#include "MeshLib/MeshEditing/ElementValueModification.h"

#include "BaseLib/FileTools.h"

#include <fstream>
#include <map>
#include <vector>

namespace FileIO
{
namespace GMSH
{

bool isGMSHMeshFile(const std::string& fname)
{
    std::ifstream input(fname.c_str());

    if (!input) {
        ERR("isGMSHMeshFile(): Could not open file {:s}.", fname);
        return false;
    }

    std::string header_first_line;
    input >> header_first_line;
    if (header_first_line.find("$MeshFormat") != std::string::npos) {
        // read version
        std::string version;
        getline(input, version);
        getline(input, version);
        INFO("isGMSHMeshFile(): Found GMSH mesh file version: {:s}.", version);
        input.close();
        return true;
    }

    return false;
}

void readNodeIDs(std::ifstream& in, unsigned n_nodes,
                 std::vector<unsigned>& node_ids,
                 std::map<unsigned, unsigned> const& id_map)
{
    unsigned idx;
    for (unsigned i = 0; i < n_nodes; i++)
    {
        in >> idx;
        node_ids.push_back(id_map.at(idx));
    }
}

std::pair<MeshLib::Element*, int> readElement(
    std::ifstream& in, std::vector<MeshLib::Node*> const& nodes,
    std::map<unsigned, unsigned> const& id_map)
{
    unsigned idx;
    unsigned type;
    unsigned n_tags;
    unsigned dummy;
    int mat_id;
    std::vector<unsigned> node_ids;

    //element format is structured like this:
    //element-id element-type n-tags physical-entity elementary entity node-ids
    in >> idx >> type >> n_tags >> dummy >> mat_id;

    switch (type)
    {
    case 1: {
        readNodeIDs(in, 2, node_ids, id_map);
        // edge_nodes array will be deleted from Line object
        auto edge_nodes = new MeshLib::Node*[2];
        edge_nodes[0] = nodes[node_ids[0]];
        edge_nodes[1] = nodes[node_ids[1]];
        return std::make_pair(new MeshLib::Line(edge_nodes), mat_id);
    }
    case 2: {
        readNodeIDs(in, 3, node_ids, id_map);
        auto tri_nodes = new MeshLib::Node*[3];
        tri_nodes[0] = nodes[node_ids[2]];
        tri_nodes[1] = nodes[node_ids[1]];
        tri_nodes[2] = nodes[node_ids[0]];
        return std::make_pair(new MeshLib::Tri(tri_nodes), mat_id);
    }
    case 3: {
        readNodeIDs(in, 4, node_ids, id_map);
        auto quad_nodes = new MeshLib::Node*[4];
        for (unsigned k(0); k < 4; k++)
        {
            quad_nodes[k] = nodes[node_ids[k]];
        }
        return std::make_pair(new MeshLib::Quad(quad_nodes), mat_id);
    }
    case 4: {
        readNodeIDs(in, 4, node_ids, id_map);
        auto tet_nodes = new MeshLib::Node*[5];
        for (unsigned k(0); k < 4; k++)
        {
            tet_nodes[k] = nodes[node_ids[k]];
        }
        return std::make_pair(new MeshLib::Tet(tet_nodes), mat_id);
    }
    case 5: {
        readNodeIDs(in, 8, node_ids, id_map);
        auto hex_nodes = new MeshLib::Node*[8];
        for (unsigned k(0); k < 8; k++)
        {
            hex_nodes[k] = nodes[node_ids[k]];
        }
        return std::make_pair(new MeshLib::Hex(hex_nodes), mat_id);
    }
    case 6: {
        readNodeIDs(in, 6, node_ids, id_map);
        auto prism_nodes = new MeshLib::Node*[6];
        for (unsigned k(0); k < 6; k++)
        {
            prism_nodes[k] = nodes[node_ids[k]];
        }
        return std::make_pair(new MeshLib::Prism(prism_nodes), mat_id);
    }
    case 7: {
        readNodeIDs(in, 5, node_ids, id_map);
        auto pyramid_nodes = new MeshLib::Node*[5];
        for (unsigned k(0); k < 5; k++)
        {
            pyramid_nodes[k] = nodes[node_ids[k]];
        }
        return std::make_pair(new MeshLib::Pyramid(pyramid_nodes), mat_id);
    }
    case 15:
        in >> dummy; // skip rest of line
        break;
    default:
        WARN("readGMSHMesh(): Unknown element type {:d}.", type);
        break;
    }
    return std::make_pair(nullptr, -1);
}

MeshLib::Mesh* readGMSHMesh(std::string const& fname)
{
    std::string line;
    std::ifstream in(fname.c_str(), std::ios::in);
    if (!in.is_open())
    {
        WARN("readGMSHMesh() - Could not open file {:s}.", fname);
        return nullptr;
    }

    getline(in, line); // $MeshFormat keyword
    if (line.find("$MeshFormat") == std::string::npos)
    {
        in.close();
        WARN ("No GMSH file format recognized.");
        return nullptr;
    }

    getline(in, line); // version-number file-type data-size
    if (line.substr(0, 3) != "2.2")
    {
        WARN("Wrong gmsh file format version '{:s}'.", line.substr(0, 3));
        return nullptr;
    }

    if (line[4] != '0')
    {
        WARN("Currently reading gmsh binary file type is not supported.");
        return nullptr;
    }
    getline(in, line); //$EndMeshFormat

    std::vector<MeshLib::Node*> nodes;
    std::vector<MeshLib::Element*> elements;
    std::vector<int> materials;
    std::map<unsigned, unsigned> id_map;
    while (line.find("$EndElements") == std::string::npos)
    {
        // Node data
        getline(in, line); //$Nodes Keywords
        if (line.find("$Nodes") != std::string::npos)
        {
            std::size_t n_nodes(0);
            long id;
            double x;
            double y;
            double z;
            in >> n_nodes >> std::ws;
            nodes.resize(n_nodes);
            for (std::size_t i = 0; i < n_nodes; i++) {
                in >> id >> x >> y >> z >> std::ws;
                id_map.insert(std::map<unsigned, unsigned>::value_type(id, i));
                nodes[i] = new MeshLib::Node(x,y,z,id);
            }
            getline(in, line); // End Node keyword $EndNodes
        }

        // Element data
        if (line.find("$Elements") != std::string::npos)
        {
            std::size_t n_elements(0);
            if (! (in >> n_elements >> std::ws)) { // number-of-elements
                ERR("Read GMSH mesh does not contain any elements");
            }
            elements.reserve(n_elements);
            materials.reserve(n_elements);
            for (std::size_t i = 0; i < n_elements; i++)
            {
                MeshLib::Element* elem(nullptr);
                int mat_id(0);
                std::tie(elem, mat_id) = readElement(in, nodes, id_map);

                if (elem) {
                    elements.push_back(elem);
                    materials.push_back(mat_id);
                }
            }
            getline(in, line); // END keyword
        }

        if (line.find("PhysicalNames") != std::string::npos)
        {
            std::size_t n_lines(0);
            in >> n_lines >> std::ws; // number-of-lines
            for (std::size_t i = 0; i < n_lines; i++)
            {
                getline(in, line);
            }
            getline(in, line); // END keyword
        }
    }
    in.close();
    if (elements.empty()) {
        for (auto& node : nodes)
        {
            delete node;
        }
        return nullptr;
    }

    MeshLib::Mesh * mesh(new MeshLib::Mesh(
        BaseLib::extractBaseNameWithoutExtension(fname), nodes, elements));

    auto* const material_ids =
        mesh->getProperties().createNewPropertyVector<int>(
            "MaterialIDs", MeshLib::MeshItemType::Cell, 1);
    if (!material_ids)
    {
        WARN("Could not create PropertyVector for MaterialIDs in Mesh.");
    }
    else
    {
        material_ids->insert(material_ids->end(), materials.cbegin(),
                             materials.cend());
    }

    MeshLib::ElementValueModification::condense(*mesh);

    INFO("\t... finished.");
    INFO("Nr. Nodes: {:d}.", nodes.size());
    INFO("Nr. Elements: {:d}.", elements.size());

    return mesh;
}

} // end namespace GMSH
} // end namespace FileIO
back to top