Merge branch 'master' of https://github.com/libigl/libigl

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/FileEncoding.h

@@ -0,0 +1,21 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Jérémie Dumas <[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_FILEENCODING_H
+#define IGL_FILEENCODING_H
+
+namespace igl
+{
+
+enum class FileEncoding {
+  Binary,
+  Ascii
+};
+
+}
+
+#endif

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/average_from_edges_onto_vertices.cpp

@@ -0,0 +1,68 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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 "average_from_edges_onto_vertices.h"
+
+template<typename DerivedF,typename DerivedE,typename DerivedoE,
+typename DeriveduE,typename DeriveduV>
+IGL_INLINE void
+igl::average_from_edges_onto_vertices(
+  const Eigen::MatrixBase<DerivedF> &F,
+  const Eigen::MatrixBase<DerivedE> &E,
+  const Eigen::MatrixBase<DerivedoE> &oE,
+  const Eigen::MatrixBase<DeriveduE> &uE,
+  Eigen::PlainObjectBase<DeriveduV> &uV)
+{
+  using Scalar = typename DeriveduE::Scalar;
+  using VecX = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+  using Int = typename DerivedF::Scalar;
+  
+  assert(E.rows()==F.rows() && "E does not match dimensions of F.");
+  assert(oE.rows()==F.rows() && "oE does not match dimensions of F.");
+  assert(E.cols()==3 && F.cols()==3 && oE.cols()==3 &&
+   "This method is for triangle meshes.");
+  
+  const Int n = F.maxCoeff()+1;
+  
+  VecX edgesPerVertex(n);
+  edgesPerVertex.setZero();
+  uV.resize(n,1);
+  uV.setZero();
+  
+  for(Eigen::Index i=0; i<F.rows(); ++i) {
+    for(int j=0; j<3; ++j) {
+      if(oE(i,j)<0) {
+        continue;
+      }
+      const Int e = E(i,j);
+      const Int vi=F(i,(j+1)%3), vj=F(i,(j+2)%3);
+      
+      //Count vertex valence
+      ++edgesPerVertex(vi);
+      ++edgesPerVertex(vj);
+      
+      //Average uE value onto vertices
+      uV(vi) += uE(e);
+      uV(vj) += uE(e);
+    }
+  }
+  
+  //Divide by valence
+  for(Int i=0; i<n; ++i) {
+    const Scalar valence = edgesPerVertex(i);
+    if(valence>0) {
+      uV(i) /= valence;
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::average_from_edges_onto_vertices<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::PartialReduxExpr<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::internal::member_norm<double>, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::PartialReduxExpr<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::internal::member_norm<double>, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::average_from_edges_onto_vertices<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::average_from_edges_onto_vertices<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/average_from_edges_onto_vertices.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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_AVERAGE_FROM_EDGES_ONTO_VERTICES_H
+#define IGL_AVERAGE_FROM_EDGES_ONTO_VERTICES_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+namespace igl
+{
+  // Move a scalar field defined on edges to vertices by averaging
+  //
+  // Input:
+  // F: triangle mesh connectivity
+  // E, oE: mapping from halfedges to edges and orientation as generated by
+  //    orient_halfedges
+  // uE: scalar field defined on edges, one per edge
+  //
+  // Output:
+  // uV: scalar field defined on vertices
+  template<typename DerivedF,typename DerivedE,typename DerivedoE,
+  typename DeriveduE,typename DeriveduV>
+  IGL_INLINE void average_from_edges_onto_vertices(
+    const Eigen::MatrixBase<DerivedF> &F,
+    const Eigen::MatrixBase<DerivedE> &E,
+    const Eigen::MatrixBase<DerivedoE> &oE,
+    const Eigen::MatrixBase<DeriveduE> &uE,
+    Eigen::PlainObjectBase<DeriveduV> &uV);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "average_from_edges_onto_vertices.cpp"
+#endif
+
+#endif

include/igl/blue_noise.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2020 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 
+//
+// 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 "blue_noise.h"
 #include "doublearea.h"
@@ -126,6 +126,7 @@ namespace igl
     std::unordered_map<BlueNoiseKeyType,int> & S,
     std::vector<int> & active)
   {
+    assert(M.count(nk));
     auto & Mvec = M.find(nk)->second;
     auto miter = Mvec.begin();
     while(miter != Mvec.end())
@@ -149,7 +150,14 @@ namespace igl
         // back)
         //std::swap(*miter,Mvec.back());
         *miter = Mvec.back();
+        bool was_last = (std::next(miter) == Mvec.end());
         Mvec.pop_back();
+        if (was_last) {
+          // popping from the vector can invalidate the iterator, if it was
+          // pointing to the last element that was popped. Alternatively,
+          // one could use indices directly...
+          miter = Mvec.end();
+        }
       }
     }
     return false;
@@ -259,10 +267,10 @@ IGL_INLINE void igl::blue_noise(
   // function of r and s; and removing the `if S=-1 checks`)
   const Scalar s = r/sqrt(3.0);
 
-  const double area = 
+  const double area =
     [&](){Eigen::VectorXd A;igl::doublearea(V,F,A);return A.array().sum()/2;}();
   // Circle packing in the plane has igl::PI*sqrt(3)/6 efficiency
-  const double expected_number_of_points = 
+  const double expected_number_of_points =
     area * (igl::PI * sqrt(3.0) / 6.0) / (igl::PI * min_r * min_r / 4.0);
 
   // Make a uniform random sampling with 30*expected_number_of_points.
@@ -272,7 +280,7 @@ IGL_INLINE void igl::blue_noise(
   igl::random_points_on_mesh(nx,V,F,XB,XFI,X);
 
   // Rescale so that s = 1
-  Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> Xs = 
+  Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> Xs =
     ((X.rowwise()-X.colwise().minCoeff())/s).template cast<int>();
   const int w = Xs.maxCoeff()+1;
   {
@@ -294,7 +302,7 @@ IGL_INLINE void igl::blue_noise(
   S.reserve(Xs.rows());
   for(int i = 0;i<Xs.rows();i++)
   {
-    BlueNoiseKeyType k = blue_noise_key(w,Xs(i,0),Xs(i,1),Xs(i,2)); 
+    BlueNoiseKeyType k = blue_noise_key(w,Xs(i,0),Xs(i,1),Xs(i,2));
     const auto Miter = M.find(k);
     if(Miter  == M.end())
     {
@@ -309,7 +317,7 @@ IGL_INLINE void igl::blue_noise(
 
   std::vector<int> active;
   // precompute r²
-  // Q: is this necessary? 
+  // Q: is this necessary?
   const double rr = r*r;
   std::vector<int> collected;
   collected.reserve(2.0*expected_number_of_points);

include/igl/cr_vector_curvature_correction.cpp

@@ -0,0 +1,189 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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 "cr_vector_curvature_correction.h"
+
+#include "orient_halfedges.h"
+#include "gaussian_curvature.h"
+
+#include "squared_edge_lengths.h"
+#include "doublearea.h"
+#include "boundary_loop.h"
+#include "internal_angles.h"
+
+#include "PI.h"
+
+
+template <typename DerivedV, typename DerivedF, typename DerivedE,
+typename DerivedOE, typename ScalarK>
+IGL_INLINE void
+igl::cr_vector_curvature_correction(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedE>& E,
+  const Eigen::MatrixBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarK>& K)
+{
+  Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic>
+  l_sq;
+  squared_edge_lengths(V, F, l_sq);
+  cr_vector_curvature_correction_intrinsic(F, l_sq, E, oE, K);
+}
+
+
+template <typename DerivedV, typename DerivedF, typename DerivedE,
+typename DerivedOE, typename ScalarK>
+IGL_INLINE void
+igl::cr_vector_curvature_correction(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedE>& E,
+  Eigen::PlainObjectBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarK>& K)
+{
+  if(E.rows()!=F.rows() || E.cols()!=F.cols() || oE.rows()!=F.rows() ||
+   oE.cols()!=F.cols()) {
+    orient_halfedges(F, E, oE);
+  }
+
+  const Eigen::PlainObjectBase<DerivedE>& cE = E;
+  const Eigen::PlainObjectBase<DerivedOE>& coE = oE;
+  cr_vector_curvature_correction(V, F, cE, coE, K);
+}
+
+
+template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
+typename DerivedOE, typename ScalarK>
+IGL_INLINE void
+igl::cr_vector_curvature_correction_intrinsic(
+  const Eigen::MatrixBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+  const Eigen::MatrixBase<DerivedE>& E,
+  const Eigen::MatrixBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarK>& K)
+{
+  Eigen::Matrix<typename DerivedL_sq::Scalar,Eigen::Dynamic,Eigen::Dynamic>
+  theta;
+  internal_angles_using_squared_edge_lengths(l_sq, theta);
+  
+  cr_vector_curvature_correction_intrinsic(F, l_sq, theta, E, oE, K);
+}
+
+
+template <typename DerivedF, typename DerivedL_sq, typename Derivedtheta,
+typename DerivedE, typename DerivedOE,
+typename ScalarK>
+IGL_INLINE void
+igl::cr_vector_curvature_correction_intrinsic(
+  const Eigen::MatrixBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+  const Eigen::MatrixBase<Derivedtheta>& theta,
+  const Eigen::MatrixBase<DerivedE>& E,
+  const Eigen::MatrixBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarK>& K)
+{
+  // Compute the angle defect kappa, set it to 0 at the boundary
+  const typename DerivedF::Scalar n = F.maxCoeff() + 1;
+  Eigen::Matrix<typename DerivedL_sq::Scalar,Eigen::Dynamic,1> kappa(n);
+  kappa.setZero();
+  for(Eigen::Index i=0; i<F.rows(); ++i) {
+    for(int j=0; j<3; ++j) {
+      kappa(F(i,j)) -= theta(i,j);
+    }
+  }
+  kappa.array() += 2 * PI;
+  std::vector<std::vector<typename DerivedF::Scalar> > b;
+  boundary_loop(F, b);
+  for(const auto& loop : b) {
+    for(auto v : loop) {
+      kappa(v) = 0;
+    }
+  }
+    
+  cr_vector_curvature_correction_intrinsic(F, l_sq, theta, kappa, E, oE, K);
+}
+
+
+template <typename DerivedF, typename DerivedL_sq, typename Derivedtheta,
+typename Derivedkappa, typename DerivedE, typename DerivedOE,
+typename ScalarK>
+IGL_INLINE void
+igl::cr_vector_curvature_correction_intrinsic(
+  const Eigen::MatrixBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+  const Eigen::MatrixBase<Derivedtheta>& theta,
+  const Eigen::MatrixBase<Derivedkappa>& kappa,
+  const Eigen::MatrixBase<DerivedE>& E,
+  const Eigen::MatrixBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarK>& K)
+{
+  assert(F.cols()==3 && "Faces have three vertices");
+  assert(E.rows()==F.rows() && E.cols()==F.cols() && oE.rows()==F.rows() &&
+   theta.rows()==F.rows() && theta.cols()==F.cols() &&
+   oE.cols()==F.cols() && "Wrong dimension in edge vectors");
+  assert(kappa.rows()==F.maxCoeff()+1 &&
+   "Wrong dimension in theta or kappa");
+  
+  const Eigen::Index m = F.rows();
+  const typename DerivedE::Scalar nE = E.maxCoeff() + 1;
+  
+  //Divide kappa by the actual angle sum to weigh consistently.
+  Derivedtheta angleSum = Derivedtheta::Zero(kappa.rows(), 1);
+  for(Eigen::Index i=0; i<F.rows(); ++i) {
+    for(int j=0; j<3; ++j) {
+      angleSum(F(i,j)) += theta(i,j);
+    }
+  }
+  const Eigen::Matrix<typename Derivedkappa::Scalar, Eigen::Dynamic, 1>
+  scaledKappa = kappa.array() / angleSum.array();
+  
+  std::vector<Eigen::Triplet<ScalarK> > tripletList;
+  tripletList.reserve(10*3*m);
+  for(Eigen::Index f=0; f<m; ++f) {
+    for(int e=0; e<3; ++e) {
+      const ScalarK eij=l_sq(f,e), ejk=l_sq(f,(e+1)%3),
+      eki=l_sq(f,(e+2)%3); //These are squared quantities.
+      const ScalarK lens = sqrt(eij*eki);
+      const ScalarK o = oE(f,e)*oE(f,(e+2)%3);
+      const typename DerivedF::Scalar i=F(f,(e+1)%3), j=F(f,(e+2)%3), k=F(f,e);
+      const ScalarK ki=scaledKappa(i)*theta(f,(e+1)%3),
+      kj=scaledKappa(j)*theta(f,(e+2)%3), kk=scaledKappa(k)*theta(f,e);
+      
+      const ScalarK costhetaidiv = (eij-ejk+eki)/(2.*lens);
+      const ScalarK sinthetaidiv = sqrt( (1.-pow(eij-ejk+eki,2)/
+        (4.*eij*eki)) );
+      
+      const ScalarK Corrijij = (ki+kj+kk);
+      tripletList.emplace_back(E(f,e), E(f,e), Corrijij);
+      tripletList.emplace_back(E(f,e)+nE, E(f,e)+nE, Corrijij);
+      
+      const ScalarK Corrijki = -o*(ki-kj-kk)*costhetaidiv;
+      tripletList.emplace_back(E(f,e), E(f,(e+2)%3), Corrijki);
+      tripletList.emplace_back(E(f,(e+2)%3), E(f,e), Corrijki);
+      tripletList.emplace_back(E(f,e)+nE, E(f,(e+2)%3)+nE, Corrijki);
+      tripletList.emplace_back(E(f,(e+2)%3)+nE, E(f,e)+nE, Corrijki);
+      
+      const ScalarK Corrijkiperp = o*(ki-kj-kk)*sinthetaidiv;
+      tripletList.emplace_back(E(f,e), E(f,(e+2)%3)+nE, Corrijkiperp);
+      tripletList.emplace_back(E(f,(e+2)%3)+nE, E(f,e), Corrijkiperp);
+      tripletList.emplace_back(E(f,e)+nE, E(f,(e+2)%3), -Corrijkiperp);
+      tripletList.emplace_back(E(f,(e+2)%3), E(f,e)+nE, -Corrijkiperp);
+    }
+  }
+  
+  K.resize(2*nE, 2*nE);
+  K.setFromTriplets(tripletList.begin(), tripletList.end());
+}
+
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::cr_vector_curvature_correction_intrinsic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cr_vector_curvature_correction<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/cr_vector_curvature_correction.h

@@ -0,0 +1,115 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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_CR_VECTOR_CURVATURE_CORRECTION_H
+#define IGL_CR_VECTOR_CURVATURE_CORRECTION_H
+
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+
+namespace igl
+{
+  // Computes the vector Crouzeix-Raviart curvature correction
+  //  term of Oded Stein, Alec Jacobson, Max Wardetzky, Eitan
+  //  Grinspun, 2020. "A Smoothness Energy without Boundary Distortion for
+  //  Curved Surfaces", but using the basis functions by Oded Stein,
+  //  Max Wardetzky, Alec Jacobson, Eitan Grinspun, 2020.
+  //  "A Simple Discretization of the Vector Dirichlet Energy"
+  //
+  // Inputs:
+  //  V, F: input mesh
+  //  E: a mapping from each halfedge to each edge, as computed with
+  //     orient_halfedges.
+  //     will be computed if not provided.
+  //  oE: the orientation of each halfedge compared to the orientation of the
+  //      actual edge, as computed with orient_halfedges.
+  //      will be computed if not provided.
+  //
+  // Outputs:
+  //  K: computed curvature correction matrix
+  //  E, oE: these are computed if they are not present, as described above
+
+  template <typename DerivedV, typename DerivedF, typename DerivedE,
+  typename DerivedOE, typename ScalarK>
+  IGL_INLINE void
+  cr_vector_curvature_correction(
+    const Eigen::MatrixBase<DerivedV>& V,
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarK>& K);
+
+  template <typename DerivedV, typename DerivedF, typename DerivedE,
+  typename DerivedOE, typename ScalarK>
+  IGL_INLINE void
+  cr_vector_curvature_correction(
+    const Eigen::MatrixBase<DerivedV>& V,
+    const Eigen::MatrixBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedE>& E,
+    Eigen::PlainObjectBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarK>& K);
+
+
+  // Version that uses intrinsic quantities as input
+  //
+  // Inputs:
+  //  F: input mesh connectivity
+  //  l_sq: squared edge lengths of each halfedge
+  //  theta: the tip angles at each halfedge
+  //  kappa: the Gaussian curvature at each vertex
+  //  E: a mapping from each halfedge to each edge.
+  //  oE: the orientation of each halfedge compared to the orientation of the
+  //      actual edge.
+  //
+  // Outputs:
+  //  K: computed curvature correction matrix
+
+  template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
+  typename DerivedOE, typename ScalarK>
+  IGL_INLINE void
+  cr_vector_curvature_correction_intrinsic(
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarK>& K);
+
+  template <typename DerivedF, typename DerivedL_sq, typename Derivedtheta,
+  typename DerivedE, typename DerivedOE,
+  typename ScalarK>
+  IGL_INLINE void
+  cr_vector_curvature_correction_intrinsic(
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+    const Eigen::MatrixBase<Derivedtheta>& theta,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarK>& K);
+
+  template <typename DerivedF, typename DerivedL_sq, typename Derivedtheta,
+  typename Derivedkappa, typename DerivedE, typename DerivedOE,
+  typename ScalarK>
+  IGL_INLINE void
+  cr_vector_curvature_correction_intrinsic(
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+    const Eigen::MatrixBase<Derivedtheta>& theta,
+    const Eigen::MatrixBase<Derivedkappa>& kappa,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarK>& K);
+}
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cr_vector_curvature_correction.cpp"
+#endif
+
+#endif

include/igl/cr_vector_laplacian.cpp

@@ -0,0 +1,128 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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 "cr_vector_laplacian.h"
+
+#include <vector>
+
+#include "orient_halfedges.h"
+
+#include "doublearea.h"
+#include "squared_edge_lengths.h"
+
+
+template <typename DerivedV, typename DerivedF, typename DerivedE,
+typename DerivedOE, typename ScalarL>
+IGL_INLINE void
+igl::cr_vector_laplacian(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedE>& E,
+  const Eigen::MatrixBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarL>& L)
+{
+  Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic>
+  l_sq;
+  squared_edge_lengths(V, F, l_sq);
+  cr_vector_laplacian_intrinsic(F, l_sq, E, oE, L);
+}
+
+
+template <typename DerivedV, typename DerivedF, typename DerivedE,
+typename DerivedOE, typename ScalarL>
+IGL_INLINE void
+igl::cr_vector_laplacian(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedE>& E,
+  Eigen::PlainObjectBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarL>& L)
+{
+  if(E.rows()!=F.rows() || E.cols()!=F.cols() || oE.rows()!=F.rows() ||
+   oE.cols()!=F.cols()) {
+    orient_halfedges(F, E, oE);
+  }
+
+  const Eigen::PlainObjectBase<DerivedE>& cE = E;
+  const Eigen::PlainObjectBase<DerivedOE>& coE = oE;
+  cr_vector_laplacian(V, F, cE, coE, L);
+}
+
+
+template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
+typename DerivedOE, typename ScalarL>
+IGL_INLINE void
+igl::cr_vector_laplacian_intrinsic(
+  const Eigen::MatrixBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+  const Eigen::MatrixBase<DerivedE>& E,
+  const Eigen::MatrixBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarL>& L)
+{
+  Eigen::Matrix<typename DerivedL_sq::Scalar, Eigen::Dynamic, Eigen::Dynamic>
+  dA;
+  DerivedL_sq l_sqrt = l_sq.array().sqrt().matrix();
+  doublearea(l_sqrt, dA);
+  cr_vector_laplacian_intrinsic(F, l_sq, dA, E, oE, L);
+}
+
+
+template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
+typename DerivedE, typename DerivedOE, typename ScalarL>
+IGL_INLINE void
+igl::cr_vector_laplacian_intrinsic(
+  const Eigen::MatrixBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+  const Eigen::MatrixBase<DeriveddA>& dA,
+  const Eigen::MatrixBase<DerivedE>& E,
+  const Eigen::MatrixBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarL>& L)
+{
+  assert(F.cols()==3 && "Faces have three vertices");
+  assert(E.rows()==F.rows() && E.cols()==F.cols() && oE.rows()==F.rows() &&
+   oE.cols()==F.cols() && "Wrong dimension in edge vectors");
+  assert(l_sq.rows()==F.rows() && l_sq.cols()==3 && "l_sq dimensions wrong");
+  assert(dA.size()==F.rows() && "dA dimensions wrong");
+
+  const Eigen::Index m = F.rows();
+  const typename DerivedE::Scalar nE = E.maxCoeff() + 1;
+
+  std::vector<Eigen::Triplet<ScalarL> > tripletList;
+  tripletList.reserve(10*3*m);
+  for(Eigen::Index f=0; f<m; ++f) {
+    for(int e=0; e<3; ++e) {
+      const ScalarL eij=l_sq(f,e), ejk=l_sq(f,(e+1)%3),
+      eki=l_sq(f,(e+2)%3); //These are squared quantities.
+      const ScalarL lens = sqrt(eij*eki);
+      const ScalarL o = oE(f,e)*oE(f,(e+2)%3);
+      
+      tripletList.emplace_back(E(f,e), E(f,e), 2./dA(f) * eij);
+      tripletList.emplace_back(E(f,e)+nE, E(f,e)+nE, 2./dA(f) * eij);
+      
+      const ScalarL Dijki = o * pow(eij-ejk+eki,2)/(2.*lens*dA(f));
+      tripletList.emplace_back(E(f,e), E(f,(e+2)%3), Dijki);
+      tripletList.emplace_back(E(f,(e+2)%3), E(f,e), Dijki);
+      tripletList.emplace_back(E(f,e)+nE, E(f,(e+2)%3)+nE, Dijki);
+      tripletList.emplace_back(E(f,(e+2)%3)+nE, E(f,e)+nE, Dijki);
+      
+      const ScalarL Dijkiperp = -o * (eij-ejk+eki)/lens;
+      tripletList.emplace_back(E(f,e), E(f,(e+2)%3)+nE, Dijkiperp);
+      tripletList.emplace_back(E(f,(e+2)%3)+nE, E(f,e), Dijkiperp);
+      tripletList.emplace_back(E(f,e)+nE, E(f,(e+2)%3), -Dijkiperp);
+      tripletList.emplace_back(E(f,(e+2)%3), E(f,e)+nE, -Dijkiperp);
+    }
+  }
+  L.resize(2*nE, 2*nE);
+  L.setFromTriplets(tripletList.begin(), tripletList.end());
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::cr_vector_laplacian<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cr_vector_laplacian_intrinsic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/cr_vector_laplacian.h

@@ -0,0 +1,99 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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_CR_VECTOR_LAPLACIAN_H
+#define IGL_CR_VECTOR_LAPLACIAN_H
+
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+
+namespace igl
+{
+  // Computes the CR vector Laplacian matrix.
+  // See Oded Stein, Max Wardetzky, Alec Jacobson, Eitan Grinspun, 2020.
+  //  "A Simple Discretization of the Vector Dirichlet Energy"
+  //
+  // Inputs:
+  //  V, F: input mesh
+  //  E: a mapping from each halfedge to each edge, as computed with
+  //     orient_halfedges.
+  //     will be computed if not provided.
+  //  oE: the orientation of each halfedge compared to the orientation of the
+  //      actual edge, as computed with orient_halfedges.
+  //      will be computed if not provided.
+  //
+  // Outputs:
+  //  L: computed Laplacian matrix
+  //  E, oE: these are computed if they are not present, as described above
+
+  template <typename DerivedV, typename DerivedF, typename DerivedE,
+  typename DerivedOE, typename ScalarL>
+  IGL_INLINE void
+  cr_vector_laplacian(
+    const Eigen::MatrixBase<DerivedV>& V,
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarL>& L);
+
+  template <typename DerivedV, typename DerivedF, typename DerivedE,
+  typename DerivedOE, typename ScalarL>
+  IGL_INLINE void
+  cr_vector_laplacian(
+    const Eigen::MatrixBase<DerivedV>& V,
+    const Eigen::MatrixBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedE>& E,
+    Eigen::PlainObjectBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarL>& L);
+
+
+  // Version that uses intrinsic quantities as input
+  //
+  // Inputs:
+  //  F: input mesh connectivity
+  //  l_sq: squared edge lengths of each halfedge
+  //  dA: double area of each face
+  //  E: a mapping from each halfedge to each edge.
+  //  oE: the orientation of each halfedge compared to the orientation of the
+  //      actual edge.
+  //
+  // Outputs:
+  //  L: computed Laplacian matrix
+
+  template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
+  typename DerivedOE, typename ScalarL>
+  IGL_INLINE void
+  cr_vector_laplacian_intrinsic(
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarL>& L);
+
+  template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
+  typename DerivedE, typename DerivedOE, typename ScalarL>
+  IGL_INLINE void
+  cr_vector_laplacian_intrinsic(
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+    const Eigen::MatrixBase<DeriveddA>& dA,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarL>& L);
+
+
+}
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cr_vector_laplacian.cpp"
+#endif
+
+#endif

include/igl/cr_vector_mass.cpp

@@ -0,0 +1,112 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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 "cr_vector_mass.h"
+
+#include <vector>
+
+#include "orient_halfedges.h"
+
+#include "doublearea.h"
+#include "squared_edge_lengths.h"
+
+
+template <typename DerivedV, typename DerivedF, typename DerivedE,
+typename DerivedOE, typename ScalarM>
+IGL_INLINE void
+igl::cr_vector_mass(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedE>& E,
+  const Eigen::MatrixBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarM>& M)
+{
+  Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic>
+  l_sq;
+  squared_edge_lengths(V, F, l_sq);
+  cr_vector_mass_intrinsic(F, l_sq, E, oE, M);
+}
+
+
+template <typename DerivedV, typename DerivedF, typename DerivedE,
+typename DerivedOE, typename ScalarM>
+IGL_INLINE void
+igl::cr_vector_mass(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedE>& E,
+  Eigen::PlainObjectBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarM>& M)
+{
+  if(E.rows()!=F.rows() || E.cols()!=F.cols() || oE.rows()!=F.rows() ||
+   oE.cols()!=F.cols()) {
+    orient_halfedges(F, E, oE);
+  }
+
+  const Eigen::PlainObjectBase<DerivedE>& cE = E;
+  const Eigen::PlainObjectBase<DerivedOE>& coE = oE;
+  cr_vector_mass(V, F, cE, coE, M);
+}
+
+
+template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
+typename DerivedOE, typename ScalarM>
+IGL_INLINE void
+igl::cr_vector_mass_intrinsic(
+  const Eigen::MatrixBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+  const Eigen::MatrixBase<DerivedE>& E,
+  const Eigen::MatrixBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarM>& M)
+{
+  Eigen::Matrix<typename DerivedL_sq::Scalar, Eigen::Dynamic, Eigen::Dynamic>
+  dA;
+  DerivedL_sq l_sqrt = l_sq.array().sqrt().matrix();
+  doublearea(l_sqrt, dA);
+  cr_vector_mass_intrinsic(F, l_sq, dA, E, oE, M);
+}
+
+
+template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
+typename DerivedE, typename DerivedOE, typename ScalarM>
+IGL_INLINE void
+igl::cr_vector_mass_intrinsic(
+ const Eigen::MatrixBase<DerivedF>& F,
+ const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+ const Eigen::MatrixBase<DeriveddA>& dA,
+ const Eigen::MatrixBase<DerivedE>& E,
+ const Eigen::MatrixBase<DerivedOE>& oE,
+ Eigen::SparseMatrix<ScalarM>& M)
+{
+  assert(F.cols()==3 && "Faces have three vertices");
+  assert(E.rows()==F.rows() && E.cols()==F.cols() && oE.rows()==F.rows() &&
+   oE.cols()==F.cols() && "Wrong dimension in edge vectors");
+
+  const Eigen::Index m = F.rows();
+  const typename DerivedE::Scalar nE = E.maxCoeff() + 1;
+
+  std::vector<Eigen::Triplet<ScalarM> > tripletList;
+  tripletList.reserve(2*3*m);
+  for(Eigen::Index f=0; f<m; ++f) {
+    for(int e=0; e<3; ++e) {
+      const typename DerivedF::Scalar v1=F(f,(e+1)%3), v2=F(f,(e+2)%3);
+      //Scaled
+      const ScalarM entry = dA(f) / 6;
+      tripletList.emplace_back(E(f,e), E(f,e), entry);
+      tripletList.emplace_back(E(f,e)+nE, E(f,e)+nE, entry);
+    }
+  }
+  M.resize(2*nE, 2*nE);
+  M.setFromTriplets(tripletList.begin(), tripletList.end());
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::cr_vector_mass<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cr_vector_mass_intrinsic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/cr_vector_mass.h

@@ -0,0 +1,100 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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_CR_VECTOR_MASS
+#define IGL_CR_VECTOR_MASS
+
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+
+namespace igl
+{
+  // Computes the CR vector mass matrix, using an arrangement of all parallel
+  //  degrees of freedom first, and all perpendicular degrees of freedom next.
+  // See Oded Stein, Max Wardetzky, Alec Jacobson, Eitan Grinspun, 2020.
+  //  "A Simple Discretization of the Vector Dirichlet Energy"
+  //
+  // Inputs:
+  //  V, F: input mesh
+  //  E: a mapping from each halfedge to each edge, as computed with
+  //     orient_halfedges.
+  //     will be computed if not provided.
+  //  oE: the orientation of each halfedge compared to the orientation of the
+  //      actual edge, as computed with orient_halfedges.
+  //      will be computed if not provided.
+  //
+  // Outputs:
+  //  M: computed mass matrix
+  //  E, oE: these are computed if they are not present, as described above
+
+  template <typename DerivedV, typename DerivedF, typename DerivedE,
+  typename DerivedOE, typename ScalarM>
+  IGL_INLINE void
+  cr_vector_mass(
+    const Eigen::MatrixBase<DerivedV>& V,
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarM>& M);
+
+  template <typename DerivedV, typename DerivedF, typename DerivedE,
+  typename DerivedOE, typename ScalarM>
+  IGL_INLINE void
+  cr_vector_mass(
+    const Eigen::MatrixBase<DerivedV>& V,
+    const Eigen::MatrixBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedE>& E,
+    Eigen::PlainObjectBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarM>& M);
+
+
+  // Version that uses intrinsic quantities as input
+  //
+  // Inputs:
+  //  F: input mesh connectivity
+  //  l_sq: squared edge lengths of each halfedge
+  //  dA: double area of each face
+  //  E: a mapping from each halfedge to each edge.
+  //  oE: the orientation of each halfedge compared to the orientation of the
+  //      actual edge.
+  //
+  // Outputs:
+  //  M: computed mass matrix
+
+  template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
+  typename DerivedOE, typename ScalarM>
+  IGL_INLINE void
+  cr_vector_mass_intrinsic(
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarM>& M);
+
+  template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
+  typename DerivedE, typename DerivedOE, typename ScalarM>
+  IGL_INLINE void
+  cr_vector_mass_intrinsic(
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+    const Eigen::MatrixBase<DeriveddA>& dA,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarM>& M);
+
+
+}
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cr_vector_mass.cpp"
+#endif
+
+#endif

include/igl/crouzeix_raviart_cotmatrix.cpp

@@ -97,4 +97,5 @@ void igl::crouzeix_raviart_cotmatrix(
 // Explicit template instantiation
 // generated by autoexplicit.sh
 template void igl::crouzeix_raviart_cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::crouzeix_raviart_cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
 #endif

include/igl/crouzeix_raviart_massmatrix.cpp

@@ -82,4 +82,5 @@ void igl::crouzeix_raviart_massmatrix(
 template void igl::crouzeix_raviart_massmatrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::crouzeix_raviart_massmatrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&);
 template void igl::crouzeix_raviart_massmatrix<float, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<float, 0, int>&);
+template void igl::crouzeix_raviart_massmatrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
 #endif

include/igl/curved_hessian_energy.cpp

@@ -0,0 +1,127 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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 "curved_hessian_energy.h"
+
+#include "orient_halfedges.h"
+#include "doublearea.h"
+#include "squared_edge_lengths.h"
+#include "cr_vector_laplacian.h"
+#include "cr_vector_mass.h"
+#include "cr_vector_curvature_correction.h"
+#include "scalar_to_cr_vector_gradient.h"
+
+
+template <typename DerivedV, typename DerivedF, typename ScalarQ>
+IGL_INLINE void
+igl::curved_hessian_energy(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedF>& F,
+  Eigen::SparseMatrix<ScalarQ>& Q)
+{
+  Eigen::MatrixXi E, oE;
+  curved_hessian_energy(V, F, E, oE, Q);
+}
+
+
+template <typename DerivedV, typename DerivedF, typename DerivedE,
+typename DerivedOE, typename ScalarQ>
+IGL_INLINE void
+igl::curved_hessian_energy(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedE>& E,
+  const Eigen::MatrixBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarQ>& Q)
+{
+  Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic>
+  l_sq;
+  squared_edge_lengths(V, F, l_sq);
+  curved_hessian_energy_intrinsic(F, l_sq, E, oE, Q);
+}
+
+
+template <typename DerivedV, typename DerivedF, typename DerivedE,
+typename DerivedOE, typename ScalarQ>
+IGL_INLINE void
+igl::curved_hessian_energy(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedE>& E,
+  Eigen::PlainObjectBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarQ>& Q)
+{
+  if(E.rows()!=F.rows() || E.cols()!=F.cols() || oE.rows()!=F.rows() ||
+   oE.cols()!=F.cols()) {
+    orient_halfedges(F, E, oE);
+  }
+
+  const Eigen::PlainObjectBase<DerivedE>& cE = E;
+  const Eigen::PlainObjectBase<DerivedOE>& coE = oE;
+  curved_hessian_energy(V, F, cE, coE, Q);
+}
+
+
+template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
+typename DerivedOE, typename ScalarQ>
+IGL_INLINE void
+igl::curved_hessian_energy_intrinsic(
+  const Eigen::MatrixBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+  const Eigen::MatrixBase<DerivedE>& E,
+  const Eigen::MatrixBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarQ>& Q)
+{
+  Eigen::Matrix<typename DerivedL_sq::Scalar, Eigen::Dynamic, Eigen::Dynamic>
+  dA;
+  Eigen::Matrix<typename DerivedL_sq::Scalar, Eigen::Dynamic, Eigen::Dynamic>
+  l_sqrt = l_sq.array().sqrt().matrix();
+  doublearea(l_sqrt, dA);
+  curved_hessian_energy_intrinsic(F, l_sq, dA, E, oE, Q);
+}
+
+
+template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
+typename DerivedE, typename DerivedOE, typename ScalarQ>
+IGL_INLINE void
+igl::curved_hessian_energy_intrinsic(
+  const Eigen::MatrixBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+  const Eigen::MatrixBase<DeriveddA>& dA,
+  const Eigen::MatrixBase<DerivedE>& E,
+  const Eigen::MatrixBase<DerivedOE>& oE,
+  Eigen::SparseMatrix<ScalarQ>& Q)
+{
+  //Matrices that need to be combined
+  Eigen::SparseMatrix<ScalarQ> M, D, L, K;
+  cr_vector_mass_intrinsic(F, l_sq, dA,  E, oE, M);
+  scalar_to_cr_vector_gradient_intrinsic(F, l_sq, dA, E, oE, D);
+  cr_vector_laplacian_intrinsic(F, l_sq, dA, E, oE, L);
+  cr_vector_curvature_correction_intrinsic(F, l_sq, E, oE, K);
+
+  //Invert M
+  std::vector<Eigen::Triplet<ScalarQ> > tripletListMi;
+  for(Eigen::Index k=0; k<M.outerSize(); ++k) {
+    for(typename Eigen::SparseMatrix<ScalarQ>::InnerIterator it(M,k);
+      it; ++it) {
+      if(it.value() > 0) {
+        tripletListMi.emplace_back(it.row(), it.col(), 1./it.value());
+      }
+    }
+  }
+  Eigen::SparseMatrix<ScalarQ> Mi(M.rows(), M.cols());
+  Mi.setFromTriplets(tripletListMi.begin(), tripletListMi.end());
+
+  //Hessian energy matrix
+  Q = D.transpose()*Mi*(L + K)*Mi*D;
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::curved_hessian_energy<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/curved_hessian_energy.h

@@ -0,0 +1,114 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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_CURVED_HESSIAN_ENERGY_H
+#define IGL_CURVED_HESSIAN_ENERGY_H
+
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+
+namespace igl
+{
+  // Computes the curved Hessian energy using the Crouzeix-Raviart
+  //  discretization.
+  // See Oded Stein, Alec Jacobson, Max Wardetzky, Eitan Grinspun, 2020.
+  //  "A Smoothness Energy without Boundary Distortion for Curved Surfaces"
+  //
+  // Inputs:
+  //  V, F: input mesh
+  //
+  // Outputs:
+  //  Q: computed Hessian energy matrix
+
+  template <typename DerivedV, typename DerivedF, typename ScalarQ>
+  IGL_INLINE void
+  curved_hessian_energy(
+    const Eigen::MatrixBase<DerivedV>& V,
+    const Eigen::MatrixBase<DerivedF>& F,
+    Eigen::SparseMatrix<ScalarQ>& Q);
+
+  // Version that exposes the edge orientation used.
+  //
+  // Inputs:
+  //  V, F: input mesh
+  //  E: a mapping from each halfedge to each edge, as computed with
+  //     orient_halfedges.
+  //     will be computed if not provided.
+  //  oE: the orientation of each halfedge compared to the orientation of the
+  //      actual edge, as computed with orient_halfedges.
+  //      will be computed if not provided.
+  //
+  // Outputs:
+  //  Q: computed Hessian energy matrix
+  //  E, oE: these are computed if they are not present, as described above
+  template <typename DerivedV, typename DerivedF, typename DerivedE,
+  typename DerivedOE, typename ScalarQ>
+  IGL_INLINE void
+  curved_hessian_energy(
+    const Eigen::MatrixBase<DerivedV>& V,
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarQ>& Q);
+
+  template <typename DerivedV, typename DerivedF, typename DerivedE,
+  typename DerivedOE, typename ScalarQ>
+  IGL_INLINE void
+  curved_hessian_energy(
+    const Eigen::MatrixBase<DerivedV>& V,
+    const Eigen::MatrixBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedE>& E,
+    Eigen::PlainObjectBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarQ>& Q);
+
+
+  // Version that uses intrinsic quantities as input
+  //
+  // Inputs:
+  //  F: input mesh connectivity
+  //  l_sq: squared edge lengths of each halfedge
+  //  dA: double area of each face
+  //  E: a mapping from each halfedge to each edge.
+  //  oE: the orientation of each halfedge compared to the orientation of the
+  //      actual edge.
+  //
+  // Outputs:
+  //  Q: computed Hessian energy matrix
+
+  template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
+  typename DerivedOE, typename ScalarQ>
+  IGL_INLINE void
+  curved_hessian_energy_intrinsic(
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarQ>& Q);
+
+  template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
+  typename DerivedE, typename DerivedOE, typename ScalarQ>
+  IGL_INLINE void
+  curved_hessian_energy_intrinsic(
+    const Eigen::MatrixBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+    const Eigen::MatrixBase<DeriveddA>& dA,
+    const Eigen::MatrixBase<DerivedE>& E,
+    const Eigen::MatrixBase<DerivedOE>& oE,
+    Eigen::SparseMatrix<ScalarQ>& Q);
+
+
+}
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "curved_hessian_energy.cpp"
+#endif
+
+#endif

include/igl/edge_midpoints.cpp

@@ -0,0 +1,48 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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 "edge_midpoints.h"
+
+template<typename DerivedV,typename DerivedF,typename DerivedE,
+typename DerivedoE, typename Derivedmps>
+IGL_INLINE void
+igl::edge_midpoints(
+	const Eigen::MatrixBase<DerivedV> &V,
+	const Eigen::MatrixBase<DerivedF> &F,
+	const Eigen::MatrixBase<DerivedE> &E,
+	const Eigen::MatrixBase<DerivedoE> &oE,
+	Eigen::PlainObjectBase<Derivedmps> &mps)
+{
+  assert(E.rows()==F.rows() && "E does not match dimensions of F.");
+  assert(oE.rows()==F.rows() && "oE does not match dimensions of F.");
+  assert(E.cols()==3 && F.cols()==3 && oE.cols()==3 &&
+    "This method is for triangle meshes.");
+  assert(F.maxCoeff()<V.rows() && "V does not seem to belong to F.");
+
+  using ScalarE = typename DerivedE::Scalar;
+  using ScalarF = typename DerivedF::Scalar;
+	
+  const ScalarE m = E.maxCoeff()+1;
+	
+  mps.resize(m, V.cols());
+  for(Eigen::Index i=0; i<F.rows(); ++i) {
+    for(int j=0; j<3; ++j) {
+      if(oE(i,j)<0) {
+        continue;
+      }
+      const ScalarE e = E(i,j);
+      const ScalarF vi=F(i,(j+1)%3), vj=F(i,(j+2)%3);
+
+      mps.row(e) = 0.5*(V.row(vi) + V.row(vj));
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::edge_midpoints<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/edge_midpoints.h

@@ -0,0 +1,38 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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_EDGE_MIDPOINTS_H
+#define IGL_EDGE_MIDPOINTS_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+namespace igl
+{
+  // Computes the midpoints of edges in a triangle mesh.
+  //
+  // Input:
+  // V, F: triangle mesh
+  // E, oE: mapping from halfedges to edges and orientation as generated by
+  //    orient_halfedges
+  //
+  // Output:
+  // mps: edge midpoints, one per edge in E
+  template<typename DerivedV,typename DerivedF,typename DerivedE,
+  typename DerivedoE, typename Derivedmps>
+  IGL_INLINE void edge_midpoints(
+    const Eigen::MatrixBase<DerivedV> &V,
+    const Eigen::MatrixBase<DerivedF> &F,
+    const Eigen::MatrixBase<DerivedE> &E,
+    const Eigen::MatrixBase<DerivedoE> &oE,
+    Eigen::PlainObjectBase<Derivedmps> &mps);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "edge_midpoints.cpp"
+#endif
+
+#endif

include/igl/edge_vectors.cpp

@@ -0,0 +1,94 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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 "edge_vectors.h"
+
+#include <Eigen/Geometry>
+#include "per_face_normals.h"
+
+#include "PI.h"
+
+
+template<typename DerivedV,typename DerivedF,typename DerivedE,
+typename DerivedoE, typename Derivedvec>
+IGL_INLINE void
+igl::edge_vectors(
+                  const Eigen::MatrixBase<DerivedV> &V,
+                  const Eigen::MatrixBase<DerivedF> &F,
+                  const Eigen::MatrixBase<DerivedE> &E,
+                  const Eigen::MatrixBase<DerivedoE> &oE,
+                  Eigen::PlainObjectBase<Derivedvec> &vec)
+{
+  Eigen::Matrix<typename Derivedvec::Scalar, Eigen::Dynamic, Eigen::Dynamic>
+  dummy;
+  edge_vectors<false>(V, F, E, oE, vec, dummy);
+}
+
+
+template<bool computePerpendicular,
+typename DerivedV,typename DerivedF,typename DerivedE,
+typename DerivedoE, typename DerivedvecParallel,
+typename DerivedvecPerpendicular>
+IGL_INLINE void
+igl::edge_vectors(
+                  const Eigen::MatrixBase<DerivedV> &V,
+                  const Eigen::MatrixBase<DerivedF> &F,
+                  const Eigen::MatrixBase<DerivedE> &E,
+                  const Eigen::MatrixBase<DerivedoE> &oE,
+                  Eigen::PlainObjectBase<DerivedvecParallel> &vecParallel,
+                  Eigen::PlainObjectBase<DerivedvecPerpendicular> &vecPerpendicular)
+{
+  using Scalar = typename DerivedvecParallel::Scalar;
+  using MatX = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
+  
+  assert(E.rows()==F.rows() && "E does not match dimensions of F.");
+  assert(oE.rows()==F.rows() && "oE does not match dimensions of F.");
+  assert(E.cols()==3 && F.cols()==3 && oE.cols()==3 &&
+         "This method is for triangle meshes.");
+  assert(F.maxCoeff()<V.rows() && "V does not seem to belong to F.");
+  
+  const typename DerivedE::Scalar m = E.maxCoeff()+1;
+  
+  //Compute edge-based normal
+  MatX N, edgeN(m, 3);
+  edgeN.setZero();
+  per_face_normals(V, F, N);
+  for(Eigen::Index i=0; i<E.rows(); ++i) {
+    for(int j=0; j<3; ++j) {
+      edgeN.row(E(i,j)) += N.row(i);
+    }
+  }
+  edgeN.rowwise().normalize();
+  
+  //Compute edge vectors
+  vecParallel.resize(m, 3);
+  if(computePerpendicular) { //This should ideally be an if constexpr
+    vecPerpendicular.resize(m, 3);
+  }
+  for(Eigen::Index i=0; i<E.rows(); ++i) {
+    for(int j=0; j<3; ++j) {
+      if(oE(i,j)<0) {
+        continue;
+      }
+      const typename DerivedE::Scalar e=E(i,j);
+      const typename DerivedF::Scalar vi=F(i,(j+1)%3), vj=F(i,(j+2)%3);
+      vecParallel.row(e) = (V.row(vj)-V.row(vi)).normalized();
+      if(computePerpendicular) { //This should ideally be an if constexpr
+        vecPerpendicular.row(e) =
+        Eigen::AngleAxis<Scalar>(0.5*PI, edgeN.row(e)) *
+        vecParallel.row(e).transpose();
+      }
+    }
+  }
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::edge_vectors<true, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/edge_vectors.h

@@ -0,0 +1,68 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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_EDGE_VECTORS_H
+#define IGL_EDGE_VECTORS_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+namespace igl
+{
+// Computes the normalized edge vectors for edges in a triangle mesh
+//
+// Input:
+// V, F: triangle mesh
+// E, oE: mapping from halfedges to edges and orientation as generated by
+//    orient_halfedges
+//
+// Output:
+// vec: normalized edge vectors for each unique edge in E, according to order
+//      of E.
+template<typename DerivedV,typename DerivedF,typename DerivedE,
+typename DerivedoE, typename Derivedvec>
+IGL_INLINE void edge_vectors(
+                             const Eigen::MatrixBase<DerivedV> &V,
+                             const Eigen::MatrixBase<DerivedF> &F,
+                             const Eigen::MatrixBase<DerivedE> &E,
+                             const Eigen::MatrixBase<DerivedoE> &oE,
+                             Eigen::PlainObjectBase<Derivedvec> &vec);
+
+
+// Computes the normalized edge vectors for edges in a triangle mesh
+//
+// Input:
+// V, F: triangle mesh
+// E, oE: mapping from halfedges to edges and orientation as generated by
+//    orient_halfedges
+// template parameter computePerpendicular: whether to compute
+//                                          vecPerpendicular or not.
+//
+// Output:
+// vecParallel: normalized edge vectors for each unique edge in E, according
+//              to order of E.
+// vecPerpendicular: tangent unit perpendicular vector to each edge, according
+//                   to orientation in oE, corresponds to each vector in
+//                   vecParallel rotated by pi/2 around an edge-based normal.
+//
+template<bool computePerpendicular=true,
+typename DerivedV,typename DerivedF,typename DerivedE,
+typename DerivedoE, typename DerivedvecParallel,
+typename DerivedvecPerpendicular>
+IGL_INLINE void edge_vectors(
+                             const Eigen::MatrixBase<DerivedV> &V,
+                             const Eigen::MatrixBase<DerivedF> &F,
+                             const Eigen::MatrixBase<DerivedE> &E,
+                             const Eigen::MatrixBase<DerivedoE> &oE,
+                             Eigen::PlainObjectBase<DerivedvecParallel> &vecParallel,
+                             Eigen::PlainObjectBase<DerivedvecPerpendicular> &vecPerpendicular);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "edge_vectors.cpp"
+#endif
+
+#endif

include/igl/internal_angles.cpp

@@ -101,4 +101,5 @@ template void igl::internal_angles<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eig
 template void igl::internal_angles<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::internal_angles<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::internal_angles_using_squared_edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::internal_angles_using_squared_edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/min_quad_with_fixed.cpp

@@ -600,4 +600,5 @@ template bool igl::min_quad_with_fixed_solve<double, Eigen::Matrix<double, -1, 1
 template bool igl::min_quad_with_fixed_solve<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template bool igl::min_quad_with_fixed_solve<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template bool igl::min_quad_with_fixed<double, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::min_quad_with_fixed<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/opengl/glfw/imgui/ImGuizmo.h

@@ -1,59 +0,0 @@
-#pragma once
-////////////////////////////////////////////////////////////////////////////////
-#include <imgui/imgui.h>
-#include <imgui/imgui_internal.h>
-#include <imguizmo/ImGuizmo.h>
-#include <Eigen/Dense>
-#include <vector>
-////////////////////////////////////////////////////////////////////////////////
-
-namespace igl
-{
-namespace opengl
-{
-namespace glfw
-{
-namespace imgui
-{
-
-void EditTransform(const float *cameraView, const float *cameraProjection, float* matrix, bool isOrthographic)
-{
-	static ImGuizmo::OPERATION mCurrentGizmoOperation(ImGuizmo::ROTATE);
-        // Alec: This could be a parameter ImGuizmo::WORLD will make senze in some
-        // circumstances. For now, 109_ImGuizmo and my other use cases make more
-        // sense in LOCAL mode.
-	static ImGuizmo::MODE mCurrentGizmoMode(ImGuizmo::LOCAL);
-
-        // Maya-style key bindings
-	if (ImGui::IsKeyPressed('W'))
-		mCurrentGizmoOperation = ImGuizmo::TRANSLATE;
-	if (ImGui::IsKeyPressed('E'))
-		mCurrentGizmoOperation = ImGuizmo::ROTATE;
-	if (ImGui::IsKeyPressed('R'))
-		mCurrentGizmoOperation = ImGuizmo::SCALE;
-
-	if (ImGui::RadioButton("Translate", mCurrentGizmoOperation == ImGuizmo::TRANSLATE))
-		mCurrentGizmoOperation = ImGuizmo::TRANSLATE;
-	ImGui::SameLine();
-	if (ImGui::RadioButton("Rotate", mCurrentGizmoOperation == ImGuizmo::ROTATE))
-		mCurrentGizmoOperation = ImGuizmo::ROTATE;
-	ImGui::SameLine();
-	if (ImGui::RadioButton("Scale", mCurrentGizmoOperation == ImGuizmo::SCALE))
-		mCurrentGizmoOperation = ImGuizmo::SCALE;
-
-	float matrixTranslation[3], matrixRotation[3], matrixScale[3];
-	ImGuizmo::DecomposeMatrixToComponents(matrix, matrixTranslation, matrixRotation, matrixScale);
-	ImGui::InputFloat3("Tr", matrixTranslation, 3);
-	ImGui::InputFloat3("Rt", matrixRotation, 3);
-	ImGui::InputFloat3("Sc", matrixScale, 3);
-	ImGuizmo::RecomposeMatrixFromComponents(matrixTranslation, matrixRotation, matrixScale, matrix);
-
-	ImGuiIO& io = ImGui::GetIO();
-	ImGuizmo::SetRect(0, 0, io.DisplaySize.x, io.DisplaySize.y);
-	ImGuizmo::Manipulate(cameraView, cameraProjection, mCurrentGizmoOperation, mCurrentGizmoMode, matrix, NULL, NULL);
-}
-
-} // end namespace
-} // end namespace
-} // end namespace
-} // end namespace

include/igl/opengl/glfw/imgui/ImGuizmoPlugin.cpp

@@ -0,0 +1,49 @@
+#include "ImGuizmoPlugin.h"
+#include <imgui/imgui.h>
+#include <imgui_impl_glfw.h>
+#include <imgui_impl_opengl3.h>
+#include <imgui_fonts_droid_sans.h>
+#include <GLFW/glfw3.h>
+
+namespace igl{ namespace opengl{ namespace glfw{ namespace imgui{
+
+IGL_INLINE void ImGuizmoPlugin::init(igl::opengl::glfw::Viewer *_viewer)
+{
+  ImGuiMenu::init(_viewer);
+}
+IGL_INLINE bool ImGuizmoPlugin::pre_draw() 
+{
+  ImGuiMenu::pre_draw();
+  ImGui::PushStyleVar(ImGuiStyleVar_WindowBorderSize, 0);
+  ImGuizmo::BeginFrame();
+  ImGui::PopStyleVar();
+  return false;
+}
+IGL_INLINE bool ImGuizmoPlugin::post_draw() 
+{
+  // Don't draw the Viewer's default menu: draw just the ImGuizmo
+  Eigen::Matrix4f view = (viewer->core().view / viewer->core().camera_zoom);
+  Eigen::Matrix4f proj = viewer->core().proj;
+  if(viewer->core().orthographic){ view(2,3) -= 1000;/* Hack depth */ }
+  // ImGuizmo doesn't like a lot of scaling in view, shift it to T
+  const float z = viewer->core().camera_base_zoom;
+  const Eigen::Matrix4f S = 
+    (Eigen::Matrix4f()<< z,0,0,0, 0,z,0,0, 0,0,z,0, 0,0,0,1).finished();
+  view = (view * S.inverse()).eval();
+  const Eigen::Matrix4f T0 = T;
+  T = (S * T).eval();
+  ImGuiIO& io = ImGui::GetIO();
+  ImGuizmo::SetRect(0, 0, io.DisplaySize.x, io.DisplaySize.y);
+  ImGuizmo::Manipulate(
+    view.data(),proj.data(),operation,ImGuizmo::LOCAL,T.data(),NULL,NULL);
+  // invert scaling that was shifted onto T
+  T = (S.inverse() * T).eval();
+  const float diff = (T-T0).array().abs().maxCoeff();
+  // Only call if actually changed; otherwise, triggers on all mouse events
+  if( diff > 1e-7) { callback(T); }
+  ImGui::Render();
+  ImGui_ImplOpenGL3_RenderDrawData(ImGui::GetDrawData());
+  return false;
+}
+
+}}}}

include/igl/opengl/glfw/imgui/ImGuizmoPlugin.h

@@ -0,0 +1,37 @@
+#ifndef IGL_OPENGL_GFLW_IMGUI_IMGUIZMOPLUGIN_H
+#define IGL_OPENGL_GFLW_IMGUI_IMGUIZMOPLUGIN_H
+#include "../../../igl_inline.h"
+#include "ImGuiMenu.h"
+#include <imgui/imgui.h>
+#include <imgui/imgui_internal.h>
+#include <imguizmo/ImGuizmo.h>
+#include <Eigen/Dense>
+
+namespace igl{ namespace opengl{ namespace glfw{ namespace imgui{
+
+class ImGuizmoPlugin : public igl::opengl::glfw::imgui::ImGuiMenu
+{
+public:
+  // callback(T) called when the stored transform T changes
+  std::function<void(const Eigen::Matrix4f &)> callback;
+  // whether rotating, translating or scaling
+  ImGuizmo::OPERATION operation;
+  // stored transformation
+  Eigen::Matrix4f T;
+  // Initilize with rotate operation on an identity transform (at origin)
+  ImGuizmoPlugin():operation(ImGuizmo::ROTATE),T(Eigen::Matrix4f::Identity()){};
+  /////////////////////////////////////////////////////////////////////////////
+  // Boilerplate
+  virtual void init(igl::opengl::glfw::Viewer *_viewer) override;
+  virtual bool pre_draw() override;
+  /////////////////////////////////////////////////////////////////////////////
+  virtual bool post_draw() override;
+};
+
+}}}}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "ImGuizmoPlugin.cpp"
+#endif
+
+#endif

include/igl/orient_halfedges.cpp

@@ -0,0 +1,48 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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 "orient_halfedges.h"
+
+#include "oriented_facets.h"
+#include "unique_simplices.h"
+
+
+template <typename DerivedF, typename DerivedE, typename DerivedOE>
+IGL_INLINE void
+igl::orient_halfedges(
+  const Eigen::MatrixBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedE>& E,
+  Eigen::PlainObjectBase<DerivedOE>& oE)
+{
+  assert(F.cols()==3 && "This only works for triangle meshes.");
+  
+  using Int = typename DerivedF::Scalar;
+  
+  const Eigen::Index m = F.rows();
+  
+  DerivedE allE, EE;
+  oriented_facets(F, allE);
+  Eigen::Matrix<Int, Eigen::Dynamic, 1> IA, IC;
+  unique_simplices(allE, EE, IA, IC);
+  
+  E.resize(m, 3);
+  oE.resize(m, 3);
+  for(Eigen::Index f=0; f<m; ++f) {
+    for(int e=0; e<3; ++e) {
+      const Int ind = f + m*e;
+      E(f,e) = IC(ind);
+      assert((EE(E(f,e),0)==allE(ind,0) || EE(E(f,e),0)==allE(ind,1)) &&
+       "Something is wrong in the edge matrix.");
+      oE(f,e) = EE(E(f,e),0)==allE(ind,0) ? 1 : -1;
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::orient_halfedges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/orient_halfedges.h

@@ -0,0 +1,41 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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_ORIENT_HALFEDGES_H
+#define IGL_ORIENT_HALFEDGES_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+
+namespace igl
+{
+  // Orients halfedges for a triangle mesh, assigning them to a unique edge.
+  //
+  // Inputs:
+  //  F: input mesh connectivity
+  //
+  // Outputs:
+  //  E: a mapping from each halfedge to each edge
+  //  oE: the orientation of each halfedge compared to the orientation of the
+  //      actual edge. Every edge appears positively oriented exactly once.
+
+  template <typename DerivedF, typename DerivedE, typename DerivedOE>
+  IGL_INLINE void
+  orient_halfedges(
+    const Eigen::MatrixBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedE>& E,
+    Eigen::PlainObjectBase<DerivedOE>& oE);
+
+}
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "orient_halfedges.cpp"
+#endif
+
+#endif

include/igl/random_points_on_mesh.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2014 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 
+//
+// 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 "random_points_on_mesh.h"
 #include "doublearea.h"
@@ -42,6 +42,7 @@ IGL_INLINE void igl::random_points_on_mesh(
   assert(R.minCoeff() >= 0);
   assert(R.maxCoeff() <= 1);
   histc(R,C,FI);
+  FI = FI.array().min(F.rows() - 1); // fix the bin when R(i) == 1 exactly
   const VectorXs S = (VectorXs::Random(n,1).array() + 1.)/2.;
   const VectorXs T = (VectorXs::Random(n,1).array() + 1.)/2.;
   B.resize(n,3);
@@ -51,9 +52,9 @@ IGL_INLINE void igl::random_points_on_mesh(
 }
 
 template <
-  typename DerivedV, 
-  typename DerivedF, 
-  typename DerivedB, 
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedB,
   typename DerivedFI,
   typename DerivedX>
 IGL_INLINE void igl::random_points_on_mesh(
@@ -70,7 +71,10 @@ IGL_INLINE void igl::random_points_on_mesh(
   {
     for(int b = 0;b<B.cols();b++)
     {
-      X.row(x) += B(x,b)*V.row(F(FI(x),b));
+        auto fi = FI(x);
+        auto f = F(fi, b);
+        auto bval = B(x, b);
+        X.row(x) += bval*V.row(f);
     }
   }
 }

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/scalar_to_cr_vector_gradient.cpp

@@ -0,0 +1,116 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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 "scalar_to_cr_vector_gradient.h"
+
+#include "orient_halfedges.h"
+
+#include "doublearea.h"
+#include "squared_edge_lengths.h"
+
+
+template <typename DerivedV, typename DerivedF, typename DerivedE,
+typename DerivedOE, typename ScalarG>
+IGL_INLINE void
+igl::scalar_to_cr_vector_gradient(
+                                  const Eigen::MatrixBase<DerivedV>& V,
+                                  const Eigen::MatrixBase<DerivedF>& F,
+                                  const Eigen::MatrixBase<DerivedE>& E,
+                                  const Eigen::MatrixBase<DerivedOE>& oE,
+                                  Eigen::SparseMatrix<ScalarG>& G)
+{
+  Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic>
+  l_sq;
+  squared_edge_lengths(V, F, l_sq);
+  scalar_to_cr_vector_gradient_intrinsic(F, l_sq, E, oE, G);
+}
+
+template <typename DerivedV, typename DerivedF, typename DerivedE,
+typename DerivedOE, typename ScalarG>
+IGL_INLINE void
+igl::scalar_to_cr_vector_gradient(
+                                  const Eigen::MatrixBase<DerivedV>& V,
+                                  const Eigen::MatrixBase<DerivedF>& F,
+                                  Eigen::PlainObjectBase<DerivedE>& E,
+                                  Eigen::PlainObjectBase<DerivedOE>& oE,
+                                  Eigen::SparseMatrix<ScalarG>& G)
+{
+  if(E.rows()!=F.rows() || E.cols()!=F.cols() || oE.rows()!=F.rows() ||
+     oE.cols()!=F.cols()) {
+    orient_halfedges(F, E, oE);
+  }
+  
+  const Eigen::PlainObjectBase<DerivedE>& cE = E;
+  const Eigen::PlainObjectBase<DerivedOE>& coE = oE;
+
+  scalar_to_cr_vector_gradient(V, F, cE, coE, G);
+}
+
+
+template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
+typename DerivedOE, typename ScalarG>
+IGL_INLINE void
+igl::scalar_to_cr_vector_gradient_intrinsic(
+                                            const Eigen::MatrixBase<DerivedF>& F,
+                                            const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+                                            const Eigen::MatrixBase<DerivedE>& E,
+                                            const Eigen::MatrixBase<DerivedOE>& oE,
+                                            Eigen::SparseMatrix<ScalarG>& G)
+{
+  Eigen::Matrix<typename DerivedL_sq::Scalar, Eigen::Dynamic, Eigen::Dynamic>
+  dA;
+  DerivedL_sq l_sqrt = l_sq.array().sqrt().matrix();
+  doublearea(l_sqrt, dA);
+  scalar_to_cr_vector_gradient_intrinsic(F, l_sq, dA, E, oE, G);
+}
+
+
+template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
+typename DerivedE, typename DerivedOE, typename ScalarG>
+IGL_INLINE void
+igl::scalar_to_cr_vector_gradient_intrinsic(
+                                            const Eigen::MatrixBase<DerivedF>& F,
+                                            const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+                                            const Eigen::MatrixBase<DeriveddA>& dA,
+                                            const Eigen::MatrixBase<DerivedE>& E,
+                                            const Eigen::MatrixBase<DerivedOE>& oE,
+                                            Eigen::SparseMatrix<ScalarG>& G)
+{
+  assert(F.cols()==3 && "Faces have three vertices");
+  assert(E.rows()==F.rows() && E.cols()==F.cols() && oE.rows()==F.rows() &&
+         oE.cols()==F.cols() && "Wrong dimension in edge vectors");
+  
+  const Eigen::Index m = F.rows();
+  const typename DerivedF::Scalar n = F.maxCoeff() + 1;
+  const typename DerivedE::Scalar nE = E.maxCoeff() + 1;
+  
+  std::vector<Eigen::Triplet<ScalarG> > tripletList;
+  tripletList.reserve(5*3*m);
+  for(Eigen::Index f=0; f<m; ++f) {
+    for(int e=0; e<3; ++e) {
+      const typename DerivedF::Scalar i=F(f,(e+1)%3), j=F(f,(e+2)%3), k=F(f,e);
+      const ScalarG o=oE(f,e),
+      eij=l_sq(f,e), ejk=l_sq(f,(e+1)%3), eki=l_sq(f,(e+2)%3); //These are squared quantities.
+      const ScalarG s_eij = sqrt(eij);
+      
+      tripletList.emplace_back(E(f,e), i, -o*dA(f)/(6.*s_eij));
+      tripletList.emplace_back(E(f,e)+nE, i, -o*(eij+ejk-eki)/(12.*s_eij));
+      tripletList.emplace_back(E(f,e), j, o*dA(f)/(6.*s_eij));
+      tripletList.emplace_back(E(f,e)+nE, j, -o*(eij-ejk+eki)/(12.*s_eij));
+      tripletList.emplace_back(E(f,e)+nE, k, o*s_eij/6.);
+    }
+  }
+  G.resize(2*nE, n);
+  G.setFromTriplets(tripletList.begin(), tripletList.end());
+}
+
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::scalar_to_cr_vector_gradient_intrinsic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/scalar_to_cr_vector_gradient.h

@@ -0,0 +1,100 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Oded Stein <[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_SCALAR_TO_CR_VECTOT_GRADIENT_H
+#define IGL_SCALAR_TO_CR_VECTOT_GRADIENT_H
+
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+
+namespace igl
+{
+  // Computes the gradient matrix with hat functions on the right, and
+  //  vector CR functions on the left.
+  // See Oded Stein, Max Wardetzky, Alec Jacobson, Eitan Grinspun, 2020.
+  //  "A Simple Discretization of the Vector Dirichlet Energy"
+  //
+  // Inputs:
+  //  V, F: input mesh
+  //  E: a mapping from each halfedge to each edge, as computed with
+  //     orient_halfedges.
+  //     will be computed if not provided.
+  //  oE: the orientation of each halfedge compared to the orientation of the
+  //      actual edge, as computed with orient_halfedges.
+  //      will be computed if not provided.
+  //
+  // Outputs:
+  //  G: computed gradient matrix
+  //  E, oE: these are computed if they are not present, as described above
+
+  template <typename DerivedV, typename DerivedF, typename DerivedE,
+  typename DerivedOE, typename ScalarG>
+  IGL_INLINE void
+  scalar_to_cr_vector_gradient(
+                               const Eigen::MatrixBase<DerivedV>& V,
+                               const Eigen::MatrixBase<DerivedF>& F,
+                               const Eigen::MatrixBase<DerivedE>& E,
+                               const Eigen::MatrixBase<DerivedOE>& oE,
+                               Eigen::SparseMatrix<ScalarG>& G);
+
+  template <typename DerivedV, typename DerivedF, typename DerivedE,
+  typename DerivedOE, typename ScalarG>
+  IGL_INLINE void
+  scalar_to_cr_vector_gradient(
+                               const Eigen::MatrixBase<DerivedV>& V,
+                               const Eigen::MatrixBase<DerivedF>& F,
+                               Eigen::PlainObjectBase<DerivedE>& E,
+                               Eigen::PlainObjectBase<DerivedOE>& oE,
+                               Eigen::SparseMatrix<ScalarG>& G);
+
+
+  // Version that uses intrinsic quantities as input
+  //
+  // Inputs:
+  //  F: input mesh connectivity
+  //  l_sq: squared edge lengths of each halfedge
+  //  dA: double area of each face
+  //  E: a mapping from each halfedge to each edge.
+  //  oE: the orientation of each halfedge compared to the orientation of the
+  //      actual edge.
+  //
+  // Outputs:
+  //  G: computed gradient matrix
+
+  template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
+  typename DerivedOE, typename ScalarG>
+  IGL_INLINE void
+  scalar_to_cr_vector_gradient_intrinsic(
+                                         const Eigen::MatrixBase<DerivedF>& F,
+                                         const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+                                         const Eigen::MatrixBase<DerivedE>& E,
+                                         const Eigen::MatrixBase<DerivedOE>& oE,
+                                         Eigen::SparseMatrix<ScalarG>& G);
+
+  template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
+  typename DerivedE, typename DerivedOE, typename ScalarG>
+  IGL_INLINE void
+  scalar_to_cr_vector_gradient_intrinsic(
+                                         const Eigen::MatrixBase<DerivedF>& F,
+                                         const Eigen::MatrixBase<DerivedL_sq>& l_sq,
+                                         const Eigen::MatrixBase<DeriveddA>& dA,
+                                         const Eigen::MatrixBase<DerivedE>& E,
+                                         const Eigen::MatrixBase<DerivedOE>& oE,
+                                         Eigen::SparseMatrix<ScalarG>& G);
+
+
+}
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "scalar_to_cr_vector_gradient.cpp"
+#endif
+
+#endif

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

include/igl/writePLY.cpp

@@ -45,7 +45,7 @@ bool writePLY(
   const std::vector<std::string> & EDheader,
 
   const std::vector<std::string> & comments,
-  bool isBinary
+  FileEncoding encoding
    )
 {
     typedef typename DerivedV::Scalar VScalar;
@@ -155,7 +155,7 @@ bool writePLY(
       file.get_comments().push_back(a);
 
     // Write a binary file
-    file.write(ply_stream, isBinary);
+    file.write(ply_stream, (encoding == FileEncoding::Binary));
 
     return true;
 }
@@ -188,12 +188,12 @@ bool writePLY(
   const std::vector<std::string> & EDheader,
 
   const std::vector<std::string> & comments,
-  bool isBinary
+  FileEncoding encoding
    )
 {
   try
   {
-    if(isBinary)
+    if(encoding == FileEncoding::Binary)
     {
       std::filebuf fb_binary;
       fb_binary.open(filename , std::ios::out | std::ios::binary);
@@ -202,7 +202,7 @@ bool writePLY(
         std::cerr << "writePLY: Error opening file " << filename << std::endl;
         return false; //throw std::runtime_error("failed to open " + filename);
       }
-      return writePLY(outstream_binary,V,F,E,N,UV,VD,VDheader,FD,FDheader,ED,EDheader,comments,isBinary);
+      return writePLY(outstream_binary,V,F,E,N,UV,VD,VDheader,FD,FDheader,ED,EDheader,comments,encoding);
     } else {
       std::filebuf fb_ascii;
       fb_ascii.open(filename, std::ios::out);
@@ -211,7 +211,7 @@ bool writePLY(
         std::cerr << "writePLY: Error opening file " << filename << std::endl;
         return false; //throw std::runtime_error("failed to open " + filename);
       }
-      return writePLY(outstream_ascii,V,F,E,N,UV,VD,VDheader,FD,FDheader,ED,EDheader,comments,isBinary);
+      return writePLY(outstream_ascii,V,F,E,N,UV,VD,VDheader,FD,FDheader,ED,EDheader,comments,encoding);
     }
   }
   catch(const std::exception& e)
@@ -234,7 +234,7 @@ bool writePLY(
   Eigen::MatrixXd _dummy;
   std::vector<std::string> _dummy_header;
 
-  return writePLY(filename,V,F,_dummy, _dummy, _dummy, _dummy, _dummy_header, _dummy, _dummy_header, _dummy, _dummy_header, _dummy_header, true);
+  return writePLY(filename,V,F,_dummy, _dummy, _dummy, _dummy, _dummy_header, _dummy, _dummy_header, _dummy, _dummy_header, _dummy_header, FileEncoding::Binary);
 }
 
 template <
@@ -252,7 +252,7 @@ bool writePLY(
   Eigen::MatrixXd _dummy;
   std::vector<std::string> _dummy_header;
 
-  return writePLY(filename,V,F,E, _dummy, _dummy, _dummy, _dummy_header, _dummy, _dummy_header, _dummy, _dummy_header, _dummy_header, true);
+  return writePLY(filename,V,F,E, _dummy, _dummy, _dummy, _dummy_header, _dummy, _dummy_header, _dummy, _dummy_header, _dummy_header, FileEncoding::Binary);
 }
 
 
@@ -273,7 +273,7 @@ bool writePLY(
   Eigen::MatrixXd _dummy;
   std::vector<std::string> _dummy_header;
 
-  return writePLY(filename,V,F,_dummy, N,UV, _dummy, _dummy_header, _dummy, _dummy_header, _dummy, _dummy_header, _dummy_header, true);
+  return writePLY(filename,V,F,_dummy, N,UV, _dummy, _dummy_header, _dummy, _dummy_header, _dummy, _dummy_header, _dummy_header, FileEncoding::Binary);
 }
 
 template <
@@ -295,7 +295,7 @@ bool writePLY(
   Eigen::MatrixXd _dummy;
   std::vector<std::string> _dummy_header;
 
-  return writePLY(filename,V,F,E, N,UV, _dummy, _dummy_header, _dummy, _dummy_header, _dummy, _dummy_header, _dummy_header, true);
+  return writePLY(filename,V,F,E, N,UV, _dummy, _dummy_header, _dummy, _dummy_header, _dummy, _dummy_header, _dummy_header, FileEncoding::Binary);
 }
 
 template <
@@ -306,14 +306,14 @@ bool writePLY(
   const std::string & filename,
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedF> & F,
-  bool force_ascii
+  FileEncoding encoding
    )
 {
   Eigen::MatrixXd _dummy(0,0);
   std::vector<std::string> _dummy_header;
 
   return writePLY(filename,V,F,_dummy, _dummy,_dummy, _dummy, _dummy_header,
-                         _dummy, _dummy_header, _dummy, _dummy_header, _dummy_header, force_ascii);
+                         _dummy, _dummy_header, _dummy, _dummy_header, _dummy_header, encoding);
 }
 
 template <
@@ -326,14 +326,14 @@ bool writePLY(
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedF> & F,
   const Eigen::MatrixBase<DerivedE> & E,
-  bool force_ascii
+  FileEncoding encoding
    )
 {
   Eigen::MatrixXd _dummy(0,0);
   std::vector<std::string> _dummy_header;
 
   return writePLY(filename,V,F,E, _dummy,_dummy, _dummy, _dummy_header,
-                         _dummy, _dummy_header, _dummy, _dummy_header, _dummy_header, force_ascii);
+                         _dummy, _dummy_header, _dummy, _dummy_header, _dummy_header, encoding);
 }
 
 
@@ -361,7 +361,7 @@ bool writePLY(
   std::vector<std::string> _dummy_header;
 
   return writePLY(filename,V,F,_dummy, N, UV, VD, VDheader,
-                         _dummy, _dummy_header, _dummy, _dummy_header, comments, true);
+                         _dummy, _dummy_header, _dummy, _dummy_header, comments, FileEncoding::Binary);
 }
 
 
@@ -391,7 +391,7 @@ bool writePLY(
   std::vector<std::string> _dummy_header;
 
   return writePLY(filename,V,F,E, N, UV, VD, VDheader,
-                         _dummy, _dummy_header, _dummy, _dummy_header, comments, true);
+                         _dummy, _dummy_header, _dummy, _dummy_header, comments, FileEncoding::Binary);
 
 }
 
@@ -404,6 +404,6 @@ bool writePLY(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 template bool igl::writePLY<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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
-template bool igl::writePLY<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, bool);
-template bool igl::writePLY<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, bool);
+template bool igl::writePLY<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, igl::FileEncoding);
+template bool igl::writePLY<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, std::vector<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::basic_string<char, std::char_traits<char>, std::allocator<char> > > > const&, igl::FileEncoding);
 #endif

include/igl/writePLY.h

@@ -1,6 +1,7 @@
 #ifndef IGL_WRITEPLY_H
 #define IGL_WRITEPLY_H
 #include <igl/igl_inline.h>
+#include <igl/FileEncoding.h>
 
 #include <string>
 #include <iostream>
@@ -16,7 +17,7 @@ namespace igl
   //   Derived from Eigen matrix parameters
   // Inputs:
   //  ply_stream  ply file output stream
-  //   V  (#V,3) matrix of vertex positions 
+  //   V  (#V,3) matrix of vertex positions
   //   F  (#F,3) list of face indices into vertex positions
   //   E  (#E,2) list of edge indices into vertex positions
   //   N  (#V,3) list of normals
@@ -28,7 +29,7 @@ namespace igl
   //   ED (#E,*) additional edge data
   //   Eheader (#E) list of edge data headers
   //   comments (*) file comments
-  //   isBinary - bool, to set binary or ascii file format
+  //   encoding - enum, to set binary or ascii file format
   // Returns true on success, false on errors
   template <
   typename DerivedV,
@@ -58,7 +59,7 @@ bool writePLY(
   const std::vector<std::string> & EDheader,
 
   const std::vector<std::string> & comments,
-  bool isBinary
+  FileEncoding encoding
    );
 
   // write triangular mesh to ply file
@@ -67,7 +68,7 @@ bool writePLY(
   //   Derived from Eigen matrix parameters
   // Inputs:
   //  filename  ply file name
-  //   V  (#V,3) matrix of vertex positions 
+  //   V  (#V,3) matrix of vertex positions
   //   F  (#F,3) list of face indices into vertex positions
   //   E  (#E,2) list of edge indices into vertex positions
   //   N  (#V,3) list of normals
@@ -79,7 +80,7 @@ bool writePLY(
   //   ED (#E,*) additional edge data
   //   Eheader (#E) list of edge data headers
   //   comments (*) file comments
-  //   isBinary - bool, to set binary or ascii file format
+  //   encoding - enum, to set binary or ascii file format
   // Returns true on success, false on errors
 template <
   typename DerivedV,
@@ -109,7 +110,7 @@ bool writePLY(
   const std::vector<std::string> & EDheader,
 
   const std::vector<std::string> & comments,
-  bool isBinary
+  FileEncoding encoding
    );
 
 template <
@@ -175,8 +176,8 @@ bool writePLY(
   const std::string & filename,
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedF> & F,
-  bool force_ascii
-   );   
+  FileEncoding encoding
+   );
 
 template <
   typename DerivedV,
@@ -188,7 +189,7 @@ bool writePLY(
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedF> & F,
   const Eigen::MatrixBase<DerivedE> & E,
-  bool force_ascii
+  FileEncoding encoding
    );
 
 template <

include/igl/writeSTL.cpp

@@ -14,11 +14,11 @@ IGL_INLINE bool igl::writeSTL(
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedF> & F,
   const Eigen::MatrixBase<DerivedN> & N,
-  const bool ascii)
+  FileEncoding encoding)
 {
   using namespace std;
   assert(N.rows() == 0 || F.rows() == N.rows());
-  if(ascii)
+  if(encoding == FileEncoding::Ascii)
   {
     FILE * stl_file = fopen(filename.c_str(),"w");
     if(stl_file == NULL)
@@ -103,21 +103,21 @@ IGL_INLINE bool igl::writeSTL(
   const std::string & filename,
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedF> & F,
-  const bool ascii)
+  FileEncoding encoding)
 {
-  return writeSTL(filename,V,F, Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic>(), ascii);
+  return writeSTL(filename,V,F, Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic>(), encoding);
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
-template bool igl::writeSTL<Eigen::Matrix<float, -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::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool);
+template bool igl::writeSTL<Eigen::Matrix<float, -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::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::FileEncoding);
 // generated by autoexplicit.sh
-template bool igl::writeSTL<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, bool);
+template bool igl::writeSTL<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, igl::FileEncoding);
 // generated by autoexplicit.sh
-template bool igl::writeSTL<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, bool);
+template bool igl::writeSTL<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, igl::FileEncoding);
 // generated by autoexplicit.sh
-template bool igl::writeSTL<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, bool);
-template bool igl::writeSTL<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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool);
-template bool igl::writeSTL<Eigen::Matrix<double, 8, 3, 0, 8, 3>, Eigen::Matrix<int, 12, 3, 0, 12, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 8, 3, 0, 8, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, 12, 3, 0, 12, 3> > const&, bool);
+template bool igl::writeSTL<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::FileEncoding);
+template bool igl::writeSTL<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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::FileEncoding);
+template bool igl::writeSTL<Eigen::Matrix<double, 8, 3, 0, 8, 3>, Eigen::Matrix<int, 12, 3, 0, 12, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 8, 3, 0, 8, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, 12, 3, 0, 12, 3> > const&, igl::FileEncoding);
 #endif

include/igl/writeSTL.h

@@ -1,13 +1,14 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 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 
+//
+// 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_WRITESTL_H
 #define IGL_WRITESTL_H
 #include "igl_inline.h"
+#include <igl/FileEncoding.h>
 
 #ifndef IGL_NO_EIGEN
 #  include <Eigen/Core>
@@ -15,7 +16,7 @@
 #include <string>
 #include <vector>
 
-namespace igl 
+namespace igl
 {
   // Write a mesh to an stl file.
   //
@@ -27,7 +28,7 @@ namespace igl
   //   V  double matrix of vertex positions  #F*3 by 3
   //   F  index matrix of triangle indices #F by 3
   //   N  double matrix of vertex positions  #F by 3
-  //   asci  write ascii file {true}
+  //   encoding enum to set file encoding (ascii by default)
   // Returns true on success, false on errors
   //
   template <typename DerivedV, typename DerivedF, typename DerivedN>
@@ -36,13 +37,13 @@ namespace igl
     const Eigen::MatrixBase<DerivedV> & V,
     const Eigen::MatrixBase<DerivedF> & F,
     const Eigen::MatrixBase<DerivedN> & N,
-    const bool ascii=true);
+    FileEncoding encoding=FileEncoding::Ascii);
   template <typename DerivedV, typename DerivedF>
   IGL_INLINE bool writeSTL(
     const std::string & filename,
     const Eigen::MatrixBase<DerivedV> & V,
     const Eigen::MatrixBase<DerivedF> & F,
-    const bool ascii=true);
+    FileEncoding encoding=FileEncoding::Ascii);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/write_triangle_mesh.cpp

@@ -21,7 +21,7 @@ IGL_INLINE bool igl::write_triangle_mesh(
   const std::string str,
   const Eigen::MatrixBase<DerivedV>& V,
   const Eigen::MatrixBase<DerivedF>& F,
-  const bool force_ascii)
+  FileEncoding encoding)
 {
   using namespace std;
   // dirname, basename, extension and filename
@@ -41,10 +41,10 @@ IGL_INLINE bool igl::write_triangle_mesh(
     return writeOFF(str,V,F);
   }else if(e == "ply")
   {
-    return writePLY(str,V,F,force_ascii);
+    return writePLY(str,V,F,encoding);
   }else if(e == "stl")
   {
-    return writeSTL(str,V,F,force_ascii);
+    return writeSTL(str,V,F,encoding);
   }else if(e == "wrl")
   {
     return writeWRL(str,V,F);
@@ -59,12 +59,12 @@ IGL_INLINE bool igl::write_triangle_mesh(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
-template bool igl::write_triangle_mesh<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, bool);
+template bool igl::write_triangle_mesh<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, igl::FileEncoding);
 // generated by autoexplicit.sh
-template bool igl::write_triangle_mesh<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, bool);
+template bool igl::write_triangle_mesh<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, igl::FileEncoding);
 // generated by autoexplicit.sh
-template bool igl::write_triangle_mesh<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, bool);
+template bool igl::write_triangle_mesh<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::FileEncoding);
 // generated by autoexplicit.sh
-template bool igl::write_triangle_mesh<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> >, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, const bool);
-template bool igl::write_triangle_mesh<Eigen::Matrix<double, 8, 3, 0, 8, 3>, Eigen::Matrix<int, 12, 3, 0, 12, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<double, 8, 3, 0, 8, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, 12, 3, 0, 12, 3> > const&, bool);
+template bool igl::write_triangle_mesh<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> >, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::FileEncoding);
+template bool igl::write_triangle_mesh<Eigen::Matrix<double, 8, 3, 0, 8, 3>, Eigen::Matrix<int, 12, 3, 0, 12, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<double, 8, 3, 0, 8, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, 12, 3, 0, 12, 3> > const&, igl::FileEncoding);
 #endif

include/igl/write_triangle_mesh.h

@@ -8,6 +8,7 @@
 #ifndef IGL_WRITE_TRIANGLE_MESH_H
 #define IGL_WRITE_TRIANGLE_MESH_H
 #include "igl_inline.h"
+#include <igl/FileEncoding.h>
 
 #include <Eigen/Core>
 #include <string>
@@ -25,14 +26,14 @@ namespace igl
   //   str  path to file
   //   V  eigen double matrix #V by 3
   //   F  eigen int matrix #F by 3
-  //   force_ascii  force ascii format even if binary is available
+  //   encoding  set file encoding (ascii or binary) when both are available
   // Returns true iff success
   template <typename DerivedV, typename DerivedF>
   IGL_INLINE bool write_triangle_mesh(
     const std::string str,
     const Eigen::MatrixBase<DerivedV>& V,
     const Eigen::MatrixBase<DerivedF>& F,
-    const bool force_ascii = true);
+    FileEncoding encoding = FileEncoding::Ascii);
 }
 
 #ifndef IGL_STATIC_LIBRARY

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/orient_halfedges.cpp

@@ -0,0 +1,102 @@
+#include <test_common.h>
+#include <igl/orient_halfedges.h>
+#include <igl/is_border_vertex.h>
+#include <igl/edges.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/triangle_triangle_adjacency.h>
+#include <igl/EPS.h>
+
+#include <vector>
+
+
+TEST_CASE("orient_halfedges: sanity checks", "[igl]")
+{
+  const auto meshes = test_common::manifold_meshes();
+
+  for(const auto& mesh : meshes) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    igl::read_triangle_mesh(test_common::data_path(mesh), V, F);
+    Eigen::MatrixXi I;
+    igl::remove_unreferenced(Eigen::MatrixXd(V), Eigen::MatrixXi(F), V, F,
+     I);
+
+    Eigen::MatrixXi TT, TTi;
+    igl::triangle_triangle_adjacency(F, TT, TTi);
+    // Fix mis-match convention
+    {
+      Eigen::PermutationMatrix<3,3> perm(3);
+      perm.indices() = Eigen::Vector3i(1,2,0);
+      TT = (TT*perm).eval();
+      TTi = (TTi*perm).eval();
+      for(int i=0;i<TTi.rows();i++) {
+        for(int j=0;j<TTi.cols();j++) {
+          TTi(i,j)=TTi(i,j)==-1?-1:(TTi(i,j)+3-1)%3;
+        }
+      }
+    }
+
+    Eigen::MatrixXi E, oE;
+    igl::orient_halfedges(F, E, oE);
+
+    const int m = E.maxCoeff()+1;
+    std::vector<bool> b(m);
+    for(int i=0; i<F.rows(); ++i) {
+      for(int j=0; j<3; ++j) {
+        b[E(i,j)] = TT(i,j)<0;
+      }
+    }
+    int nb = 0;
+    for(int i=0; i<b.size(); ++i) {
+      if(b[i]) {
+        ++nb;
+      }
+    }
+
+
+    //Perform a variety of sanity checks.
+
+    //Correct number of edges
+    Eigen::MatrixXi sanityEdges;
+    igl::edges(F, sanityEdges);
+    REQUIRE(m == sanityEdges.rows());
+
+    //All border halfedges edges have orientation 1 only, all others
+    // orientations 1 and -1, so oE must sum to the number of border
+    // edges.
+    REQUIRE(nb == oE.array().sum());
+
+    //Every border halfedge has orientation 1. Every other edge appeats
+    // with orientation 1 and -1.
+    std::vector<int> appeared1(m,0), appearedm1(m,0);
+    for(int i=0; i<F.rows(); ++i) {
+      for(int j=0; j<3; ++j) {
+        if(oE(i,j)==1) {
+          ++appeared1[E(i,j)];
+        } else if(oE(i,j)==-1) {
+          ++appearedm1[E(i,j)];
+        } else {
+          REQUIRE(false); //Only 1 and -1 should occur.
+        }
+      }
+    }
+    for(int i=0; i<m; ++i) {
+      REQUIRE(appeared1[i]==1);
+      if(b[i]) {
+        REQUIRE(appearedm1[i]==0);
+      } else {
+        REQUIRE(appearedm1[i]==1);
+      }
+    }
+
+    //Two opposite halfedges always map to the same edge
+    for(int i=0; i<F.rows(); ++i) {
+      for(int j=0; j<3; ++j) {
+        if(TT(i,j)>=0) {
+          REQUIRE(E(i,j) == E(TT(i,j),TTi(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));
+    }
+}

tests/include/igl/writePLY.cpp

@@ -10,7 +10,7 @@ TEST_CASE("writePLY: bunny.ply", "[igl]")
     std::ifstream f(test_common::data_path("bunny.ply"));
     REQUIRE (f.good());
     f.close();
-    
+
     Eigen::MatrixXd V1,N1,UV1,VD1,FD1,ED1;
     std::vector<std::string> Vheader1,Fheader1,Eheader1,comments1;
 
@@ -19,7 +19,7 @@ TEST_CASE("writePLY: bunny.ply", "[igl]")
     // load test data first
     REQUIRE (igl::readPLY(test_common::data_path("bunny.ply"), V1, F1, E1, N1, UV1, VD1,Vheader1, FD1,Fheader1, ED1,Eheader1,comments1));
 
-    // add more data 
+    // add more data
     Vheader1.push_back("dummy_data");
     Eigen::VectorXd dummy_data(V1.rows());
     for(size_t i=0;i<V1.rows();++i)
@@ -29,7 +29,7 @@ TEST_CASE("writePLY: bunny.ply", "[igl]")
 
 
     // test that saving preserves all the data, including new data column
-    REQUIRE (igl::writePLY("test_bunny.ply", V1, F1, E1, N1, UV1, VD2, Vheader1, FD1,Fheader1, ED1, Eheader1, comments1, true));
+    REQUIRE (igl::writePLY("test_bunny.ply", V1, F1, E1, N1, UV1, VD2, Vheader1, FD1,Fheader1, ED1, Eheader1, comments1, igl::FileEncoding::Binary));
 
     Eigen::MatrixXd V,N,UV,VD,FD,ED;
     Eigen::MatrixXi F,E;
@@ -85,7 +85,7 @@ TEST_CASE("writePLY: bunny.ply float", "[igl]")
     std::ifstream f(test_common::data_path("bunny.ply"));
     REQUIRE (f.good());
     f.close();
-    
+
     Eigen::MatrixXf V1,N1,UV1,VD1,FD1,ED1;
     std::vector<std::string> Vheader1,Fheader1,Eheader1,comments1;
 
@@ -94,7 +94,7 @@ TEST_CASE("writePLY: bunny.ply float", "[igl]")
     // load test data first
     REQUIRE (igl::readPLY(test_common::data_path("bunny.ply"), V1, F1, E1, N1, UV1, VD1,Vheader1, FD1,Fheader1, ED1,Eheader1,comments1));
 
-    // add more data 
+    // add more data
     Vheader1.push_back("dummy_data");
     Eigen::VectorXf dummy_data(V1.rows());
     for(size_t i=0;i<V1.rows();++i)
@@ -104,7 +104,7 @@ TEST_CASE("writePLY: bunny.ply float", "[igl]")
 
 
     // test that saving preserves all the data, including new data column
-    REQUIRE (igl::writePLY("test_bunny.ply", V1, F1, E1, N1, UV1, VD2, Vheader1, FD1,Fheader1, ED1, Eheader1, comments1, true));
+    REQUIRE (igl::writePLY("test_bunny.ply", V1, F1, E1, N1, UV1, VD2, Vheader1, FD1,Fheader1, ED1, Eheader1, comments1, igl::FileEncoding::Binary));
 
     Eigen::MatrixXf V,N,UV,VD,FD,ED;
     Eigen::MatrixXi F,E;

tutorial/109_ImGuizmo/main.cpp

@@ -1,101 +1,53 @@
+#include <igl/read_triangle_mesh.h>
 #include <igl/opengl/glfw/Viewer.h>
-#include <igl/opengl/glfw/imgui/ImGuiMenu.h>
-#include <igl/opengl/glfw/imgui/ImGuizmo.h>
-#include <imgui_impl_glfw.h>
-#include <imgui_impl_opengl3.h>
+#include <igl/opengl/glfw/imgui/ImGuizmoPlugin.h>
 #include <GLFW/glfw3.h>
 #include <imgui/imgui.h>
-#include "tutorial_shared_path.h"
+#include <imgui/imgui_internal.h>
+#include <imgui_impl_glfw.h>
+#include <imgui_impl_opengl3.h>
+#include <imguizmo/ImGuizmo.h>
 
-class SlicingPlugin : public igl::opengl::glfw::imgui::ImGuiMenu
+int main(int argc, char *argv[])
 {
-  igl::opengl::ViewerData data;
-  // Hard Coded Plane Mesh
-  const Eigen::MatrixXd OV = (Eigen::MatrixXd(4, 3) <<
-    -0.5, -0.5, 0.0,
-    -0.5,  0.5, 0.0,
-    0.5,  0.5, 0.0,
-    0.5, -0.5, 0.0).finished();
-  const Eigen::MatrixXi OF = (Eigen::MatrixXi(2, 3) <<
-    0, 2, 1,
-    0, 3, 2).finished();
-
-  virtual void init(igl::opengl::glfw::Viewer *_viewer) override
+  // Load a mesh from file
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  igl::read_triangle_mesh(argc>1?argv[1]: TUTORIAL_SHARED_PATH "/cow.off",V,F);
+  // Set up viewer
+  igl::opengl::glfw::Viewer vr;
+  vr.data().set_mesh(V,F);
+  // Custom menu
+  igl::opengl::glfw::imgui::ImGuizmoPlugin plugin;
+  vr.plugins.push_back(&plugin);
+  // Initialize ImGuizmo at mesh centroid
+  plugin.T.block(0,3,3,1) = 
+    0.5*(V.colwise().maxCoeff() + V.colwise().minCoeff()).transpose().cast<float>();
+  // Update can be applied relative to this remembered initial transform
+  const Eigen::Matrix4f T0 = plugin.T;
+  // Attach callback to apply imguizmo's transform to mesh
+  plugin.callback = [&](const Eigen::Matrix4f & T)
   {
-    ImGuiMenu::init(_viewer);
-    // Load slicing plane into viewer (a square mesh)
-    data.set_mesh(SlicingPlugin::OV, SlicingPlugin::OF);
-    data.set_face_based(true);
-    data.set_colors(Eigen::RowVector4d(224, 86, 253, 128)/255.0);
-    data.show_lines = false;
-  }
-
-  virtual bool pre_draw() override {
-    ImGuiMenu::pre_draw();
-    ImGui::PushStyleVar(ImGuiStyleVar_WindowBorderSize, 0);
-    ImGuizmo::BeginFrame();
-    ImGui::PopStyleVar();
-    return false;
-  }
-
-  virtual bool post_draw() override {
-    viewer->core().draw(data);
-    ImGuiMenu::post_draw();
-    return false;
-  }
-
-  virtual void draw_custom_window() override {
-    float menu_width = 180.f * menu_scaling();
-    ImGui::SetNextWindowPos(ImVec2(menu_width, 0.0f), ImGuiCond_FirstUseEver);
-    ImGui::SetNextWindowSize(ImVec2(0.0f, 0.0f), ImGuiCond_FirstUseEver);
-    ImGui::SetNextWindowSizeConstraints(ImVec2(1.5f * menu_width, -1.0f), ImVec2(1.5f * menu_width, -1.0f));
-    bool _guizmo_menu_visible = true;
-    ImGui::Begin(
-      "ImGuizmo Tools", &_guizmo_menu_visible,
-      ImGuiWindowFlags_NoSavedSettings
-      | ImGuiWindowFlags_AlwaysAutoResize
-    );
-
-    draw_imguizmo_menu();
-
-    ImGui::End();
-  }
-
-  virtual void draw_imguizmo_menu() {
-    static Eigen::Matrix4f matrix = Eigen::Matrix4f::Identity();
-    const float z = viewer->core().camera_base_zoom;
-    Eigen::Affine3f rescale = Eigen::Scaling(0.5f * z)
-      * Eigen::Translation3f(viewer->core().camera_base_translation);
-    Eigen::Affine3f view = Eigen::Affine3f( viewer->core().view * 1./viewer->core().camera_zoom ) * rescale.inverse();
-    Eigen::Matrix4f proj = viewer->core().proj;
-    if(viewer->core().orthographic)
+    const Eigen::Matrix4d TT = (T*T0.inverse()).cast<double>().transpose();
+    vr.data().set_vertices(
+      (V.rowwise().homogeneous()*TT).rowwise().hnormalized());
+    vr.data().compute_normals();
+  };
+  // Maya-style keyboard shortcuts for operation
+  vr.callback_key_pressed = [&](decltype(vr) &,unsigned int key, int mod)
+  {
+    switch(key)
     {
-      view(2,3) -= 50; // Adjust depth for view transform
+      case 'W': case 'w': plugin.operation = ImGuizmo::TRANSLATE; return true;
+      case 'E': case 'e': plugin.operation = ImGuizmo::ROTATE;    return true;
+      case 'R': case 'r': plugin.operation = ImGuizmo::SCALE;     return true;
     }
-
-    igl::opengl::glfw::imgui::EditTransform(view.matrix().data(), proj.data(), matrix.data(), viewer->core().orthographic);
-
-    // Transform the slicing plane according to
-    // ImGuizmo tool manipulations in the viewer
-    Eigen::Affine3f model = Eigen::Affine3f(rescale.inverse()*matrix);
-    Eigen::MatrixXd V = (SlicingPlugin::OV.rowwise().homogeneous()
-      * model.matrix().cast<double>().transpose()).rowwise().hnormalized();
-
-    if (data.V.rows() == V.rows())
-      data.set_vertices(V);
-  }
-
-};
-
-int main(int argc, char *argv[])
-{
-  // Plot the mesh
-  igl::opengl::glfw::Viewer viewer;
-  viewer.load_mesh_from_file(argc>1?argv[1]: TUTORIAL_SHARED_PATH "/cow.off");
-
-  // Custom menu
-  SlicingPlugin menu;
-  viewer.plugins.push_back(&menu);
-
-  viewer.launch();
+    return false;
+  };
+  std::cout<<R"(
+W,w  Switch to translate operation
+E,e  Switch to rotate operation
+R,r  Switch to scale operation
+)";
+  vr.launch();
 }

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/712_DataSmoothing/main.cpp

@@ -1,5 +1,6 @@
 #include <igl/read_triangle_mesh.h>
 #include <igl/hessian_energy.h>
+#include <igl/curved_hessian_energy.h>
 #include <igl/massmatrix.h>
 #include <igl/cotmatrix.h>
 #include <igl/isolines_map.h>
@@ -7,6 +8,7 @@
 #include <igl/vertex_components.h>
 #include <igl/remove_unreferenced.h>
 #include <igl/opengl/glfw/Viewer.h>
+#include <igl/heat_geodesics.h>
 
 #include <Eigen/Core>
 #include <Eigen/SparseCholesky>
@@ -22,86 +24,160 @@
 
 int main(int argc, char * argv[])
 {
-    typedef Eigen::SparseMatrix<double> SparseMat;
-
-    //Read our mesh
-    Eigen::MatrixXd V;
-    Eigen::MatrixXi F;
-    if(!igl::read_triangle_mesh(
-        argc>1?argv[1]: TUTORIAL_SHARED_PATH "/beetle.off",V,F)) {
-        std::cout << "Failed to load mesh." << std::endl;
+  typedef Eigen::SparseMatrix<double> SparseMat;
+  srand(57);
+  
+  //Read our mesh
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  if(!igl::read_triangle_mesh(
+                              argc>1?argv[1]: TUTORIAL_SHARED_PATH "/beetle.off",V,F)) {
+    std::cout << "Failed to load mesh." << std::endl;
+  }
+  
+  //Constructing an exact function to smooth
+  igl::HeatGeodesicsData<double> hgData;
+  igl::heat_geodesics_precompute(V, F, hgData);
+  Eigen::VectorXd heatDist;
+  Eigen::VectorXi gamma(1); gamma << 1947; //1631;
+  igl::heat_geodesics_solve(hgData, gamma, heatDist);
+  Eigen::VectorXd zexact =
+  0.1*(heatDist.array() + (-heatDist.maxCoeff())).pow(2)
+  + 3*V.block(0,1,V.rows(),1).array().cos();
+  
+  //Make the exact function noisy
+  const double s = 0.1*(zexact.maxCoeff() - zexact.minCoeff());
+  Eigen::VectorXd znoisy = zexact + s*Eigen::VectorXd::Random(zexact.size());
+  
+  //Constructing the squared Laplacian and squared Hessian energy
+  SparseMat L, M;
+  igl::cotmatrix(V, F, L);
+  igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_BARYCENTRIC, M);
+  Eigen::SimplicialLDLT<SparseMat> solver(M);
+  SparseMat MinvL = solver.solve(L);
+  SparseMat QL = L.transpose()*MinvL;
+  SparseMat QH;
+  igl::hessian_energy(V, F, QH);
+  SparseMat QcH;
+  igl::curved_hessian_energy(V, F, QcH);
+  
+  //Solve to find Laplacian-smoothed Hessian-smoothed, and
+  // curved-Hessian-smoothed solutions
+  const double al = 3e-7;
+  Eigen::SimplicialLDLT<SparseMat> lapSolver(al*QL + (1.-al)*M);
+  Eigen::VectorXd zl = lapSolver.solve(al*M*znoisy);
+  const double ah = 2e-7;
+  Eigen::SimplicialLDLT<SparseMat> hessSolver(ah*QH + (1.-ah)*M);
+  Eigen::VectorXd zh = hessSolver.solve(ah*M*znoisy);
+  const double ach = 3e-7;
+  Eigen::SimplicialLDLT<SparseMat> curvedHessSolver(al*QcH + (1.-ach)*M);
+  Eigen::VectorXd zch = curvedHessSolver.solve(ach*M*znoisy);
+  
+  //Viewer that shows all functions: zexact, znoisy, zl, zh
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(V,F);
+  viewer.data().show_lines = false;
+  viewer.callback_key_down =
+  [&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
+  {
+    //Graduate result to show isolines, then compute color matrix
+    const Eigen::VectorXd* z;
+    switch(key) {
+      case '1':
+        z = &zexact;
+        break;
+      case '2':
+        z = &znoisy;
+        break;
+      case '3':
+        z = &zl;
+        break;
+      case '4':
+        z = &zh;
+        break;
+      case '5':
+        z = &zch;
+        break;
+      default:
+        return false;
     }
-
-    //Constructing an exact function to smooth
-    Eigen::VectorXd zexact = V.block(0,2,V.rows(),1).array()
-        + 0.5*V.block(0,1,V.rows(),1).array()
-        + V.block(0,1,V.rows(),1).array().pow(2)
-        + V.block(0,2,V.rows(),1).array().pow(3);
-
-    //Make the exact function noisy
-    srand(5);
-    const double s = 0.2*(zexact.maxCoeff() - zexact.minCoeff());
-    Eigen::VectorXd znoisy = zexact + s*Eigen::VectorXd::Random(zexact.size());
-
-    //Constructing the squared Laplacian and squared Hessian energy
-    SparseMat L, M;
-    igl::cotmatrix(V, F, L);
-    igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_BARYCENTRIC, M);
-    Eigen::SimplicialLDLT<SparseMat> solver(M);
-    SparseMat MinvL = solver.solve(L);
-    SparseMat QL = L.transpose()*MinvL;
-    SparseMat QH;
-    igl::hessian_energy(V, F, QH);
-
-    //Solve to find Laplacian-smoothed and Hessian-smoothed solutions
-    const double al = 8e-4;
-    Eigen::SimplicialLDLT<SparseMat> lapSolver(al*QL + (1.-al)*M);
-    Eigen::VectorXd zl = lapSolver.solve(al*M*znoisy);
-    const double ah = 5e-6;
-    Eigen::SimplicialLDLT<SparseMat> hessSolver(ah*QH + (1.-ah)*M);
-    Eigen::VectorXd zh = hessSolver.solve(ah*M*znoisy);
-
-    //Viewer that shows all functions: zexact, znoisy, zl, zh
-    igl::opengl::glfw::Viewer viewer;
-    viewer.data().set_mesh(V,F);
-    viewer.data().show_lines = false;
-    viewer.callback_key_down =
-      [&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
-      {
-        //Graduate result to show isolines, then compute color matrix
-        const Eigen::VectorXd* z;
-        switch(key) {
-            case '1':
-                z = &zexact;
-                break;
-            case '2':
-                z = &znoisy;
-                break;
-            case '3':
-                z = &zl;
-                break;
-            case '4':
-                z = &zh;
-                break;
-            default:
-                return false;
-        }
-        viewer.data().set_data(*z);
-        return true;
-    };
-    std::cout << R"(Usage:
+    viewer.data().set_data(*z);
+    return true;
+  };
+  std::cout << R"(Smoothing a noisy function.
+Usage:
 1  Show original function
 2  Show noisy function
 3  Biharmonic smoothing (zero Neumann boundary)
-4  Biharmonic smoothing (natural Hessian boundary)
+4  Biharmonic smoothing (natural planar Hessian boundary)
+5  Biharmonic smoothing (natural curved Hessian boundary)
 
 )";
-    Eigen::MatrixXd CM;
-    igl::parula(Eigen::VectorXd::LinSpaced(21,0,1).eval(),false,CM);
-    igl::isolines_map(Eigen::MatrixXd(CM),CM);
-    viewer.data().set_colormap(CM);
-    viewer.data().set_data(znoisy);
-    viewer.launch();
+  Eigen::MatrixXd CM;
+  igl::parula(Eigen::VectorXd::LinSpaced(21,0,1).eval(),false,CM);
+  igl::isolines_map(Eigen::MatrixXd(CM),CM);
+  viewer.data().set_colormap(CM);
+  viewer.data().set_data(znoisy);
+  viewer.launch();
+  
+  
+  //Constructing a step function to smooth
+  Eigen::VectorXd zstep = Eigen::VectorXd::Zero(V.rows());
+  for(int i=0; i<V.rows(); ++i) {
+    zstep(i) = V(i,2)<-0.25 ? 1. : (V(i,2)>0.31 ? 2. : 0);
+  }
+  
+  //Smooth that function
+  const double sl = 2e-5;
+  Eigen::SimplicialLDLT<SparseMat> stepLapSolver(sl*QL + (1.-sl)*M);
+  Eigen::VectorXd stepzl = stepLapSolver.solve(al*M*zstep);
+  const double sh = 6e-6;
+  Eigen::SimplicialLDLT<SparseMat> stepHessSolver(sh*QH + (1.-sh)*M);
+  Eigen::VectorXd stepzh = stepHessSolver.solve(ah*M*zstep);
+  const double sch = 2e-5;
+  Eigen::SimplicialLDLT<SparseMat> stepCurvedHessSolver(sl*QcH + (1.-sch)*M);
+  Eigen::VectorXd stepzch = stepCurvedHessSolver.solve(ach*M*zstep);
+  
+  //Display functions
+  igl::opengl::glfw::Viewer viewer2;
+  viewer2.data().set_mesh(V,F);
+  viewer2.data().show_lines = false;
+  viewer2.callback_key_down =
+  [&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
+  {
+    //Graduate result to show isolines, then compute color matrix
+    const Eigen::VectorXd* z;
+    switch(key) {
+      case '1':
+        z = &zstep;
+        break;
+      case '2':
+        z = &stepzl;
+        break;
+      case '3':
+        z = &stepzh;
+        break;
+      case '4':
+        z = &stepzch;
+        break;
+      default:
+        return false;
+    }
+    viewer.data().set_data(*z);
+    return true;
+  };
+  std::cout << R"(Smoothing a step function.
+Usage:
+1  Show step function
+2  Biharmonic smoothing (zero Neumann boundary)
+3  Biharmonic smoothing (natural planar Hessian boundary)
+4  Biharmonic smoothing (natural curved Hessian boundary)
 
-    return 0;
+)";
+  
+  viewer2.data().set_colormap(CM);
+  viewer2.data().set_data(zstep);
+  viewer2.launch();
+  
+  return 0;
 }

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/721_VectorFieldSmoothing/CMakeLists.txt

@@ -0,0 +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)

tutorial/721_VectorFieldSmoothing/main.cpp

@@ -0,0 +1,147 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/parula.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/per_face_normals.h>
+#include <igl/orient_halfedges.h>
+#include <igl/cr_vector_laplacian.h>
+#include <igl/cr_vector_mass.h>
+#include <igl/edge_midpoints.h>
+#include <igl/edge_vectors.h>
+#include <igl/average_from_edges_onto_vertices.h>
+#include <igl/PI.h>
+
+#include <Eigen/Core>
+#include <Eigen/SparseCholesky>
+#include <Eigen/Geometry>
+
+#include <iostream>
+#include <set>
+#include <limits>
+#include <stdlib.h>
+
+#include "tutorial_shared_path.h"
+
+
+
+int main(int argc, char * argv[])
+{
+  typedef Eigen::SparseMatrix<double> SparseMat;
+  
+  //Constants used for smoothing
+  const double howMuchToSmoothBy = 1e-1;
+  const int howManySmoothingInterations = 50;
+  
+  //Read our mesh
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  if(!igl::read_triangle_mesh
+     (argc>1?argv[1]: TUTORIAL_SHARED_PATH "/elephant.obj",V,F)) {
+    std::cout << "Failed to load mesh." << std::endl;
+  }
+  
+  //Orient edges for plotting
+  Eigen::MatrixXi E, oE;
+  igl::orient_halfedges(F, E, oE);
+  
+  //Compute edge midpoints & edge vectors
+  Eigen::MatrixXd edgeMps, parVec, perpVec;
+  igl::edge_midpoints(V, F, E, oE, edgeMps);
+  igl::edge_vectors(V, F, E, oE, parVec, perpVec);
+  
+  //Constructing a function to add noise to
+  const auto zraw_function = [] (const Eigen::Vector3d& x) {
+    return Eigen::Vector3d(0.2*x(1) + cos(2*x(1)+0.2),
+                           0.5*x(0) + 0.15,
+                           0.3*cos(0.2+igl::PI*x(2)));
+  };
+  
+  Eigen::VectorXd zraw(2*edgeMps.rows());
+  for(int i=0; i<edgeMps.rows(); ++i) {
+    const Eigen::Vector3d f = zraw_function(edgeMps.row(i));
+    zraw(i) = f.dot(parVec.row(i));
+    zraw(i+edgeMps.rows()) = f.dot(perpVec.row(i));
+  }
+  
+  //Add noise
+  srand(71);
+  const double l = 15;
+  Eigen::VectorXd znoisy = zraw + l*Eigen::VectorXd::Random(zraw.size());
+  
+  //Denoise function using the vector Dirichlet energy
+  Eigen::VectorXd zsmoothed = znoisy;
+  for(int i=0; i<howManySmoothingInterations; ++i) {
+    //Compute Laplacian and mass matrix
+    SparseMat L, M;
+    igl::cr_vector_mass(V, F, E, oE, M);
+    igl::cr_vector_laplacian(V, F, E, oE, L);
+    
+    //Implicit step
+    Eigen::SimplicialLDLT<SparseMat> rhsSolver(M + howMuchToSmoothBy*L);
+    zsmoothed = rhsSolver.solve(M*zsmoothed);
+  }
+  
+  //Convert vector fields for plotting
+  const auto cr_result_to_vecs_and_colors = [&]
+  (const Eigen::VectorXd& z, Eigen::MatrixXd& vecs, Eigen::MatrixXd& colors) {
+    vecs.resize(edgeMps.rows(), 3);
+    for(int i=0; i<edgeMps.rows(); ++i) {
+      vecs.row(i) = z(i)*parVec.row(i)
+      + z(i+edgeMps.rows())*perpVec.row(i);
+    }
+    igl::average_from_edges_onto_vertices
+    (F, E, oE, vecs.rowwise().norm(), colors);
+  };
+  Eigen::MatrixXd noisyvecs, noisycolors, smoothedvecs, smoothedcolors,
+  rawvecs, rawcolors;
+  cr_result_to_vecs_and_colors(znoisy, noisyvecs, noisycolors);
+  cr_result_to_vecs_and_colors(zsmoothed, smoothedvecs, smoothedcolors);
+  cr_result_to_vecs_and_colors(zraw, rawvecs, rawcolors);
+  
+  
+  //Viewer that shows noisy and denoised functions
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(V,F);
+  viewer.data().show_lines = false;
+  viewer.callback_key_down =
+  [&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
+  {
+    const Eigen::MatrixXd *vecs, *colors;
+    switch(key) {
+      case '1':
+        vecs = &rawvecs;
+        colors = &rawcolors;
+        break;
+      case '2':
+        vecs = &noisyvecs;
+        colors = &noisycolors;
+        break;
+      case '3':
+        vecs = &smoothedvecs;
+        colors = &smoothedcolors;
+        break;
+      default:
+        return false;
+    }
+    viewer.data().set_data(*colors);
+    viewer.data().clear_edges();
+    const double s = 0.08; //How much to scale vectors during plotting
+    viewer.data().add_edges(edgeMps, edgeMps + s*(*vecs),
+                            Eigen::RowVector3d(0.1, 0.1, 0.1));
+    return true;
+  };
+  std::cout << R"(Usage:
+1  Show raw function
+2  Show noisy function
+3  Show smoothed function
+
+)";
+  Eigen::MatrixXd CM;
+  igl::parula(Eigen::VectorXd::LinSpaced
+              (500,znoisy.minCoeff(), znoisy.maxCoeff()).eval(), true, CM);
+  viewer.data().set_colormap(CM);
+  viewer.callback_key_down(viewer, '1', 0);
+  viewer.launch();
+  
+  return 0;
+}

tutorial/722_VectorParallelTransport/CMakeLists.txt

@@ -0,0 +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)

tutorial/722_VectorParallelTransport/main.cpp

@@ -0,0 +1,150 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/parula.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/per_face_normals.h>
+#include <igl/orient_halfedges.h>
+#include <igl/cr_vector_laplacian.h>
+#include <igl/cr_vector_mass.h>
+#include <igl/crouzeix_raviart_cotmatrix.h>
+#include <igl/crouzeix_raviart_massmatrix.h>
+#include <igl/edge_midpoints.h>
+#include <igl/edge_vectors.h>
+#include <igl/average_from_edges_onto_vertices.h>
+#include <igl/min_quad_with_fixed.h>
+#include <igl/heat_geodesics.h>
+
+#include <Eigen/Core>
+#include <Eigen/SparseCholesky>
+#include <Eigen/Geometry>
+
+#include <iostream>
+#include <set>
+#include <limits>
+#include <stdlib.h>
+
+#include "tutorial_shared_path.h"
+
+
+
+int main(int argc, char * argv[])
+{
+  typedef Eigen::SparseMatrix<double> SparseMat;
+  typedef Eigen::Matrix<double, 1, 1> Vector1d;
+  typedef Eigen::Matrix<int, 1, 1> Vector1i;
+  
+  //Constants used for smoothing
+  const double howMuchToSmoothBy = 1e-1;
+  const int howManySmoothingInterations = 50;
+  
+  //Read our mesh
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  if(!igl::read_triangle_mesh
+     (argc>1?argv[1]: TUTORIAL_SHARED_PATH "/cheburashka.off",V,F)) {
+    std::cout << "Failed to load mesh." << std::endl;
+  }
+  
+  //Compute vector Laplacian and mass matrix
+  Eigen::MatrixXi E, oE;//Compute Laplacian and mass matrix
+  SparseMat vecL, vecM;
+  igl::cr_vector_mass(V, F, E, oE, vecM);
+  igl::cr_vector_laplacian(V, F, E, oE, vecL);
+  const int m = vecL.rows()/2; //The number of edges in the mesh
+  
+  //Convert the E / oE matrix format to list of edges / EMAP format required
+  // by the functions constructing scalar Crouzeix-Raviart functions
+  Eigen::MatrixXi Elist(m,2), EMAP(3*F.rows(),1);
+  for(int i=0; i<F.rows(); ++i) {
+    for(int j=0; j<3; ++j) {
+      const int e = E(i,j);
+      EMAP(i+j*F.rows()) = e;
+      if(oE(i,j)>0) {
+        Elist.row(e) << F(i, (j+1)%3), F(i, (j+2)%3);
+      }
+    }
+  }
+  SparseMat scalarL, scalarM;
+  igl::crouzeix_raviart_massmatrix(V, F, Elist, EMAP, scalarM);
+  igl::crouzeix_raviart_cotmatrix(V, F, Elist, EMAP, scalarL);
+  
+  //Compute edge midpoints & edge vectors
+  Eigen::MatrixXd edgeMps, parVec, perpVec;
+  igl::edge_midpoints(V, F, E, oE, edgeMps);
+  igl::edge_vectors(V, F, E, oE, parVec, perpVec);
+  
+  //Perform the vector heat method
+  const int initialIndex = 14319;
+  const double initialPara=0.95, initialPerp=0.08;
+  const double t = 0.01;
+  
+  SparseMat Aeq;
+  Eigen::VectorXd Beq;
+  Eigen::VectorXi known = Eigen::Vector2i(initialIndex, initialIndex+m);
+  Eigen::VectorXd knownVals = Eigen::Vector2d(initialPara, initialPerp);
+  Eigen::VectorXd Y0 = Eigen::VectorXd::Zero(2*m), Yt;
+  Y0(initialIndex) = initialPara; Y0(initialIndex+m) = initialPerp;
+  igl::min_quad_with_fixed
+  (SparseMat(vecM+t*vecL), Eigen::VectorXd(-vecM*Y0), known, knownVals,
+   Aeq, Beq, false, Yt);
+  
+  Eigen::VectorXd u0 = Eigen::VectorXd::Zero(m), ut;
+  u0(initialIndex) = sqrt(initialPara*initialPara + initialPerp*initialPerp);
+  Eigen::VectorXi knownScal = Vector1i(initialIndex);
+  Eigen::VectorXd knownScalVals = Vector1d(u0(initialIndex));
+  igl::min_quad_with_fixed
+  (SparseMat(scalarM+t*scalarL), Eigen::VectorXd(-scalarM*u0), knownScal,
+   knownScalVals, Aeq, Beq, false, ut);
+  
+  Eigen::VectorXd phi0 = Eigen::VectorXd::Zero(m), phit;
+  phi0(initialIndex) = 1;
+  Eigen::VectorXd knownScalValsPhi = Vector1d(1);
+  igl::min_quad_with_fixed
+  (SparseMat(scalarM+t*scalarL), Eigen::VectorXd(-scalarM*phi0), knownScal,
+   knownScalValsPhi, Aeq, Beq, false, phit);
+  
+  Eigen::ArrayXd Xtfactor = ut.array() /
+  (phit.array() * (Yt.array().segment(0,m)*Yt.array().segment(0,m)
+                   + Yt.array().segment(m,m)*Yt.array().segment(m,m)).sqrt());
+  Eigen::VectorXd Xt(2*m);
+  Xt.segment(0,m) = Xtfactor * Yt.segment(0,m).array();
+  Xt.segment(m,m) = Xtfactor * Yt.segment(m,m).array();
+  
+  
+  //Compute scalar heat colors
+  igl::HeatGeodesicsData<double> hgData;
+  igl::heat_geodesics_precompute(V, F, hgData);
+  Eigen::VectorXd heatColor;
+  Eigen::VectorXi gamma = Elist.row(initialIndex);
+  igl::heat_geodesics_solve(hgData, gamma, heatColor);
+  
+  
+  //Convert vector field for plotting
+  Eigen::MatrixXd vecs(m, 3);
+  for(int i=0; i<edgeMps.rows(); ++i) {
+    vecs.row(i) = Xt(i)*parVec.row(i) + Xt(i+edgeMps.rows())*perpVec.row(i);
+  }
+  
+  
+  //Viewer that shows parallel transported vector
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(V,F);
+  viewer.data().show_lines = false;
+  viewer.data().set_data(heatColor.maxCoeff()-heatColor.array(), //invert colormap
+                         igl::COLOR_MAP_TYPE_VIRIDIS);
+  const double s = 0.012; //How much to scale vectors during plotting
+  Eigen::MatrixXd vecColors(m, 3);
+  for(int i=0; i<m; ++i) {
+    vecColors.row(i) << 0.1, 0.1, 0.1;
+  }
+  vecColors.row(initialIndex) << 0.9, 0.1, 0.1;
+  viewer.data().add_edges(edgeMps, edgeMps + s*vecs, vecColors);
+  
+  std::cout << R"(The red vector is parallel transported to every point on the surface.
+The surface is shaded by geodesic distance from the red vector.
+)"
+  << std::endl;
+  viewer.launch();
+  
+  return 0;
+}

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
@@ -166,4 +169,6 @@ if(TUTORIALS_CHAPTER7)
   add_subdirectory("718_IterativeClosestPoint")
   add_subdirectory("719_ExplodedView")
   add_subdirectory("720_BlueNoise")
+  add_subdirectory("721_VectorFieldSmoothing")
+  add_subdirectory("722_VectorParallelTransport")
 endif()