Merge branch 'master' of https://boogie.inm.ras.ru/terekhov/INMOST-developement

.gitlab-ci.yml

+# use the official gcc image, based on debian
+# can use versions as well, like gcc:5.2
+image: gcc
+
+stages:
+  - build
+  - test
+
+build_debug:
+  stage: build
+  script:
+    - mkdir build_debug
+    - cd build_debug
+    - cmake -DCOMPILE_TESTS=ON -DUSE_OMP=ON -DCMAKE_CXX_FLAGS="-O0 -g" -DCMAKE_C_FLAGS="-O0 -g" ..
+    - make VERBOSE=1
+  artifacts:
+    paths:
+    - build_debug/
+
+build_opt:
+  stage: build
+  script:
+    - mkdir build_opt
+    - cd build_opt
+    - cmake -DCOMPILE_TESTS=ON -DCOMPILE_EXAMPLES=ON -DUSE_OMP=ON -DCMAKE_CXX_FLAGS="-Ofast -march=native" -DCMAKE_C_FLAGS="-Ofast -march=native" ..
+    - make VERBOSE=1
+  artifacts:
+    paths:
+    - build_opt/
+
+test_debug:
+  stage: test
+  script:
+    - cd build_debug
+    - ctest --output-on-failure
+  dependencies:
+    - build_debug
+
+test_opt:
+  stage: test
+  script:
+    - cd build_opt
+    - ctest --output-on-failure
+  dependencies:
+    - build_opt

CMakeLists.txt

@@ -16,6 +16,7 @@ option(USE_MPI_P2P "Use MPI point to point functionality, may be faster with har
 option(USE_MPI_FILE "Use MPI extension to work with files, may save a lot of memory" ON)
 option(USE_MPI2 "Use MPI-2 extensions, useful if your MPI library warns you to use new functions" ON)
 option(USE_OMP "Compile with OpenMP support (experimental)" OFF)
+option(USE_OPENCL "Use OpenCL where possible (experimental)" OFF)
 
 option(USE_MESH "Compile mesh capabilities" ON)
 option(USE_SOLVER "Compile solver capabilities" ON)
@@ -33,7 +34,6 @@ option(USE_SOLVER_MONDRIAAN "Use Mondriaan for matrix reordering" OFF)
 option(USE_SOLVER_PETSC "Use PETSc solvers" OFF)
 option(USE_SOLVER_TRILINOS "Use Trilinos solvers" OFF)
 option(USE_SOLVER_SUPERLU "Use SuperLU solver" OFF)
-#option(USE_AUTODIFF_OPENCL "Use OpenCL for automatic differentiation (under work)" OFF)
 #option(USE_AUTODIFF_ASMJIT "Use AsmJit for automatic differentiation" OFF)
 #option(USE_AUTODIFF_EXPRESSION_TEMPLATES "Use c++ expression templates for automatic differentiation" OFF)
 
@@ -206,13 +206,13 @@ if(USE_SOLVER_SUPERLU)
 	endif()
 endif()
 
-if(USE_AUTODIFF_OPENCL)
+if(USE_OPENCL)
 	find_package(OpenCL)
 	if(OPENCL_FOUND)
 		include_directories(${OPENCL_INCLUDE_DIRS})
-		link_directories(${PETSC_LIBRARY_DIRS})
+		link_directories(${OpenCL_LIBRARY})
 	else()
-		set(USE_AUTODIFF_OPENCL OFF CACHE BOOL "Use OpenCL for automatic differentiation" FORCE)
+		set(USE_OPENCL OFF CACHE BOOL "Use OpenCL where possible (experimental)" FORCE)
 		message("OpenCL not found")
 	endif()
 endif()

Examples/ADFVDiscr/main.cpp

+#include "inmost.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
 #include <string.h>
 
-#include "inmost.h"
 using namespace INMOST;
 
 #ifndef M_PI
@@ -90,7 +90,7 @@ int main(int argc,char ** argv)
 		{ // currently only non-distributed meshes are supported by Inner_RCM partitioner
 			ttt = Timer();
 			Partitioner * p = new Partitioner(m);
-			p->SetMethod(Partitioner::Inner_RCM,Partitioner::Partition); // Specify the partitioner
+			p->SetMethod(Partitioner::INNER_KMEANS,Partitioner::Partition); // Specify the partitioner
 			p->Evaluate(); // Compute the partitioner and store new processor ID in the mesh
 			delete p;
 			BARRIER;

Examples/ADMFD/main.cpp

+#include "inmost.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
 #include <string.h>
 
-#include "inmost.h"
 using namespace INMOST;
 
 #ifndef M_PI
@@ -125,7 +125,7 @@ int main(int argc,char ** argv)
 			{ // Compute mesh partitioning
 				ttt = Timer();
 				Partitioner p(m); //Create Partitioning object
-				p.SetMethod(Partitioner::Inner_RCM,repartition ? Partitioner::Repartition : Partitioner::Partition); // Specify the partitioner
+				p.SetMethod(Partitioner::INNER_KMEANS,repartition ? Partitioner::Repartition : Partitioner::Partition); // Specify the partitioner
 				p.Evaluate(); // Compute the partitioner and store new processor ID in the mesh
 				BARRIER
 				if( m->GetProcessorRank() == 0 ) std::cout << "Evaluate: " << Timer()-ttt << std::endl;

Examples/AdaptiveMesh/amesh.cpp

@@ -332,13 +332,15 @@ namespace INMOST
 		//free created tag
 		DeleteTag(indicator,FACE|EDGE);
 		//11. Restore parallel connectivity, global ids
-		//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
 		ApplyModification();
 		//14. Done
 		EndModification();
+		//ExchangeData(hanging_nodes,CELL | FACE,0);
 		//reorder element's data to free up space
 		ReorderEmpty(CELL|FACE|EDGE|NODE);
 		//return number of refined cells
@@ -640,12 +642,14 @@ namespace INMOST
 		//free created tag
 		DeleteTag(indicator,FACE|EDGE);
 		//todo:
-		//ResolveModification();
+		ResolveShared(true);
+		ResolveModification();
 		//todo:
 		//let the user update their data
 		ApplyModification();
 		//done
 		EndModification();
+		//ExchangeData(hanging_nodes,CELL | FACE,0);
 		//cleanup null links to hanging nodes
 		for(ElementType etype = FACE; etype <= CELL; etype = NextElementType(etype))
 		{
@@ -666,4 +670,4 @@ namespace INMOST
 		call_counter++;
 		return ret != 0;
 	}
-}
\ No newline at end of file
+}

Examples/AdaptiveMesh/main.cpp

@@ -9,6 +9,7 @@ int main(int argc, char ** argv)
 	if( argc > 1 )
 	{
 		AdaptiveMesh m;
+		m.SetCommunicator(INMOST_MPI_COMM_WORLD);
 		m.Load(argv[1]);
 		//m.SetTopologyCheck(NEED_TEST_CLOSURE);
 		//m.SetTopologyCheck(PROHIBIT_MULTILINE);
@@ -29,6 +30,8 @@ int main(int argc, char ** argv)
 				if( it->Coords()[d] < cmin[d] ) cmin[d] = it->Coords()[d];
 			}
 		}
+		m.AggregateMax(cmax,3);
+		m.AggregateMin(cmin,3);
 		r = 1;
 		for(int d = 0; d < 3; ++d)
 		{
@@ -47,18 +50,20 @@ int main(int argc, char ** argv)
 				numref = 0;
 				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( m.GetLevel(it->self()) < max_levels )
 				{
-					it->Centroid(xyz);
+					it->Barycenter(xyz);
 					//refine a circle
 					q = 0;
 					for(int d = 0; d < 3; ++d)
 						q += (xyz[d]-cnt[d])*(xyz[d]-cnt[d]);
 					q = sqrt(q);
 					if( q < r*(k+1) && q > r*k)
+					//if( q < 0.02 )
 					{
 						indicator[it->self()] = 1;
 						numref++;
 					}
 				}
+				numref = m.Integrate(numref);
 				if( numref )
 				{
 					std::cout << "k " << k << " refcnt " << refcnt << " " <<  r*k << " < r < " << r*(k+1) << " numref " << numref << " cells " << m.NumberOfCells() << std::endl;
@@ -89,7 +94,7 @@ int main(int argc, char ** argv)
 				numref = 0;
 				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( m.GetLevel(it->self()) > 0 )
 				{
-					it->Centroid(xyz);
+					it->Barycenter(xyz);
 					//refine a circle
 					q = 0;
 					for(int d = 0; d < 3; ++d)
@@ -101,6 +106,7 @@ int main(int argc, char ** argv)
 						numref++;
 					}
 				}
+				numref = m.Integrate(numref);
 				if( numref )
 				{
 					std::cout << "k " << k << " crscnt " << refcnt << " " << r*k << " < r < " << r*(k+1) << " numcrs " << numref << " cells " << m.NumberOfCells() <<  std::endl;
@@ -131,11 +137,19 @@ int main(int argc, char ** argv)
 			
 			{
 				std::stringstream file;
-				file << "step_" << k << ".vtk";
+				file << "step_" << k << ".pvtk";
 				m.Save(file.str());
 			}
 			
 			{
+				TagInteger tag_owner = m.CreateTag("OWN",DATA_INTEGER,CELL,NONE,1);
+				TagInteger tag_owner0 = m.GetTag("OWNER_PROCESSOR");
+				TagInteger tag_stat = m.CreateTag("STAT",DATA_INTEGER,CELL,NONE,1);
+				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
+				{
+					tag_owner[*it] = tag_owner0[*it];
+					tag_stat[*it] = it->GetStatus();
+				}
 				std::stringstream file;
 				file << "step_" << k << ".pmf";
 				m.Save(file.str());
@@ -143,4 +157,6 @@ int main(int argc, char ** argv)
 		}
 	}
 	else std::cout << "Usage: " << argv[0] << " mesh_file [max_levels=2]" << std::endl;
-}
\ No newline at end of file
+	
+	Mesh::Finalize();
+}

Examples/FVDiscr/main.cpp

+#include "inmost.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
 #include <string.h>
 
-#include "inmost.h"
 using namespace INMOST;
 
 #ifndef M_PI
@@ -92,7 +92,7 @@ int main(int argc,char ** argv)
     { // currently only non-distributed meshes are supported by Inner_RCM partitioner
 			ttt = Timer();
 			Partitioner * p = new Partitioner(m);
-			p->SetMethod(Partitioner::Inner_RCM,Partitioner::Partition); // Specify the partitioner
+			p->SetMethod(Partitioner::INNER_KMEANS,Partitioner::Partition); // Specify the partitioner
 			p->Evaluate(); // Compute the partitioner and store new processor ID in the mesh
 			delete p;
 			BARRIER
@@ -327,4 +327,4 @@ int main(int argc,char ** argv)
 #endif
 	Solver::Finalize(); // Finalize solver and close MPI activity
 	return 0;
-}
\ No newline at end of file
+}

Examples/GridGen/main.cpp

+#include "inmost.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
 #include <string.h>
 
-#include "inmost.h"
 using namespace INMOST;
 
 #define MESH_SIZE 32
@@ -227,7 +227,7 @@ Mesh * ParallelGenerator(INMOST_MPI_Comm comm, int ng, int nx, int ny, int nz)
 				{
 					CreateCubeElement(m,verts);
 				}
-				else if ((i + j) % 2 == 0) // Create two prism cells
+				else if ((i + j) % 2 == 0 || ng == 6) // Create two prism cells
 				{
 					CreateNWPrismElements(m,verts);
 				}
@@ -327,4 +327,4 @@ int main(int argc, char *argv[])
 	delete mesh;
 	Mesh::Finalize();
 	return 0;
-}
\ No newline at end of file
+}

Examples/GridTools/CMakeLists.txt

@@ -16,6 +16,8 @@ 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_library(FractureLib fracture.cpp fracture.h)
 
 target_link_libraries(FixFaults inmost)
@@ -185,6 +187,26 @@ if(USE_MPI)
 endif(USE_MPI)
 install(TARGETS MeshDifference EXPORT inmost-targets RUNTIME DESTINATION bin)
 
+target_link_libraries(Kmeans inmost)
+if(USE_MPI)
+  message("linking Kmeans with MPI")
+  target_link_libraries(Kmeans ${MPI_LIBRARIES})
+  if(MPI_LINK_FLAGS)
+    set_target_properties(Kmeans PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif()
+endif(USE_MPI)
+install(TARGETS Kmeans EXPORT inmost-targets RUNTIME DESTINATION bin)
+
+target_link_libraries(Agglomerate inmost)
+if(USE_MPI)
+  message("linking Agglomerate with MPI")
+  target_link_libraries(Agglomerate ${MPI_LIBRARIES})
+  if(MPI_LINK_FLAGS)
+    set_target_properties(Agglomerate PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif()
+endif(USE_MPI)
+install(TARGETS Agglomerate EXPORT inmost-targets RUNTIME DESTINATION bin)
+
 
 set_property(TARGET FractureLib PROPERTY PUBLIC_HEADER "${CMAKE_CURRENT_SOURCE_DIR}/fracture.h")
 

Examples/GridTools/agglomerate.cpp

+#include "inmost.h"
+
+using namespace INMOST;
+
+//todo: want to separate all faces into those that share edges
+//see todo: inside text
+
+
+int main(int argc, char ** argv)
+{
+	if (argc < 2)
+	{
+		std::cout << "Usage: " << argv[0] << " mesh [tag=MATERIAL] [mesh_out=grid.vtk]" << std::endl;
+		return -1;
+	}
+	
+	Mesh::Initialize(&argc,&argv);
+	
+	std::string tag_name = "MATERIAL";
+	if( argc > 2) tag_name = std::string(argv[2]);
+	
+	std::string grid_out = "grid.vtk";
+	if (argc > 3) grid_out = std::string(argv[3]);
+
+
+	Mesh m;
+	m.SetCommunicator(INMOST_MPI_COMM_WORLD);
+	m.Load(argv[1]);
+	
+	if( !m.HaveTag(tag_name) )
+	{
+		std::cout << "mesh does not have tag named " << tag_name << std::endl;
+		exit(-1);
+	}
+	
+	TagInteger mat = m.GetTag(tag_name);
+
+	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;
+
+	std::map< int,ElementArray<Cell> > mat_cells;
+	
+	for(int q = 0; q < m.CellLastLocalID(); ++q) if( m.isValidCell(q) )
+	{
+		Cell n = m.CellByLocalID(q);
+		mat_cells[ mat[n] ].push_back(n);
+	}
+	
+	//Unite cells
+	//todo: split cells that are not connected by face
+	std::cout << "Unite cells" << std::endl;
+	int nunited = 0;
+	for(std::map< int,ElementArray<Cell> >::iterator it = mat_cells.begin();
+		it != mat_cells.end(); ++it)
+	{
+		if( it->second.size() > 1 )
+		{
+			mat[Cell::UniteCells(it->second,0)] = it->first;
+			nunited++;
+		}
+	}
+	std::cout << "united: " << nunited << std::endl;
+	/*
+	std::cout << "Unite faces" << std::endl;
+	
+	std::map< std::pair<int,int>, ElementArray<Face> > faces;
+	
+	nunited = 0;
+	
+	
+	
+	for(int q = 0; q < m.FaceLastLocalID(); ++q) if( m.isValidFace(q) )
+	{
+		Face n = m.FaceByLocalID(q);
+		std::pair<int,int> f_cells;
+		int bc = n.BackCell().LocalID();;
+		int fc = n.FrontCell().isValid() ? n.FrontCell().LocalID() : -1;
+		f_cells.first = std::min(bc,fc);
+		f_cells.second = std::max(bc,fc);
+		faces[f_cells].push_back(n);
+	}
+	
+	std::cout << "computed faces" << std::endl;
+	
+	m.BeginModification();
+	
+	 
+	 //todo: split faces that are not connected by edge
+	for(std::map< std::pair<int,int>, ElementArray<Face> >::iterator it = faces.begin();
+		it != faces.end(); ++it)
+	{
+		if( it->second.size() > 1 )
+		{
+			std::cout << it->first.first << "<->" << it->first.second << " faces " << it->second.size() << std::endl;
+			Face::UniteFaces(it->second,0);
+			nunited++;
+		}
+	}
+	
+	m.EndModification();
+	
+	std::cout << "united: " << nunited << 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);
+	
+	Mesh::Finalize();
+	return 0;
+}

Examples/GridTools/collapse_degenerate.cpp

-#include <stdio.h>
 #include "inmost.h"
+#include <stdio.h>
 
 
 using namespace INMOST;
@@ -87,4 +87,4 @@ int main(int argc, char ** argv)
 	}
 
 	return 0;
-}
\ No newline at end of file
+}

Examples/GridTools/convert.cpp

-#include <stdio.h>
 #include "inmost.h"
+#include <stdio.h>
 
 
 using namespace INMOST;

Examples/GridTools/difference_map.cpp

+#include "inmost.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <time.h>
-#include "inmost.h"
 
 using namespace INMOST;
 

Examples/GridTools/difference_same.cpp

+#include "inmost.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <time.h>
-#include "inmost.h"
 
 using namespace INMOST;
 

Examples/GridTools/dual.cpp

-#include <stdio.h>
 #include "inmost.h"
+#include <stdio.h>
 
 
 using namespace INMOST;
@@ -273,4 +273,4 @@ int main(int argc, char ** argv)
 	}
 
 	return 0;
-}
\ No newline at end of file
+}

Examples/GridTools/move.cpp

-#include <stdio.h>
 #include "inmost.h"
+#include <stdio.h>
 
 
 using namespace INMOST;

Examples/GridTools/scale.cpp

-#include <stdio.h>
 #include "inmost.h"
+#include <stdio.h>
 
 
 using namespace INMOST;

Examples/GridTools/sector.cpp

-#include <stdio.h>
 #include "inmost.h"
+#include <stdio.h>
 
 
 using namespace INMOST;
@@ -10,7 +10,7 @@ int main(int argc, char ** argv)
 {
 	if( argc < 2 ) 
 	{
-		printf("Usage: %s input_mesh [rotation_angle=0 degrees] [output_mesh]\n",argv[0]);
+		printf("Usage: %s input_mesh [rotation_angle=0 degrees] [refine_boundary=1] [output_mesh]\n",argv[0]);
 		return -1;
 	}
 
@@ -19,7 +19,9 @@ int main(int argc, char ** argv)
 	
 	const double pi = 3.1415926536;
 	double theta = 0, ct, st;
+	int refine = 1;
 	if( argc > 2 ) theta = atof(argv[2])/180.0*pi;
+	if( argc > 3 ) refine = atoi(argv[3]);
 	ct = cos(theta);
 	st = sin(theta);
 	
@@ -39,6 +41,16 @@ int main(int argc, char ** argv)
 	std::cout << "x: " << xmin << ":" << xmax << std::endl;
 	std::cout << "y: " << ymin << ":" << ymax << std::endl;
 	
+	if( refine )
+	{
+		for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
+		{
+			double a = (n->Coords()[0]-xmin)/(xmax-xmin);
+			a = sin(3.14159265359*a/2.0);
+			n->Coords()[0] = (xmax-xmin)*a + xmin;
+		}
+	}
+	
 	for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
 	{
 		double x = n->Coords()[0], mx = x;
@@ -66,10 +78,12 @@ int main(int argc, char ** argv)
 	std::cout << "x: " << xmin << ":" << xmax << std::endl;
 	std::cout << "y: " << ymin << ":" << ymax << std::endl;
 	
-	if( argc > 3 )
+	
+	
+	if( argc > 4 )
 	{
-		std::cout << "Save to " << argv[3] << std::endl;
-		m.Save(argv[3]);
+		std::cout << "Save to " << argv[4] << std::endl;
+		m.Save(argv[4]);
 	}
 	else
 	{

Examples/GridTools/sew.cpp

-#include <stdio.h>
 #include "inmost.h"
+#include <stdio.h>
 
 
 using namespace INMOST;

Examples/GridTools/slice_func.cpp

@@ -4,6 +4,7 @@ using namespace INMOST;
 
 double func(double x, double y, double z, int n)
 {
+            
 	if( n == 0 )
 		return sqrt(x*x+y*y)-0.25;
 	else if( n == 1 )
@@ -38,7 +39,7 @@ double func(double x, double y, double z, int n)
 		return fmin( fmin(x-Lx,Rx-x), fmin(y-Ly, Ry-y) );
 		*/
 		
-		double Lx = 0., Rx = 2, Ly = 0.0, Ry = 4.5;
+            double Lx = 0., Rx = 20, Ly = -2.5, Ry = 17.5;
 	    if (x > Rx)
 	    {
 	        if (y > Ry) return -sqrt( (x-Rx)*(x-Rx) + (y-Ry)*(y-Ry) );
@@ -184,12 +185,21 @@ int main(int argc, char ** argv)
 	//m.RemTopologyCheck(THROW_EXCEPTION);
 	TagInteger material = m.CreateTag("MATERIAL",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
 	Tag sliced = m.CreateTag("SLICED",DATA_BULK,FACE|EDGE|NODE,FACE|EDGE|NODE,1);
+        TagReal level = m.CreateTag("level",DATA_REAL,NODE,NONE,1);
 
 	std::cout << "Cells: " << m.NumberOfCells() << std::endl;
 	std::cout << "Faces: " << m.NumberOfFaces() << std::endl;
 	
 	MarkerType original = m.CreateMarker();
-	for(Mesh::iteratorNode it = m.BeginNode(); it != m.EndNode(); ++it) it->SetMarker(original);
+#pragma omp parallel for
+        for(int k = 0; k < m.NodeLastLocalID(); ++k) if(m.isValidNode(k))
+        {
+            Node it = m.NodeByLocalID(k);
+            it->SetMarker(original);
+            material[it] = 3;
+            Storage::real_array c0 = it->Coords();
+            level[it] = func(c0[0],c0[1],c0[2],type);
+        }
 	for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it) it->SetMarker(original);
 
 	MarkerType slice = m.CreateMarker();
@@ -202,8 +212,8 @@ int main(int argc, char ** argv)
 		Storage::real_array c1 = it->getEnd()->Coords();
 		double r0 = func(c0[0],c0[1],c0[2],type);
 		double r1 = func(c1[0],c1[1],c1[2],type);
-		material[it->getBeg()] = (r0 <= 0)? 0 : 1;
-		material[it->getEnd()] = (r1 <= 0)? 0 : 1;
+                if (material[it->getBeg()] == 3) material[it->getBeg()] = (r0 <= 0)? 0 : 1;
+                if (material[it->getEnd()] == 3) material[it->getEnd()] = (r1 <= 0)? 0 : 1;
 		//std::cout << "r0 " << r0 << " r1 " << r1 << std::endl;
 		if( (r0*r1 < -1.0e-12) || (fabs(r0*r1) < 1.0e-12 && ((fabs(r0) < 1.0e-6) ^ (fabs(r1) < 1.0e-6))) )
 		{
@@ -322,6 +332,7 @@ int main(int argc, char ** argv)
 		{
 			material[*it] = 2;
 			std::cout << "oops, materials for edge nodes were not split, 0: " << mat[0] << ", 1: " << mat[1] << ", 2: " << mat[2] << std::endl;
+                        
 		}
 		else if( mat[0] != 0 ) material[*it] = 0;
 		else if( mat[1] != 0 ) material[*it] = 1;
@@ -359,24 +370,14 @@ int main(int argc, char ** argv)
 				//std::cout << "Sliced nodes " << nodes.size() << " total nodes: ";
 				nodes = it->getNodes();
 				//std::cout << nodes.size() << std::endl;
-				for(int q = 0; q < (int)nodes.size(); ++q)
-				{
-					Storage::real_array c0 = nodes[q].Coords();
-					double r0 = func(c0[0],c0[1],c0[2],type);
-					//std::cout << "NODE:" << nodes[q].LocalID() << " " << (nodes[q]->GetMarker(slice)?"":"not ") << "slice " << (nodes[q]->GetMarker(original)?"":"not ") << "original r="<<r0 << " material " << material[nodes[q]] << std::endl;
-				}
-				ElementArray<Edge> fedges = it->getEdges();
-				for(int q = 0; q < (int)fedges.size(); ++q)
-				{
-					//std::cout << "EDGE:" << fedges[q].LocalID() << " " << fedges[q].getBeg().LocalID() << "<->" << fedges[q].getEnd().LocalID() << " " << (fedges[q]->GetMarker(slice)?"":"not ") << "slice material " << material[fedges[q]] << std::endl;
-				}
 
 				
 				double c0[3],c1[3],pc0[3],pc1[3],p[3];
 				it->Centroid(c0);
 				double r0 = func(c0[0],c0[1],c0[2],type);
 				int material0 = (r0 <= 0)? 0 : 1;
-				bool s0 = false;
+                                
+                                bool s0 = false;
 				if( fabs(r0) < 1.0e-6 ) s0 = true;
 				
 				//std::cout << "Centernode " << (s0?"sliced":"") << " r=" << r0 << std::endl;
@@ -631,7 +632,7 @@ int main(int argc, char ** argv)
 			if( simple )
 			{
 				//gather loop by loop and create faces
-				int loop_cnt = 0;
+				//int loop_cnt = 0;
 				//order edges
 				ElementArray<Edge> order_edges(&m);
 				
@@ -748,7 +749,7 @@ int main(int argc, char ** argv)
 				double c0[3],c1[3],pc0[3],pc1[3],p[3];
 				it->Centroid(c0);
 				double r0 = func(c0[0],c0[1],c0[2],type);
-				int material0 = (r0 <= 0)? 0 : 1;
+
 				bool s0 = false;
 				if( fabs(r0) < 1.0e-6 ) s0 = true;
 				
@@ -800,22 +801,26 @@ int main(int argc, char ** argv)
 							l0 = sqrt(l0);
 							l1 = sqrt(l1);
 							l = l0+l1;
-							if( l0 < 5.0e-2*l ) //edge goes through centernode
+							if( l0 < 1.0e-3*l ) //edge goes through centernode
 							{
 								if( !centernode.isValid() )
 									centernode = m.CreateNode(c0);
 								cutnodes[q] = centernode;
 								//std::cout << "selected centernode " << cutnodes[q].LocalID() << std::endl;
 							}
-							else if( l1 > 5.0e-2*l )
+							else if( l1 > 1.0e-3*l )
 							{