Merge branch 'master' of https://github.com/INMOST-DEV/INMOST

Conflicts:
	Examples/Octree/main_test.cpp

Examples/AdaptiveGrid/amesh.cpp

-#include "agrid.h"
-#include <math.h>
-
-
-Storage::integer min_nonzero(Storage::integer a, Storage::integer b)
-{
-	Storage::integer get_max = std::max(a,b);
-	Storage::integer get_min = std::min(a,b);
-	if( get_min == 0 ) return get_max; else return get_min;
-}
-
-void AdaptiveGrid
-
-void AdaptiveGrid::Refine(const TagInteger & indicator)
-{
-	int schedule_counter = 1;
-	bool scheduled = true;
-	//0. Extend indicator for edges and faces
-	indicator = CreateTag(indicator.GetTagName(),DATA_INTEGER,FACE|EDGE,NONE,1);
-	while(scheduled)
-	{
-		//1.Communicate indicator - it may be not synced
-		ExchangeTag(indicator,0,CELL);
-		//2.Propogate indicator down to the faces,edges
-		//  select latest schedule for them
-		//
-		//possible problem - a cell needs all it's elements to be splitted in the same schedule
-		//so the cell must mark all it's elements without hanging nodes into the same schedule,
-		//unless there is an even earlier cell
-		//
-		//should we select earliest possible for elements without hanging nodes and
-		//latest possible for elements with hanging nodes?
-		//
-		//with loop over cells we would mark elements without hanging nodes of the cells
-		//in current schedule
-		for(ElementType etype = FACE; etype >= EDGE; etype = PrevElementType(etype))
-		{
-#if defined(USE_OMP)
-#pragma omp parallel for
-#endif
-			for(Storage::integer it = 0; it < LastLocalID(etype); ++it) if( isValidElement(etype,it) )
-			{
-				Element e = ElementByLocalID(etype,it);
-				ElementArray<Element> adj = f.getAdjElements(NextElementType(etype));
-				//here latest schedule is selected
-				if( e->nbAdjElements(NODE,hanging_node,0) > 0 ) //latest possible schedule
-				{
-					for(ElementArray<Element>::size_type kt = 0; kt < adj.size(); ++kt)
-						if( indicator[adj[kt]] )
-							indicator[f] = indicator[f] ? std::min(indicator[f],indicator[adj[kt]]) : indicator[adj[kt]];
-				}
-				else //earliest possible schedule
-				{
-					for(ElementArray<Element>::size_type kt = 0; kt < adj.size(); ++kt)
-						indicator[f] = std::max(indicator[f],indicator[adj[kt]]);
-				}
-				
-			}
-		}
-		//3.Communicate indicator on faces and edges
-		ExchangeTag(indicator,0,FACE|EDGE);
-		//4.Check for each cell without indicator if there is
-		//  any hanging node with adjacent in a need to refine,
-		//  schedule for refinement earlier.
-		scheduled = false;
-#if defined(USE_OMP)
-#pragma omp parallel for
-#endif
-		for(Storage::integer it = 0; it < CellLastLocalID(); ++it) if( isValidCell(it) )
-		{
-			Cell c = CellByLocalID(it);
-			if( indicator[c] == 0 )
-			{
-				bool scheduled_c = false;
-				//retrive only hanging nodes, this may be empty
-				ElementArray<Node> hanging = c->getNodes(hanging_node,0);
-				for(ElementArray<Node>::size_type kt = 0; kt < hanging.size() && !scheduled_c; ++kt)
-				{
-					//adjacent edges may be scheduled for refinement
-					ElementArray<Edge> adj = hanging[kt].getEdges();
-					for(ElementArray<Edge>::size_type lt = 0; lt < adj.size() && !scheduled_c; ++lt)
-						if( indicator[adj[lt]] != 0 )
-						{
-							indicator[c] = schedulde_counter+1;
-							scheduled = scheduled_c = true;
-						}
-				}
-			}
-		}
-		//5.Go back to 1 until no new elements scheduled
-		if( scheduled ) schedulde_counter++;
-	}
-	//5.Refine
-	BeginModification();
-	while(scheduled_counter)
-	{
-		scheduled_counter--;
-	}
-	ResolveModification();
-	//Let the user update their data
-	ApplyModification();
-	EndModification();
-}
-
-void AdaptiveGrid::Coarse(const TagInteger & indicator)
-{
-}

Examples/AdaptiveGrid/amesh.h

-#ifndef _AGRID_H
-#define _AGRID_H
-#include "inmost.h"
-
-
-class AdaptiveMesh : public Mesh
-{
-	MarkerType hanging_node;
-	
-	void SplitFace();
-public:
-	AdaptiveMesh() : Mesh() {}
-	~AdaptiveMesh() {}
-	/// Indicator must be 1 on cells to be refined
-	/// and 0 on all other cells
-	void Refine(const TagInteger & indicator);
-	void Coarse(const TagInteger & indicator);
-};
-
-#endif //_AGRID_H

Examples/AdaptiveMesh/CMakeLists.txt

+project(AdaptiveMesh)
+
+
+add_library(AdaptiveMeshLib amesh.cpp amesh.h)
+add_executable(AdaptiveMesh main.cpp)
+add_executable(AdaptiveMeshCube main_cube.cpp)
+add_executable(AdaptiveMeshSphere main_sphere.cpp)
+
+target_link_libraries(AdaptiveMesh inmost AdaptiveMeshLib)
+target_link_libraries(AdaptiveMeshCube inmost AdaptiveMeshLib)
+target_link_libraries(AdaptiveMeshSphere inmost AdaptiveMeshLib)
+
+if(USE_MPI)
+  message("linking AdaptiveMesh with MPI")
+  target_link_libraries(AdaptiveMesh ${MPI_LIBRARIES}) 
+  target_link_libraries(AdaptiveMeshCube ${MPI_LIBRARIES}) 
+  target_link_libraries(AdaptiveMeshSphere ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(AdaptiveMesh PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+    set_target_properties(AdaptiveMeshCube PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+    set_target_properties(AdaptiveMeshSphere PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+
+set_property(TARGET AdaptiveMeshLib PROPERTY PUBLIC_HEADER "${CMAKE_CURRENT_SOURCE_DIR}/amesh.h")
+
+install(TARGETS AdaptiveMeshLib EXPORT inmost-targets 
+        LIBRARY DESTINATION lib 
+        ARCHIVE DESTINATION lib
+        PUBLIC_HEADER DESTINATION include)
+install(TARGETS AdaptiveMesh EXPORT inmost-targets RUNTIME DESTINATION bin)
+install(TARGETS AdaptiveMeshCube EXPORT inmost-targets RUNTIME DESTINATION bin)
+install(TARGETS AdaptiveMeshSphere EXPORT inmost-targets RUNTIME DESTINATION bin)
\ No newline at end of file

Examples/AdaptiveMesh/amesh.h

+#ifndef _AMESH_H
+#define _AMESH_H
+#include "inmost.h"
+
+namespace INMOST
+{
+	class AdaptiveMesh : public Mesh
+	{
+		ElementSet root; //< Root set that links all the other sets for coarsements
+		TagInteger level; //< Refinement level of the cell
+		TagReference parent_set; //<Link to the set that contains an element.
+		TagReferenceArray hanging_nodes; //< Link to current hanging nodes of the cell.
+		/// Prepare sets for coarsements.
+		/// Do not do this in constructor, since mesh may contain no cells.
+		void PrepareSet();
+	public:
+		Storage::integer GetLevel(const Storage & e) {return level[e];}
+		AdaptiveMesh();
+		~AdaptiveMesh();
+		/// Indicator must be 1 on cells to be refined
+		/// and 0 on all other cells
+		bool Refine(TagInteger & indicator);
+		bool Coarse(TagInteger & indicator);
+		/// Delete all data related to mesh refinement-coarsement.
+		void ClearData();
+	};
+}
+
+#endif //_AMESH_H

Examples/AdaptiveMesh/main.cpp

+#include "amesh.h"
+
+using namespace INMOST;
+
+int main(int argc, char ** argv)
+{
+	Mesh::Initialize(&argc,&argv);
+	
+	if( argc > 1 )
+	{
+		AdaptiveMesh m;
+		m.Load(argv[1]);
+		//m.SetTopologyCheck(NEED_TEST_CLOSURE);
+		//m.SetTopologyCheck(PROHIBIT_MULTILINE);
+		//m.SetTopologyCheck(PROHIBIT_MULTIPOLYGON);
+		TagInteger indicator = m.CreateTag("INDICATOR",DATA_INTEGER,CELL,NONE,1);
+		int max_levels = 2;
+		if( argc > 2 ) max_levels = atoi(argv[2]);
+		//bounding box around mesh
+		Storage::real cmax[3] = {-1.0e20,-1.0e20,-1.0e20};
+		Storage::real cmin[3] = {+1.0e20,+1.0e20,+1.0e20};
+		Storage::real xyz[3], r, q, cnt[3];
+		//find bounding box around mesh
+		for(Mesh::iteratorNode it = m.BeginNode(); it != m.EndNode(); ++it)
+		{
+			for(int d = 0; d < 3; ++d)
+			{
+				if( it->Coords()[d] > cmax[d] ) cmax[d] = it->Coords()[d];
+				if( it->Coords()[d] < cmin[d] ) cmin[d] = it->Coords()[d];
+			}
+		}
+		r = 1;
+		for(int d = 0; d < 3; ++d)
+		{
+			r *= cmax[d]-cmin[d];
+			cnt[d] = (cmax[d]+cmin[d])*0.5;
+		}
+		r = pow(r,1.0/3.0)/20.0;
+		
+		for(int k = 0; k < 16; ++k)
+		{
+			
+			int numref;
+			int refcnt = 0;
+			do
+			{
+				numref = 0;
+				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( m.GetLevel(it->self()) < max_levels )
+				{
+					it->Centroid(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)
+					{
+						indicator[it->self()] = 1;
+						numref++;
+					}
+				}
+				if( numref )
+				{
+					std::cout << "k " << k << " refcnt " << refcnt << " " <<  r*k << " < r < " << r*(k+1) << " numref " << numref << " cells " << m.NumberOfCells() << std::endl;
+					/*
+					{
+						std::stringstream file;
+						file << "indicator_" << k << "_" << refcnt << "r.pmf";
+						m.Save(file.str());
+					}
+					*/
+					if( !m.Refine(indicator) ) break;
+					
+					{
+						std::stringstream file;
+						file << "dump_" << k << "_" << refcnt << "r.pmf";
+						m.Save(file.str());
+					}
+					
+					//cleanup indicator
+					for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) indicator[it->self()] = 0;
+				}
+				refcnt++;
+			}
+			while(numref);
+			refcnt = 0;
+			do
+			{
+				numref = 0;
+				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( m.GetLevel(it->self()) > 0 )
+				{
+					it->Centroid(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)
+					{
+						indicator[it->self()] = 1;
+						numref++;
+					}
+				}
+				if( numref )
+				{
+					std::cout << "k " << k << " crscnt " << refcnt << " " << r*k << " < r < " << r*(k+1) << " numcrs " << numref << " cells " << m.NumberOfCells() <<  std::endl;
+					/*
+					{
+						std::stringstream file;
+						file << "indicator_" << k << "_" << refcnt << "c.pmf";
+						m.Save(file.str());
+					}
+					 */
+					if( !m.Coarse(indicator) ) break;
+					
+					{
+						std::stringstream file;
+						file << "dump_" << k << "_" << refcnt << "c.pmf";
+						m.Save(file.str());
+					}
+					
+					
+					
+					//cleanup indicator
+					for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) indicator[it->self()] = 0;
+				}
+				//return -1;
+				refcnt++;
+			}
+			while(numref);
+			
+			{
+				std::stringstream file;
+				file << "step_" << k << ".vtk";
+				m.Save(file.str());
+			}
+			
+			{
+				std::stringstream file;
+				file << "step_" << k << ".pmf";
+				m.Save(file.str());
+			}
+		}
+	}
+	else std::cout << "Usage: " << argv[0] << " mesh_file [max_levels=2]" << std::endl;
+}
\ No newline at end of file

Examples/AdaptiveMesh/main_cube.cpp

+#include "amesh.h"
+
+using namespace INMOST;
+
+int main(int argc, char ** argv)
+{
+	Mesh::Initialize(&argc,&argv);
+	
+	if( argc > 1 )
+	{
+		AdaptiveMesh m;
+		m.Load(argv[1]);
+		TagInteger indicator = m.CreateTag("INDICATOR",DATA_INTEGER,CELL,NONE,1);
+		
+		int max_levels = 2;
+		if( argc > 2 ) max_levels = atoi(argv[2]);
+		
+		int numref;
+		do
+		{
+			numref = 0;
+			for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
+				if( m.GetLevel(it->self()) < max_levels )
+				{
+					double x[3];
+					it->Centroid(x);
+					if( x[0] > 0.3 && x[0] < 0.7 && x[1] > 0.3 && x[1] < 0.7 && x[2] > 0.3 && x[2] < 0.7)
+					{
+						indicator[it->self()] = 1;
+						numref++;
+					}
+				}
+			if( numref )
+			{
+				if( !m.Refine(indicator) ) break;
+				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) indicator[it->self()] = 0;
+			}
+		}
+		while(numref);
+		std::string file = "out.vtk";
+		if( argc > 3 ) file = std::string(argv[3]);
+		m.Save(file);
+	}
+	else std::cout << "Usage: " << argv[0] << " mesh_file [max_levels=2] [mesh_out=out.vtk]" << std::endl;
+}
\ No newline at end of file

Examples/AdaptiveMesh/main_sphere.cpp

+#include "amesh.h"
+
+using namespace INMOST;
+
+int main(int argc, char ** argv)
+{
+	Mesh::Initialize(&argc,&argv);
+	
+	if( argc > 1 )
+	{
+		AdaptiveMesh m;
+		m.Load(argv[1]);
+		TagInteger indicator = m.CreateTag("INDICATOR",DATA_INTEGER,CELL,NONE,1);
+		
+		int max_levels = 2;
+		if( argc > 2 ) max_levels = atoi(argv[2]);
+		
+		int numref;
+		do
+		{
+			numref = 0;
+			for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
+				if( m.GetLevel(it->self()) < max_levels )
+				{
+					double x[3];
+					it->Centroid(x);
+					if( sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5)+(x[2]-0.5)*(x[2]-0.5)) < 0.25 )
+					{
+						indicator[it->self()] = 1;
+						numref++;
+					}
+				}
+			if( numref )
+			{
+				if( !m.Refine(indicator) ) break;
+				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) indicator[it->self()] = 0;
+			}
+		}
+		while(numref);
+		std::string file = "out.vtk";
+		if( argc > 3 ) file = std::string(argv[3]);
+		m.Save(file);
+	}
+	else std::cout << "Usage: " << argv[0] << " mesh_file [max_levels=2] [mesh_out=out.vtk]" << std::endl;
+}
\ No newline at end of file

Examples/CMakeLists.txt

 if(USE_MESH)
+add_subdirectory(AdaptiveMesh)
 add_subdirectory(DrawGrid)
 add_subdirectory(OldDrawGrid)
 add_subdirectory(GridGen)
@@ -14,6 +15,7 @@ endif(USE_SOLVER AND USE_MESH)
 if(USE_AUTODIFF AND USE_SOLVER AND USE_MESH)
 add_subdirectory(ADFVDiscr)
 add_subdirectory(ADMFD)
+add_subdirectory(MFD-ES)
 endif(USE_AUTODIFF AND USE_SOLVER AND USE_MESH)
 #add_subdirectory(OctreeCutcell)
 if(USE_MESH AND USE_PARTITIONER)

Examples/MFD-ES/CMakeLists.txt

+project(MFDES)
+set(SOURCE main.cpp)
+
+add_executable(MFDES ${SOURCE})
+
+target_link_libraries(MFDES inmost)
+
+if(USE_SOLVER)
+  if(USE_SOLVER_ANI)
+    message("linking MFDES with ani3d and BLAS")
+    target_link_libraries(MFDES ani3d ${BLAS_LIBRARIES})
+    if(BLAS_LINKER_FLAGS)
+      set_target_properties(MFDES PROPERTIES LINK_FLAGS "${BLAS_LINKER_FLAGS}")
+    endif()
+  endif()
+  if(USE_SOLVER_PETSC)
+    message("linking MFDES with PETSc")
+    target_link_libraries(MFDES ${PETSC_LIBRARIES})
+  endif()
+  if(USE_SOLVER_TRILINOS)
+    message("linking MFDES with Trilinos")
+    target_link_libraries(MFDES ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
+  endif()
+  if(USE_SOLVER_METIS)
+    message("linking MFDES with Metis")
+    target_link_libraries(MFDES ${METIS_LIBRARIES})
+  endif()
+  if(USE_SOLVER_MONDRIAAN)
+    message("linking MFDES with Mondriaan")
+    target_link_libraries(MFDES ${MONDRIAAN_LIBRARIES})
+  endif()
+  if(USE_SOLVER_SUPERLU)
+    message("linking MFDES with SuperLU")
+    target_link_libraries(MFDES ${SUPERLU_LIBRARIES})
+  endif()
+endif()
+
+
+if(USE_PARTITIONER)
+  if(USE_PARTITIONER_ZOLTAN)
+    message("linking MFDES with Zoltan")
+    target_link_libraries(MFDES ${ZOLTAN_LIBRARIES})
+  endif()
+  if(USE_PARTITIONER_PARMETIS)
+    message("linking MFDES with ParMETIS")
+    target_link_libraries(MFDES ${PARMETIS_LIBRARIES})
+  endif()
+endif()
+
+if(USE_MPI)
+  message("linking MFDES with MPI")
+  target_link_libraries(MFDES ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(MFDES PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+
+install(TARGETS MFDES EXPORT inmost-targets RUNTIME DESTINATION bin)

Examples/Octree/CMakeLists.txt

@@ -48,6 +48,7 @@ if(OPENGL_FOUND)
 	endif(GLUT_FOUND)
 else(OPENGL_FOUND)
 	message("OpenGL not found for ${vtarget}")
+	
 endif(OPENGL_FOUND)
 
 if(USE_PARTITIONER_ZOLTAN)
@@ -61,4 +62,6 @@ endif()
 
 endforeach(vtarget ${TARGETS})
 
-
+if(NOT GLUT_FOUND)
+	set_target_properties(Octree PROPERTIES EXCLUDE_FROM_ALL 1 EXCLUDE_FROM_DEFAULT_BUILD 1)
+endif(NOT GLUT_FOUND)

Examples/OldDrawGrid/main.cpp

@@ -106,7 +106,7 @@ void printtext(const char * fmt, ... )
 {
 	
 	unsigned int i;
-	char stext[4096];
+	char stext[1048576];
 	va_list ap;
 	if ( fmt == NULL ) return;
 	va_start(ap,fmt);
@@ -4102,13 +4102,14 @@ double display_elem_info(Element e, double top, double left, double interval)
 	top -= interval;
 	glRasterPos2f(left,top);
 	printtext("%s %d",ElementTypeName(e->GetElementType()),e->LocalID());
+	if( e->GetElementType() == ESET ) printtext(" size %d",e->getAsSet().Size());
 	glColor3f(0.2,0.2,0.2);
 	for(Mesh::iteratorTag t = mesh->BeginTag(); t != mesh->EndTag(); ++t) if( t->isDefined(e->GetElementType()) )
 	{
 		if( e->HaveData(*t) )
 		{
-			char str[4096];
-			char temp[4096];
+			char str[1048576];
+			char temp[1048576];
 			str[0] = '\0';
 			int dsize;
 			switch(t->GetDataType())
@@ -4475,7 +4476,7 @@ void draw_screen()
 		for(int k = 0; k < (int)orphans.size(); ++k)
 			DrawElement(orphans[k],color_t(1,0,1),color_t(0,1,1),color_t(0,0,1));
 
-	if( disp_e.isValid() )
+	if( disp_e.isValid() && disp_e.GetElementType() != ESET)
 		DrawElement(disp_e,color_t(1,1,0),color_t(1,0,0),color_t(0,0,1));
 
 
@@ -4588,6 +4589,7 @@ void draw_screen()
 							else if ( stype == "edge" ) visualization_type = EDGE;
 							else if ( stype == "face" ) visualization_type = FACE;
 							else if ( stype == "cell" ) visualization_type = CELL;
+							else if ( stype == "eset" ) visualization_type = ESET;
 							if( visualization_type == NONE )
 								printf("unknown element type %s\n",typen);
 							if( !mesh->isValidElement(visualization_type,comp) )
@@ -4596,9 +4598,12 @@ void draw_screen()
 							if( mesh->isValidElement(visualization_type,comp) )
 							{
 								disp_e = mesh->ElementByLocalID(visualization_type,comp);
-								disp_e->Centroid(shift);
-								for(int r = 0; r < 3; ++r)
-									shift[r] = -shift[r];
+								if( disp_e.GetElementType() != ESET )
+								{
+									disp_e->Centroid(shift);
+									for(int r = 0; r < 3; ++r)
+										shift[r] = -shift[r];
+								}
 							}
 						}
 

Source/Autodiff/autodiff.cpp

@@ -54,34 +54,74 @@ namespace INMOST
 #endif //USE_MESH
 	
 #if defined(USE_MESH)
-	Automatizator::Automatizator() :first_num(0), last_num(0) {}
+	Automatizator::Automatizator(const Automatizator & b) : name(b.name+"_copy")
+	{
+		std::vector<INMOST_DATA_ENUM_TYPE> regs = b.ListRegisteredTags();
+		for(std::vector<INMOST_DATA_ENUM_TYPE>::iterator kt = regs.begin(); kt != regs.end(); ++kt)
+			RegisterTag(b.GetValueTag(*kt),b.GetElementType(*kt),b.GetMask(*kt));
+		if( b.last_num != 0 ) EnumerateTags();
+	}
+	Automatizator & Automatizator::operator =(Automatizator const & b)
+	{