Some algorithms and bugfixes

Fix a bug requesting back/front cell of a face during mesh modification

Comment unfinished parts of inmost in AdaptiveMesh

Update elasticity example in Examples/ADMFD

Add algorithms to split faces and to glue faces of the mesh (for tests
with MFD elasticty)

Examples/AdaptiveMesh/amesh.cpp

@@ -332,8 +332,8 @@ namespace INMOST
 		//free created tag
 		DeleteTag(indicator,FACE|EDGE);
 		//11. Restore parallel connectivity, global ids
-		ResolveShared(true);
-		ResolveModification();
+		//ResolveShared(true);
+		//ResolveModification();
 		//12. Let the user update their data
 		//todo: call back function for New() cells
 		//13. Delete old elements of the mesh
@@ -642,8 +642,8 @@ namespace INMOST
 		//free created tag
 		DeleteTag(indicator,FACE|EDGE);
 		//todo:
-		ResolveShared(true);
-		ResolveModification();
+		//ResolveShared(true);
+		//ResolveModification();
 		//todo:
 		//let the user update their data
 		ApplyModification();

Examples/AdaptiveMesh/main.cpp

@@ -15,6 +15,14 @@ int main(int argc, char ** argv)
 		//m.SetTopologyCheck(PROHIBIT_MULTILINE);
 		//m.SetTopologyCheck(PROHIBIT_MULTIPOLYGON);
 		TagInteger indicator = m.CreateTag("INDICATOR",DATA_INTEGER,CELL,NONE,1);
+		/*
+		for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
+			indicator[*it] = 1;
+		if( !m.Refine(indicator) ) std::cout << "refine failed" << std::endl;
+		else std::cout << "refine ok!" << std::endl;
+		m.Save("grid.pmf");
+		return 0;
+		*/
 		int max_levels = 2;
 		if( argc > 2 ) max_levels = atoi(argv[2]);
 		//bounding box around mesh
@@ -159,4 +167,5 @@ int main(int argc, char ** argv)
 	else std::cout << "Usage: " << argv[0] << " mesh_file [max_levels=2]" << std::endl;
 	
 	Mesh::Finalize();
+	return 0;
 }

Examples/GridTools/CMakeLists.txt

@@ -25,8 +25,31 @@ 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_library(FractureLib fracture.cpp fracture.h)
+add_executable(sinusoidal_mesh sinusoidal_grid.cpp)
+add_executable(split_faces split_faces.cpp)
+add_executable(glue_faces glue_faces.cpp)
+add_library(FractureLib fracture.cpp fracture.h)
+
+target_link_libraries(split_faces inmost)
+if(USE_MPI)
+  message("linking split_faces with MPI")
+  target_link_libraries(split_faces ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    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)
+
+target_link_libraries(glue_faces inmost)
+if(USE_MPI)
+  message("linking glue_faces with MPI")
+  target_link_libraries(glue_faces ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(glue_faces PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+install(TARGETS glue_faces EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
+
 
 target_link_libraries(FixFaults inmost)
 if(USE_MPI)

Examples/GridTools/glue_faces.cpp

+#include "inmost.h"
+#include <stdio.h>
+#include <deque>
+//#include "tetgen/tetgen.h"
+using namespace INMOST;
+typedef Storage::real real;
+typedef Storage::real_array real_array;
+
+
+
+int main(int argc, char ** argv)
+{
+	if( argc < 2 ) 
+	{
+		printf("Usage: %s input_mesh [output_mesh]\n",argv[0]);
+		return -1;
+	}
+
+	Mesh m;
+	m.SetFileOption("ECL_CURVILINEAR","FALSE");
+	m.SetFileOption("ECL_SPLIT_GLUED","TRUE");
+	m.Load(argv[1]);
+	Tag sliced;
+	if( m.HaveTag("SLICED") ) sliced = m.GetTag("SLICED");
+	
+	ElementArray<Cell> unite(&m,2);
+	
+	m.BeginModification();
+	
+	//MarkerType mrk = m.CreateMarker();
+	//Face last, next, test(&m,ComposeHandle(FACE,382));
+	for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it)
+	{
+		//std::cout << "test face " << test.LocalID() << " back " << test.BackCell().LocalID() << " front " << test.FrontCell().LocalID() << std::endl;
+		/*
+		int a = last.isValid() ? last.LocalID() : -1;
+		if( last.isValid() && next.isValid() )
+		{
+			//std::cout << "unite " << last.BackCell().LocalID();
+			//std::cout << " (" << Element::GeometricTypeName(last.BackCell().GetGeometricType()) << ":" << last.BackCell().nbAdjElements(FACE) << ")";
+			//std::cout << " and " << last.FrontCell().LocalID();
+			//std::cout << " (" << Element::GeometricTypeName(last.FrontCell().GetGeometricType()) << ":" << last.FrontCell().nbAdjElements(FACE) << ")";
+			//std::cout << " for " << last.LocalID();
+			//std::cout << " next " << next.LocalID() << std::endl;
+			unite[0] = last->BackCell();
+			unite[1] = last->FrontCell();
+			Cell c = Cell::UniteCells(unite,0);
+			c.SetMarker(mrk);
+			last = next;
+			next = InvalidFace();
+		}
+		*/
+		if( !it->Boundary() && (!sliced.isValid() || !it->HaveData(sliced)) )
+		{
+			//if( !it->BackCell().GetMarker(mrk) && !it->FrontCell().GetMarker(mrk) )
+			if( !it->BackCell().New() && !it->FrontCell().New() )
+			{
+				
+				//std::cout << "candidate " << it->LocalID();
+				//std::cout << " back " << it->BackCell()->LocalID();
+				//std::cout << " (" << Element::GeometricTypeName(it->BackCell().GetGeometricType()) <<  ":" << it->BackCell().nbAdjElements(FACE) <<")";
+				//std::cout << " front " << it->FrontCell()->LocalID();
+				//std::cout << " (" << Element::GeometricTypeName(it->FrontCell().GetGeometricType()) <<  ":" << it->BackCell().nbAdjElements(FACE) <<")";
+				//std::cout << std::endl;
+				/*
+				it->BackCell().SetMarker(mrk);
+				it->FrontCell().SetMarker(mrk);
+				if( !last.isValid() )
+					last = it->self();
+				else
+					next = it->self();
+				 */
+				//break;
+				
+				unite[0] = it->BackCell();
+				unite[1] = it->FrontCell();
+				Cell c = Cell::UniteCells(unite,0);
+				
+			}
+		}
+	}
+	
+	
+	m.EndModification();
+	
+	TagInteger tetra = m.CreateTag("TET",DATA_INTEGER,CELL,NONE,1);
+	
+	int ntet = 0;
+	for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
+		if( it->GetGeometricType() == Element::Tet )
+		{
+			tetra[*it] = 1;
+			ElementArray<Cell> cells = it->BridgeAdjacencies2Cell(FACE);
+			for(int k = 0; k < cells.size(); ++k)
+				if( cells[k]->GetGeometricType() == Element::Tet )
+					tetra[*it] = 2;
+			//std::cout << "Cell:" << it->LocalID() << " " << it->nbAdjElements(FACE) << std::endl;
+			ntet ++;
+		}
+	
+	std::cout << "tetrahedras: " << ntet << std::endl;
+	
+							  
+	if( argc > 2 )
+		m.Save(argv[2]);
+	else
+		m.Save("out.vtk");
+
+	return 0;
+}

Examples/GridTools/split_faces.cpp

+#include "inmost.h"
+#include <stdio.h>
+#include <deque>
+//#include "tetgen/tetgen.h"
+using namespace INMOST;
+typedef Storage::real real;
+typedef Storage::real_array real_array;
+
+
+
+int main(int argc, char ** argv)
+{
+	if( argc < 2 ) 
+	{
+		printf("Usage: %s input_mesh [output_mesh] [special_triangle_split=1]\n",argv[0]);
+		return -1;
+	}
+	
+	int special_triangle_split = 1;
+	if( argc > 3 ) special_triangle_split = atoi(argv[3]);
+
+	Mesh m;
+	m.SetFileOption("ECL_CURVILINEAR","FALSE");
+	m.SetFileOption("ECL_SPLIT_GLUED","TRUE");
+	m.Load(argv[1]);
+	double cnt[3];
+	ElementArray<Node> edge_nodes(&m,2);
+	MarkerType new_node = m.CreateMarker();
+	m.BeginModification();
+	for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it) if( !it->New() )
+	{
+		for(int k = 0; k < 3; ++k)
+			cnt[k] = (it->getBeg().Coords()[k]+it->getEnd().Coords()[k])*0.5;
+		Node n = m.CreateNode(cnt);
+		n.SetMarker(new_node);
+		Edge::SplitEdge(it->self(),ElementArray<Node>(&m,1,n.GetHandle()),0);
+	}
+	for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) if( !it->New() )
+	{
+		ElementArray<Node> new_nodes = it->getNodes(new_node);
+		ElementArray<Edge> new_edges(&m,new_nodes.size());
+		if( new_nodes.size() == 3 )
+		{
+			for(int k = 0; k < 3; ++k)
+			{
+				edge_nodes[0] = new_nodes[k];
+				edge_nodes[1] = new_nodes[(k+1)%3];
+				new_edges[k] = m.CreateEdge(edge_nodes).first;
+			}
+		}
+		else
+		{
+			cnt[0] = cnt[1] = cnt[2] = 0;
+			for(int k = 0; k < new_nodes.size(); ++k)
+			{
+				cnt[0] += new_nodes[k].Coords()[0];
+				cnt[1] += new_nodes[k].Coords()[1];
+				cnt[2] += new_nodes[k].Coords()[2];
+			}
+			cnt[0] /= new_nodes.size();
+			cnt[1] /= new_nodes.size();
+			cnt[2] /= new_nodes.size();
+			Node n = m.CreateNode(cnt);
+			edge_nodes[0] = n;
+			for(int k = 0; k < new_nodes.size(); ++k)
+			{
+				edge_nodes[1] = new_nodes[k];
+				new_edges[k] = m.CreateEdge(edge_nodes).first;
+			}
+		}
+		Face::SplitFace(it->self(),new_edges,0);
+	}
+	
+	m.EndModification();
+	m.ReleaseMarker(new_node,NODE);
+							  
+							  
+	if( argc > 2 )
+		m.Save(argv[2]);
+	else
+		m.Save("out.vtk");
+
+	return 0;
+}

Source/Mesh/face.cpp

@@ -48,7 +48,8 @@ namespace INMOST
 				enumerator i = ENUMUNDEF;
 				MarkerType hm = m->HideMarker();
 				i = m->getNext(hc.data(),static_cast<enumerator>(hc.size()),i,hm); //found first
-				i = m->getNext(hc.data(),static_cast<enumerator>(hc.size()),i,hm); //found second
+				if( i != hc.size() )
+					i = m->getNext(hc.data(),static_cast<enumerator>(hc.size()),i,hm); //found second
 				if( i != hc.size() ) return Cell(m,hc[i]);
 			}
 			return Cell(m,InvalidHandle());
@@ -106,13 +107,16 @@ namespace INMOST
 				enumerator i = ENUMUNDEF;
 				MarkerType hm = m->HideMarker();
 				i = m->getNext(lc.data(),static_cast<enumerator>(lc.size()),i,hm);
-				i = m->getNext(lc.data(),static_cast<enumerator>(lc.size()),i,hm);
-				if( i != static_cast<enumerator>(lc.size()) ) 
+				if( i != static_cast<enumerator>(lc.size()) )
 				{
-					adj_type const & llc = m->LowConn(lc[i]);
-					enumerator j = ENUMUNDEF;
-					j = m->getNext(llc.data(),static_cast<enumerator>(llc.size()),j,hm);
-					if( j != static_cast<enumerator>(llc.size()) ) return Node(m,llc[j]);
+					i = m->getNext(lc.data(),static_cast<enumerator>(lc.size()),i,hm);
+					if( i != static_cast<enumerator>(lc.size()) )
+					{
+						adj_type const & llc = m->LowConn(lc[i]);
+						enumerator j = ENUMUNDEF;
+						j = m->getNext(llc.data(),static_cast<enumerator>(llc.size()),j,hm);
+						if( j != static_cast<enumerator>(llc.size()) ) return Node(m,llc[j]);
+					}
 				}
 			}
 			return Node(m,InvalidHandle());

Source/Mesh/modify.cpp

@@ -247,6 +247,7 @@ namespace INMOST
 		//delete inner faces
 		for(dynarray<HandleType,64>::size_type j = 0; j < inner_faces.size(); j++)
 		{
+			//std::cout << "delete face " << GetHandleID(inner_faces[j]) << std::endl;
 			if( m->GetMarker(inner_faces[j],rem) )
 			{
 				m->RemMarker(inner_faces[j],rem);
@@ -1865,7 +1866,7 @@ namespace INMOST
 	{
 		Mesh * m = GetMeshLink();
 		
-		if( m->HaveGeometricData(ORIENTATION,FACE) )
+		//if( m->HaveGeometricData(ORIENTATION,FACE) )
 		{
 			//retrive faces
 			MarkerType hm = m->HideMarker();
@@ -1883,7 +1884,7 @@ namespace INMOST
 						std::swap(hc[k1],hc[k2]);
 						//hc[k2] = GetHandle(); //cannot use the cell because virtualization table is already destroyed and FixNormalOrientation will do bad things
 						//hc.resize(1); //just remove element, we will do this anyway later
-						Face(m,lc[it])->FixNormalOrientation(); //restore orientation
+						if( m->HaveGeometricData(ORIENTATION,FACE) ) Face(m,lc[it])->FixNormalOrientation(); //restore orientation
 					}
 				}
 			}