Added improved .msh reader/writer + tests + example (#1578)

* Added improved .msh reader/writer + tests + example,
with fixes proposed by @qnzhou
put append_mat_to_vec into igl::internal

* Update tutorial number.

* Changed tutorial code to use example from https://github.com/libigl/libigl-tutorial-data/pull/1 and flat shading

* Remove commented code.

* Rename tutorial number + update tutorial data.

Co-authored-by: Jérémie Dumas <[email protected]>
Co-authored-by: Jérémie Dumas <[email protected]>

cmake/LibiglDownloadExternal.cmake

@@ -184,7 +184,7 @@ function(igl_download_test_data)
 	igl_download_project_aux(test_data
 		"${LIBIGL_EXTERNAL}/../tests/data"
 		GIT_REPOSITORY https://github.com/libigl/libigl-tests-data
-		GIT_TAG        b5dddf45e329af685cd107e38770a28cfc18eb15
+		GIT_TAG        19cedf96d70702d8b3a83eb27934780c542356fe
 	)
 endfunction()
 
@@ -193,7 +193,7 @@ function(igl_download_tutorial_data)
 	igl_download_project_aux(tutorial_data
 		"${LIBIGL_EXTERNAL}/../tutorial/data"
 		GIT_REPOSITORY https://github.com/libigl/libigl-tutorial-data
-		GIT_TAG        fb5fa00bc4ede64b36002d703ce541552370b6e9
+		GIT_TAG        37d4e836054c9c2d2125a817c489ed8e07cd56fc
 	)
 endfunction()
 

include/igl/MshLoader.cpp

@@ -0,0 +1,497 @@
+// based on MSH reader from PyMesh 
+
+// Copyright (c) 2015 Qingnan Zhou <[email protected]>           
+// Copyright (C) 2020 Vladimir Fonov <[email protected]> 
+//
+// This Source Code Form is subject to the terms of the Mozilla 
+// Public License v. 2.0. If a copy of the MPL was not distributed 
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 
+
+#include "MshLoader.h"
+
+#include <cassert>
+#include <iostream>
+#include <sstream>
+#include <vector>
+
+#include <string.h>
+
+namespace igl {
+    // helper function
+    void inline _msh_eat_white_space(std::ifstream& fin) {
+        char next = fin.peek();
+        while (next == '\n' || next == ' ' || next == '\t' || next == '\r') {
+            fin.get();
+            next = fin.peek();
+        }
+    }
+}
+
+IGL_INLINE igl::MshLoader::MshLoader(const std::string &filename) {
+    std::ifstream fin(filename, std::ios::in | std::ios::binary);
+
+    if (!fin.is_open()) {
+        std::stringstream err_msg;
+        err_msg << "failed to open file \"" << filename << "\"";
+        throw std::ios_base::failure(err_msg.str());
+    }
+    // Parse header
+    std::string buf;
+    double version;
+    int type;
+    fin >> buf;
+    if (buf != "$MeshFormat") { throw std::runtime_error("Unexpected .msh format"); }
+
+    fin >> version >> type >> m_data_size;
+    m_binary = (type == 1);
+    if(version>2.2 || version<2.0)
+    {
+        // probably unsupported version
+        std::stringstream err_msg;
+        err_msg << "Error: Unsupported file version:" << version << std::endl;
+        throw std::runtime_error(err_msg.str());
+
+    }
+    // Some sanity check.
+    if (m_data_size != 8) {
+        std::stringstream err_msg;
+        err_msg << "Error: data size must be 8 bytes." << std::endl;
+        throw std::runtime_error(err_msg.str());
+    }
+    if (sizeof(int) != 4) {
+        std::stringstream err_msg;
+        err_msg << "Error: code must be compiled with int size 4 bytes." << std::endl;
+        throw std::runtime_error(err_msg.str());
+    }
+
+    // Read in extra info from binary header.
+    if (m_binary) {
+        int one;
+        igl::_msh_eat_white_space(fin);
+        fin.read(reinterpret_cast<char*>(&one), sizeof(int));
+        if (one != 1) {
+            std::stringstream err_msg;
+                err_msg << "Binary msh file " << filename
+                << " is saved with different endianness than this machine."
+                << std::endl;
+            throw std::runtime_error(err_msg.str());
+        }
+    }
+
+    fin >> buf;
+    if (buf != "$EndMeshFormat") 
+    { 
+        std::stringstream err_msg;
+        err_msg << "Unexpected contents in the file header." << std::endl;
+        throw std::runtime_error(err_msg.str());
+    }
+
+    while (!fin.eof()) {
+        buf.clear();
+        fin >> buf;
+        if (buf == "$Nodes") {
+            parse_nodes(fin);
+            fin >> buf;
+            if (buf != "$EndNodes") { throw std::runtime_error("Unexpected tag"); }
+        } else if (buf == "$Elements") {
+            parse_elements(fin);
+            fin >> buf;
+            if (buf != "$EndElements") { throw std::runtime_error("Unexpected tag"); }
+        } else if (buf == "$NodeData") {
+            parse_node_field(fin);
+            fin >> buf;
+            if (buf != "$EndNodeData") { throw std::runtime_error("Unexpected tag"); }
+        } else if (buf == "$ElementData") {
+            parse_element_field(fin);
+            fin >> buf;
+            if (buf != "$EndElementData") { throw std::runtime_error("Unexpected tag"); }
+        } else if (fin.eof()) {
+            break;
+        } else {
+            parse_unknown_field(fin, buf);
+        }
+    }
+    fin.close();
+}
+
+IGL_INLINE void igl::MshLoader::parse_nodes(std::ifstream& fin) {
+    size_t num_nodes;
+    fin >> num_nodes;
+    m_nodes.resize(num_nodes*3);
+
+    if (m_binary) {
+		size_t stride = (4+3*m_data_size);
+        size_t num_bytes = stride * num_nodes;
+        char* data = new char[num_bytes];
+        igl::_msh_eat_white_space(fin);
+        fin.read(data, num_bytes);
+
+        for (size_t i=0; i<num_nodes; i++) {
+            int node_idx;
+			memcpy(&node_idx, data+i*stride, sizeof(int));
+			node_idx-=1;
+			// directly move into vector storage
+			// this works only when m_data_size==sizeof(Float)==sizeof(double)
+			memcpy(&m_nodes[node_idx*3], data+i*stride + 4, m_data_size*3);
+        }
+        delete [] data;
+    } else {
+        int node_idx;
+        for (size_t i=0; i<num_nodes; i++) {
+            fin >> node_idx;
+            node_idx -= 1;
+            // here it's 3D node explicitly
+            fin >> m_nodes[node_idx*3]
+                >> m_nodes[node_idx*3+1]
+                >> m_nodes[node_idx*3+2];
+        }
+    }
+}
+
+IGL_INLINE void igl::MshLoader::parse_elements(std::ifstream& fin) {
+    m_elements_tags.resize(2); //hardcoded to have 2 tags
+    size_t num_elements;
+    fin >> num_elements;
+
+    size_t nodes_per_element;
+
+    if (m_binary) {
+        igl::_msh_eat_white_space(fin);
+        int elem_read = 0;
+        while (elem_read < num_elements) {
+            // Parse element header.
+            int elem_type, num_elems, num_tags;
+            fin.read((char*)&elem_type, sizeof(int));
+            fin.read((char*)&num_elems, sizeof(int));
+            fin.read((char*)&num_tags,  sizeof(int));
+            nodes_per_element = num_nodes_per_elem_type(elem_type);
+
+            // store node info
+            for (size_t i=0; i<num_elems; i++) {
+                int elem_idx;
+
+                // all elements in the segment share the same elem_type and number of nodes per element
+                m_elements_types.push_back(elem_type);
+                m_elements_lengths.push_back(nodes_per_element);
+
+                fin.read((char*)&elem_idx, sizeof(int));
+                elem_idx -= 1;
+                m_elements_ids.push_back(elem_idx);
+
+                // read first two tags
+                for (size_t j=0; j<num_tags; j++) {
+                    int tag;
+                    fin.read((char*)&tag, sizeof(int));
+                    if(j<2) m_elements_tags[j].push_back(tag);
+                }
+
+                for (size_t j=num_tags; j<2; j++) 
+                    m_elements_tags[j].push_back(-1); // fill up tags if less then 2
+
+                m_elements_nodes_idx.push_back(m_elements.size());
+                // Element values.
+                for (size_t j=0; j<nodes_per_element; j++) {
+                    int idx;
+                    fin.read((char*)&idx, sizeof(int));
+                    
+                    m_elements.push_back(idx-1);
+                }
+            }
+            elem_read += num_elems;
+        }
+    } else {
+        for (size_t i=0; i<num_elements; i++) {
+            // Parse per element header
+            int elem_num, elem_type, num_tags;
+            fin >> elem_num >> elem_type >> num_tags;
+
+            // read tags.
+            for (size_t j=0; j<num_tags; j++) {
+                int tag;
+                fin >> tag;
+                if(j<2) m_elements_tags[j].push_back(tag);
+            }
+            for (size_t j=num_tags; j<2; j++) 
+                m_elements_tags[j].push_back(-1); // fill up tags if less then 2
+            
+            nodes_per_element = num_nodes_per_elem_type(elem_type);
+            m_elements_types.push_back(elem_type);
+            m_elements_lengths.push_back(nodes_per_element);
+
+            elem_num -= 1;
+            m_elements_ids.push_back(elem_num);
+            m_elements_nodes_idx.push_back(m_elements.size());
+            // Parse node idx.
+            for (size_t j=0; j<nodes_per_element; j++) {
+                int idx;
+                fin >> idx;
+                m_elements.push_back(idx-1); // msh index starts from 1.
+            }
+        }
+    }
+    // debug
+    assert(m_elements_types.size()   == m_elements_ids.size());
+    assert(m_elements_tags[0].size() == m_elements_ids.size());
+    assert(m_elements_tags[1].size() == m_elements_ids.size());
+    assert(m_elements_lengths.size() == m_elements_ids.size());
+}
+
+IGL_INLINE void igl::MshLoader::parse_node_field( std::ifstream& fin ) {
+    size_t num_string_tags;
+    size_t num_real_tags;
+    size_t num_int_tags;
+
+    fin >> num_string_tags;
+    std::vector<std::string> str_tags(num_string_tags);
+
+    for (size_t i=0; i<num_string_tags; i++) {
+        igl::_msh_eat_white_space(fin);
+        if (fin.peek() == '\"') {
+            // Handle field name between quotes.
+            char buf[128];
+            fin.get(); // remove the quote at the beginning.
+            fin.getline(buf, 128, '\"');
+            str_tags[i] = std::string(buf);
+        } else {
+            fin >> str_tags[i];
+        }
+    }
+
+    fin >> num_real_tags;
+    std::vector<Float> real_tags(num_real_tags);
+    for (size_t i=0; i<num_real_tags; i++)
+        fin >> real_tags[i];
+
+    fin >> num_int_tags;
+    std::vector<int> int_tags(num_int_tags);
+    for (size_t i=0; i<num_int_tags; i++)
+        fin >> int_tags[i];
+
+    if (num_string_tags <= 0 || num_int_tags <= 2) {
+        throw std::runtime_error("Unexpected number of field tags");
+    }
+    std::string fieldname = str_tags[0];
+    int num_components    = int_tags[1];
+    int num_entries       = int_tags[2];
+
+    std::vector<Float> field( num_entries*num_components );
+
+    if (m_binary) {
+        size_t num_bytes = (num_components * m_data_size + 4) * num_entries;
+        char* data = new char[num_bytes];
+        igl::_msh_eat_white_space(fin);
+        fin.read(data, num_bytes);
+        for (size_t i=0; i<num_entries; i++) {
+			int node_idx;
+			memcpy(&node_idx,&data[i*(4+num_components*m_data_size)],4);
+			
+            if(node_idx<1) throw std::runtime_error("Negative or zero index");
+            node_idx -= 1;
+			
+            if(node_idx>=num_entries) throw std::runtime_error("Index too big");
+            size_t base_idx = i*(4+num_components*m_data_size) + 4;
+            // TODO: make this work when m_data_size != sizeof(double) ?
+			memcpy(&field[node_idx*num_components], &data[base_idx], num_components*m_data_size);
+        }
+        delete [] data;
+    } else {
+        int node_idx;
+        for (size_t i=0; i<num_entries; i++) {
+            fin >> node_idx;
+            node_idx -= 1;
+            for (size_t j=0; j<num_components; j++) {
+                fin >> field[node_idx*num_components+j];
+            }
+        }
+    }
+    
+    m_node_fields_names.push_back(fieldname);
+    m_node_fields.push_back(field);
+    m_node_fields_components.push_back(num_components);
+}
+
+IGL_INLINE void igl::MshLoader::parse_element_field(std::ifstream& fin) {
+    size_t num_string_tags;
+    size_t num_real_tags;
+    size_t num_int_tags;
+
+    fin >> num_string_tags;
+    std::vector<std::string> str_tags(num_string_tags);
+    for (size_t i=0; i<num_string_tags; i++) {
+        igl::_msh_eat_white_space(fin);
+        if (fin.peek() == '\"') {
+            // Handle field name between quoates.
+            char buf[128];
+            fin.get(); // remove the quote at the beginning.
+            fin.getline(buf, 128, '\"');
+            str_tags[i] = buf;
+        } else {
+            fin >> str_tags[i];
+        }
+    }
+
+    fin >> num_real_tags;
+    std::vector<Float> real_tags(num_real_tags);
+    for (size_t i=0; i<num_real_tags; i++)
+        fin >> real_tags[i];
+
+    fin >> num_int_tags;
+    std::vector<int> int_tags(num_int_tags);
+    for (size_t i=0; i<num_int_tags; i++)
+        fin >> int_tags[i];
+
+    if (num_string_tags <= 0 || num_int_tags <= 2) {
+        throw std::runtime_error("Invalid file format");
+    }
+    std::string fieldname = str_tags[0];
+    int num_components = int_tags[1];
+    int num_entries = int_tags[2];
+    std::vector<Float> field(num_entries*num_components);
+
+    if (m_binary) {
+        size_t num_bytes = (num_components * m_data_size + 4) * num_entries;
+        char* data = new char[num_bytes];
+        igl::_msh_eat_white_space(fin);
+        fin.read(data, num_bytes);
+        for (int i=0; i<num_entries; i++) {
+			int elem_idx;
+			// works with sizeof(int)==4
+			memcpy(&elem_idx, &data[i*(4+num_components*m_data_size)],4);
+            elem_idx -= 1;
+			
+			// directly copy data into vector storage space
+			memcpy(&field[elem_idx*num_components], &data[i*(4+num_components*m_data_size) + 4], m_data_size*num_components);
+        }
+        delete [] data;
+    } else {
+        int elem_idx;
+        for (size_t i=0; i<num_entries; i++) {
+            fin >> elem_idx;
+            elem_idx -= 1;
+            for (size_t j=0; j<num_components; j++) {
+                fin >> field[elem_idx*num_components+j];
+            }
+        }
+    }
+    m_element_fields_names.push_back(fieldname);
+    m_element_fields.push_back(field);
+    m_element_fields_components.push_back(num_components);
+}
+
+IGL_INLINE void igl::MshLoader::parse_unknown_field(std::ifstream& fin,
+        const std::string& fieldname) {
+    std::cerr << "Warning: \"" << fieldname << "\" not supported yet.  Ignored." << std::endl;
+    std::string endmark = fieldname.substr(0,1) + "End"
+        + fieldname.substr(1,fieldname.size()-1);
+
+    std::string buf("");
+    while (buf != endmark && !fin.eof()) {
+        fin >> buf;
+    }
+}
+
+IGL_INLINE int igl::MshLoader::num_nodes_per_elem_type(int elem_type) {
+    int nodes_per_element = 0;
+    switch (elem_type) {
+        case ELEMENT_LINE:         // 2-node line
+            nodes_per_element = 2; 
+            break;
+        case ELEMENT_TRI:
+            nodes_per_element = 3; // 3-node triangle
+            break;
+        case ELEMENT_QUAD:
+            nodes_per_element = 4; // 5-node quad
+            break;
+        case ELEMENT_TET:
+            nodes_per_element = 4; // 4-node tetrahedra
+            break;
+        case ELEMENT_HEX:          // 8-node hexahedron
+            nodes_per_element = 8; 
+            break;
+        case ELEMENT_PRISM:        // 6-node prism
+            nodes_per_element = 6; 
+            break;
+        case ELEMENT_LINE_2ND_ORDER: 
+            nodes_per_element = 3;
+            break;
+        case ELEMENT_TRI_2ND_ORDER: 
+            nodes_per_element = 6;
+            break;
+        case ELEMENT_QUAD_2ND_ORDER: 
+            nodes_per_element = 9;
+            break;
+        case ELEMENT_TET_2ND_ORDER: 
+            nodes_per_element = 10;
+            break;
+        case ELEMENT_HEX_2ND_ORDER: 
+            nodes_per_element = 27;
+            break;
+        case ELEMENT_PRISM_2ND_ORDER: 
+            nodes_per_element = 18;
+            break;
+        case ELEMENT_PYRAMID_2ND_ORDER: 
+            nodes_per_element = 14;
+            break;
+        case ELEMENT_POINT:        // 1-node point
+            nodes_per_element = 1; 
+            break;
+        default:
+            std::stringstream err_msg;
+                err_msg << "Element type (" << elem_type << ") is not supported yet."
+                << std::endl;
+            throw std::runtime_error(err_msg.str());
+    }
+    return nodes_per_element;
+}
+
+
+IGL_INLINE bool igl::MshLoader::is_element_map_identity() const
+{
+    for(int i=0;i<m_elements_ids.size();i++) {
+        int id=m_elements_ids[i];
+        if (id!=i) return false;
+    }
+    return true;
+}
+
+
+IGL_INLINE void igl::MshLoader::index_structures(int tag_column)
+{
+    //cleanup
+    m_structure_index.clear();
+    m_structures.clear();
+    m_structure_length.clear();
+
+    //index structure tags
+    for(auto i=0; i != m_elements_tags[tag_column].size(); ++i )
+    {
+        m_structure_index.insert(
+            std::pair<msh_struct,int>(
+                msh_struct( m_elements_tags[tag_column][i], 
+                            m_elements_types[i]), i)
+            );
+    }
+
+    // identify unique structures 
+    std::vector<StructIndex::value_type> _unique_structs;
+    std::unique_copy(std::begin(m_structure_index), 
+                     std::end(m_structure_index), 
+                     std::back_inserter(_unique_structs),
+                     [](const StructIndex::value_type &c1, const StructIndex::value_type &c2)
+                     { return c1.first == c2.first; });
+
+    std::for_each( _unique_structs.begin(), _unique_structs.end(), 
+        [this](const StructIndex::value_type &n){ this->m_structures.push_back(n.first); });
+
+    for(auto t = m_structures.begin(); t != m_structures.end(); ++t)
+    {
+        // identify all elements corresponding to this tag
+        auto structure_range = m_structure_index.equal_range( *t );
+        int cnt=0;
+
+        for(auto i=structure_range.first; i!=structure_range.second; i++)
+            cnt++;
+
+        m_structure_length.insert( std::pair<msh_struct,int>( *t, cnt));
+    }
+}

include/igl/MshLoader.h

@@ -0,0 +1,190 @@
+// based on MSH reader from PyMesh 
+
+// Copyright (c) 2015 Qingnan Zhou <[email protected]>           
+// Copyright (C) 2020 Vladimir Fonov <[email protected]> 
+//
+// This Source Code Form is subject to the terms of the Mozilla 
+// Public License v. 2.0. If a copy of the MPL was not distributed 
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 
+#ifndef IGL_MSH_LOADER_H
+#define IGL_MSH_LOADER_H
+#include "igl_inline.h"
+
+#include <fstream>
+#include <map>
+#include <string>
+#include <vector>
+#include <algorithm>
+
+namespace igl {
+
+// Class for loading information from .msh file
+// depends only on c++stl library
+class MshLoader {
+    public:
+
+        struct msh_struct {
+            int tag,el_type;
+            msh_struct(int _tag=0,int _type=0):
+                tag(_tag),el_type(_type){}
+            bool operator== (const msh_struct& a) const {
+                return this->tag==a.tag && 
+                       this->el_type==a.el_type;
+            }
+
+            bool operator< (const msh_struct& a) const {
+                return (this->tag*100+this->el_type) < 
+                       (a.tag*100+a.el_type);
+            }
+        };
+
+        typedef double Float;
+        
+        typedef std::vector<int>      IndexVector;
+        typedef std::vector<int>      IntVector;
+        typedef std::vector<Float>    FloatVector;
+        typedef std::vector<FloatVector> FloatField;
+        typedef std::vector<IntVector> IntField;
+        typedef std::vector<std::string> FieldNames;
+        typedef std::multimap<msh_struct,int> StructIndex;
+        typedef std::vector<msh_struct> StructVector;
+
+        enum {ELEMENT_LINE=1, ELEMENT_TRI=2, ELEMENT_QUAD=3, 
+              ELEMENT_TET=4,  ELEMENT_HEX=5, ELEMENT_PRISM=6,
+              ELEMENT_PYRAMID=7,
+              // 2nd order elements
+              ELEMENT_LINE_2ND_ORDER=8, ELEMENT_TRI_2ND_ORDER=9, 
+              ELEMENT_QUAD_2ND_ORDER=10,ELEMENT_TET_2ND_ORDER=11, 
+              ELEMENT_HEX_2ND_ORDER=12, ELEMENT_PRISM_2ND_ORDER=13, 
+              ELEMENT_PYRAMID_2ND_ORDER=14,
+              // other elements
+              ELEMENT_POINT=15 };
+    public:
+        MshLoader(const std::string &filename);
+
+    public:
+
+        // get nodes , x,y,z sequentially
+        const FloatVector& get_nodes()    const { return m_nodes; } 
+        // get elements , identifying nodes that create an element
+        // variable length per element
+        const IndexVector& get_elements() const { return m_elements; }
+
+        // get element types 
+        const IntVector& get_elements_types() const { return m_elements_types; }
+        // get element lengths
+        const IntVector& get_elements_lengths() const { return m_elements_lengths; }
+        // get element tags ( physical (0) and elementary (1) )
+        const IntField&  get_elements_tags() const { return m_elements_tags; }
+        // get element IDs
+        const IntVector& get_elements_ids() const { return m_elements_ids; }
+
+        // get reverse index from node to element
+        const IndexVector& get_elements_nodes_idx() const { return m_elements_nodes_idx; }
+
+        // get fields assigned per node, all fields and components sequentially
+        const FloatField& get_node_fields() const { return m_node_fields;}
+        // get node field names, 
+        const FieldNames& get_node_fields_names() const { return m_node_fields_names;}
+        // get number of node field components
+        const IntVector&  get_node_fields_components() const {return m_node_fields_components;}
+
+        int get_node_field_components(size_t c)  const 
+        {
+            return m_node_fields_components[c];
+        }
+
+        // get fields assigned per element, all fields and components sequentially
+        const FloatField& get_element_fields() const { return m_element_fields;}
+        // get element field names
+        const FieldNames& get_element_fields_names() const { return m_element_fields_names;}
+        // get number of element field components
+        const IntVector&  get_element_fields_components() const {return m_element_fields_components;}
+
+        int get_element_field_components(size_t c)  const {
+            return m_element_fields_components[c];
+        }
+        // check if field is present at node level
+        bool is_node_field(const std::string& fieldname)  const {
+            return (std::find(std::begin(m_node_fields_names),
+                              std::end(m_node_fields_names),
+                              fieldname) != std::end(m_node_fields_names) );
+        }
+        // check if field is present at element level
+        bool is_element_field(const std::string& fieldname) const {
+            return (std::find(std::begin(m_element_fields_names),
+                              std::end(m_element_fields_names),
+                              fieldname) != std::end(m_node_fields_names) );
+        }
+
+        // check if all elements have ids assigned sequentially
+        bool is_element_map_identity() const ;
+
+        // create tag index
+        // tag_column: ( physical (0) or elementary (1) ) specifying which tag to use
+        void index_structures(int tag_column); 
+
+        // get tag index, call index_structure_tags first
+        const StructIndex& get_structure_index() const 
+        {
+            return m_structure_index;
+        }
+
+        // get size of a structure identified by tag and element type
+        const StructIndex& get_structure_length() const 
+        {
+            return m_structure_length;
+        }
+
+        //! get list of structures
+        const StructVector& get_structures() const 
+        {
+            return m_structures;
+        }
+        
+    public:
+        // helper function, calculate number of nodes associated with an element
+        static int num_nodes_per_elem_type(int elem_type);
+
+    private:
+        void parse_nodes(std::ifstream& fin);
+        void parse_elements(std::ifstream& fin);
+        void parse_node_field(std::ifstream& fin);
+        void parse_element_field(std::ifstream& fin);
+        void parse_unknown_field(std::ifstream& fin,
+                const std::string& fieldname);
+
+    private:
+        bool   m_binary;
+        size_t m_data_size;
+        
+        FloatVector m_nodes;    // len x 3 vector 
+
+        IndexVector m_elements; // linear array for nodes corresponding to each element 
+        IndexVector m_elements_nodes_idx; // element indexes  
+
+        IntVector   m_elements_ids;     // element id's 
+        IntVector   m_elements_types;   // Element types 
+        IntVector   m_elements_lengths; // Element lengths 
+        IntField    m_elements_tags;    // Element tags, currently 2xtags per element 
+
+        FloatField  m_node_fields;      // Float field defined at each node 
+        IntVector   m_node_fields_components; // Number of components for node field 
+        FieldNames  m_node_fields_names; // Node field name 
+
+        FloatField  m_element_fields;    // Float field defined at each element 
+        IntVector   m_element_fields_components; // Number of components for element field 
+        FieldNames  m_element_fields_names; // Element field name 
+
+        StructIndex  m_structure_index; // index tag ids  
+        StructVector m_structures;  // unique structures
+        StructIndex  m_structure_length; // length of structures with consistent element type
+};
+
+} //igl
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "MshLoader.cpp"
+#endif
+
+#endif //IGL_MSH_LOADER_H

include/igl/MshSaver.cpp

@@ -0,0 +1,347 @@
+// based on MSH writer from PyMesh 
+
+// Copyright (c) 2015 Qingnan Zhou <[email protected]>           
+// Copyright (C) 2020 Vladimir Fonov <[email protected]> 
+//
+// This Source Code Form is subject to the terms of the Mozilla 
+// Public License v. 2.0. If a copy of the MPL was not distributed 
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 
+
+#include "MshSaver.h"
+
+#include <cassert>
+#include <iostream>
+#include <sstream>
+#include <exception>
+
+
+IGL_INLINE igl::MshSaver::MshSaver(const std::string& filename, bool binary) :
+    m_binary(binary), m_num_nodes(0), m_num_elements(0) {
+        if (!m_binary) {
+            fout.open(filename.c_str(), std::fstream::out);
+        } else {
+            fout.open(filename.c_str(), std::fstream::binary);
+        }
+        if (!fout) {
+            std::stringstream err_msg;
+            err_msg << "Error opening " << filename << " to write msh file." << std::endl;
+            throw std::ios_base::failure(err_msg.str());
+        }
+}
+
+IGL_INLINE igl::MshSaver::~MshSaver() {
+    fout.close();
+}
+
+IGL_INLINE void igl::MshSaver::save_mesh(
+    const FloatVector& nodes, 
+    const IndexVector& elements, 
+    const IntVector& element_lengths,
+    const IntVector& element_types,
+    const IntVector& element_tags
+     ) {
+
+    save_header();
+
+    save_nodes(nodes);
+
+    save_elements(elements, element_lengths, element_types, element_tags );
+}
+
+IGL_INLINE void igl::MshSaver::save_header() {
+    if (!m_binary) {
+        fout << "$MeshFormat" << std::endl;
+        fout << "2.2 0 " << sizeof(double) << std::endl;
+        fout << "$EndMeshFormat" << std::endl;
+        fout.precision(17);
+    } else {
+        fout << "$MeshFormat" << std::endl;
+        fout << "2.2 1 " << sizeof(double) << std::endl;
+        int one = 1;
+        fout.write((char*)&one, sizeof(int));
+        fout << "\n$EndMeshFormat" << std::endl;
+    }
+    fout.flush();
+}
+
+IGL_INLINE void igl::MshSaver::save_nodes(const FloatVector& nodes) {
+    // Save nodes.
+    // 3D hadrcoded
+    m_num_nodes = nodes.size() / 3;
+    fout << "$Nodes" << std::endl;
+    fout << m_num_nodes << std::endl;
+    if (!m_binary) {
+        for (size_t i=0; i<nodes.size(); i+=3) {
+            //const VectorF& v = nodes.segment(i,m_dim);
+            int node_idx = i/3 + 1;
+            fout << node_idx << " " << nodes[i] << " " << nodes[i+1] << " " << nodes[i+2] << std::endl;
+        }
+    } else {
+        for (size_t i=0; i<nodes.size(); i+=3) {
+            //const VectorF& v = nodes.segment(i,m_dim);
+            int node_idx = i/3 + 1;
+            fout.write((const char*)&node_idx, sizeof(int));
+            fout.write((const char*)&nodes[i], sizeof(Float)*3);
+        }
+    }
+    fout << "$EndNodes" << std::endl;
+    fout.flush();
+}
+
+IGL_INLINE void igl::MshSaver::save_elements(const IndexVector& elements, 
+            const IntVector& element_lengths,
+            const IntVector& element_types,
+            const IntVector& element_tags) 
+    {
+
+    m_num_elements = element_tags.size();
+    assert(element_lengths.size() == element_types.size() );
+    assert(element_lengths.size() == element_tags.size() );
+    // TODO: sum up all lengths
+    // Save elements.
+    // node inxes are 1-based
+    fout << "$Elements" << std::endl;
+    fout << m_num_elements << std::endl;
+
+    if (m_num_elements > 0) {
+        //int elem_type = el_type;
+        int num_elems = m_num_elements;
+        //int tags = 0;
+        if (!m_binary) {
+            size_t el_ptr=0;
+            for (size_t i=0;i<m_num_elements;++i) {
+                
+                int elem_num = (int) i + 1;
+                ///VectorI elem = elements.segment(i, nodes_per_element) + VectorI::Ones(nodes_per_element);
+                // hardcoded: duplicate tags (I don't know why)
+                fout << elem_num << " " << element_types[i] << " " << 2 << " "<< element_tags[i] << " "<< element_tags[i] << " ";
+                for (size_t j=0; j<element_lengths[i]; j++) {
+                    fout << elements[el_ptr + j] + 1 << " ";
+                }
+                fout << std::endl;
+                el_ptr+=element_lengths[i];
+            }
+        } else {
+            size_t el_ptr=0,i=0;
+            while(i<m_num_elements) {
+
+                // write elements in consistent chunks
+                // TODO: refactor this code to be able to specify different elements
+                // more effeciently 
+
+                int elem_type=-1;
+                int elem_len=-1;
+                size_t j=i;
+                for(;j<m_num_elements;++j)
+                {
+                    if( elem_type==-1 ) 
+                    {
+                        elem_type=element_types[j];
+                        elem_len=element_lengths[j];
+                    } else if( elem_type!=element_types[j] || 
+                               elem_len!=element_lengths[j]) {
+                        break; // found the edge of the segment
+                    }
+                }
+
+                //hardcoded: 2 tags
+                int num_elems=j-i, num_tags=2;
+
+                fout.write((const char*)& elem_type, sizeof(int));
+                fout.write((const char*)& num_elems, sizeof(int));
+                fout.write((const char*)& num_tags,  sizeof(int));
+
+                for(int k=0;k<num_elems; ++k,++i){
+                    int elem_num = (int )i + 1;
+                    fout.write((const char*)&elem_num, sizeof(int));
+
+                    // HACK: hardcoded 2 tags
+                    fout.write((const char*)& element_tags[i], sizeof(int));
+                    fout.write((const char*)& element_tags[i], sizeof(int));
+
+                    for (size_t e=0; e<elem_len; e++) {
+                        int _elem = static_cast<int>( elements[el_ptr + e] )+1;
+                        fout.write((const char*)&_elem, sizeof(int));
+                    }
+                    el_ptr+=elem_len;
+                }
+            }
+        }
+    }
+    fout << "$EndElements" << std::endl;
+    fout.flush();
+}
+
+IGL_INLINE void igl::MshSaver::save_scalar_field(const std::string& fieldname, const FloatVector& field) {
+    assert(field.size() == m_num_nodes);
+    fout << "$NodeData" << std::endl;
+    fout << "1" << std::endl; // num string tags.
+    fout << "\"" << fieldname << "\"" << std::endl;
+    fout << "1" << std::endl; // num real tags.
+    fout << "0.0" << std::endl; // time value.
+    fout << "3" << std::endl; // num int tags.
+    fout << "0" << std::endl; // the time step
+    fout << "1" << std::endl; // 1-component scalar field.
+    fout << m_num_nodes << std::endl; // number of nodes
+
+    if (m_binary) {
+        for (size_t i=0; i<m_num_nodes; i++) {
+            int node_idx = i+1;
+            fout.write((char*)&node_idx, sizeof(int));
+            fout.write((char*)&field[i], sizeof(Float));
+        }
+    } else {
+        for (size_t i=0; i<m_num_nodes; i++) {
+            int node_idx = i+1;
+            fout << node_idx << " " << field[i] << std::endl;
+        }
+    }
+    fout << "$EndNodeData" << std::endl;
+    fout.flush();
+}
+
+IGL_INLINE void igl::MshSaver::save_vector_field(const std::string& fieldname, const FloatVector& field) {
+    assert(field.size() == 3 * m_num_nodes);
+
+    fout << "$NodeData" << std::endl;
+    fout << "1" << std::endl; // num string tags.
+    fout << "\"" << fieldname << "\"" << std::endl;
+    fout << "1" << std::endl; // num real tags.
+    fout << "0.0" << std::endl; // time value.
+    fout << "3" << std::endl; // num int tags.
+    fout << "0" << std::endl; // the time step
+    fout << "3" << std::endl; // 3-component vector field.
+    fout << m_num_nodes << std::endl; // number of nodes
+
+    const Float zero = 0.0;
+    if (m_binary) {
+        for (size_t i=0; i<m_num_nodes; i++) {
+            int node_idx = i+1;
+            fout.write((const char*)&node_idx, sizeof(int));
+            fout.write((const char*)&field[i*3], sizeof(Float)*3);
+        }
+    } else {
+        for (size_t i=0; i<m_num_nodes; i++) {
+            int node_idx = i+1;
+                fout << node_idx
+                    << " " << field[i*3]
+                    << " " << field[i*3+1]
+                    << " " << field[i*3+2]
+                    << std::endl;
+        }
+    }
+    fout << "$EndNodeData" << std::endl;
+    fout.flush();
+}
+
+IGL_INLINE void igl::MshSaver::save_elem_scalar_field(const std::string& fieldname, const FloatVector& field) {
+    assert(field.size() == m_num_elements);
+    fout << "$ElementData" << std::endl;
+    fout << 1 << std::endl; // num string tags.
+    fout << "\"" << fieldname << "\"" << std::endl;
+    fout << "1" << std::endl; // num real tags.
+    fout << "0.0" << std::endl; // time value.
+    fout << "3" << std::endl; // num int tags.
+    fout << "0" << std::endl; // the time step
+    fout << "1" << std::endl; // 1-component scalar field.
+    fout << m_num_elements << std::endl; // number of elements
+
+    if (m_binary) {
+        for (size_t i=0; i<m_num_elements; i++) {
+            int elem_idx = i+1;
+            fout.write((const char*)&elem_idx, sizeof(int));
+            fout.write((const char*)&field[i], sizeof(Float));
+        }
+    } else {
+        for (size_t i=0; i<m_num_elements; i++) {
+            int elem_idx = i+1;
+            fout << elem_idx << " " << field[i] << std::endl;
+        }
+    }
+
+    fout << "$EndElementData" << std::endl;
+    fout.flush();
+}
+
+IGL_INLINE void igl::MshSaver::save_elem_vector_field(const std::string& fieldname, const FloatVector& field) {
+    assert(field.size() == m_num_elements * 3);
+    fout << "$ElementData" << std::endl;
+    fout << 1 << std::endl; // num string tags.
+    fout << "\"" << fieldname << "\"" << std::endl;
+    fout << "1" << std::endl; // num real tags.
+    fout << "0.0" << std::endl; // time value.
+    fout << "3" << std::endl; // num int tags.
+    fout << "0" << std::endl; // the time step
+    fout << "3" << std::endl; // 3-component vector field.
+    fout << m_num_elements << std::endl; // number of elements
+
+    const Float zero = 0.0;
+    if (m_binary) {
+        for (size_t i=0; i<m_num_elements; ++i) {
+            int elem_idx = i+1;
+            fout.write((const char*)&elem_idx, sizeof(int));
+            fout.write((const char*)&field[i*3], sizeof(Float) * 3);
+        }
+    } else {
+        for (size_t i=0; i<m_num_elements; ++i) {
+            int elem_idx = i+1;
+            fout << elem_idx
+                << " " << field[i*3]
+                << " " << field[i*3+1]
+                << " " << field[i*3+2]
+                << std::endl;
+        }
+    }
+
+    fout << "$EndElementData" << std::endl;
+    fout.flush();
+}
+
+
+IGL_INLINE void igl::MshSaver::save_elem_tensor_field(const std::string& fieldname, const FloatVector& field) {
+    assert(field.size() == m_num_elements * 3 * (3 + 1) / 2);
+    fout << "$ElementData" << std::endl;
+    fout << 1 << std::endl; // num string tags.
+    fout << "\"" << fieldname << "\"" << std::endl;
+    fout << "1" << std::endl; // num real tags.
+    fout << "0.0" << std::endl; // time value.
+    fout << "3" << std::endl; // num int tags.
+    fout << "0" << std::endl; // the time step
+    fout << "9" << std::endl; // 9-component tensor field.
+    fout << m_num_elements << std::endl; // number of elements
+
+    const Float zero = 0.0;
+    
+    if (m_binary) {
+        for (size_t i=0; i<m_num_elements; i++) {
+            int elem_idx = i+1;
+            fout.write((char*)&elem_idx, sizeof(int));
+            //const VectorF& val = field.segment(i*6, 6);
+            const Float* val = &field[i*6];
+            Float tensor[9] = {
+                val[0], val[5], val[4],
+                val[5], val[1], val[3],
+                val[4], val[3], val[2] };
+            fout.write((char*)tensor, sizeof(Float) * 9);
+        }
+    } else {
+        for (size_t i=0; i<m_num_elements; i++) {
+            int elem_idx = i+1;
+            const Float* val = &field[i*6];
+            fout << elem_idx
+                << " " << val[0]
+                << " " << val[5]
+                << " " << val[4]
+                << " " << val[5]
+                << " " << val[1]
+                << " " << val[3]
+                << " " << val[4]
+                << " " << val[3]
+                << " " << val[2]
+                << std::endl;
+        }
+    }
+
+    fout << "$EndElementData" << std::endl;
+    fout.flush();
+}

include/igl/MshSaver.h

@@ -0,0 +1,84 @@
+// based on MSH writer from PyMesh 
+
+// Copyright (c) 2015 Qingnan Zhou <[email protected]>           
+// Copyright (C) 2020 Vladimir Fonov <[email protected]> 
+//
+// This Source Code Form is subject to the terms of the Mozilla 
+// Public License v. 2.0. If a copy of the MPL was not distributed 
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 
+#ifndef IGL_MSH_SAVER_H
+#define IGL_MSH_SAVER_H
+#include "igl_inline.h"
+
+#include <fstream>
+#include <string>
+#include <vector>
+
+namespace igl {
+
+// Class for dumping information to .msh file
+// depends only on c++stl library
+// current implementation works only with 3D information
+class MshSaver {
+    public:
+        typedef double Float;
+
+        typedef std::vector<int>         IndexVector;
+        typedef std::vector<int>         IntVector;
+        typedef std::vector<Float>       FloatVector;
+        typedef std::vector<FloatVector> FloatField;
+        typedef std::vector<IntVector>   IntField;
+        typedef std::vector<std::string> FieldNames;
+
+        MshSaver(const std::string& filename, bool binary=true);
+        ~MshSaver();
+
+    public:
+        // Only these element types are supported right now
+        enum {ELEMENT_LINE=1, ELEMENT_TRI=2, ELEMENT_QUAD=3, 
+              ELEMENT_TET=4,  ELEMENT_HEX=5, ELEMENT_PRISM=6 };
+
+    public:
+        // save mesh geometry
+        void save_mesh(
+            const FloatVector& nodes, 
+            const IndexVector& elements, 
+            const IntVector& element_lengths,
+            const IntVector& element_type,
+            const IntVector& element_tags );
+        
+        // save additional fields associated with the mesh
+
+        // add node scalar field
+        void save_scalar_field(const std::string& fieldname, const FloatVector& field);
+        // add node vectot field
+        void save_vector_field(const std::string& fieldname, const FloatVector& field);
+        // add element scalar field
+        void save_elem_scalar_field(const std::string& fieldname, const FloatVector& field);
+        // add element vector field
+        void save_elem_vector_field(const std::string& fieldname, const FloatVector& field);
+        // add element tensor field
+        void save_elem_tensor_field(const std::string& fieldname, const FloatVector& field);
+
+    protected:
+        void save_header();
+        void save_nodes(const FloatVector& nodes);
+        void save_elements(const IndexVector& elements, 
+            const IntVector& element_lengths,
+            const IntVector& element_type,
+            const IntVector& element_tags);
+
+    private:
+        bool m_binary;
+        size_t m_num_nodes;
+        size_t m_num_elements;
+
+        std::ofstream fout;
+};
+} //igl
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "MshSaver.cpp"
+#endif
+
+#endif //MSH_SAVER_H

include/igl/readMSH.cpp

@@ -1,514 +1,211 @@
+/* high level interface for MshLoader */
+
+/* Copyright (C) 2020 Vladimir Fonov <[email protected]> */
+/*
+/* This Source Code Form is subject to the terms of the Mozilla */
+/* Public License v. 2.0. If a copy of the MPL was not distributed */
+/* with this file, You can obtain one at http://mozilla.org/MPL/2.0/. */
 
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2018 Alec Jacobson <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
-// obtain one at http://mozilla.org/MPL/2.0/.
-//
 #include "readMSH.h"
+#include "MshLoader.h"
 #include <iostream>
-#include <sstream>
-#include <fstream>
-#include <vector>
-#include <map>
 
-template <
-  typename DerivedV,
-  typename DerivedT>
-IGL_INLINE bool igl::readMSH(
-  const std::string & filename,
-  Eigen::PlainObjectBase<DerivedV> & V,
-  Eigen::PlainObjectBase<DerivedT> & T)
+IGL_INLINE  bool  igl::readMSH(const std::string &msh,
+            Eigen::MatrixXd &X,
+            Eigen::MatrixXi &Tri,
+            Eigen::MatrixXi &Tet,
+            Eigen::VectorXi &TriTag,
+            Eigen::VectorXi &TetTag,
+            std::vector<std::string>     &XFields,
+            std::vector<Eigen::MatrixXd> &XF,
+            std::vector<std::string>     &EFields,
+            std::vector<Eigen::MatrixXd> &TriF,
+            std::vector<Eigen::MatrixXd> &TetF
+            )
 {
-  // https://github.com/Yixin-Hu/TetWild/blob/master/pymesh/MshSaver.cpp
-  // Original copyright: /* This file is part of PyMesh. Copyright (c) 2015 by Qingnan Zhou */
-  typedef typename DerivedV::Scalar Float;
-  typedef Eigen::Matrix<Float,Eigen::Dynamic,1> VectorF;
-  typedef Eigen::Matrix<int,Eigen::Dynamic,1> VectorI;
-  typedef std::map<std::string, VectorF> FieldMap;
-  typedef std::vector<std::string> FieldNames;
-  VectorF m_nodes;
-  VectorI m_elements;
-  FieldMap m_node_fields;
-  FieldMap m_element_fields;
-
-  bool m_binary;
-  size_t m_data_size;
-  size_t m_nodes_per_element;
-  size_t m_element_type;
-  std::ifstream fin(filename.c_str(), std::ios::in | std::ios::binary);
-  if (!fin.is_open())
-  {
-    std::stringstream err_msg;
-    err_msg << "failed to open file \"" << filename << "\"";
-    return false;
-  }
-  // Parse header
-  std::string buf;
-  double version;
-  int type;
-  fin >> buf;
-  const auto invalid_format = []()->bool
-  {
-    assert(false && "Invalid format");
-    return false;
-  };
-  const auto not_implemented = []()->bool
-  {
-    assert(false && "Not implemented");
-    return false;
-  };
-  if (buf != "$MeshFormat") { return invalid_format(); }
-
-  fin >> version >> type >> m_data_size;
-  m_binary = (type == 1);
-
-  // Some sanity check.
-  if (m_data_size != 8) {
-      std::cerr << "Error: data size must be 8 bytes." << std::endl;
-      return not_implemented();
-  }
-  if (sizeof(int) != 4) {
-      std::cerr << "Error: code must be compiled with int size 4 bytes." << std::endl;
-      return not_implemented();
-  }
-  const auto eat_white_space = [](std::ifstream& fin)
-  {
-    char next = fin.peek();
-    while (next == '\n' || next == ' ' || next == '\t' || next == '\r') 
+    try 
     {
-      fin.get();
-      next = fin.peek();
-    }
-  };
-
-  // Read in extra info from binary header.
-  if (m_binary) {
-      int one;
-      eat_white_space(fin);
-      fin.read(reinterpret_cast<char*>(&one), sizeof(int));
-      if (one != 1) {
-          std::cerr << "Warning: binary msh file " << filename
-              << " is saved with different endianness than this machine."
-              << std::endl;
-          return not_implemented();
-      }
-  }
-
-  fin >> buf;
-  if (buf != "$EndMeshFormat") { return not_implemented(); }
-
-  const auto num_nodes_per_elem_type = [](int elem_type)->int 
-  {
-    size_t nodes_per_element = 0;
-    switch (elem_type) {
-        case 2:
-            nodes_per_element = 3; // Triangle
-            break;
-        case 3:
-            nodes_per_element = 4; // Quad
-            break;
-        case 4:
-            nodes_per_element = 4; // Tet
-            break;
-        case 5:
-            nodes_per_element = 8; // hexahedron
-            break;
-        default:
-            assert(false && "not implemented");
-            nodes_per_element = -1;
-            break;
-    }
-    return nodes_per_element;
-  };
+        igl::MshLoader _loader(msh);
+        const int USETAG = 1;
 
-  const auto parse_nodes = [&](std::ifstream& fin) 
-  {
-    size_t num_nodes;
-    fin >> num_nodes;
-    m_nodes.resize(num_nodes*3);
+        #ifndef NDEBUG
+        std::cout<<"readMSH:Total number of nodes:" << _loader.get_nodes().size()<<std::endl;   
+        std::cout<<"readMSH:Total number of elements:" << _loader.get_elements().size()<<std::endl;
 
-    if (m_binary) {
-        size_t num_bytes = (4+3*m_data_size) * num_nodes;
-        char* data = new char[num_bytes];
-        eat_white_space(fin);
-        fin.read(data, num_bytes);
-
-        for (size_t i=0; i<num_nodes; i++) {
-            int node_idx          = *reinterpret_cast<int*>  (&data[i*(4+3*m_data_size)]) - 1;
-            m_nodes[node_idx*3]   = *reinterpret_cast<Float*>(&data[i*(4+3*m_data_size) + 4]);
-            m_nodes[node_idx*3+1] = *reinterpret_cast<Float*>(&data[i*(4+3*m_data_size) + 4 + m_data_size]);
-            m_nodes[node_idx*3+2] = *reinterpret_cast<Float*>(&data[i*(4+3*m_data_size) + 4 + 2*m_data_size]);
+        std::cout<<"readMSH:Node fields:" << std::endl;
+        for(auto i=std::begin(_loader.get_node_fields_names()); i!=std::end(_loader.get_node_fields_names()); i++)
+        {
+            std::cout << i->c_str() << ":" << _loader.get_node_fields()[i-std::begin(_loader.get_node_fields_names())].size() << std::endl;
         }
-
-        delete [] data;
-    } else {
-        int node_idx;
-        for (size_t i=0; i<num_nodes; i++) {
-            fin >> node_idx;
-            node_idx -= 1;
-            fin >> m_nodes[node_idx*3]
-                >> m_nodes[node_idx*3+1]
-                >> m_nodes[node_idx*3+2];
+        
+        std::cout << "readMSH:Element fields:" << std::endl;
+        for(auto i=std::begin(_loader.get_element_fields_names()); i!=std::end(_loader.get_element_fields_names()); i++)
+        {
+            std::cout << i->c_str() << ":" << _loader.get_element_fields()[i-std::begin(_loader.get_element_fields_names())].size() << std::endl;
         }
-    }
-  };
-
-  const auto parse_elements = [&](std::ifstream& fin) 
-  {
-    size_t num_elements;
-    fin >> num_elements;
-
-    // Tmp storage of elements;
-    std::vector<int> triangle_element_idx;
-    std::vector<int> triangle_elements;
-    std::vector<int> quad_element_idx;
-    std::vector<int> quad_elements;
-    std::vector<int> tet_element_idx;
-    std::vector<int> tet_elements;
-    std::vector<int> hex_element_idx;
-    std::vector<int> hex_elements;
-
-    auto get_element_storage = [&](int elem_type) -> std::vector<int>* {
-        switch (elem_type) {
-            default:
-                assert(false && "Unsupported element type encountered");
-            case 2:
-                return &triangle_elements;
-            case 3:
-                return &quad_elements;
-            case 4:
-                return &tet_elements;
-            case 5:
-                return &hex_elements;
-        };
-    };
-
-    auto get_element_idx_storage = [&](int elem_type) -> std::vector<int>* {
-        switch (elem_type) {
-            default:
-                assert(false && "Unsupported element type encountered");
-            case 2:
-                return &triangle_element_idx;
-            case 3:
-                return &quad_element_idx;
-            case 4:
-                return &tet_element_idx;
-            case 5:
-                return &hex_element_idx;
-        };
-    };
-
-    size_t nodes_per_element;
-    int glob_elem_type = -1;
-
 
-  if (m_binary) 
-  {
-    eat_white_space(fin);
-    int elem_read = 0;
-    while (elem_read < num_elements) {
-        // Parse element header.
-        int elem_type, num_elems, num_tags;
-        fin.read((char*)&elem_type, sizeof(int));
-        fin.read((char*)&num_elems, sizeof(int));
-        fin.read((char*)&num_tags, sizeof(int));
-        nodes_per_element = num_nodes_per_elem_type(elem_type);
-        std::vector<int>& elements = *get_element_storage(elem_type);
-        std::vector<int>& element_idx = *get_element_idx_storage(elem_type);
-
-        for (size_t i=0; i<num_elems; i++) {
-            int elem_idx;
-            fin.read((char*)&elem_idx, sizeof(int));
-            elem_idx -= 1;
-            element_idx.push_back(elem_idx);
-
-            // Eat up tags.
-            for (size_t j=0; j<num_tags; j++) {
-                int tag;
-                fin.read((char*)&tag, sizeof(int));
-            }
-
-            // Element values.
-            for (size_t j=0; j<nodes_per_element; j++) {
-                int idx;
-                fin.read((char*)&idx, sizeof(int));
-                elements.push_back(idx-1);
-            }
+        if(_loader.is_element_map_identity())
+            std::cout<<"readMSH:Element ids map is identity"<<std::endl;
+        else
+            std::cout<<"readMSH:Element ids map is NOT identity"<<std::endl;
+        
+        #endif
+        
+        // convert nodes
+        // hadrcoded for 3D 
+        Eigen::Map< const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > 
+            node_map( _loader.get_nodes().data(), _loader.get_nodes().size()/3, 3 );
+
+        X = node_map;
+        XFields = _loader.get_element_fields_names();
+        XF.resize(_loader.get_node_fields().size());
+        XFields = _loader.get_node_fields_names();
+        for(size_t i=0;i<_loader.get_node_fields().size();++i)
+        {
+            Eigen::Map< const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > 
+                field_map( _loader.get_node_fields()[i].data(), 
+                        _loader.get_node_fields()[i].size()/_loader.get_node_fields_components()[i], 
+                        _loader.get_node_fields_components()[i] );
+            XF[i] = field_map;
         }
 
-        elem_read += num_elems;
-    }
-  } else 
-  {
-        for (size_t i=0; i<num_elements; i++) {
-            // Parse per element header
-            int elem_num, elem_type, num_tags;
-            fin >> elem_num >> elem_type >> num_tags;
-            for (size_t j=0; j<num_tags; j++) {
-                int tag;
-                fin >> tag;
-            }
-            nodes_per_element = num_nodes_per_elem_type(elem_type);
-            std::vector<int>& elements = *get_element_storage(elem_type);
-            std::vector<int>& element_idx = *get_element_idx_storage(elem_type);
-
-            elem_num -= 1;
-            element_idx.push_back(elem_num);
-
-            // Parse node idx.
-            for (size_t j=0; j<nodes_per_element; j++) {
-                int idx;
-                fin >> idx;
-                elements.push_back(idx-1); // msh index starts from 1.
-            }
-        }
-    }
-
-    auto copy_to_array = [&](
-            const std::vector<int>& elements,
-            const int nodes_per_element) {
-        const size_t num_elements = elements.size() / nodes_per_element;
-        if (elements.size() % nodes_per_element != 0) {
-            assert(false && "parsing element failed");
-            return;
+        // calculate number of elements 
+        std::map<int,int> element_counts;
+       
+        for(auto i:_loader.get_elements_types())
+        {
+            auto j=element_counts.insert({i,1});
+            if(!j.second) (*j.first).second+=1;
         }
-        m_elements.resize(elements.size());
-        std::copy(elements.begin(), elements.end(), m_elements.data());
-        m_nodes_per_element = nodes_per_element;
-    };
-
-    if (!tet_elements.empty()) {
-        copy_to_array(tet_elements, 4);
-        m_element_type = 4;
-    } else if (!hex_elements.empty()) {
-        copy_to_array(hex_elements, 8);
-        m_element_type = 5;
-    } else if (!triangle_elements.empty()) {
-        copy_to_array(triangle_elements, 3);
-        m_element_type = 2;
-    } else if (!quad_elements.empty()) {
-        copy_to_array(quad_elements, 4);
-        m_element_type = 3;
-    } else {
-        // 0 elements, use triangle by default.
-        m_element_type = 2;
-    }
-  };
-  const auto parse_element_field = [&](std::ifstream& fin) 
-  {
-    size_t num_string_tags;
-    size_t num_real_tags;
-    size_t num_int_tags;
-
-    fin >> num_string_tags;
-    std::string* str_tags = new std::string[num_string_tags];
-    for (size_t i=0; i<num_string_tags; i++) {
-        eat_white_space(fin);
-        if (fin.peek() == '\"') {
-            // Handle field name between quoates.
-            char buf[128];
-            fin.get(); // remove the quote at the beginning.
-            fin.getline(buf, 128, '\"');
-            str_tags[i] = std::string(buf);
-        } else {
-            fin >> str_tags[i];
+        #ifndef NDEBUG
+        std::cout<<"ReadMSH: elements found"<<std::endl;
+        for(auto i:element_counts)
+            std::cout<<"\t"<<i.first<<":"<<i.second<<std::endl;
+        #endif
+        int n_tri_el=0;
+        int n_tet_el=0;
+        
+        auto n_tri_el_=element_counts.find(igl::MshLoader::ELEMENT_TRI);
+        auto n_tet_el_=element_counts.find(igl::MshLoader::ELEMENT_TET);
+        if(n_tri_el_!=std::end(element_counts)) 
+            n_tri_el=n_tri_el_->second;
+        if(n_tet_el_!=std::end(element_counts)) 
+            n_tet_el=n_tet_el_->second;
+
+        Tri.resize(n_tri_el,3);
+        Tet.resize(n_tet_el,4);
+        TriTag.resize(n_tri_el);
+        TetTag.resize(n_tet_el);
+        size_t el_start = 0;
+        TriF.resize(_loader.get_element_fields().size());
+        TetF.resize(_loader.get_element_fields().size());
+        for(size_t i=0;i<_loader.get_element_fields().size();++i)
+        {
+            TriF[i].resize(n_tri_el,_loader.get_element_fields_components()[i]);
+            TetF[i].resize(n_tet_el,_loader.get_element_fields_components()[i]);
         }
-    }
-
-    fin >> num_real_tags;
-    Float* real_tags = new Float[num_real_tags];
-    for (size_t i=0; i<num_real_tags; i++)
-        fin >> real_tags[i];
-
-    fin >> num_int_tags;
-    int* int_tags = new int[num_int_tags];
-    for (size_t i=0; i<num_int_tags; i++)
-        fin >> int_tags[i];
-
-    if (num_string_tags <= 0 || num_int_tags <= 2) {
-        delete[] str_tags;
-        delete[] real_tags;
-        delete[] int_tags;
-        assert(false && "invalid format");
-        return;
-    }
-    std::string fieldname = str_tags[0];
-    int num_components = int_tags[1];
-    int num_entries = int_tags[2];
-    VectorF field(num_entries * num_components);
-
-    delete [] str_tags;
-    delete [] real_tags;
-    delete [] int_tags;
-
-    if (m_binary) {
-        size_t num_bytes = (num_components * m_data_size + 4) * num_entries;
-        char* data = new char[num_bytes];
-        eat_white_space(fin);
-        fin.read(data, num_bytes);
-        for (size_t i=0; i<num_entries; i++) {
-            int elem_idx = *reinterpret_cast<int*>(&data[i*(4+num_components*m_data_size)]);
-            elem_idx -= 1;
-            size_t base_idx = i*(4+num_components*m_data_size) + 4;
-            for (size_t j=0; j<num_components; j++) {
-                field[elem_idx * num_components + j] = *reinterpret_cast<Float*>(&data[base_idx+j*m_data_size]);
-            }
-        }
-        delete [] data;
-    } else {
-        int elem_idx;
-        for (size_t i=0; i<num_entries; i++) {
-            fin >> elem_idx;
-            elem_idx -= 1;
-            for (size_t j=0; j<num_components; j++) {
-                fin >> field[elem_idx * num_components + j];
+        EFields = _loader.get_element_fields_names();
+        int i_tri = 0;
+        int i_tet = 0;
+
+        for(size_t i=0;i<_loader.get_elements_lengths().size();++i)
+        {
+            if(_loader.get_elements_types()[i]==MshLoader::ELEMENT_TRI )
+            {
+                assert(_loader.get_elements_lengths()[i]==3);
+
+                Tri(i_tri, 0) = _loader.get_elements()[el_start  ];
+                Tri(i_tri, 1) = _loader.get_elements()[el_start+1];
+                Tri(i_tri, 2) = _loader.get_elements()[el_start+2];
+
+                TriTag(i_tri) = _loader.get_elements_tags()[1][i];
+
+                for(size_t j=0;j<_loader.get_element_fields().size();++j)
+                    for(size_t k=0;k<_loader.get_element_fields_components()[j];++k)
+                        TriF[j](i_tri,k) = _loader.get_element_fields()[j][_loader.get_element_fields_components()[j]*i+k];
+
+                ++i_tri;
+            } else if(_loader.get_elements_types()[i]==MshLoader::ELEMENT_TET ) {
+                assert(_loader.get_elements_lengths()[i]==4);
+
+                Tet(i_tet, 0) = _loader.get_elements()[el_start  ];
+                Tet(i_tet, 1) = _loader.get_elements()[el_start+1];
+                Tet(i_tet, 2) = _loader.get_elements()[el_start+2];
+                Tet(i_tet, 3) = _loader.get_elements()[el_start+3];
+
+                TetTag(i_tet) = _loader.get_elements_tags()[USETAG][i];
+
+                for(size_t j=0;j<_loader.get_element_fields().size();++j)
+                    for(size_t k=0;k<_loader.get_element_fields_components()[j];++k)
+                        TetF[j](i_tet,k) = _loader.get_element_fields()[j][_loader.get_element_fields_components()[j]*i+k];
+                
+                ++i_tet;
+            } else {
+                // else: it's unsupported type of the element, ignore for now
+                std::cerr<<"readMSH: unsupported element type: "<<_loader.get_elements_types()[i] << 
+                           ", length: "<< _loader.get_elements_lengths()[i] <<std::endl;
             }
-        }
-    }
-
-    m_element_fields[fieldname] = field;
-  };
 
-  const auto parse_node_field = [&](std::ifstream& fin) 
-  {
-    size_t num_string_tags;
-    size_t num_real_tags;
-    size_t num_int_tags;
-
-    fin >> num_string_tags;
-    std::string* str_tags = new std::string[num_string_tags];
-    for (size_t i=0; i<num_string_tags; i++) {
-        eat_white_space(fin);
-        if (fin.peek() == '\"') {
-            // Handle field name between quoates.
-            char buf[128];
-            fin.get(); // remove the quote at the beginning.
-            fin.getline(buf, 128, '\"');
-            str_tags[i] = std::string(buf);
-        } else {
-            fin >> str_tags[i];
+            el_start += _loader.get_elements_lengths()[i];
         }
-    }
-
-    fin >> num_real_tags;
-    Float* real_tags = new Float[num_real_tags];
-    for (size_t i=0; i<num_real_tags; i++)
-        fin >> real_tags[i];
-
-    fin >> num_int_tags;
-    int* int_tags = new int[num_int_tags];
-    for (size_t i=0; i<num_int_tags; i++)
-        fin >> int_tags[i];
-
-    if (num_string_tags <= 0 || num_int_tags <= 2) {
-        delete[] str_tags;
-        delete[] real_tags;
-        delete[] int_tags;
-        assert(false && "invalid format");
-        return;
-    }
-    std::string fieldname = str_tags[0];
-    int num_components = int_tags[1];
-    int num_entries = int_tags[2];
-    VectorF field(num_entries * num_components);
 
-    delete [] str_tags;
-    delete [] real_tags;
-    delete [] int_tags;
-
-    if (m_binary) {
-        size_t num_bytes = (num_components * m_data_size + 4) * num_entries;
-        char* data = new char[num_bytes];
-        eat_white_space(fin);
-        fin.read(data, num_bytes);
-        for (size_t i=0; i<num_entries; i++) {
-            int node_idx = *reinterpret_cast<int*>(&data[i*(4+num_components*m_data_size)]);
-            node_idx -= 1;
-            size_t base_idx = i*(4+num_components*m_data_size) + 4;
-            for (size_t j=0; j<num_components; j++) {
-                field[node_idx * num_components + j] = *reinterpret_cast<Float*>(&data[base_idx+j*m_data_size]);
-            }
-        }
-        delete [] data;
-    } else {
-        int node_idx;
-        for (size_t i=0; i<num_entries; i++) {
-            fin >> node_idx;
-            node_idx -= 1;
-            for (size_t j=0; j<num_components; j++) {
-                fin >> field[node_idx * num_components + j];
-            }
-        }
+        assert(i_tet == n_tet_el);
+        assert(i_tri == n_tri_el);
+    } catch(const std::exception& e) {
+        std::cerr << e.what() << std::endl;
+        return false;
     }
+    return true;
+}
 
-    m_node_fields[fieldname] = field;
-  };
-  const auto parse_unknown_field = [](std::ifstream& fin,
-        const std::string& fieldname) 
-  {
-    std::cerr << "Warning: \"" << fieldname << "\" not supported yet.  Ignored." << std::endl;
-    std::string endmark = fieldname.substr(0,1) + "End"
-        + fieldname.substr(1,fieldname.size()-1);
+IGL_INLINE bool igl::readMSH(const std::string &msh,
+                Eigen::MatrixXd &X,
+                Eigen::MatrixXi &Tri,
+                Eigen::MatrixXi &Tet,
+                Eigen::VectorXi &TriTag,
+                Eigen::VectorXi &TetTag
+                )
+{
+    std::vector<std::string>     XFields;
+    std::vector<Eigen::MatrixXd> XF;
+    std::vector<std::string>     EFields;
+    std::vector<Eigen::MatrixXd> TriF;
+    std::vector<Eigen::MatrixXd> TetF;
+    return igl::readMSH(msh,X,Tri,Tet,TriTag,TetTag,XFields,XF,EFields,TriF,TetF);
+}
 
-    std::string buf("");
-    while (buf != endmark && !fin.eof()) {
-        fin >> buf;
-    }
-  };
+IGL_INLINE bool igl::readMSH(const std::string &msh,
+                Eigen::MatrixXd &X,
+                Eigen::MatrixXi &Tri,
+                Eigen::VectorXi &TriTag
+                )
+{
+    Eigen::MatrixXi Tet;
+    Eigen::VectorXi TetTag;
 
+    std::vector<std::string>     XFields;
+    std::vector<Eigen::MatrixXd> XF;
+    std::vector<std::string>     EFields;
+    std::vector<Eigen::MatrixXd> TriF;
+    std::vector<Eigen::MatrixXd> TetF;
 
-  while (!fin.eof()) {
-      buf.clear();
-      fin >> buf;
-      if (buf == "$Nodes") {
-          parse_nodes(fin);
-          fin >> buf;
-          if (buf != "$EndNodes") { return invalid_format(); }
-      } else if (buf == "$Elements") {
-          parse_elements(fin);
-          fin >> buf;
-          if (buf != "$EndElements") { return invalid_format(); }
-      } else if (buf == "$NodeData") {
-          parse_node_field(fin);
-          fin >> buf;
-          if (buf != "$EndNodeData") { return invalid_format(); }
-      } else if (buf == "$ElementData") {
-          parse_element_field(fin);
-          fin >> buf;
-          if (buf != "$EndElementData") { return invalid_format(); }
-      } else if (fin.eof()) {
-          break;
-      } else {
-          parse_unknown_field(fin, buf);
-      }
-  }
-  fin.close();
-  V.resize(m_nodes.rows()/3,3);
-  for (int i = 0; i < m_nodes.rows() / 3; i++) 
-  {
-    for (int j = 0; j < 3; j++)
-    {
-      V(i,j) = m_nodes(i * 3 + j);
-    }
-  }
-  int ss = num_nodes_per_elem_type(m_element_type);
-  T.resize(m_elements.rows()/ss,ss);
-  for (int i = 0; i < m_elements.rows() / ss; i++) 
-  {
-    for (int j = 0; j < ss; j++)
-    {
-      T(i, j) = m_elements(i * ss + j);
-    }
-  }
-  return true;
+    return igl::readMSH(msh,X,Tri,Tet,TriTag,TetTag,XFields,XF,EFields,TriF,TetF);
 }
 
+IGL_INLINE bool igl::readMSH(const std::string &msh,
+                Eigen::MatrixXd &X,
+                Eigen::MatrixXi &Tri
+                )
+{
+    Eigen::MatrixXi Tet;
+    Eigen::VectorXi TetTag;
+    Eigen::VectorXi TriTag;
+
+    std::vector<std::string>     XFields;
+    std::vector<Eigen::MatrixXd> XF;
+    std::vector<std::string>     EFields;
+    std::vector<Eigen::MatrixXd> TriF;
+    std::vector<Eigen::MatrixXd> TetF;
 
-#ifdef IGL_STATIC_LIBRARY
-// Explicit template instantiation
-// generated by autoexplicit.sh
-template bool igl::readMSH<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-#endif
+    return igl::readMSH(msh,X,Tri,Tet,TriTag,TetTag,XFields,XF,EFields,TriF,TetF);
+}

include/igl/readMSH.h

@@ -1,38 +1,103 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2018 Alec Jacobson <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
-// obtain one at http://mozilla.org/MPL/2.0/.
-#ifndef IGL_READMSH_H
-#define IGL_READMSH_H
+/* high level interface for MshLoader.h/.cpp */
+
+/* Copyright (C) 2020 Vladimir Fonov <[email protected]> */
+/*
+/* This Source Code Form is subject to the terms of the Mozilla */
+/* Public License v. 2.0. If a copy of the MPL was not distributed */
+/* with this file, You can obtain one at http://mozilla.org/MPL/2.0/. */
+#ifndef IGL_READ_MSH_H
+#define IGL_READ_MSH_H
 #include "igl_inline.h"
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 #include <string>
+#include <vector>
+
 
-namespace igl
+namespace igl 
 {
-  // Read a mesh (e.g., tet mesh) from a gmsh .msh file
-  // 
-  // Inputs:
-  //   filename  path to .msh file
-  // Outputs:
-  //    V  #V by 3 list of 3D mesh vertex positions
-  //    T  #T by ss list of 3D ss-element indices into V (e.g., ss=4 for tets)
-  // Returns true on success
-  template <
-    typename DerivedV,
-    typename DerivedT>
-  IGL_INLINE bool readMSH(
-    const std::string & filename,
-    Eigen::PlainObjectBase<DerivedV> & V,
-    Eigen::PlainObjectBase<DerivedT> & T);
+    // read triangle surface mesh and tetrahedral volume mesh from .msh file
+    // Inputs:
+    //   msh - file name
+    // Outputs: 
+    //   X  eigen double matrix of vertex positions  #X by 3
+    //   Tri  #Tri eigen integer matrix of triangular faces indices into vertex positions
+    //   Tet  #Tet eigen integer matrix of tetrahedral indices into vertex positions
+    //   TriTag #Tri eigen integer vector of tags associated with surface faces
+    //   TetTag #Tet eigen integer vector of tags associated with volume elements
+    //   XFields #XFields list of strings with field names associated with nodes
+    //   XF      #XFields list of eigen double matrices, fields associated with nodes 
+    //   EFields #EFields list of strings with field names associated with elements
+    //   TriF    #EFields list of eigen double matrices, fields associated with surface elements
+    //   TetF    #EFields list of eigen double matrices, fields associated with volume elements
+    // Known bugs: 
+    //     only version 2.2 of .msh file is supported (gmsh 3.X)
+    //     only triangle surface elements and tetrahedral volumetric elements are supported
+    //     only 3D information is supported
+    //     only the 1st tag per element is returned (physical) 
+    //     same element fields are expected to be associated with surface elements and volumetric elements
+    IGL_INLINE bool readMSH(const std::string &msh,
+                Eigen::MatrixXd &X,
+                Eigen::MatrixXi &Tri,
+                Eigen::MatrixXi &Tet,
+                Eigen::VectorXi &TriTag,
+                Eigen::VectorXi &TetTag,
+                std::vector<std::string> &XFields,
+                std::vector<Eigen::MatrixXd> &XF,
+                std::vector<std::string> &EFields,
+                std::vector<Eigen::MatrixXd> &TriF,
+                std::vector<Eigen::MatrixXd> &TetF
+                );
+
+    // read triangle surface mesh and tetrahedral volume mesh from .msh file
+    // ignoring any fields
+    // Inputs:
+    //   msh - file name
+    // Outputs: 
+    //   X  eigen double matrix of vertex positions  #X by 3
+    //   Tri  #Tri eigen integer matrix of triangular faces indices into vertex positions
+    //   Tet  #Tet eigen integer matrix of tetrahedral indices into vertex positions
+    //   TriTag #Tri eigen integer vector of tags associated with surface faces
+    //   TetTag #Tet eigen integer vector of tags associated with volume elements
+    IGL_INLINE bool readMSH(const std::string &msh,
+                Eigen::MatrixXd &X,
+                Eigen::MatrixXi &Tri,
+                Eigen::MatrixXi &Tet,
+                Eigen::VectorXi &TriTag,
+                Eigen::VectorXi &TetTag
+                );
+    
+    // read triangle surface mesh and tetrahedral volume mesh from .msh file
+    // ignoring any fields, and any volumetric elements
+    // Inputs:
+    //   msh - file name
+    // Outputs: 
+    //   X  eigen double matrix of vertex positions  #X by 3
+    //   Tri  #Tri eigen integer matrix of triangular faces indices into vertex positions
+    //   TriTag #Tri eigen integer vector of tags associated with surface faces
+    IGL_INLINE bool readMSH(const std::string &msh,
+                Eigen::MatrixXd &X,
+                Eigen::MatrixXi &Tri,
+                Eigen::VectorXi &TriTag
+                );
+    
+    // read triangle surface mesh and tetrahedral volume mesh from .msh file
+    // ignoring any fields, and any volumetric elements and tags
+    // Inputs:
+    //   msh - file name
+    // Outputs: 
+    //   X  eigen double matrix of vertex positions  #X by 3
+    //   Tri  #Tri eigen integer matrix of triangular faces indices into vertex positions
+    IGL_INLINE bool readMSH(const std::string &msh,
+                Eigen::MatrixXd &X,
+                Eigen::MatrixXi &Tri
+                );
+
 }
 
 
 #ifndef IGL_STATIC_LIBRARY
 #  include "readMSH.cpp"
 #endif
-#endif
+
+#endif //IGL_READ_MSH_H

include/igl/writeMSH.cpp

@@ -0,0 +1,149 @@
+/* high level interface for MshSaver*/
+
+/* Copyright (C) 2020 Vladimir Fonov <[email protected]> */
+/*
+/* This Source Code Form is subject to the terms of the Mozilla */
+/* Public License v. 2.0. If a copy of the MPL was not distributed */
+/* with this file, You can obtain one at http://mozilla.org/MPL/2.0/. */
+
+#include "writeMSH.h"
+#include "MshSaver.h"
+#include "MshLoader.h"
+#include <iostream>
+
+
+namespace igl {
+namespace internal {
+
+    // helper function, appends contents of Eigen matrix to an std::vector, in RowMajor fashion
+    template <typename T, typename Derived>
+    void append_mat_to_vec(std::vector<T> &vec, const Eigen::PlainObjectBase<Derived> & mat)
+    {
+        size_t st = vec.size();
+        vec.resize(st + mat.size());
+
+        Eigen::Map< Eigen::Matrix<typename Derived::Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >
+            _map_vec( reinterpret_cast<T *>( vec.data() + st ), mat.rows(), mat.cols() );
+        _map_vec = mat;
+    }
+
+}
+}
+
+IGL_INLINE bool igl::writeMSH(
+             const std::string    &msh,
+             const Eigen::MatrixXd &X,
+             const Eigen::MatrixXi &Tri,
+             const Eigen::MatrixXi &Tet,
+             const Eigen::MatrixXi &TriTag,
+             const Eigen::MatrixXi &TetTag,
+             const std::vector<std::string>     &XFields,
+             const std::vector<Eigen::MatrixXd> &XF,
+             const std::vector<std::string>     &EFields,
+             const std::vector<Eigen::MatrixXd> &TriF,
+             const std::vector<Eigen::MatrixXd> &TetF
+             )
+{
+    using namespace internal;
+
+    try
+    {
+        // error checks
+        if(!XFields.empty())
+        {
+            if(XFields.size()!=XF.size())
+                throw std::invalid_argument("Vertex field count mismatch");
+            for(int i=0;i<XFields.size();++i)
+                if(XF[i].rows()!=X.rows())
+                    throw std::invalid_argument("Vertex field size mismatch");
+        }
+
+        if(!EFields.empty())
+        {
+            if(EFields.size()!=TriF.size())
+                throw std::invalid_argument("Triangle field count mismatch");
+            if(EFields.size()!=TetF.size())
+                throw std::invalid_argument("Tetrahedra field count mismatch");
+
+            for(int i=0;i<EFields.size();++i)
+            {
+                if(TriF[i].rows()!=Tri.rows())
+                    throw std::invalid_argument("Triangle field size mismatch");
+                if(TetF[i].rows()!=Tet.rows())
+                    throw std::invalid_argument("Tetrahedra field size mismatch");
+            }
+        }
+
+        // this is not the most optimal , it would be faster to modify RRMshSaver to work with Eiged data types
+        std::vector<double> _X;
+        append_mat_to_vec(_X, X);
+
+        std::vector<int> _Tri_Tet;
+        append_mat_to_vec( _Tri_Tet, Tri);
+        append_mat_to_vec( _Tri_Tet, Tet);
+
+        std::vector<int> _Tri_Tet_len(Tri.rows(), 3); //each is 3 elements long
+        _Tri_Tet_len.insert(_Tri_Tet_len.end(), Tet.rows(), 4);
+
+        std::vector<int> _Tri_Tet_type(Tri.rows(), MshLoader::ELEMENT_TRI);
+        _Tri_Tet_type.insert(_Tri_Tet_type.end(), Tet.rows(), MshLoader::ELEMENT_TET);
+
+        std::vector<int> _Tri_Tet_tag;
+        append_mat_to_vec(_Tri_Tet_tag, TriTag);
+        append_mat_to_vec(_Tri_Tet_tag, TetTag);
+
+
+        igl::MshSaver msh_saver(msh, true);
+        msh_saver.save_mesh( _X,
+            _Tri_Tet,
+            _Tri_Tet_len,
+            _Tri_Tet_type,
+            _Tri_Tet_tag);
+
+        // append vertex data
+        for(size_t i=0;i<XFields.size();++i)
+        {
+            assert(X.rows()==XF[i].rows());
+
+            std::vector<double> _XF;
+            append_mat_to_vec(_XF, XF[i]);
+
+            if(XF[i].cols() == 1)
+                msh_saver.save_scalar_field(XFields[i], _XF );
+            else if(XF[i].cols() == 3)
+                msh_saver.save_vector_field(XFields[i], _XF );
+            else
+            {
+                throw std::invalid_argument("unsupported vertex field dimensionality");
+            }
+        }
+
+        // append node data
+        for(size_t i=0; i<EFields.size(); ++i)
+        {
+            assert(TriF[i].cols() == TetF[i].cols());
+            assert(TriF[i].rows() == Tri.rows());
+            assert(TetF[i].rows() == Tet.rows());
+
+            std::vector<double> _EF;
+            append_mat_to_vec(_EF, TriF[i]);
+            append_mat_to_vec(_EF, TetF[i]);
+
+            assert(_EF.size() == (TriF[i].size()+TetF[i].size()));
+
+            if( TriF[i].cols() == 1 )
+                msh_saver.save_elem_scalar_field(EFields[i], _EF );
+            else if( TriF[i].cols() == 3 )
+                msh_saver.save_elem_vector_field(EFields[i], _EF );
+            else
+            {
+                throw std::invalid_argument("unsupported node field dimensionality");
+            }
+        }
+    } catch(const std::exception& e) {
+        std::cerr << e.what() << std::endl;
+        return false;
+    }
+    return true;
+}
+

include/igl/writeMSH.h

@@ -0,0 +1,57 @@
+/* high level interface for MshSaver */
+
+/* Copyright (C) 2020 Vladimir Fonov <[email protected]> */
+/*
+/* This Source Code Form is subject to the terms of the Mozilla */
+/* Public License v. 2.0. If a copy of the MPL was not distributed */
+/* with this file, You can obtain one at http://mozilla.org/MPL/2.0/. */
+#ifndef IGL_WRITE_MSH_H
+#define IGL_WRITE_MSH_H
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+#include <string>
+#include <vector>
+
+namespace igl
+{
+    // write triangle surface mesh and tetrahedral volume mesh to .msh file
+    // Inputs:
+    //   msh - file name
+    //   X  eigen double matrix of vertex positions  #X by 3
+    //   Tri  #Tri eigen integer matrix of triangular faces indices into vertex positions
+    //   Tet  #Tet eigen integer matrix of tetrahedral indices into vertex positions
+    //   TriTag #Tri eigen integer vector of tags associated with surface faces
+    //   TetTag #Tet eigen integer vector of tags associated with volume elements
+    //   XFields #XFields list of strings with field names associated with nodes
+    //   XF      #XFields list of eigen double matrices, fields associated with nodes 
+    //   EFields #EFields list of strings with field names associated with elements
+    //   TriF    #EFields list of eigen double matrices, fields associated with surface elements
+    //   TetF    #EFields list of eigen double matrices, fields associated with volume elements
+    // Known bugs:
+    //     files are always stored in binary format
+    //     file format is 2.2
+    //     only triangle surface elements and tetrahedral volumetric elements are supported
+    //     only 3D information is supported
+    //     the tag id is duplicated for physical (0) and elementary (1)
+    //     same element fields are expected to be associated with surface elements and volumetric elements
+    IGL_INLINE bool writeMSH(const std::string   &msh,
+                const Eigen::MatrixXd &X,
+                const Eigen::MatrixXi &Tri,
+                const Eigen::MatrixXi &Tet,
+                const Eigen::MatrixXi &TriTag,
+                const Eigen::MatrixXi &TetTag,
+                const std::vector<std::string>     &XFields,
+                const std::vector<Eigen::MatrixXd> &XF,
+                const std::vector<std::string>     &EFields,
+                const std::vector<Eigen::MatrixXd> &TriF,
+                const std::vector<Eigen::MatrixXd> &TetF
+                );
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "writeMSH.cpp"
+#endif
+
+#endif //IGL_WRITE_MSH_H

tests/include/igl/MshLoader.cpp

@@ -0,0 +1,59 @@
+#include <test_common.h>
+
+#include <catch2/catch.hpp>
+#include <igl/MshLoader.h>
+
+
+TEST_CASE("MshLoader","[igl]")
+{
+    igl::MshLoader msh_loader(test_common::data_path("sphere_lowres_TMS_1-0001_Magstim_70mm_Fig8_nii_scalar.msh"));
+
+    //
+    // size of nodes - 3xNodes (x,y,z)
+    //
+    REQUIRE( msh_loader.get_nodes().size()    == (398+4506)*3 );
+
+    // 
+    // size of elements: number of triangles x 3 
+    //                   number of tetrahedra x 4
+    // 
+    REQUIRE( msh_loader.get_elements().size()         == (8988*3 + 25937*4) );
+    // all these should have the same size
+    REQUIRE( msh_loader.get_elements_lengths().size() == (8988 + 25937) );
+    REQUIRE( msh_loader.get_elements_ids().size()     == (8988 + 25937) );
+    REQUIRE( msh_loader.get_elements_types().size()   == (8988 + 25937) );
+
+    // the test file should have identity map
+    REQUIRE( msh_loader.is_element_map_identity() );
+
+    // there should be 12 tags: 6 surface structures and 6 volume structures 
+    // index 1st structure tag: Volumes
+    msh_loader.index_structures(1); 
+
+    REQUIRE( msh_loader.is_element_map_identity() );
+    REQUIRE( msh_loader.get_structures().size() == 6+6 );
+
+    // verify that all triangle elements have length of 3 and 
+    // tetrahedral elements have length 4
+
+    for(auto t =msh_loader.get_structures().begin(); 
+             t!=msh_loader.get_structures().end(); ++t )
+    {
+        REQUIRE(( t->el_type == igl::MshLoader::ELEMENT_TRI || 
+                  t->el_type == igl::MshLoader::ELEMENT_TET ));
+
+        auto element_range = msh_loader.get_structure_index().equal_range(*t);
+
+        for(auto e=element_range.first; e != element_range.second; ++e)
+        { 
+            if( t->el_type == igl::MshLoader::ELEMENT_TRI )
+            {
+                REQUIRE( msh_loader.get_elements_lengths()[e->second] == 3 );
+            } else {
+                // t->first.el_type==MshLoader::ELEMENT_TET
+                REQUIRE( msh_loader.get_elements_lengths()[e->second] == 4 );
+            }
+        }
+    }
+}
+

tests/include/igl/MshSaver.cpp

@@ -0,0 +1,72 @@
+#include <test_common.h>
+
+#include <catch2/catch.hpp>
+#include <igl/MshLoader.h>
+#include <igl/MshSaver.h>
+
+TEST_CASE("MshSaver","[igl]")
+{
+    igl::MshLoader msh_loader1(test_common::data_path("sphere_lowres_TMS_1-0001_Magstim_70mm_Fig8_nii_scalar.msh"));
+    REQUIRE( msh_loader1.is_element_map_identity());
+    
+    igl::MshSaver msh_saver("test_binary_sphere.msh",true);
+    msh_saver.save_mesh( msh_loader1.get_nodes(), 
+        msh_loader1.get_elements(),
+        msh_loader1.get_elements_lengths(), 
+        msh_loader1.get_elements_types(),
+        msh_loader1.get_elements_tags()[1]);
+
+    for(size_t i=0;i<msh_loader1.get_node_fields_names().size();++i)
+    {
+        if(msh_loader1.get_node_fields_components()[i]==3)
+            msh_saver.save_vector_field(msh_loader1.get_node_fields_names()[i],msh_loader1.get_node_fields()[i]);
+        else
+            msh_saver.save_scalar_field(msh_loader1.get_node_fields_names()[i],msh_loader1.get_node_fields()[i]);
+    }
+
+    for(size_t i=0;i<msh_loader1.get_element_fields_names().size();++i)
+    {
+        if(msh_loader1.get_element_fields_components()[i]==3)
+            msh_saver.save_elem_vector_field(msh_loader1.get_element_fields_names()[i],msh_loader1.get_element_fields()[i]);
+        else
+            msh_saver.save_elem_scalar_field(msh_loader1.get_element_fields_names()[i],msh_loader1.get_element_fields()[i]);
+    }
+
+    igl::MshLoader msh_loader("test_binary_sphere.msh");
+
+    for(size_t i=0;i<msh_loader1.get_elements().size();++i)
+        REQUIRE(msh_loader.get_elements()[i] == msh_loader1.get_elements()[i]);
+
+    for(size_t i=0;i<msh_loader1.get_elements_lengths().size();++i)
+        REQUIRE(msh_loader.get_elements_lengths()[i] == msh_loader1.get_elements_lengths()[i]);
+
+    for(size_t i=0;i<msh_loader1.get_elements_types().size();++i)
+        REQUIRE(msh_loader.get_elements_types()[i] == msh_loader1.get_elements_types()[i]);
+
+    for(size_t j=0;j<2;++j)
+        for(size_t i=0;i<msh_loader1.get_elements_tags()[j].size();++i)
+            REQUIRE(msh_loader.get_elements_tags()[j][i] == msh_loader1.get_elements_tags()[j][i]);
+
+    REQUIRE(msh_loader.get_node_fields_names().size() == msh_loader1.get_node_fields_names().size());
+
+    for(size_t i=0;i<msh_loader1.get_node_fields_names().size();++i)
+    {
+        REQUIRE(msh_loader.get_node_fields_names()[i] == msh_loader1.get_node_fields_names()[i]);
+        REQUIRE(msh_loader.get_node_fields_components()[i] == msh_loader.get_node_fields_components()[i]);
+
+        for(size_t j=0;j<msh_loader1.get_node_fields()[i].size();++j)
+            REQUIRE(msh_loader1.get_node_fields()[i][j] == msh_loader.get_node_fields()[i][j]);
+    }
+
+    REQUIRE(msh_loader.get_element_fields_names().size() == msh_loader1.get_element_fields_names().size());
+    for(size_t i=0;i<msh_loader1.get_element_fields_names().size();++i)
+    {
+        REQUIRE(msh_loader.get_element_fields_names()[i] == msh_loader1.get_element_fields_names()[i]);
+        REQUIRE(msh_loader.get_element_fields_components()[i] == msh_loader1.get_element_fields_components()[i]);
+
+        for(size_t j=0;j<msh_loader1.get_element_fields()[i].size();++j)
+            REQUIRE(msh_loader1.get_element_fields()[i][j] == msh_loader.get_element_fields()[i][j]);
+    }
+
+}
+

tests/include/igl/readMSH.cpp

@@ -0,0 +1,72 @@
+#include <test_common.h>
+
+#include <catch2/catch.hpp>
+
+#include <igl/readMSH.h>
+
+#include <set>
+
+
+TEST_CASE("readMSH","[igl]")
+{
+    Eigen::MatrixXd X;
+    Eigen::MatrixXi Tri;
+    Eigen::MatrixXi Tet;
+    Eigen::VectorXi TriTag;
+    Eigen::VectorXi TetTag;
+
+    std::vector<std::string> XFields;
+    std::vector<std::string> EFields;
+
+    std::vector<Eigen::MatrixXd> XF;
+    std::vector<Eigen::MatrixXd> TriF;
+    std::vector<Eigen::MatrixXd> TetF;
+
+    REQUIRE(igl::readMSH(test_common::data_path("sphere_lowres_TMS_1-0001_Magstim_70mm_Fig8_nii_scalar.msh"), 
+        X, Tri, Tet, TriTag, TetTag, XFields, XF, EFields, TriF, TetF));
+
+    REQUIRE(X.cols() == 3);
+    REQUIRE(X.rows() == (398+4506));
+
+    REQUIRE(Tri.cols() == 3);
+    REQUIRE(Tri.rows() == 8988);
+    REQUIRE(TriTag.rows() == 8988);
+
+    // determine all tags
+    std::set<int> tri_tags_unique;
+    for(size_t i=0; i<TriTag.rows(); ++i) tri_tags_unique.insert(TriTag(i));
+    REQUIRE(tri_tags_unique.size()==6);
+
+    // make sure we have tags 1001-1006
+    for(int i=1;i<6;++i)
+        REQUIRE(tri_tags_unique.find(i+1000)!=std::end(tri_tags_unique));
+
+    REQUIRE(Tet.cols() == 4);
+    REQUIRE(Tet.rows() == 25937);
+    REQUIRE(TetTag.rows() == 25937);
+    // determine all tags
+    std::set<int> tet_tags_unique;
+    for(size_t i=0; i<TetTag.rows(); ++i) tet_tags_unique.insert(TetTag(i));
+    REQUIRE(tet_tags_unique.size()==6);
+
+    // make sure we have tags 1-6
+    for(int i=1;i<6;++i)
+        REQUIRE(tet_tags_unique.find(i)!=std::end(tet_tags_unique));
+
+    REQUIRE(XFields.size()==0);
+    REQUIRE(EFields.size()==1);
+
+    REQUIRE(EFields[0]=="normE");
+
+    //make sure field sizes are correct
+    REQUIRE(XF.size()==0);
+    REQUIRE(TriF.size()==1);
+    REQUIRE(TetF.size()==1);
+
+    // normE , scalar field
+    REQUIRE(TriF[0].cols()==1);
+    REQUIRE(TriF[0].rows()==8988);
+    REQUIRE(TetF[0].cols()==1);
+    REQUIRE(TetF[0].rows()==25937);
+}
+

tests/include/igl/writeMSH.cpp

@@ -0,0 +1,122 @@
+#include <test_common.h>
+#include <catch2/catch.hpp>
+
+#include <igl/readMSH.h>
+#include <igl/writeMSH.h>
+
+#include <set>
+
+TEST_CASE("writeMSH","[igl]")
+{
+    Eigen::MatrixXd X;
+    Eigen::MatrixXi Tri;
+    Eigen::MatrixXi Tet;
+    Eigen::VectorXi TriTag;
+    Eigen::VectorXi TetTag;
+
+    std::vector<std::string> XFields;
+    std::vector<std::string> EFields;
+
+    std::vector<Eigen::MatrixXd> XF;
+    std::vector<Eigen::MatrixXd> TriF;
+    std::vector<Eigen::MatrixXd> TetF;
+
+    // load source data
+    REQUIRE(igl::readMSH(test_common::data_path("sphere_lowres_TMS_1-0001_Magstim_70mm_Fig8_nii_scalar.msh"), 
+        X, Tri, Tet, TriTag, TetTag, XFields, XF, EFields, TriF, TetF));
+
+    // save check data
+    REQUIRE(igl::writeMSH("test_binary_sphere_2.msh",
+        X, Tri, Tet, TriTag, TetTag, XFields, XF, EFields, TriF, TetF));
+
+    // load again
+    Eigen::MatrixXd _X;
+    Eigen::MatrixXi _Tri;
+    Eigen::MatrixXi _Tet;
+    Eigen::VectorXi _TriTag;
+    Eigen::VectorXi _TetTag;
+
+    std::vector<std::string> _XFields;
+    std::vector<std::string> _EFields;
+
+    std::vector<Eigen::MatrixXd> _XF;
+    std::vector<Eigen::MatrixXd> _TriF;
+    std::vector<Eigen::MatrixXd> _TetF;
+
+    REQUIRE(igl::readMSH("test_binary_sphere_2.msh", 
+        _X, _Tri, _Tet, _TriTag, _TetTag, _XFields, _XF, _EFields, _TriF, _TetF));
+
+    REQUIRE(_X.size()   == X.size());
+    REQUIRE(_Tri.size() == Tri.size());
+    REQUIRE(_Tet.size() == Tet.size());
+    REQUIRE(_TriTag.size() == TriTag.size());
+    REQUIRE(_TetTag.size() == TetTag.size());
+    REQUIRE(_XFields.size() == XFields.size());
+    REQUIRE(_XF.size() == XF.size());
+    REQUIRE(_EFields.size() == EFields.size());
+    REQUIRE(_TriF.size() == TriF.size());
+    REQUIRE(_TetF.size() == TetF.size());
+
+    REQUIRE(_X.cols() == 3);
+    REQUIRE(_X.rows() == (398+4506));
+
+    REQUIRE(_Tri.cols() == 3);
+    REQUIRE(_Tri.rows() == 8988);
+    REQUIRE(_TriTag.rows() == 8988);
+
+    REQUIRE(_Tet.cols() == 4);
+    REQUIRE(_Tet.rows() == 25937);
+    REQUIRE(_TetTag.rows() == 25937);
+
+    //make sure field sizes are correct
+    REQUIRE(_XF.size()==0);
+    REQUIRE(_TriF.size()==1);
+    REQUIRE(_TetF.size()==1);
+
+    // normE , scalar field
+    REQUIRE(_TriF[0].cols()==1);
+    REQUIRE(_TriF[0].rows()==8988);
+
+    REQUIRE(_TetF[0].cols()==1);
+    REQUIRE(_TetF[0].rows()==25937);
+
+    // check the contents too
+    for(size_t i=0;i<X.rows();++i)
+        for(size_t j=0;j<X.cols();++j)
+            REQUIRE(_X(i,j) == X(i,j));
+
+    for(size_t i=0;i<Tri.rows();++i)
+        for(size_t j=0;j<Tri.cols();++j)
+            REQUIRE(_Tri(i,j) == Tri(i,j));
+
+    for(size_t i=0;i<Tet.rows();++i)
+        for(size_t j=0;j<Tet.cols();++j)
+            REQUIRE(_Tet(i,j) == Tet(i,j));
+
+    for(size_t i=0;i<XFields.size();++i)
+    {
+        REQUIRE(XFields[i]==_XFields[i]);
+        REQUIRE(XF[i].rows()==_XF[i].rows());
+        REQUIRE(XF[i].cols()==_XF[i].cols());
+        for(size_t j=0;j<XF[i].rows();++j)
+            for(size_t k=0;k<XF[i].cols();++k)
+                REQUIRE(XF[i](j,k)==_XF[i](j,k));
+    }
+
+    for(size_t i=0;i<EFields.size();++i)
+    {
+        REQUIRE(EFields[i]==_EFields[i]);
+        REQUIRE(TriF[i].rows()==_TriF[i].rows());
+        REQUIRE(TriF[i].cols()==_TriF[i].cols());
+
+        for(size_t j=0;j<TriF[i].rows();++j)
+            for(size_t k=0;k<TriF[i].cols();++k)
+                REQUIRE(TriF[i](j,k)==_TriF[i](j,k));
+        
+        REQUIRE(TetF[i].rows()==_TetF[i].rows());
+        REQUIRE(TetF[i].cols()==_TetF[i].cols());
+        for(size_t j=0;j<TetF[i].rows();++j)
+            for(size_t k=0;k<TetF[i].cols();++k)
+                REQUIRE(TetF[i](j,k)==_TetF[i](j,k));
+    }
+}

tutorial/110_MshView/CMakeLists.txt

@@ -0,0 +1,5 @@
+get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+project(${PROJECT_NAME})
+
+add_executable(${PROJECT_NAME} main.cpp)
+target_link_libraries(${PROJECT_NAME} igl::core igl::opengl igl::opengl_glfw tutorials)

tutorial/110_MshView/main.cpp

@@ -0,0 +1,111 @@
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/barycenter.h>
+#include <igl/colormap.h>
+
+#include <igl/readMSH.h>
+#include <igl/readMESH.h>
+
+
+Eigen::MatrixXd X,B;
+Eigen::MatrixXi Tri;
+Eigen::MatrixXi Tet;
+Eigen::VectorXi TriTag;
+Eigen::VectorXi TetTag;
+
+Eigen::VectorXd D;
+
+std::vector<std::string> XFields;
+std::vector<std::string> EFields;
+
+std::vector<Eigen::MatrixXd> XF;
+std::vector<Eigen::MatrixXd> TriF;
+std::vector<Eigen::MatrixXd> TetF;
+
+
+// This function is called every time a keyboard button is pressed
+bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  if (key >= '1' && key <= '9')
+  {
+    double t = double((key - '1')+1) / 9.0;
+
+    VectorXd v = B.col(2).array() - B.col(2).minCoeff();
+    v /= v.col(0).maxCoeff();
+
+    vector<int> s;
+
+    for (unsigned i=0; i<v.size();++i)
+      if (v(i) < t && v(i)>(t-0.1)) // select a thick slab
+        s.push_back(i);
+
+    MatrixXd V_temp(s.size()*4,3);
+    MatrixXi F_temp(s.size()*4,3);
+    VectorXd D_temp(s.size()*4);
+
+    for (unsigned i=0; i<s.size();++i)
+    {
+      V_temp.row(i*4+0) = X.row(Tet(s[i],0));
+      V_temp.row(i*4+1) = X.row(Tet(s[i],1));
+      V_temp.row(i*4+2) = X.row(Tet(s[i],2));
+      V_temp.row(i*4+3) = X.row(Tet(s[i],3));
+
+      F_temp.row(i*4+0) << (i*4)+0, (i*4)+1, (i*4)+3;
+      F_temp.row(i*4+1) << (i*4)+0, (i*4)+2, (i*4)+1;
+      F_temp.row(i*4+2) << (i*4)+3, (i*4)+2, (i*4)+0;
+      F_temp.row(i*4+3) << (i*4)+1, (i*4)+2, (i*4)+3;
+
+      D_temp(i*4+0) = D(s[i]);
+      D_temp(i*4+1) = D(s[i]);
+      D_temp(i*4+2) = D(s[i]);
+      D_temp(i*4+3) = D(s[i]);
+    }
+
+    viewer.data().clear();
+    viewer.data().set_mesh(V_temp, F_temp);
+
+    Eigen::MatrixXd C;
+    igl::colormap(igl::COLOR_MAP_TYPE_VIRIDIS, D_temp, true, C);
+    viewer.data().set_face_based(true);
+    viewer.data().set_colors(C);
+
+  }
+
+  return false;
+}
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+
+  igl::readMSH(argc > 1 ? argv[1] : TUTORIAL_SHARED_PATH "/hand.msh", X, Tri, Tet, TriTag, TetTag, XFields, XF, EFields, TriF, TetF);
+
+  for(auto i:EFields)
+    std::cout<<i<<"\t";
+  std::cout<<std::endl;
+
+  // search for a predefined field name "E"
+  for(int i=0;i<EFields.size();++i)
+  {
+    if(EFields[i]=="E")
+      D = TetF[i].rowwise().norm(); // take a row-wise norm
+  }
+  std::cout<<"D:"<<D.rows()<<"x"<<D.cols()<<std::endl;
+
+  // generate fake data
+  if(D.rows()==0)
+    D = TetTag.cast<double>();
+
+  // Compute barycenters
+  igl::barycenter(X, Tet, B);
+
+  // Plot the generated mesh
+  igl::opengl::glfw::Viewer viewer;
+
+  viewer.callback_key_down = &key_down;
+  key_down(viewer,'5',0);
+  viewer.launch();
+}

tutorial/720_BlueNoise/CMakeLists.txt

@@ -1,5 +1,5 @@
 get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
 project(${PROJECT_NAME})
 
-add_executable(${PROJECT_NAME}_bin main.cpp)
-target_link_libraries(${PROJECT_NAME}_bin igl::core igl::opengl igl::opengl_glfw tutorials)
+add_executable(${PROJECT_NAME} main.cpp)
+target_link_libraries(${PROJECT_NAME} igl::core igl::opengl igl::opengl_glfw tutorials)

tutorial/CMakeLists.txt

@@ -58,10 +58,13 @@ if(TUTORIALS_CHAPTER1)
   if(LIBIGL_WITH_OPENGL_GLFW_IMGUI)
     add_subdirectory("105_Overlays")
     add_subdirectory("106_ViewerMenu")
-    add_subdirectory("109_ImGuizmo")
   endif()
   add_subdirectory("107_MultipleMeshes")
   add_subdirectory("108_MultipleViews")
+  if(LIBIGL_WITH_OPENGL_GLFW_IMGUI)
+    add_subdirectory("109_ImGuizmo")
+  endif()
+  add_subdirectory("110_MshView")
 endif()
 
 # Chapter 2