Commit de8259b3 authored by Kirill Terekhov's avatar Kirill Terekhov

Some fixes in edge collapse algorithm, some more fixes will follow

parent 060d0851
Pipeline #168 failed with stages
in 10 minutes and 6 seconds
......@@ -2,8 +2,9 @@ project(GridGen)
add_executable(FixFaults fix_faults.cpp)
add_executable(FixTiny fix_tiny.cpp)
add_executable(FixTinyCollapse fix_tiny_collapse.cpp)
add_executable(UniteFaces unite_faces.cpp)
add_executable(Dual dual.cpp)
add_executable(Dual dual.cpp)
add_executable(Tetra tetra.cpp)
add_executable(Slice slice.cpp)
add_executable(SliceFunc slice_func.cpp)
......@@ -16,21 +17,32 @@ add_executable(Scale scale.cpp)
add_executable(Move move.cpp)
add_executable(Convert convert.cpp)
add_executable(SameMeshDifference difference_same.cpp)
add_executable(MeshDifference difference_map.cpp)
add_executable(Kmeans kmeans.cpp)
add_executable(Agglomerate agglomerate.cpp)
add_executable(acute_mesh acute_grid.cpp)
add_executable(disturbed_mesh disturbed_grid.cpp)
add_executable(hex_mesh hex_grid.cpp)
add_executable(kershaw_mesh kershaw_grid.cpp)
add_executable(nonconvex_mesh nonconvex_grid.cpp)
add_executable(shestakov_mesh shestakov_grid.cpp)
add_executable(sinusoidal_mesh sinusoidal_grid.cpp)
add_executable(split_faces split_faces.cpp)
add_executable(glue_faces glue_faces.cpp)
add_executable(MeshDifference difference_map.cpp)
add_executable(Kmeans kmeans.cpp)
add_executable(Agglomerate agglomerate.cpp)
add_executable(acute_mesh acute_grid.cpp)
add_executable(disturbed_mesh disturbed_grid.cpp)
add_executable(hex_mesh hex_grid.cpp)
add_executable(kershaw_mesh kershaw_grid.cpp)
add_executable(nonconvex_mesh nonconvex_grid.cpp)
add_executable(shestakov_mesh shestakov_grid.cpp)
add_executable(sinusoidal_mesh sinusoidal_grid.cpp)
add_executable(split_faces split_faces.cpp)
add_executable(glue_faces glue_faces.cpp)
add_executable(check_collapse check_collapse.cpp)
add_library(FractureLib fracture.cpp fracture.h)
add_library(FractureLib fracture.cpp fracture.h)
target_link_libraries(FixTinyCollapse inmost)
if(USE_MPI)
message("linking FixTinyCollapse with MPI")
target_link_libraries(FixTinyCollapse ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(FixTinyCollapse PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS FixTinyCollapse EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(split_faces inmost)
if(USE_MPI)
message("linking split_faces with MPI")
......@@ -39,8 +51,8 @@ if(USE_MPI)
set_target_properties(split_faces PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS split_faces EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS split_faces EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(check_collapse inmost)
if(USE_MPI)
message("linking check_collapse with MPI")
......@@ -51,7 +63,7 @@ if(USE_MPI)
endif(USE_MPI)
install(TARGETS check_collapse EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(glue_faces inmost)
if(USE_MPI)
message("linking glue_faces with MPI")
......@@ -91,7 +103,8 @@ if(USE_MPI)
set_target_properties(UniteFaces PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS UniteFaces EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS UniteFaces EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Dual inmost)
......@@ -102,8 +115,8 @@ if(USE_MPI)
set_target_properties(Dual PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS Dual EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS Dual EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Tetra inmost)
if(USE_MPI)
message("linking Tetra with MPI")
......@@ -112,7 +125,9 @@ if(USE_MPI)
set_target_properties(Tetra PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS Tetra EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS Tetra EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Slice inmost)
......@@ -123,7 +138,8 @@ if(USE_MPI)
set_target_properties(Slice PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS Slice EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS Slice EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(SliceFunc inmost)
......@@ -134,7 +150,8 @@ if(USE_MPI)
set_target_properties(SliceFunc PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS SliceFunc EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS SliceFunc EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(SplitNonplanar inmost)
......@@ -145,7 +162,8 @@ if(USE_MPI)
set_target_properties(SplitNonplanar PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS SplitNonplanar EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS SplitNonplanar EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(CollapseDegenerate inmost)
......@@ -156,7 +174,8 @@ if(USE_MPI)
set_target_properties(CollapseDegenerate PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS CollapseDegenerate EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS CollapseDegenerate EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Bnd2Stl inmost)
......@@ -167,7 +186,8 @@ if(USE_MPI)
set_target_properties(Bnd2Stl PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS Bnd2Stl EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS Bnd2Stl EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Sector inmost)
if(USE_MPI)
......@@ -177,7 +197,8 @@ if(USE_MPI)
set_target_properties(Sector PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS Sector EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS Sector EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Sew inmost)
if(USE_MPI)
......@@ -187,7 +208,8 @@ if(USE_MPI)
set_target_properties(Sew PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS Sew EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS Sew EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Scale inmost)
if(USE_MPI)
......@@ -197,7 +219,8 @@ if(USE_MPI)
set_target_properties(Scale PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS Scale EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS Scale EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Move inmost)
......@@ -208,7 +231,8 @@ if(USE_MPI)
set_target_properties(Move PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS Move EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS Move EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Convert inmost)
if(USE_MPI)
......@@ -218,7 +242,8 @@ if(USE_MPI)
set_target_properties(Convert PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS Convert EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS Convert EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(SameMeshDifference inmost)
if(USE_MPI)
......@@ -228,7 +253,8 @@ if(USE_MPI)
set_target_properties(SameMeshDifference PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS SameMeshDifference EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS SameMeshDifference EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(MeshDifference inmost)
if(USE_MPI)
......@@ -239,7 +265,7 @@ if(USE_MPI)
endif()
endif(USE_MPI)
install(TARGETS MeshDifference EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Kmeans inmost)
if(USE_MPI)
message("linking Kmeans with MPI")
......@@ -249,7 +275,7 @@ if(USE_MPI)
endif()
endif(USE_MPI)
install(TARGETS Kmeans EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Agglomerate inmost)
if(USE_MPI)
message("linking Agglomerate with MPI")
......@@ -259,7 +285,7 @@ if(USE_MPI)
endif()
endif(USE_MPI)
install(TARGETS Agglomerate EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(acute_mesh inmost)
if(USE_MPI)
message("linking acute_mesh with MPI")
......@@ -269,7 +295,7 @@ if(USE_MPI)
endif()
endif(USE_MPI)
install(TARGETS acute_mesh EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(disturbed_mesh inmost)
if(USE_MPI)
message("linking disturbed_mesh with MPI")
......@@ -279,7 +305,7 @@ if(USE_MPI)
endif()
endif(USE_MPI)
install(TARGETS disturbed_mesh EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(hex_mesh inmost)
if(USE_MPI)
message("linking hex_mesh with MPI")
......@@ -288,8 +314,9 @@ if(USE_MPI)
set_target_properties(hex_mesh PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS hex_mesh EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS hex_mesh EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(kershaw_mesh inmost)
if(USE_MPI)
message("linking kershaw_mesh with MPI")
......@@ -299,7 +326,7 @@ if(USE_MPI)
endif()
endif(USE_MPI)
install(TARGETS kershaw_mesh EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(nonconvex_mesh inmost)
if(USE_MPI)
message("linking nonconvex_mesh with MPI")
......@@ -309,7 +336,7 @@ if(USE_MPI)
endif()
endif(USE_MPI)
install(TARGETS nonconvex_mesh EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(shestakov_mesh inmost)
if(USE_MPI)
message("linking shestakov_mesh with MPI")
......@@ -319,7 +346,7 @@ if(USE_MPI)
endif()
endif(USE_MPI)
install(TARGETS shestakov_mesh EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(sinusoidal_mesh inmost)
if(USE_MPI)
message("linking sinusoidal_mesh with MPI")
......@@ -328,7 +355,8 @@ if(USE_MPI)
set_target_properties(sinusoidal_mesh PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS sinusoidal_mesh EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS sinusoidal_mesh EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
set_property(TARGET FractureLib PROPERTY PUBLIC_HEADER "${CMAKE_CURRENT_SOURCE_DIR}/fracture.h")
......
#include "inmost.h"
using namespace INMOST;
int main(int argc, char ** argv)
{
if (argc < 2)
{
std::cout << "Usage: " << argv[0] << " mesh [mesh_out=grid.vtk]" << std::endl;
return -1;
}
std::string grid_out = "grid.vtk";
if (argc > 2) grid_out = std::string(argv[2]);
Mesh m;
m.Load(argv[1]);
Mesh::GeomParam table;
table[MEASURE] = CELL;
m.PrepareGeometricData(table);
std::cout << "Start:" << std::endl;
std::cout << "Cells: " << m.NumberOfCells() << std::endl;
std::cout << "Faces: " << m.NumberOfFaces() << std::endl;
std::cout << "Edges: " << m.NumberOfEdges() << std::endl;
double tt = Timer();
MarkerType collapse = m.CreateMarker();
int ncells = 0;
#pragma omp parallel for reduction(+:ncells)
for(int k = 0; k < m.CellLastLocalID(); ++k) if( m.isValidCell(k) )
{
Cell c = m.CellByLocalID(k);
ElementArray<Cell> around = c.BridgeAdjacencies2Cell(FACE);
double vol0 = c.Volume(), vol = vol0;
for(unsigned l = 0; l < around.size(); ++l)
vol += around[l].Volume();
if( vol0 / vol < 0.03 )
{
c.SetMarker(collapse);
ncells++;
}
}
std::cout << "collapse cells: " << ncells << std::endl;
int nedges = 0;
#pragma omp parallel for reduction(+:nedges)
for(int k = 0; k < m.EdgeLastLocalID(); ++k) if( m.isValidEdge(k) )
{
Edge e = m.EdgeByLocalID(k);
ElementArray<Cell> around = e.getCells();
for(unsigned l = 0; l < around.size(); ++l)
if( around[l].GetMarker(collapse) )
{
e.SetMarker(collapse);
nedges++;
break;
}
}
std::cout << "collapse edges: " << nedges << std::endl;
int ncollapsed = 0;
for(int k = 0; k < m.EdgeLastLocalID(); ++k) if( m.isValidEdge(k) )
{
Edge e = m.EdgeByLocalID(k);
if( e.GetMarker(collapse) )
{
Edge::CollapseEdge(e,0);
ncollapsed++;
if( ncollapsed % 100 == 0 )
m.Save("collapsed.vtk");
if( !Element::CheckConnectivity(&m) )
{
std::cout << "connectivity is bad" << std::endl;
throw -1;
}
}
}
std::cout << "Time to fix tiny:" << Timer() - tt << std::endl;
std::cout << "Cells: " << m.NumberOfCells() << std::endl;
std::cout << "Faces: " << m.NumberOfFaces() << std::endl;
std::cout << "Edges: " << m.NumberOfEdges() << std::endl;
m.Save(grid_out);
return 0;
}
......@@ -1084,6 +1084,7 @@ namespace INMOST
{
Node n = InvalidNode();
Mesh & m = *e.GetMeshLink();
e.PrintElementConnectivity();
ElementArray<Cell> cells = e.getCells();
Element::adj_type & faces = m.HighConn(e.GetHandle());
std::cout << "number of faces " << faces.size() << std::endl;
......@@ -1121,6 +1122,7 @@ namespace INMOST
coords[k] = 0.5*(e.getBeg().Coords()[k] + e.getEnd().Coords()[k]);
n = m.CreateNode(coords.data());
Element::adj_type & n_edges = m.HighConn(n.GetHandle());
Element::adj_type & n_cells = m.LowConn(n.GetHandle());
dynarray<HandleType,128> all_cells, all_faces, all_edges;
HandleType nn[2] = {e.getBeg().GetHandle(),e.getEnd().GetHandle()};
//Collect all cells, faces and edges
......@@ -1130,6 +1132,7 @@ namespace INMOST
{
std::cout << "node " << GetHandleID(nn[k]) << std::endl;
Element::adj_type & cells = m.LowConn(nn[k]);
std::cout << "cells: " << cells.size() << std::endl;
for(unsigned l = 0; l < cells.size(); ++l)
{
if( !m.GetMarker(cells[l],unique) )
......@@ -1141,7 +1144,7 @@ namespace INMOST
Element::adj_type & cell_nodes = m.HighConn(cells[l]);
for(unsigned j = 0; j < cell_nodes.size(); ++j)
{
std::cout << " cell " << GetHandleID(cells[l]) << " node " << GetHandleID(cell_nodes[j]) << std::endl;
//std::cout << " cell " << GetHandleID(cells[l]) << " node " << GetHandleID(cell_nodes[j]) << std::endl;
if( cell_nodes[j] == nn[k] )
{
std::cout << "replace node " << GetHandleID(nn[k]) << " in cell " << GetHandleID(cells[l]) << " with " << n.LocalID() << std::endl;
......@@ -1151,6 +1154,7 @@ namespace INMOST
}
cells.clear();
Element::adj_type & edges = m.HighConn(nn[k]);
std::cout << "edges: " << edges.size() << std::endl;
for(unsigned l = 0; l < edges.size(); ++l) if( edges[l] != e.GetHandle() )
{
if( !m.GetMarker(edges[l],unique) )
......@@ -1163,13 +1167,13 @@ namespace INMOST
std::cout << "visit edge " << GetHandleID(edges[l]) << std::endl;
if( nodes[0] == nn[k] )
{
std::cout << "(0) replace node " << GetHandleID(nn[k]) << " in edge " << GetHandleID(edges[l]) << std::endl;
std::cout << "replace node " << GetHandleID(nn[k]) << " in edge " << GetHandleID(edges[l]) << std::endl;
nodes[0] = n.GetHandle();
n_edges.push_back(edges[l]);
}
if( nodes[1] == nn[k] )
{
std::cout << "(1) replace node " << GetHandleID(nn[k]) << " in edge " << GetHandleID(edges[l]) << std::endl;
std::cout << "replace node " << GetHandleID(nn[k]) << " in edge " << GetHandleID(edges[l]) << std::endl;
nodes[1] = n.GetHandle();
n_edges.push_back(edges[l]);
}
......@@ -1182,29 +1186,34 @@ namespace INMOST
}
edges.clear();
//delete node
std::cout << "before number of faces " << faces.size() << std::endl;
std::cout << "delete node " << GetHandleID(nn[k]) << std::endl;
m.Delete(nn[k]);
std::cout << "after number of faces " << faces.size() << std::endl;
}
m.RemMarkerArray(all_cells.data(),all_cells.size(),unique);
m.RemMarkerArray(all_faces.data(),all_faces.size(),unique);
m.RemMarkerArray(all_edges.data(),all_edges.size(),unique);
m.ReleaseMarker(unique);
//m.ReleaseMarker(unique);
//Check faces to be deleted
//also should erase edge
MarkerType del = m.CreateMarker();
std::cout << "number of faces " << faces.size() << std::endl;
for(unsigned k = 0; k < faces.size(); ++k)
{
std::cout << "check face " << GetHandleID(faces[k]) << std::endl;
std::cout << "check face " << GetHandleID(faces[k]) << " " << Element::GeometricTypeName(m.GetGeometricType(faces[k])) << std::endl;
//two edges will match, should replace them with one
if( m.GetGeometricType(faces[k]) == Element::Tri )
Element::adj_type & edges = m.LowConn(faces[k]);
if( edges.size() == 3 )
{
std::cout << "collapse triangle " << GetHandleID(faces[k]) << std::endl;
Element::adj_type & edges = m.LowConn(faces[k]);
std::cout << "edges: ";
for(unsigned l = 0; l < 3; ++l)
std::cout << GetHandleID(edges[l]) << " ";
std::cout << std::endl;
for(unsigned l = 0; l < 3; ++l)
if( edges[l] == e.GetHandle() )
{
std::cout << "found collapsed edge at " << l << std::endl;
//delete e2 keep e1
HandleType e1 = edges[(l+1)%3], e2 = edges[(l+2)%3];
if( m.GetMarker(e1,del_protect) )
......@@ -1212,9 +1221,18 @@ namespace INMOST
std::swap(e1,e2);
assert(!m.GetMarker(e1,del_protect));
}
std::cout << "delete edge " << GetHandleID(e1) << " keep " << GetHandleID(e2) << std::endl;
if( e1 == e2 ) throw -1;
//reconnect adjacent faces to e1
Element::adj_type & edge_faces = m.HighConn(e2);
Element::adj_type & edge_faces_e1 = m.HighConn(e1);
m.SetMarkerArray(edge_faces_e1.data(),edge_faces_e1.size(),unique);
std::cout << "was faces: ";
for(unsigned r = 0; r < edge_faces_e1.size(); ++r) std::cout << GetHandleID(edge_faces_e1[r]) << " ";
std::cout << std::endl;
//replace link from e2 to e1 in all of the faces of the e2
std::cout << "replace edge in faces connected to " << GetHandleID(e2) << std::endl;
for(unsigned r = 0; r < edge_faces.size(); ++r)
if( edge_faces[r] != faces[k] )
{
......@@ -1222,12 +1240,28 @@ namespace INMOST
for(unsigned j = 0; j < face_edges.size(); ++j)
if( face_edges[j] == e2 )
{
std::cout << "replace edge in face " << GetHandleID(edge_faces[r]) << std::endl;
face_edges[j] = e1;
std::cout << "add face " << GetHandleID(edge_faces[r]) << " to edge " << GetHandleID(e1) << std::endl;
if( !m.GetMarker(edge_faces[r],unique) )
edge_faces_e1.push_back(edge_faces[r]);
break;
}
}
//cleanup links to faces in current edge
m.RemMarkerArray(edge_faces_e1.data(),edge_faces_e1.size(),unique);
//cleanup links to faces in e2
edge_faces.clear();
//erase link to face in e1
for(unsigned r = 0; r < edge_faces_e1.size(); ++r)
if( edge_faces_e1[r] == faces[k] )
{
std::cout << "erase face " << GetHandleID(faces[k]) << " from edge " << GetHandleID(e1) << std::endl;
edge_faces_e1.erase(edge_faces_e1.begin()+r);
break;
}
std::cout << "now faces: ";
for(unsigned r = 0; r < edge_faces_e1.size(); ++r) std::cout << GetHandleID(edge_faces_e1[r]) << " ";
std::cout << std::endl;
//delete this edge
m.SetMarker(e2,del);
break;
......@@ -1256,66 +1290,169 @@ namespace INMOST
{
std::cout << "erase edge in face " << GetHandleID(faces[k]) << std::endl;
//remove link to edge
Element::adj_type & face_edges = m.LowConn(faces[k]);
for(unsigned l = 0; l < face_edges.size(); ++l)
if( face_edges[l] == e.GetHandle() )
for(unsigned l = 0; l < edges.size(); ++l)
if( edges[l] == e.GetHandle() )
{
std::cout << "edge erased at " << l << std::endl;
face_edges.erase(face_edges.begin()+l);
edges.erase(edges.begin()+l);
break;
}
}
}
//check cells to be deleted
//also correct two links to the same node
//correct two links to the same node
for(unsigned k = 0; k < cells.size(); ++k)
{
if( m.GetGeometricType(cells.at(k)) == Element::Tet )
{
std::cout << "collapse tetrahedron " << GetHandleID(cells.at(k)) << std::endl;
//there are 2 faces now that should become just one
Element::adj_type & cell_faces = m.LowConn(cells.at(k));
assert(cell_faces.size() == 2);
HandleType f1 = cell_faces[0], f2 = cell_faces[1];
//reconnect adjacent cells to f1
Element::adj_type & face_cells = m.HighConn(f2);
//replace link from f2 to f1 in all of the cells of the f2
for(unsigned r = 0; r < face_cells.size(); ++r)
if( face_cells[r] != cells.at(k) )
std::cout << "erase node in cell " << GetHandleID(cells.at(k)) << std::endl;
//there are two subsequent links to the same node
Element::adj_type & cell_nodes = m.HighConn(cells.at(k));
int cnt = 0;
for(unsigned l = 0; l < cell_nodes.size(); ++l)
if( cell_nodes[l] == n.GetHandle() ) cnt++;
if( cnt == 2 )
{
for(unsigned l = 0; l < cell_nodes.size(); ++l)
if( cell_nodes[l] == n.GetHandle() )
{
Element::adj_type & cell_faces2 = m.LowConn(face_cells[r]);
for(unsigned j = 0; j < cell_faces2.size(); ++j)
if( cell_faces2[j] == f2 )
cell_faces2[j] = f1;
std::cout << "node erased at " << l << std::endl;
cell_nodes.erase(cell_nodes.begin()+l);
break;
}
face_cells.clear();
m.SetMarker(f2,del);
m.SetMarker(cells.at(k),del);
}
else
else
{
std::cout << "erase node in cell " << GetHandleID(cells.at(k)) << std::endl;
//there are two subsequent links to the same node
Element::adj_type & cell_nodes = m.HighConn(cells.at(k));
int cnt = 0;
for(unsigned l = 0; l < cell_nodes.size(); ++l)
if( cell_nodes[l] == n.GetHandle() ) cnt++;
if( cnt == 2 )
std::cout << "node encountered " << cnt << std::endl;
assert(cnt < 2);
}
}
//connect new node to remaining cells
for(int k = 0; k < all_cells.size(); ++k)
{
if( !m.GetMarker(all_cells[k],del) )
n_cells.push_back(all_cells[k]);
}
//delete degenerate cells
for(unsigned k = 0; k < cells.size(); ++k)
{
std::cout << "check cell " << GetHandleID(cells.at(k)) << " " << Element::GeometricTypeName(m.GetGeometricType(cells.at(k))) << std::endl;
Element::adj_type & cell_faces = m.LowConn(cells.at(k));
Element::adj_type & cell_nodes = m.HighConn(cells.at(k));
std::cout << "faces: ";
for(unsigned l = 0; l < cell_faces.size(); ++l) std::cout << GetHandleID(cell_faces[l]) << " ";
std::cout << std::endl;
cells[k].PrintElementConnectivity();
bool collapse_cell = false;
if( cell_faces.size() < 4 )
collapse_cell = true;
else //advanced check for degeneracy
{
for(unsigned l = 0; l < cell_faces.size(); ++l)
if( m.LowConn(cell_faces[l]).size() == cell_nodes.size() ) // number of edges in one of the faces is equal to all nodes
{
collapse_cell = true;
break;
}
}
if( collapse_cell )
{
std::cout << "collapse cell " << GetHandleID(cells.at(k)) << std::endl;
std::cout << "faces: ";
for( unsigned l = 0; l < cell_faces.size(); ++l) std::cout << GetHandleID(cell_faces[l]) << " ";
std::cout << std::endl;
for( unsigned l = 0; l < cell_faces.size(); ++l)
{
Face(&m,cell_faces[l]).PrintElementConnectivity();
}
//first case - there are 2 remaining faces
if( cell_faces.size() == 2 )
{
for(unsigned l = 0; l < cell_nodes.size(); ++l)
if( cell_nodes[l] == n.GetHandle() )
//there are 2 faces now that should become just one
assert(cell_faces.size() == 2);
HandleType f1 = cell_faces[0], f2 = cell_faces[1];
//reconnect adjacent cells to f1
Element::adj_type & face_cells = m.HighConn(f2);
Element::adj_type & face_cells_f1 = m.HighConn(f1);
std::cout << "was cells: ";
for(unsigned r = 0; r < face_cells_f1.size(); ++r) std::cout << GetHandleID(face_cells_f1[r]) << " ";
std::cout << std::endl;
//replace link from f2 to f1 in all of the cells of the f2
for(unsigned r = 0; r < face_cells.size(); ++r)
if( face_cells[r] != cells.at(k) )
{
std::cout << "node erased at " << l << std::endl;
cell_nodes.erase(cell_nodes.begin()+l);
Element::adj_type & cell_faces2 = m.LowConn(face_cells[r]);
for(unsigned j = 0; j < cell_faces2.size(); ++j)
if( cell_faces2[j] == f2 )
{
cell_faces2[j] = f1;
face_cells_f1.push_back(face_cells[r]);
//assert(face_cells_f1.size() <= 2); //maybe have to delete current cell?
}
}
std::cout << "now cells: ";
for(unsigned r = 0; r < face_cells_f1.size(); ++r) std::cout << GetHandleID(face_cells_f1[r]) << " ";
std::cout << std::endl;
face_cells.clear();
m.SetMarker(f2,del);
}
else //second case - replace face that holds all the nodes by multiple other faces
{
HandleType fdel = InvalidHandle();
for(unsigned l = 0; l < cell_faces.size(); ++l)
if( m.LowConn(cell_faces[l]).size() == cell_nodes.size() ) // number of edges in one of the faces is equal to all nodes
{
fdel = cell_faces[l];
break;
}
assert(fdel != InvalidHandle());