sizing field + bg mesh added (6fae0833) · Commits · Intelligent 3D Graphics / Fast Tetrahedral Meshing in the Wild

src/MeshImprovement.cpp

+165 −2

Original line number	Diff line number	Diff line
		@@ -241,6 +241,12 @@ void floatTetWild::optimization(const std::vector<Vector3> &input_vertices, cons
		v.sizing_scalar = 1;
		}
		operation(input_vertices, input_faces, input_tags, is_face_inserted, mesh, tree, std::array<int, 5>({{0, 1, 0, 0, 0}}));


		///apply sizing field
		if(mesh.params.background_mesh != ""){
		apply_sizingfield(mesh, tree);
		}
		}

		void floatTetWild::cleanup_empty_slots(Mesh &mesh, double percentage) {
		@@ -297,8 +303,7 @@ void floatTetWild::cleanup_empty_slots(Mesh &mesh, double percentage) {
		cout<<mesh.tets.size()<<endl;
		}

		void floatTetWild::operation(const std::vector<Vector3> &input_vertices, const std::vector<Vector3i> &input_faces, const std::vector<int> &input_tags, std::vector<bool> &is_face_inserted,
		Mesh &mesh, AABBWrapper& tree, const std::array<int, 5> &ops) {
		void floatTetWild::operation(Mesh &mesh, AABBWrapper& tree, const std::array<int, 4> &ops){
		igl::Timer igl_timer;
		int v_num, t_num;
		double max_energy, avg_energy;
		@@ -378,6 +383,91 @@ void floatTetWild::operation(const std::vector<Vector3> &input_vertices, const s
		stats().record(StateInfo::smoothing_id, time, v_num, t_num, max_energy, avg_energy);
		output_info(mesh, tree);
		}
		}

		void floatTetWild::operation(const std::vector<Vector3> &input_vertices, const std::vector<Vector3i> &input_faces, const std::vector<int> &input_tags, std::vector<bool> &is_face_inserted,
		Mesh &mesh, AABBWrapper& tree, const std::array<int, 5> &ops) {
		operation(mesh, tree, {{ops[0], ops[1], ops[2], ops[3]}});

		igl::Timer igl_timer;
		// int v_num, t_num;
		// double max_energy, avg_energy;
		// double time;
		//
		// for (int i = 0; i < ops[0]; i++) {
		// igl_timer.start();
		// cout << "edge splitting..." << endl;
		// untangle(mesh);
		// edge_splitting(mesh, tree);
		// time = igl_timer.getElapsedTime();
		// cout << "edge splitting done!" << endl;
		// cout << "time = " << time << "s" << endl;
		// v_num = mesh.get_v_num();
		// t_num = mesh.get_t_num();
		// get_max_avg_energy(mesh, max_energy, avg_energy);
		// cout << "#v = " << v_num << endl;
		// cout << "#t = " << t_num << endl;
		// cout << "max_energy = " << max_energy << endl;
		// cout << "avg_energy = " << avg_energy << endl;
		// stats().record(StateInfo::splitting_id, time, v_num, t_num, max_energy, avg_energy);
		// output_info(mesh, tree);
		// }
		//
		// for (int i = 0; i < ops[1]; i++) {
		// igl_timer.start();
		// cout << "edge collapsing..." << endl;
		// untangle(mesh);
		// edge_collapsing(mesh, tree);
		// time = igl_timer.getElapsedTime();
		// cout << "edge collapsing done!" << endl;
		// cout << "time = " << time << "s" << endl;
		// v_num = mesh.get_v_num();
		// t_num = mesh.get_t_num();
		// get_max_avg_energy(mesh, max_energy, avg_energy);
		// cout << "#v = " << v_num << endl;
		// cout << "#t = " << t_num << endl;
		// cout << "max_energy = " << max_energy << endl;
		// cout << "avg_energy = " << avg_energy << endl;
		// stats().record(StateInfo::collapsing_id, time, v_num, t_num, max_energy, avg_energy);
		// output_info(mesh, tree);
		// }
		//
		// for (int i = 0; i < ops[2]; i++) {
		// igl_timer.start();
		// cout << "edge swapping..." << endl;
		// untangle(mesh);
		// edge_swapping(mesh);
		// time = igl_timer.getElapsedTime();
		// cout << "edge swapping done!" << endl;
		// cout << "time = " << time << "s" << endl;
		// v_num = mesh.get_v_num();
		// t_num = mesh.get_t_num();
		// get_max_avg_energy(mesh, max_energy, avg_energy);
		// cout << "#v = " << v_num << endl;
		// cout << "#t = " << t_num << endl;
		// cout << "max_energy = " << max_energy << endl;
		// cout << "avg_energy = " << avg_energy << endl;
		// stats().record(StateInfo::swapping_id, time, v_num, t_num, max_energy, avg_energy);
		// output_info(mesh, tree);
		// }
		//
		// for (int i = 0; i < ops[3]; i++) {
		// igl_timer.start();
		// cout << "vertex smoothing..." << endl;
		// vertex_smoothing(mesh, tree);
		// time = igl_timer.getElapsedTime();
		// cout << "vertex smoothing done!" << endl;
		// cout << "time = " << time << "s" << endl;
		// v_num = mesh.get_v_num();
		// t_num = mesh.get_t_num();
		// get_max_avg_energy(mesh, max_energy, avg_energy);
		// cout << "#v = " << v_num << endl;
		// cout << "#t = " << t_num << endl;
		// cout << "max_energy = " << max_energy << endl;
		// cout << "avg_energy = " << avg_energy << endl;
		// stats().record(StateInfo::smoothing_id, time, v_num, t_num, max_energy, avg_energy);
		// output_info(mesh, tree);
		// }

		if (!mesh.is_input_all_inserted) {
		pausee();
		@@ -1102,6 +1192,79 @@ void floatTetWild::output_surface(Mesh& mesh, const std::string& filename) {
		igl::writeSTL(filename + ".stl", Eigen::MatrixXd(V_sf), Eigen::MatrixXi(F_sf));
		}

		#include <floattetwild/MshLoader.h>
		#include <geogram/mesh/mesh_AABB.h>
		void floatTetWild::apply_sizingfield(Mesh& mesh, AABBWrapper& tree) {
		auto &tet_vertices = mesh.tet_vertices;
		auto &tets = mesh.tets;

		PyMesh::MshLoader mshLoader(mesh.params.background_mesh);
		Eigen::VectorXd V_in = mshLoader.get_nodes();
		Eigen::VectorXi T_in = mshLoader.get_elements();
		Eigen::VectorXd values = mshLoader.get_node_field("values");
		if (V_in.rows() == 0 \|\| T_in.rows() == 0 \|\| values.rows() == 0)
		return;

		logger().debug("Applying sizing field...");

		GEO::Mesh bg_mesh;
		bg_mesh.vertices.clear();
		bg_mesh.vertices.create_vertices((int) V_in.rows() / 3);
		for (int i = 0; i < V_in.rows() / 3; i++) {
		GEO::vec3 &p = bg_mesh.vertices.point(i);
		for (int j = 0; j < 3; j++)
		p[j] = V_in(i * 3 + j);
		}
		bg_mesh.cells.clear();
		bg_mesh.cells.create_tets((int) T_in.rows() / 4);
		for (int i = 0; i < T_in.rows() / 4; i++) {
		for (int j = 0; j < 4; j++)
		bg_mesh.cells.set_vertex(i, j, T_in(i * 4 + j));
		}

		GEO::MeshCellsAABB bg_aabb(bg_mesh, false);
		for (auto &p: tet_vertices) {
		if (p.is_removed)
		continue;

		p.sizing_scalar = 1;//reset scalar
		GEO::vec3 geo_p(p.pos[0], p.pos[1], p.pos[2]);
		int bg_t_id = bg_aabb.containing_tet(geo_p);
		if (bg_t_id == GEO::MeshCellsAABB::NO_TET)
		continue;

		//compute barycenter
		std::array<Vector3, 4> vs;
		for (int j = 0; j < 4; j++) {
		vs[j] = Vector3(V_in(T_in(bg_t_id * 4 + j) * 3), V_in(T_in(bg_t_id * 4 + j) * 3 + 1),
		V_in(T_in(bg_t_id * 4 + j) * 3 + 2));
		}
		double value = 0;
		for (int j = 0; j < 4; j++) {
		Vector3 n = ((vs[(j + 1) % 4] - vs[j]).cross(vs[(j + 2) % 4] - vs[j])).normalized();
		double d = (vs[(j + 3) % 4] - vs[j]).dot(n);
		if(d == 0)
		continue;
		double weight = abs((p.pos - vs[j]).dot(n) / d);
		value += weight * values(T_in(bg_t_id * 4 + (j + 3) % 4));
		}
		p.sizing_scalar = value / mesh.params.ideal_edge_length;
		// cout<<p.sizing_scalar<<endl;
		}

		int num_tets = mesh.get_t_num();
		for (int i = 0; i < 20; i++) {
		operation(mesh, tree);
		double tmp_num_tets = mesh.get_t_num();
		double max_energy = mesh.get_max_energy();
		cout<<"/////////"<<i<<" "<<max_energy<<endl;
		if ((tmp_num_tets - num_tets) / num_tets < 0.02
		&& max_energy < mesh.params.stop_energy) //refinement and quality enough
		break;
		num_tets = tmp_num_tets;
		}
		}

		#include <floattetwild/bfs_orient.h>
		#include <igl/unique_rows.h>
		#include <igl/remove_duplicate_vertices.h>

src/MeshImprovement.h

+3 −0

Original line number	Diff line number	Diff line
		@@ -20,6 +20,7 @@ namespace floatTetWild {
		void cleanup_empty_slots(Mesh &mesh, double percentage = 0.7);
		void operation(const std::vector<Vector3> &input_vertices, const std::vector<Vector3i> &input_faces, const std::vector<int> &input_tags, std::vector<bool> &is_face_inserted,
		Mesh &mesh, AABBWrapper& tree, const std::array<int, 5> &ops = {{1, 1, 1, 1, 1}});
		void operation(Mesh &mesh, AABBWrapper& tree, const std::array<int, 4> &ops = {{1, 1, 1, 1}});
		bool update_scaling_field(Mesh &mesh, Scalar max_energy);

		int get_max_p(const Mesh &mesh);
		@@ -37,6 +38,8 @@ namespace floatTetWild {
		void manifold_edges(Mesh& mesh);
		void manifold_vertices(Mesh& mesh);

		void apply_sizingfield(Mesh& mesh, AABBWrapper& tree);

		void output_info(Mesh& mesh, const AABBWrapper& tree);
		void check_envelope(Mesh& mesh, const AABBWrapper& tree);
		void output_surface(Mesh& mesh, const std::string& filename);

src/Parameters.h

+2 −0

Original line number	Diff line number	Diff line
		@@ -26,6 +26,8 @@ class Parameters
		std::string tag_path = "";
		std::string postfix = "";

		std::string background_mesh = "";

		std::string envelope_log = "";
		std::string envelope_log_csv = "";

src/external/CMakeLists.txt

+2 −0

Original line number	Diff line number	Diff line
		@@ -9,6 +9,8 @@ set(SOURCES

		MshSaver.h
		MshSaver.cpp
		MshLoader.h
		MshLoader.cpp
		Exception.h

		triangle_triangle_intersection.cpp

src/external/MshLoader.cpp

0 → 100644

+454 −0

Original line number	Diff line number	Diff line
		/* This file is part of PyMesh. Copyright (c) 2015 by Qingnan Zhou */
		#include "MshLoader.h"

		#include <iostream>
		#include <sstream>
		#include <vector>

		#include "Exception.h"

		using namespace PyMesh;

		MshLoader::MshLoader(const std::string& filename) {
		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;
		throw IOError(err_msg.str());
		}
		// Parse header
		std::string buf;
		double version;
		int type;
		fin >> buf;
		if (buf != "$MeshFormat") { throw 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;
		throw NOT_IMPLEMENTED;
		}
		if (sizeof(int) != 4) {
		std::cerr << "Error: code must be compiled with int size 4 bytes." << std::endl;
		throw NOT_IMPLEMENTED;
		}

		// 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;
		throw NOT_IMPLEMENTED;
		}
		}

		fin >> buf;
		if (buf != "$EndMeshFormat") { throw NOT_IMPLEMENTED; }

		while (!fin.eof()) {
		buf.clear();
		fin >> buf;
		if (buf == "$Nodes") {
		parse_nodes(fin);
		fin >> buf;
		if (buf != "$EndNodes") { throw INVALID_FORMAT; }
		} else if (buf == "$Elements") {
		parse_elements(fin);
		fin >> buf;
		if (buf != "$EndElements") { throw INVALID_FORMAT; }
		} else if (buf == "$NodeData") {
		parse_node_field(fin);
		fin >> buf;
		if (buf != "$EndNodeData") { throw INVALID_FORMAT; }
		} else if (buf == "$ElementData") {
		parse_element_field(fin);
		fin >> buf;
		if (buf != "$EndElementData") { throw INVALID_FORMAT; }
		} else if (fin.eof()) {
		break;
		} else {
		parse_unknown_field(fin, buf);
		}
		}
		fin.close();
		}

		MshLoader::FieldNames MshLoader::get_node_field_names() const {
		FieldNames result;

		for (FieldMap::const_iterator itr = m_node_fields.begin();
		itr != m_node_fields.end(); itr++) {
		result.push_back(itr->first);
		}
		return result;
		}

		MshLoader::FieldNames MshLoader::get_element_field_names() const {
		FieldNames result;

		for (FieldMap::const_iterator itr = m_element_fields.begin();
		itr != m_element_fields.end(); itr++) {
		result.push_back(itr->first);
		}
		return result;
		}

		void MshLoader::parse_nodes(std::ifstream& fin) {
		size_t num_nodes;
		fin >> num_nodes;
		m_nodes.resize(num_nodes*3);

		if (m_binary) {
		size_t num_bytes = (4+3m_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+3m_data_size)]) - 1;
		m_nodes[node_idx3] = reinterpret_cast<Float>(&data[i(4+3*m_data_size) + 4]);
		m_nodes[node_idx3+1] = reinterpret_cast<Float>(&data[i(4+3*m_data_size) + 4 + m_data_size]);
		m_nodes[node_idx3+2] = reinterpret_cast<Float>(&data[i(4+3m_data_size) + 4 + 2m_data_size]);
		}

		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];
		}
		}
		if (!m_nodes.allFinite()) {
		throw IOError("NaN or Inf detected in input file.");
		// return;
		}
		}

		void MshLoader::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) {
		case 2:
		return &triangle_elements;
		case 3:
		return &quad_elements;
		case 4:
		return &tet_elements;
		case 5:
		return &hex_elements;
		default:
		throw IOError("Unsupported element type encountered");
		};
		};

		auto get_element_idx_storage = [&](int elem_type) -> std::vector<int>* {
		switch (elem_type) {
		case 2:
		return &triangle_element_idx;
		case 3:
		return &quad_element_idx;
		case 4:
		return &tet_element_idx;
		case 5:
		return &hex_element_idx;
		default:
		throw IOError("Unsupported element type encountered");
		};
		};

		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);
		}
		}

		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) {
		throw IOError("Parsing elements failed!");
		}
		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;
		}
		}

		void 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::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];
		}
		}

		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) { throw INVALID_FORMAT; }
		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_componentsm_data_size)]);
		node_idx -= 1;
		size_t base_idx = i(4+num_componentsm_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];
		}
		}
		}

		m_node_fields[fieldname] = field;
		}

		void 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::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];
		}
		}

		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) { throw INVALID_FORMAT; }
		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_componentsm_data_size)]);
		elem_idx -= 1;
		size_t base_idx = i(4+num_componentsm_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];
		}
		}
		}

		m_element_fields[fieldname] = field;
		}

		void 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;
		}
		}

		void MshLoader::eat_white_space(std::ifstream& fin) {
		char next = fin.peek();
		while (next == '\n' \|\| next == ' ' \|\| next == '\t' \|\| next == '\r') {
		fin.get();
		next = fin.peek();
		}
		}

		int MshLoader::num_nodes_per_elem_type(int elem_type) {
		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:
		std::cerr << "Warning: Element type (" << elem_type << ") is not supported yet."
		<< std::endl;
		throw NOT_IMPLEMENTED;
		}
		return nodes_per_element;
		}
		No newline at end of file