src/parsec/disk-image/parsec/parsec-benchmark/pkgs/apps/facesim/src/Public_Library/Grids/TETRAHEDRON_MESH.h - public/gem5-resources - Git at Google

 //#####################################################################
 // Copyright 2002-2004, Zhaosheng Bao, Robert Bridson, Ron Fedkiw, Geoffrey Irving, Sergey Koltakov, Neil Molino, Joseph Teran.
 // This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
 //#####################################################################
 // Class TETRAHEDRON_MESH
 //#####################################################################
 #ifndef __TETRAHEDRON_MESH__
 #define __TETRAHEDRON_MESH__

 #include "../Arrays/LIST_ARRAYS.h"
 #include "../Matrices_And_Vectors/VECTOR_2D.h"
 namespace PhysBAM
 {

 class SEGMENT_MESH;
 class TRIANGLE_MESH;
 template<class T> class LIST_ARRAY;

 class TETRAHEDRON_MESH
 {
 public:
 	int number_nodes; // number of nodes in the mesh
 	LIST_ARRAYS<int> tetrahedrons; // array of 4 indices for each tetrahedron - tetrahedron(j,i) is j'th index in tetrahedron i
 	LIST_ARRAY<LIST_ARRAY<int> >* neighbor_nodes; // for each node, list of neighboring nodes
 	LIST_ARRAY<LIST_ARRAY<int> >* incident_tetrahedrons; // for each node, list of neighboring tetrahedrons that contain it
 	LIST_ARRAY<LIST_ARRAY<int> >* adjacent_tetrahedrons; // for each tetrahedron, list of (up to 4) adjacent neighboring tetrahedrons
 	LIST_ARRAY<LIST_ARRAY<int> >* neighbor_tetrahedrons; // for each tetrahedron, list of tetrahedrons sharing at least one of its node
 	SEGMENT_MESH* segment_mesh; // segment mesh consisting of all the edges
 	TRIANGLE_MESH* triangle_mesh; // triangle mesh consisting of all the triangles
 	LIST_ARRAYS<int>* tetrahedron_edges; // array of 6 indices for each tetrahedron - edges(j,i) is edge j in tetrahedron i
 	LIST_ARRAYS<int>* tetrahedron_faces; // array of 4 indices for each tetrahedron - faces(j,i) is the face opposite vertex j in tetrahedron i
 	TRIANGLE_MESH* boundary_mesh;
 	LIST_ARRAY<bool>* node_on_boundary; // node_on_boundary(i) is true if node i is on the boundary
 	LIST_ARRAY<int>* boundary_nodes; // array containing indices of all particles on the boundary
 	LIST_ARRAY<LIST_ARRAY<int> >* edge_tetrahedrons; // for each edge, list of tets containing that edge
 	LIST_ARRAYS<int>* triangle_tetrahedrons; // for each face, list of incident tets

 	TETRAHEDRON_MESH() // simplest constructor - null mesh
 		: number_nodes (0), tetrahedrons (4, 0), neighbor_nodes (0), incident_tetrahedrons (0), adjacent_tetrahedrons (0), neighbor_tetrahedrons (0), segment_mesh (0), triangle_mesh (0),
 		  tetrahedron_edges (0), tetrahedron_faces (0), boundary_mesh (0), node_on_boundary (0), boundary_nodes (0), edge_tetrahedrons (0), triangle_tetrahedrons (0)
 	{}

 	TETRAHEDRON_MESH (const int m, const int n, const int p, bool octahedron = true)
 		: number_nodes (0), neighbor_nodes (0), incident_tetrahedrons (0), adjacent_tetrahedrons (0), neighbor_tetrahedrons (0), segment_mesh (0), triangle_mesh (0), tetrahedron_edges (0),
 		  tetrahedron_faces (0), boundary_mesh (0), node_on_boundary (0), boundary_nodes (0), edge_tetrahedrons (0), triangle_tetrahedrons (0)
 	{
 		if (octahedron) Initialize_Octahedron_Mesh (m, n, p);
 		else Initialize_Cube_Mesh (m, n, p);
 	}

 	TETRAHEDRON_MESH (const LIST_ARRAYS<int>& tetrahedron_list)
 		: number_nodes (0), neighbor_nodes (0), incident_tetrahedrons (0), adjacent_tetrahedrons (0), neighbor_tetrahedrons (0), segment_mesh (0), triangle_mesh (0), tetrahedron_edges (0),
 		  tetrahedron_faces (0), boundary_mesh (0), node_on_boundary (0), boundary_nodes (0), edge_tetrahedrons (0), triangle_tetrahedrons (0)
 	{
 		Initialize_Tetrahedron_Mesh (tetrahedron_list);
 	}

 	TETRAHEDRON_MESH (const TETRAHEDRON_MESH& tetrahedron_mesh)
 		: number_nodes (0), tetrahedrons(), neighbor_nodes (0), incident_tetrahedrons (0), adjacent_tetrahedrons (0), neighbor_tetrahedrons (0), segment_mesh (0), triangle_mesh (0),
 		  tetrahedron_edges (0), tetrahedron_faces (0), boundary_mesh (0), node_on_boundary (0), boundary_nodes (0), edge_tetrahedrons (0), triangle_tetrahedrons (0)
 	{
 		Initialize_Tetrahedron_Mesh (tetrahedron_mesh);
 	}

 	~TETRAHEDRON_MESH();

 	void Clean_Up_Memory()
 	{
 		number_nodes = 0;
 		tetrahedrons.Exact_Resize_Array (4, 0);
 		Delete_Auxiliary_Structures();
 	}

 	void Refresh_Auxiliary_Structures()
 	{
 		if (neighbor_nodes) Initialize_Neighbor_Nodes();

 		if (incident_tetrahedrons) Initialize_Incident_Tetrahedrons();

 		if (adjacent_tetrahedrons) Initialize_Adjacent_Tetrahedrons();

 		if (neighbor_tetrahedrons) Initialize_Neighbor_Tetrahedrons();

 		if (segment_mesh) Initialize_Segment_Mesh();

 		if (triangle_mesh) Initialize_Triangle_Mesh();

 		if (tetrahedron_edges) Initialize_Tetrahedron_Edges();

 		if (tetrahedron_faces) Initialize_Tetrahedron_Faces();

 		if (boundary_mesh) Initialize_Boundary_Mesh();

 		if (node_on_boundary) Initialize_Node_On_Boundary();

 		if (boundary_nodes) Initialize_Boundary_Nodes();

 		if (edge_tetrahedrons) Initialize_Edge_Tetrahedrons();

 		if (triangle_tetrahedrons) Initialize_Triangle_Tetrahedrons();
 	}

 	void Initialize_Octahedron_Mesh (const int m, const int n, const int p)
 	{
 		Clean_Up_Memory();
 		number_nodes = m * n * p + (m + 1) * (n + 1) * (p + 1);
 		tetrahedrons.Exact_Resize_Array (4, 4 * (m - 1) *n * p + 4 * m * (n - 1) *p + 4 * m * n * (p - 1));
 		int t = 0, i, j, k;

 		for (i = 1; i <= m; i++) for (j = 1; j <= n; j++) for (k = 1; k <= p - 1; k++) // loop over k-oriented edges in inner cube mesh
 				{
 					tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j, k + 1, m, n, p), Half_Lattice (i - 1, j - 1, k, m, n, p), Half_Lattice (i, j - 1, k, m, n, p));
 					tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j, k + 1, m, n, p), Half_Lattice (i, j - 1, k, m, n, p), Half_Lattice (i, j, k, m, n, p));
 					tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j, k + 1, m, n, p), Half_Lattice (i, j, k, m, n, p), Half_Lattice (i - 1, j, k, m, n, p));
 					tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j, k + 1, m, n, p), Half_Lattice (i - 1, j, k, m, n, p), Half_Lattice (i - 1, j - 1, k, m, n, p));
 				}

 		for (i = 1; i <= m; i++) for (j = 1; j <= n - 1; j++) for (k = 1; k <= p; k++) // loop over j-oriented edge in inner cube mesh
 				{
 					tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j + 1, k, m, n, p), Half_Lattice (i, j, k - 1, m, n, p), Half_Lattice (i - 1, j, k - 1, m, n, p));
 					tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j + 1, k, m, n, p), Half_Lattice (i, j, k, m, n, p), Half_Lattice (i, j, k - 1, m, n, p));
 					tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j + 1, k, m, n, p), Half_Lattice (i - 1, j, k, m, n, p), Half_Lattice (i, j, k, m, n, p));
 					tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j + 1, k, m, n, p), Half_Lattice (i - 1, j, k - 1, m, n, p), Half_Lattice (i - 1, j, k, m, n, p));
 				}

 		for (i = 1; i <= m - 1; i++) for (j = 1; j <= n; j++) for (k = 1; k <= p; k++) // loop over i-oriented edge in inner cube mesh
 				{
 					tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i + 1, j, k, m, n, p), Half_Lattice (i, j - 1, k - 1, m, n, p), Half_Lattice (i, j, k - 1, m, n, p));
 					tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i + 1, j, k, m, n, p), Half_Lattice (i, j, k - 1, m, n, p), Half_Lattice (i, j, k, m, n, p));
 					tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i + 1, j, k, m, n, p), Half_Lattice (i, j, k, m, n, p), Half_Lattice (i, j - 1, k, m, n, p));
 					tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i + 1, j, k, m, n, p), Half_Lattice (i, j - 1, k, m, n, p), Half_Lattice (i, j - 1, k - 1, m, n, p));
 				}
 	}

 	void Initialize_Cube_Mesh (const int m, const int n, const int p) // 5 tetrahedrons per cube
 	{
 		Clean_Up_Memory();
 		number_nodes = m * n * p;
 		tetrahedrons.Exact_Resize_Array (4, 5 * (m - 1) * (n - 1) * (p - 1));
 		int t = 0;

 		for (int i = 1; i <= m - 1; i++) for (int j = 1; j <= n - 1; j++) for (int k = 1; k <= p - 1; k++)
 				{
 					if ( (i + j + k) % 2 == 0)
 					{
 						tetrahedrons.Set (++t, i + m * (j - 1) + m * n * (k - 1), i + 1 + m * (j - 1) + m * n * (k - 1), i + m * j + m * n * (k - 1), i + m * (j - 1) + m * n * k);
 						tetrahedrons.Set (++t, i + 1 + m * (j - 1) + m * n * (k - 1), i + 1 + m * (j - 1) + m * n * k, i + 1 + m * j + m * n * k, i + m * (j - 1) + m * n * k);
 						tetrahedrons.Set (++t, i + m * j + m * n * (k - 1), i + 1 + m * j + m * n * (k - 1), i + 1 + m * j + m * n * k, i + 1 + m * (j - 1) + m * n * (k - 1));
 						tetrahedrons.Set (++t, i + m * j + m * n * k, i + 1 + m * j + m * n * k, i + m * (j - 1) + m * n * k, i + m * j + m * n * (k - 1));
 						tetrahedrons.Set (++t, i + 1 + m * (j - 1) + m * n * (k - 1), i + m * (j - 1) + m * n * k, i + 1 + m * j + m * n * k, i + m * j + m * n * (k - 1));
 					}
 					else
 					{
 						tetrahedrons.Set (++t, i + m * (j - 1) + m * n * (k - 1), i + 1 + m * (j - 1) + m * n * (k - 1), i + 1 + m * j + m * n * (k - 1), i + 1 + m * (j - 1) + m * n * k);
 						tetrahedrons.Set (++t, i + m * (j - 1) + m * n * (k - 1), i + m * j + m * n * (k - 1), i + m * j + m * n * k, i + 1 + m * j + m * n * (k - 1));
 						tetrahedrons.Set (++t, i + m * j + m * n * k, i + 1 + m * (j - 1) + m * n * k, i + m * (j - 1) + m * n * k, i + m * (j - 1) + m * n * (k - 1));
 						tetrahedrons.Set (++t, i + m * j + m * n * k, i + 1 + m * j + m * n * k, i + 1 + m * (j - 1) + m * n * k, i + 1 + m * j + m * n * (k - 1));
 						tetrahedrons.Set (++t, i + m * j + m * n * k, i + m * (j - 1) + m * n * (k - 1), i + 1 + m * j + m * n * (k - 1), i + 1 + m * (j - 1) + m * n * k);
 					}
 				}
 	}

 	void Initialize_Tetrahedron_Mesh (const LIST_ARRAYS<int>& tetrahedron_list) // construct a mesh given a list of tetrahedrons
 	{
 		Clean_Up_Memory();
 		tetrahedrons.Exact_Resize_Array (4, tetrahedron_list.m);
 		LIST_ARRAYS<int>::copy (tetrahedron_list, tetrahedrons);

 		for (int t = 1; t <= tetrahedrons.m; t++)
 		{
 			int i, j, k, l;
 			tetrahedrons.Get (t, i, j, k, l);
 			number_nodes = max (number_nodes, i, j, k, l);
 		}
 	}

 	void Initialize_Tetrahedron_Mesh (const TETRAHEDRON_MESH& tetrahedron_mesh) // works with the copy constructor
 	{
 		Initialize_Tetrahedron_Mesh (tetrahedron_mesh.tetrahedrons);
 	}

 	bool Node_In_Tetrahedron (const int node, const int tetrahedron) const
 	{
 		int i, j, k, l;
 		tetrahedrons.Get (tetrahedron, i, j, k, l);
 		return i == node || j == node || k == node || l == node;
 	}

 	bool Triangle_In_Tetrahedron (const int node1, const int node2, const int node3, const int tetrahedron) const
 	{
 		return Node_In_Tetrahedron (node1, tetrahedron) && Node_In_Tetrahedron (node2, tetrahedron) && Node_In_Tetrahedron (node3, tetrahedron);
 	}

 	bool Face_Neighbors (const int tetrahedron1, const int tetrahedron2) const
 	{
 		int i, j, k, l;
 		tetrahedrons.Get (tetrahedron1, i, j, k, l);
 		return Triangle_In_Tetrahedron (j, k, l, tetrahedron2) || Triangle_In_Tetrahedron (i, k, l, tetrahedron2) ||
 		       Triangle_In_Tetrahedron (i, j, l, tetrahedron2) || Triangle_In_Tetrahedron (i, j, k, tetrahedron2);
 	}

 	void Other_Three_Nodes (const int node, const int tetrahedron, int& other_node1, int& other_node2, int& other_node3) const
 	{
 		assert (Node_In_Tetrahedron (node, tetrahedron));
 		int i, j, k, l;
 		tetrahedrons.Get (tetrahedron, i, j, k, l);

 		if (i == node)
 		{
 			other_node1 = j;
 			other_node2 = k;
 			other_node3 = l;
 		}
 		else if (j == node)
 		{
 			other_node1 = i;
 			other_node2 = l;
 			other_node3 = k;
 		}
 		else if (k == node)
 		{
 			other_node1 = l;
 			other_node2 = i;
 			other_node3 = j;
 		}
 		else if (l == node)
 		{
 			other_node1 = k;
 			other_node2 = j;
 			other_node3 = i;
 		}
 	}

 	void Other_Two_Nodes (const int node1, const int node2, const int tetrahedron, int& other_node1, int& other_node2) const
 	{
 		assert (Node_In_Tetrahedron (node1, tetrahedron) && Node_In_Tetrahedron (node2, tetrahedron));
 		int i, j, k, l;
 		tetrahedrons.Get (tetrahedron, i, j, k, l);

 		if (i == node1)
 		{
 			if (j == node2)
 			{
 				other_node1 = k;
 				other_node2 = l;
 			}
 			else if (k == node2)
 			{
 				other_node1 = l;
 				other_node2 = j;
 			}
 			else
 			{
 				other_node1 = j;
 				other_node2 = k;
 			}
 		}
 		else if (j == node1)
 		{
 			if (i == node2)
 			{
 				other_node1 = l;
 				other_node2 = k;
 			}
 			else if (k == node2)
 			{
 				other_node1 = i;
 				other_node2 = l;
 			}
 			else
 			{
 				other_node1 = k;
 				other_node2 = i;
 			}
 		}
 		else if (k == node1)
 		{
 			if (i == node2)
 			{
 				other_node1 = j;
 				other_node2 = l;
 			}
 			else if (l == node2)
 			{
 				other_node1 = i;
 				other_node2 = j;
 			}
 			else
 			{
 				other_node1 = l;
 				other_node2 = i;
 			}
 		}
 		else if (l == node1)
 		{
 			if (i == node2)
 			{
 				other_node1 = k;
 				other_node2 = j;
 			}
 			else if (k == node2)
 			{
 				other_node1 = j;
 				other_node2 = i;
 			}
 			else
 			{
 				other_node1 = i;
 				other_node2 = k;
 			}
 		}
 	}

 private:
 	int Lattice (const int i, const int j, const int k, const int m, const int n, const int p) const
 	{
 		return i + m * (j - 1) + m * n * (k - 1);
 	}

 	int Half_Lattice (const int i, const int j, const int k, const int m, const int n, const int p) const
 	{
 		return m * n * p + i + 1 + (m + 1) * j + (m + 1) * (n + 1) * k;
 	}

 public:
 	template <class RW>
 	void Read (std::istream& input_stream)
 	{
 		Clean_Up_Memory();
 		Read_Binary<RW> (input_stream, number_nodes);
 		tetrahedrons.template Read<RW> (input_stream);
 	}

 	template <class RW>
 	void Write (std::ostream& output_stream) const
 	{
 		Write_Binary<RW> (output_stream, number_nodes);
 		tetrahedrons.template Write<RW> (output_stream);
 	}

 //#####################################################################
 	void Delete_Auxiliary_Structures();
 	int Tetrahedron (const int node1, const int node2, const int node3, const int node4) const;
 	void Initialize_Neighbor_Nodes();
 	void Initialize_Incident_Tetrahedrons();
 	void Initialize_Adjacent_Tetrahedrons();
 private: // helper function for Initialize_Adjacent_Tetrahedrons
 	void Find_And_Append_Adjacent_Tetrahedrons (const int tetrahedron, const int node1, const int node2, const int node3);
 	int Number_Of_Tetrahedrons_Across_Face (const int tetrahedron, const int node1, const int node2, const int node3) const;
 public:
 	void Initialize_Neighbor_Tetrahedrons();
 	void Initialize_Segment_Mesh();
 	void Initialize_Triangle_Mesh();
 	void Initialize_Tetrahedron_Edges();
 	void Initialize_Tetrahedron_Faces();
 	void Initialize_Boundary_Mesh();
 	void Initialize_Node_On_Boundary();
 	void Initialize_Boundary_Nodes();
 	void Initialize_Edge_Tetrahedrons();
 	void Initialize_Triangle_Tetrahedrons();
 	int Delete_Tetrahedrons_With_Missing_Nodes(); // returns the number deleted
 	void Delete_Tetrahedrons (const LIST_ARRAY<int>& deletion_list);
 	int Number_Of_Boundary_Tetrahedrons();
 	int Number_Of_Interior_Tetrahedrons();
 	int Number_Of_Nodes_With_Minimum_Valence();
 	int Minimum_Valence (int* index = 0);
 	bool Edge_Neighbors (const int tet1, const int tet2) const;
 	int Number_Of_Edge_Neighbors (const int segment) const;
 	void Initialize_Segment_Mesh_Of_Subset (SEGMENT_MESH& segment_mesh_of_subset, const LIST_ARRAY<bool>& subset) const;
 	void Initialize_Boundary_Mesh_Of_Subset (TRIANGLE_MESH& boundary_mesh_of_subset, const LIST_ARRAY<bool>& subset);
 	void Update_Adjacent_Tetrahedrons_From_Incident_Tetrahedrons (const int node);
 	void Update_Neighbor_Nodes_From_Incident_Tetrahedrons (const int node);
 	void Mark_Face_Connected_Component_Incident_On_A_Node (const int node, const int tetrahedron_index_in_incident_tetrahedrons, ARRAY<bool>& marked) const;
 	int Tetrahedrons_On_Edge (const int node1, const int node2, LIST_ARRAY<int>* tetrahedrons_on_edge) const;
 	int Tetrahedrons_Across_Face (const int tetrahedron, const int node1, const int node2, const int node3, LIST_ARRAY<int>& tetrahedrons_across_face) const;
 	void Identify_Face_Connected_Components (ARRAY<int>& label);
 	void Identify_Edge_Connected_Components (ARRAY<int>& label);
 //#####################################################################
 };
 }
 #endif
	//#####################################################################
	// Copyright 2002-2004, Zhaosheng Bao, Robert Bridson, Ron Fedkiw, Geoffrey Irving, Sergey Koltakov, Neil Molino, Joseph Teran.
	// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
	//#####################################################################
	// Class TETRAHEDRON_MESH
	//#####################################################################
	#ifndef __TETRAHEDRON_MESH__
	#define __TETRAHEDRON_MESH__

	#include "../Arrays/LIST_ARRAYS.h"
	#include "../Matrices_And_Vectors/VECTOR_2D.h"
	namespace PhysBAM
	{

	class SEGMENT_MESH;
	class TRIANGLE_MESH;
	template<class T> class LIST_ARRAY;

	class TETRAHEDRON_MESH
	{
	public:
	int number_nodes; // number of nodes in the mesh
	LIST_ARRAYS<int> tetrahedrons; // array of 4 indices for each tetrahedron - tetrahedron(j,i) is j'th index in tetrahedron i
	LIST_ARRAY<LIST_ARRAY<int> >* neighbor_nodes; // for each node, list of neighboring nodes
	LIST_ARRAY<LIST_ARRAY<int> >* incident_tetrahedrons; // for each node, list of neighboring tetrahedrons that contain it
	LIST_ARRAY<LIST_ARRAY<int> >* adjacent_tetrahedrons; // for each tetrahedron, list of (up to 4) adjacent neighboring tetrahedrons
	LIST_ARRAY<LIST_ARRAY<int> >* neighbor_tetrahedrons; // for each tetrahedron, list of tetrahedrons sharing at least one of its node
	SEGMENT_MESH* segment_mesh; // segment mesh consisting of all the edges
	TRIANGLE_MESH* triangle_mesh; // triangle mesh consisting of all the triangles
	LIST_ARRAYS<int>* tetrahedron_edges; // array of 6 indices for each tetrahedron - edges(j,i) is edge j in tetrahedron i
	LIST_ARRAYS<int>* tetrahedron_faces; // array of 4 indices for each tetrahedron - faces(j,i) is the face opposite vertex j in tetrahedron i
	TRIANGLE_MESH* boundary_mesh;
	LIST_ARRAY<bool>* node_on_boundary; // node_on_boundary(i) is true if node i is on the boundary
	LIST_ARRAY<int>* boundary_nodes; // array containing indices of all particles on the boundary
	LIST_ARRAY<LIST_ARRAY<int> >* edge_tetrahedrons; // for each edge, list of tets containing that edge
	LIST_ARRAYS<int>* triangle_tetrahedrons; // for each face, list of incident tets

	TETRAHEDRON_MESH() // simplest constructor - null mesh
	: number_nodes (0), tetrahedrons (4, 0), neighbor_nodes (0), incident_tetrahedrons (0), adjacent_tetrahedrons (0), neighbor_tetrahedrons (0), segment_mesh (0), triangle_mesh (0),
	tetrahedron_edges (0), tetrahedron_faces (0), boundary_mesh (0), node_on_boundary (0), boundary_nodes (0), edge_tetrahedrons (0), triangle_tetrahedrons (0)
	{}

	TETRAHEDRON_MESH (const int m, const int n, const int p, bool octahedron = true)
	: number_nodes (0), neighbor_nodes (0), incident_tetrahedrons (0), adjacent_tetrahedrons (0), neighbor_tetrahedrons (0), segment_mesh (0), triangle_mesh (0), tetrahedron_edges (0),
	tetrahedron_faces (0), boundary_mesh (0), node_on_boundary (0), boundary_nodes (0), edge_tetrahedrons (0), triangle_tetrahedrons (0)
	{
	if (octahedron) Initialize_Octahedron_Mesh (m, n, p);
	else Initialize_Cube_Mesh (m, n, p);
	}

	TETRAHEDRON_MESH (const LIST_ARRAYS<int>& tetrahedron_list)
	: number_nodes (0), neighbor_nodes (0), incident_tetrahedrons (0), adjacent_tetrahedrons (0), neighbor_tetrahedrons (0), segment_mesh (0), triangle_mesh (0), tetrahedron_edges (0),
	tetrahedron_faces (0), boundary_mesh (0), node_on_boundary (0), boundary_nodes (0), edge_tetrahedrons (0), triangle_tetrahedrons (0)
	{
	Initialize_Tetrahedron_Mesh (tetrahedron_list);
	}

	TETRAHEDRON_MESH (const TETRAHEDRON_MESH& tetrahedron_mesh)
	: number_nodes (0), tetrahedrons(), neighbor_nodes (0), incident_tetrahedrons (0), adjacent_tetrahedrons (0), neighbor_tetrahedrons (0), segment_mesh (0), triangle_mesh (0),
	tetrahedron_edges (0), tetrahedron_faces (0), boundary_mesh (0), node_on_boundary (0), boundary_nodes (0), edge_tetrahedrons (0), triangle_tetrahedrons (0)
	{
	Initialize_Tetrahedron_Mesh (tetrahedron_mesh);
	}

	~TETRAHEDRON_MESH();

	void Clean_Up_Memory()
	{
	number_nodes = 0;
	tetrahedrons.Exact_Resize_Array (4, 0);
	Delete_Auxiliary_Structures();
	}

	void Refresh_Auxiliary_Structures()
	{
	if (neighbor_nodes) Initialize_Neighbor_Nodes();

	if (incident_tetrahedrons) Initialize_Incident_Tetrahedrons();

	if (adjacent_tetrahedrons) Initialize_Adjacent_Tetrahedrons();

	if (neighbor_tetrahedrons) Initialize_Neighbor_Tetrahedrons();

	if (segment_mesh) Initialize_Segment_Mesh();

	if (triangle_mesh) Initialize_Triangle_Mesh();

	if (tetrahedron_edges) Initialize_Tetrahedron_Edges();

	if (tetrahedron_faces) Initialize_Tetrahedron_Faces();

	if (boundary_mesh) Initialize_Boundary_Mesh();

	if (node_on_boundary) Initialize_Node_On_Boundary();

	if (boundary_nodes) Initialize_Boundary_Nodes();

	if (edge_tetrahedrons) Initialize_Edge_Tetrahedrons();

	if (triangle_tetrahedrons) Initialize_Triangle_Tetrahedrons();
	}

	void Initialize_Octahedron_Mesh (const int m, const int n, const int p)
	{
	Clean_Up_Memory();
	number_nodes = m * n * p + (m + 1) * (n + 1) * (p + 1);
	tetrahedrons.Exact_Resize_Array (4, 4 * (m - 1) n p + 4 * m * (n - 1) p + 4 m * n * (p - 1));
	int t = 0, i, j, k;

	for (i = 1; i <= m; i++) for (j = 1; j <= n; j++) for (k = 1; k <= p - 1; k++) // loop over k-oriented edges in inner cube mesh
	{
	tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j, k + 1, m, n, p), Half_Lattice (i - 1, j - 1, k, m, n, p), Half_Lattice (i, j - 1, k, m, n, p));
	tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j, k + 1, m, n, p), Half_Lattice (i, j - 1, k, m, n, p), Half_Lattice (i, j, k, m, n, p));
	tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j, k + 1, m, n, p), Half_Lattice (i, j, k, m, n, p), Half_Lattice (i - 1, j, k, m, n, p));
	tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j, k + 1, m, n, p), Half_Lattice (i - 1, j, k, m, n, p), Half_Lattice (i - 1, j - 1, k, m, n, p));
	}

	for (i = 1; i <= m; i++) for (j = 1; j <= n - 1; j++) for (k = 1; k <= p; k++) // loop over j-oriented edge in inner cube mesh
	{
	tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j + 1, k, m, n, p), Half_Lattice (i, j, k - 1, m, n, p), Half_Lattice (i - 1, j, k - 1, m, n, p));
	tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j + 1, k, m, n, p), Half_Lattice (i, j, k, m, n, p), Half_Lattice (i, j, k - 1, m, n, p));
	tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j + 1, k, m, n, p), Half_Lattice (i - 1, j, k, m, n, p), Half_Lattice (i, j, k, m, n, p));
	tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i, j + 1, k, m, n, p), Half_Lattice (i - 1, j, k - 1, m, n, p), Half_Lattice (i - 1, j, k, m, n, p));
	}

	for (i = 1; i <= m - 1; i++) for (j = 1; j <= n; j++) for (k = 1; k <= p; k++) // loop over i-oriented edge in inner cube mesh
	{
	tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i + 1, j, k, m, n, p), Half_Lattice (i, j - 1, k - 1, m, n, p), Half_Lattice (i, j, k - 1, m, n, p));
	tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i + 1, j, k, m, n, p), Half_Lattice (i, j, k - 1, m, n, p), Half_Lattice (i, j, k, m, n, p));
	tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i + 1, j, k, m, n, p), Half_Lattice (i, j, k, m, n, p), Half_Lattice (i, j - 1, k, m, n, p));
	tetrahedrons.Set (++t, Lattice (i, j, k, m, n, p), Lattice (i + 1, j, k, m, n, p), Half_Lattice (i, j - 1, k, m, n, p), Half_Lattice (i, j - 1, k - 1, m, n, p));
	}
	}

	void Initialize_Cube_Mesh (const int m, const int n, const int p) // 5 tetrahedrons per cube
	{
	Clean_Up_Memory();
	number_nodes = m * n * p;
	tetrahedrons.Exact_Resize_Array (4, 5 * (m - 1) * (n - 1) * (p - 1));
	int t = 0;

	for (int i = 1; i <= m - 1; i++) for (int j = 1; j <= n - 1; j++) for (int k = 1; k <= p - 1; k++)
	{
	if ( (i + j + k) % 2 == 0)
	{
	tetrahedrons.Set (++t, i + m * (j - 1) + m * n * (k - 1), i + 1 + m * (j - 1) + m * n * (k - 1), i + m * j + m * n * (k - 1), i + m * (j - 1) + m * n * k);
	tetrahedrons.Set (++t, i + 1 + m * (j - 1) + m * n * (k - 1), i + 1 + m * (j - 1) + m * n * k, i + 1 + m * j + m * n * k, i + m * (j - 1) + m * n * k);
	tetrahedrons.Set (++t, i + m * j + m * n * (k - 1), i + 1 + m * j + m * n * (k - 1), i + 1 + m * j + m * n * k, i + 1 + m * (j - 1) + m * n * (k - 1));
	tetrahedrons.Set (++t, i + m * j + m * n * k, i + 1 + m * j + m * n * k, i + m * (j - 1) + m * n * k, i + m * j + m * n * (k - 1));
	tetrahedrons.Set (++t, i + 1 + m * (j - 1) + m * n * (k - 1), i + m * (j - 1) + m * n * k, i + 1 + m * j + m * n * k, i + m * j + m * n * (k - 1));
	}
	else
	{
	tetrahedrons.Set (++t, i + m * (j - 1) + m * n * (k - 1), i + 1 + m * (j - 1) + m * n * (k - 1), i + 1 + m * j + m * n * (k - 1), i + 1 + m * (j - 1) + m * n * k);
	tetrahedrons.Set (++t, i + m * (j - 1) + m * n * (k - 1), i + m * j + m * n * (k - 1), i + m * j + m * n * k, i + 1 + m * j + m * n * (k - 1));
	tetrahedrons.Set (++t, i + m * j + m * n * k, i + 1 + m * (j - 1) + m * n * k, i + m * (j - 1) + m * n * k, i + m * (j - 1) + m * n * (k - 1));
	tetrahedrons.Set (++t, i + m * j + m * n * k, i + 1 + m * j + m * n * k, i + 1 + m * (j - 1) + m * n * k, i + 1 + m * j + m * n * (k - 1));
	tetrahedrons.Set (++t, i + m * j + m * n * k, i + m * (j - 1) + m * n * (k - 1), i + 1 + m * j + m * n * (k - 1), i + 1 + m * (j - 1) + m * n * k);
	}
	}
	}

	void Initialize_Tetrahedron_Mesh (const LIST_ARRAYS<int>& tetrahedron_list) // construct a mesh given a list of tetrahedrons
	{
	Clean_Up_Memory();
	tetrahedrons.Exact_Resize_Array (4, tetrahedron_list.m);
	LIST_ARRAYS<int>::copy (tetrahedron_list, tetrahedrons);

	for (int t = 1; t <= tetrahedrons.m; t++)
	{
	int i, j, k, l;
	tetrahedrons.Get (t, i, j, k, l);
	number_nodes = max (number_nodes, i, j, k, l);
	}
	}

	void Initialize_Tetrahedron_Mesh (const TETRAHEDRON_MESH& tetrahedron_mesh) // works with the copy constructor
	{
	Initialize_Tetrahedron_Mesh (tetrahedron_mesh.tetrahedrons);
	}

	bool Node_In_Tetrahedron (const int node, const int tetrahedron) const
	{
	int i, j, k, l;
	tetrahedrons.Get (tetrahedron, i, j, k, l);
	return i == node \|\| j == node \|\| k == node \|\| l == node;
	}

	bool Triangle_In_Tetrahedron (const int node1, const int node2, const int node3, const int tetrahedron) const
	{
	return Node_In_Tetrahedron (node1, tetrahedron) && Node_In_Tetrahedron (node2, tetrahedron) && Node_In_Tetrahedron (node3, tetrahedron);
	}

	bool Face_Neighbors (const int tetrahedron1, const int tetrahedron2) const
	{
	int i, j, k, l;
	tetrahedrons.Get (tetrahedron1, i, j, k, l);
	return Triangle_In_Tetrahedron (j, k, l, tetrahedron2) \|\| Triangle_In_Tetrahedron (i, k, l, tetrahedron2) \|\|
	Triangle_In_Tetrahedron (i, j, l, tetrahedron2) \|\| Triangle_In_Tetrahedron (i, j, k, tetrahedron2);
	}

	void Other_Three_Nodes (const int node, const int tetrahedron, int& other_node1, int& other_node2, int& other_node3) const
	{
	assert (Node_In_Tetrahedron (node, tetrahedron));
	int i, j, k, l;
	tetrahedrons.Get (tetrahedron, i, j, k, l);

	if (i == node)
	{
	other_node1 = j;
	other_node2 = k;
	other_node3 = l;
	}
	else if (j == node)
	{
	other_node1 = i;
	other_node2 = l;
	other_node3 = k;
	}
	else if (k == node)
	{
	other_node1 = l;
	other_node2 = i;
	other_node3 = j;
	}
	else if (l == node)
	{
	other_node1 = k;
	other_node2 = j;
	other_node3 = i;
	}
	}

	void Other_Two_Nodes (const int node1, const int node2, const int tetrahedron, int& other_node1, int& other_node2) const
	{
	assert (Node_In_Tetrahedron (node1, tetrahedron) && Node_In_Tetrahedron (node2, tetrahedron));
	int i, j, k, l;
	tetrahedrons.Get (tetrahedron, i, j, k, l);

	if (i == node1)
	{
	if (j == node2)
	{
	other_node1 = k;
	other_node2 = l;
	}
	else if (k == node2)
	{
	other_node1 = l;
	other_node2 = j;
	}
	else
	{
	other_node1 = j;
	other_node2 = k;
	}
	}
	else if (j == node1)
	{
	if (i == node2)
	{
	other_node1 = l;
	other_node2 = k;
	}
	else if (k == node2)
	{
	other_node1 = i;
	other_node2 = l;
	}
	else
	{
	other_node1 = k;
	other_node2 = i;
	}
	}
	else if (k == node1)
	{
	if (i == node2)
	{
	other_node1 = j;
	other_node2 = l;
	}
	else if (l == node2)
	{
	other_node1 = i;
	other_node2 = j;
	}
	else
	{
	other_node1 = l;
	other_node2 = i;
	}
	}
	else if (l == node1)
	{
	if (i == node2)
	{
	other_node1 = k;
	other_node2 = j;
	}
	else if (k == node2)
	{
	other_node1 = j;
	other_node2 = i;
	}
	else
	{
	other_node1 = i;
	other_node2 = k;
	}
	}
	}

	private:
	int Lattice (const int i, const int j, const int k, const int m, const int n, const int p) const
	{
	return i + m * (j - 1) + m * n * (k - 1);
	}

	int Half_Lattice (const int i, const int j, const int k, const int m, const int n, const int p) const
	{
	return m * n * p + i + 1 + (m + 1) * j + (m + 1) * (n + 1) * k;
	}

	public:
	template <class RW>
	void Read (std::istream& input_stream)
	{
	Clean_Up_Memory();
	Read_Binary<RW> (input_stream, number_nodes);
	tetrahedrons.template Read<RW> (input_stream);
	}

	template <class RW>
	void Write (std::ostream& output_stream) const
	{
	Write_Binary<RW> (output_stream, number_nodes);
	tetrahedrons.template Write<RW> (output_stream);
	}

	//#####################################################################
	void Delete_Auxiliary_Structures();
	int Tetrahedron (const int node1, const int node2, const int node3, const int node4) const;
	void Initialize_Neighbor_Nodes();
	void Initialize_Incident_Tetrahedrons();
	void Initialize_Adjacent_Tetrahedrons();
	private: // helper function for Initialize_Adjacent_Tetrahedrons
	void Find_And_Append_Adjacent_Tetrahedrons (const int tetrahedron, const int node1, const int node2, const int node3);
	int Number_Of_Tetrahedrons_Across_Face (const int tetrahedron, const int node1, const int node2, const int node3) const;
	public:
	void Initialize_Neighbor_Tetrahedrons();
	void Initialize_Segment_Mesh();
	void Initialize_Triangle_Mesh();
	void Initialize_Tetrahedron_Edges();
	void Initialize_Tetrahedron_Faces();
	void Initialize_Boundary_Mesh();
	void Initialize_Node_On_Boundary();
	void Initialize_Boundary_Nodes();
	void Initialize_Edge_Tetrahedrons();
	void Initialize_Triangle_Tetrahedrons();
	int Delete_Tetrahedrons_With_Missing_Nodes(); // returns the number deleted
	void Delete_Tetrahedrons (const LIST_ARRAY<int>& deletion_list);
	int Number_Of_Boundary_Tetrahedrons();
	int Number_Of_Interior_Tetrahedrons();
	int Number_Of_Nodes_With_Minimum_Valence();
	int Minimum_Valence (int* index = 0);
	bool Edge_Neighbors (const int tet1, const int tet2) const;
	int Number_Of_Edge_Neighbors (const int segment) const;
	void Initialize_Segment_Mesh_Of_Subset (SEGMENT_MESH& segment_mesh_of_subset, const LIST_ARRAY<bool>& subset) const;
	void Initialize_Boundary_Mesh_Of_Subset (TRIANGLE_MESH& boundary_mesh_of_subset, const LIST_ARRAY<bool>& subset);
	void Update_Adjacent_Tetrahedrons_From_Incident_Tetrahedrons (const int node);
	void Update_Neighbor_Nodes_From_Incident_Tetrahedrons (const int node);
	void Mark_Face_Connected_Component_Incident_On_A_Node (const int node, const int tetrahedron_index_in_incident_tetrahedrons, ARRAY<bool>& marked) const;
	int Tetrahedrons_On_Edge (const int node1, const int node2, LIST_ARRAY<int>* tetrahedrons_on_edge) const;
	int Tetrahedrons_Across_Face (const int tetrahedron, const int node1, const int node2, const int node3, LIST_ARRAY<int>& tetrahedrons_across_face) const;
	void Identify_Face_Connected_Components (ARRAY<int>& label);
	void Identify_Edge_Connected_Components (ARRAY<int>& label);
	//#####################################################################
	};
	}
	#endif