| //##################################################################### |
| // 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 |