Improvements

Corrected iteration over elements and calculation of number of elements
during modification of the mesh.

Corrected that saving meshes with infinity in data can produce bad vtk
files.

Examples/OctreeCutcell/main.cpp

@@ -6,7 +6,7 @@
 
 // in obj.cpp z coordinate is not normalized from 0 to 1!
 // may be a problem, lines 791-794 
-#include "../../inmost.h"
+#include "inmost.h"
 #include "octgrid.h"
 #include "my_glut.h"
 #include "rotate.h"
@@ -43,8 +43,16 @@ double gleft, gright, gbottom, gtop, gnear, gfar, zoom = 1;
 bool transparent = false;
 
 
+//#define SPHERE_MESH
+#define OBJ_READERS
+#define INIT_PERM
+
+
 void transformation(double xyz[3]) 
 {
+#if defined(SPHERE_MESH)
+	return;
+#endif
 	make_proj.Project(xyz);
 	/*
     double tmp[3];
@@ -64,11 +72,21 @@ void transformation(double xyz[3])
 
 mat_ret_type get_material_types(double xyz[3])
 {
+#if defined(SPHERE_MESH)
+	mat_ret_type ret;
+	const double eps = 1.0e-5;
+	double sphere = 0.5 - sqrt((xyz[0]-0.5)*(xyz[0]-0.5) + (xyz[1]-0.5)*(xyz[1]-0.5) + (xyz[2]-0.5)*(xyz[2]-0.5));
+	//double sphere = 1.0 - sqrt((xyz[0])*(xyz[0]) + (xyz[1])*(xyz[1]) + (xyz[2])*(xyz[2]));
+	if( sphere < eps ) ret.push_back(-1);
+	if( sphere > -eps ) ret.push_back(1);
+	return ret;
+#else
 	mat_ret_type ret(get_type.Size());
 	ret.resize(get_type.Layers(xyz,&ret[0]));
 	std::sort(ret.begin(),ret.end());
 	ret.resize(std::unique(ret.begin(),ret.end())-ret.begin());
 	return ret;
+#endif
 }
 
 
@@ -333,7 +351,9 @@ void cell_to_INMOST(struct grid * g, int cell, Cell r)
 		else throw -1;
 		*/
 	}
+#if defined(INIT_PERM)
 	fill_K(center,&r->RealArrayDF(g->Kvec)[0],&r->RealArrayDF(g->K)[0]);
+#endif
 }
 
 void init_mesh(struct grid * g)
@@ -368,6 +388,13 @@ std::map<Tag,Storage::real> cell_small_unite(ElementArray<Cell> & unite)
 	return ret;
 }
 
+void onclose()
+{
+	gridDelete(&thegrid);
+	DestroyObj();
+}
+
+
 #if defined( __GRAPHICS__)
 int show_region = 0;
 int width = 800, height = 600;
@@ -701,11 +728,6 @@ void reshape(int w, int h)
 
 int s = 0;
 
-void onclose()
-{
-	gridDelete(&thegrid);
-	DestroyObj();
-}
 
 void keyboard(unsigned char key, int x, int y)
 {
@@ -828,7 +850,7 @@ int global_test_number = -1;
 int main(int argc, char ** argv)
 {
 	int i;
-	int n[3] = {8,8,10};
+	int n[3] = {8,8,8};
 	last_time = Timer();
 	
 	if( argc > 3 )
@@ -838,6 +860,7 @@ int main(int argc, char ** argv)
 		n[2] = atoi(argv[3]);
 	}
 	
+#if defined(OBJ_READERS)
 	make_proj.SetGridStep(n[2]);
 	
 	InitObj();
@@ -872,7 +895,7 @@ int main(int argc, char ** argv)
 		*/
 		get_type.ReadLayers(layers);
 	}
-	
+#endif
 	
 	
 	thegrid.transformation = transformation;
@@ -947,15 +970,49 @@ int main(int argc, char ** argv)
 	DestroyObj();
 	return 0;
 	*/
+#if defined(INIT_PERM)
 	for(Mesh::iteratorFace it = thegrid.mesh->BeginFace(); it != thegrid.mesh->EndFace(); ++it)
 	{
 		Storage::real cnt[3];
 		it->Centroid(cnt);
 		fill_K(cnt,&it->RealArrayDF(thegrid.Kvec)[0],&it->RealArrayDF(thegrid.K)[0]);
 	}
+#endif
 	
-	
-	
+#if defined(SPHERE_MESH)
+	thegrid.mesh->BeginModification();
+	int cnt = 0;
+	for(Mesh::iteratorCell it = thegrid.mesh->BeginCell(); it != thegrid.mesh->EndCell(); ++it)
+		if( it->Integer(thegrid.cell_material) < 0 )
+		{
+			it->Delete();
+			cnt++;
+		}
+
+	std::cout << "deleted CELL " << cnt << std::endl;
+	{
+		ElementType etype = CELL;
+		do
+		{
+			cnt = 0;
+			etype = PrevElementType(etype);
+			for(Mesh::iteratorElement it = thegrid.mesh->BeginElement(etype); it != thegrid.mesh->EndElement(); ++it)
+				if( it->nbAdjElements(NextElementType(etype)) == 0 )
+				{
+					it->Delete();
+					cnt++;
+				}
+
+			std::cout << "deleted " << ElementTypeName(etype) << " " << cnt << std::endl;
+		}
+		while(etype != NODE);
+	}
+	thegrid.mesh->Save("sphere.vtk");
+	thegrid.mesh->SwapModification();
+	thegrid.mesh->EndModification();
+	onclose();
+	return 0;
+#endif
 
 
 
@@ -1014,4 +1071,5 @@ int main(int argc, char ** argv)
 	glutMainLoop();
 #endif
 	//~ std::cout << "Hello!" << std::endl;
+	onclose();
 }

Examples/OctreeCutcell/octgrid.cpp

@@ -755,7 +755,7 @@ public:
 	incident_matrix(InputIterator beg, InputIterator end, unsigned num_inner)
 	: head_column(beg,end), min_loop()
 	{
-		isInputForwardIterators<T,InputIterator>();
+		//isInputForwardIterators<T,InputIterator>();
 		if( !head_column.empty() )
 		{
 			Mesh * m = head_column[0]->GetMeshLink();

Examples/OctreeCutcell/octgrid.h

 #ifndef _OCTGRID_H
 #define _OCTGRID_H
-#include "../../inmost.h"
+#include "inmost.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>

Source/Headers/inmost_mesh.h

@@ -1220,7 +1220,7 @@ namespace INMOST
                                         senum;
 	private:
     std::string                         name;
-		Storage::real                       epsilon;
+		real                                epsilon;
 		empty_container                     empty_space[6];
 		empty_container                     empty_links[6];
 		links_container                     links[6];
@@ -1229,13 +1229,15 @@ namespace INMOST
 		Tag                                 tag_low_conn;
 		Tag                                 tag_high_conn;
 		Tag                                 tag_markers;
-    Tag *                               tag_private_markers;
+		Tag *                               tag_private_markers;
 		Tag                                 tag_geom_type;
 		Tag                                 tag_setname;
 		Tag                                 tag_setcomparator;
 		MeshState                           m_state;
 		integer                             dim;
 		HandleType                          last_created;
+		integer								hidden_count[6];
+		integer								hidden_count_zero[6];
 	private:
 		ElementType                         have_global_id;
 		INMOST_DATA_BIG_ENUM_TYPE           parallel_mesh_unique_id;
@@ -2832,18 +2834,24 @@ namespace INMOST
 		__INLINE integer                  FaceLastLocalID    () const {return static_cast<integer>(links[2].size());}
 		__INLINE integer                  CellLastLocalID    () const {return static_cast<integer>(links[3].size());}
 		__INLINE integer                  EsetLastLocalID    () const {return static_cast<integer>(links[4].size());}
+
+		integer                           LastLocalIDNum     (integer n) const {assert(n >= 0 && n < 6); return n == 5 ? 1 : static_cast<integer>(links[n].size());}
+
 		integer                           NextLocalID        (ElementType etype, integer lid) const {integer q = ElementNum(etype); ++lid; while(lid < static_cast<integer>(links[q].size()) && links[q][lid] == -1) ++lid; return lid;}
 		integer                           PrevLocalID        (ElementType etype, integer lid) const {integer q = ElementNum(etype); --lid; while(lid > 0 && links[q][lid] == -1) --lid; return lid;}
 		integer                           FirstLocalID       (ElementType etype) const;
-		integer                           LastLocalID        (integer n) const {assert(n >= 0 && n < 6); return n == 5 ? 1 : static_cast<integer>(links[n].size());}
-		integer                           LastLocalID        (ElementType etype) const {assert(OneType(etype)); return LastLocalID(ElementNum(etype));}
+		integer                           LastLocalID        (ElementType etype) const {assert(OneType(etype)); return LastLocalIDNum(ElementNum(etype));}
 
+		integer                           NextLocalIDIter    (ElementType etype, integer lid) const;
+		integer                           PrevLocalIDIter    (ElementType etype, integer lid) const;
+		integer                           FirstLocalIDIter   (ElementType etype) const;
+		integer                           LastLocalIDIter    (ElementType etype) const;
 
-		__INLINE integer                  NumberOfSets       () const { return static_cast<integer>(links[4].size() - empty_links[4].size()); }
-		__INLINE integer                  NumberOfCells      () const { return static_cast<integer>(links[3].size() - empty_links[3].size());}
-		__INLINE integer                  NumberOfFaces      () const { return static_cast<integer>(links[2].size() - empty_links[2].size()); }
-		__INLINE integer                  NumberOfEdges      () const { return static_cast<integer>(links[1].size() - empty_links[1].size()); }
-		__INLINE integer                  NumberOfNodes      () const { return static_cast<integer>(links[0].size() - empty_links[0].size()); }
+		__INLINE integer                  NumberOfSets       () const { return static_cast<integer>(links[4].size() - empty_links[4].size()) - hidden_count[4]; }
+		__INLINE integer                  NumberOfCells      () const { return static_cast<integer>(links[3].size() - empty_links[3].size()) - hidden_count[3];}
+		__INLINE integer                  NumberOfFaces      () const { return static_cast<integer>(links[2].size() - empty_links[2].size()) - hidden_count[2]; }
+		__INLINE integer                  NumberOfEdges      () const { return static_cast<integer>(links[1].size() - empty_links[1].size()) - hidden_count[1]; }
+		__INLINE integer                  NumberOfNodes      () const { return static_cast<integer>(links[0].size() - empty_links[0].size()) - hidden_count[0]; }
 		__INLINE integer                  NumberOfElements   () const { return NumberOfCells() + NumberOfFaces() + NumberOfEdges() + NumberOfNodes(); }
 		__INLINE integer                  NumberOfAll        () const { return NumberOfSets() + NumberOfElements(); }
 		integer                           NumberOf           (ElementType t) const;

Source/IO/mesh_vtk_file.cpp

@@ -26,6 +26,10 @@
 
 #define R_QUIT		  100
 
+int isnan(double x) { return x != x; }
+//int isinf(double x) { return !isnan(x) && isnan(x - x); }
+int isinf(double x) { return fabs(x) > DBL_MAX; }
+int isbad(double x) { return isnan(x) || isinf(x); }
 
 template<typename T>
 void ReadCoords(FILE * f,INMOST_DATA_REAL_TYPE c[3])
@@ -276,7 +280,7 @@ safe_output:
 									for(unsigned int m = 0; m < comps; m++) 
 									{
 										double val = static_cast<double>(arr[m]);
-										fprintf(f,"%14e ",val != val ? -0.9999E30 : val);
+										fprintf(f,"%14e ",isbad(val) ? -0.9999E30 : val);
 									}
 									fprintf(f,"\n");
 								}
@@ -295,7 +299,7 @@ safe_output:
 									for(unsigned int m = 0; m < comps; m++) 
 									{
 										double val = static_cast<double>(arr[m].GetValue());
-										fprintf(f,"%14e ",val != val ? -0.9999E30 : val);
+										fprintf(f,"%14e ",isbad(val) ? -0.9999E30 : val);
 									}
 									fprintf(f,"\n");
 								}
@@ -345,7 +349,7 @@ safe_output:
               || tags[i].GetDataType() == DATA_VARIABLE
 #endif
               ) type_str = "double";
-            fprintf(f,"SCALARS %s %s %d\n",tags[i].GetTagName().c_str(),type_str.c_str(),comps);
+						fprintf(f,"SCALARS %s %s %d\n",tags[i].GetTagName().c_str(),type_str.c_str(),comps);
 						fprintf(f,"LOOKUP_TABLE default\n");
 						for(Mesh::iteratorNode it = BeginNode(); it != EndNode(); it++)
 						{
@@ -357,7 +361,7 @@ safe_output:
 									for(unsigned int m = 0; m < comps; m++) 
 									{
 										double val = static_cast<double>(arr[m]);
-										fprintf(f,"%14e ",(val != val ? -0.9999E30 : val));
+										fprintf(f,"%14e ",(isbad(val) ? -0.9999E30 : val));
 									}
 									fprintf(f,"\n");
 								}
@@ -376,7 +380,7 @@ safe_output:
 									for(unsigned int m = 0; m < comps; m++) 
 									{
 										double val = static_cast<double>(arr[m].GetValue());
-										fprintf(f,"%14e ",(val != val ? -0.9999E30 : val));
+										fprintf(f,"%14e ",(isbad(val) ? -0.9999E30 : val));
 									}
 									fprintf(f,"\n");
 								}

Source/Mesh/iterator.cpp

@@ -41,7 +41,7 @@ namespace INMOST
 		++id;
 		while( num < 5 )
 		{
-			while(id < static_cast<integer>(links[num].size()) && links[num][id] == -1) ++id;
+			while(id < static_cast<integer>(links[num].size()) && (links[num][id] == -1)) ++id;
 			if( id == static_cast<integer>(links[num].size()) ) 
 			{
 				id = 0;
@@ -65,7 +65,7 @@ namespace INMOST
 		}
 		while( num >= 0 )
 		{
-			while(id >= 0 && links[num][id] == -1) --id;
+			while(id >= 0 && (links[num][id] == -1)) --id;
 			if( id == -1 ) 
 			{
 				--num;
@@ -83,7 +83,7 @@ namespace INMOST
 		++id;
 		while( num < 5 )
 		{
-			while(id < static_cast<integer>(links[num].size()) && links[num][id] == -1) ++id;
+			while(id < static_cast<integer>(links[num].size()) && (links[num][id] == -1)) ++id;
 			if( id == static_cast<integer>(links[num].size()) ) 
 			{
 				bool stop = true;
@@ -131,7 +131,7 @@ namespace INMOST
 		}
 		while( num >= 0 )
 		{
-			while(id >= 0 && links[num][id] == -1) --id;
+			while(id >= 0 && (links[num][id] == -1)) --id;
 			if( id == -1 ) 
 			{
 				bool stop = true;
@@ -158,7 +158,7 @@ namespace INMOST
 		assert(OneType(etype)); 
 		integer ret = 0, n = ElementNum(etype); 
 		if(n == 5) return 0; 
-		while(ret < static_cast<integer>(links[n].size()) && links[n][ret] == -1) ++ret; 
+		while(ret < static_cast<integer>(links[n].size()) && (links[n][ret] == -1)) ++ret; 
 		return ret;
 	}
 
@@ -166,7 +166,7 @@ namespace INMOST
 	{
 		integer ret = 0;
 		for(int m = 0; m < 5; m++) if( (1 << m) & t )
-			ret += static_cast<integer>(links[m].size() - empty_links[m].size());
+			ret += static_cast<integer>(links[m].size() - empty_links[m].size()) - hidden_count[m];
 		if( t & MESH ) ret++;
 		return ret;
 	}
@@ -176,7 +176,7 @@ namespace INMOST
 	{
 		while( etype != NONE ) 
 		{
-			lid = m->NextLocalID(etype,lid);
+			lid = m->NextLocalIDIter(etype,lid);
 			if( lid == m->LastLocalID(etype) )
 			{
 				etype = GetNextType(etype,types);
@@ -193,12 +193,12 @@ namespace INMOST
 	{
 		while( etype != NONE ) 
 		{
-			lid = m->PrevLocalID(etype,lid);
+			lid = m->PrevLocalIDIter(etype,lid);
 			if( lid == -1 )
 			{
 				etype = GetPrevType(etype,types);
 				if( etype ) 
-					lid = m->LastLocalID(etype);
+					lid = m->LastLocalIDIter(etype);
 				else
 				{
 					lid = -1;
@@ -236,6 +236,36 @@ namespace INMOST
 		printf("Number: %10d CurrentType %x types %x\n",lid,etype,types);
 	}
 
+	Storage::integer Mesh::NextLocalIDIter(ElementType etype, integer lid) const 
+	{
+		integer q = ElementNum(etype); 
+		++lid; 
+		while(lid < static_cast<integer>(links[q].size()) && (links[q][lid] == -1|| Hidden(ComposeHandleNum(q,lid)))) ++lid; 
+		return lid;
+	}
+
+	Storage::integer Mesh::PrevLocalIDIter(ElementType etype, integer lid) const 
+	{
+		integer q = ElementNum(etype); 
+		--lid; 
+		while(lid > 0 && (links[q][lid] == -1 || Hidden(ComposeHandleNum(q,lid)))) --lid; 
+		return lid;
+	}
+	Storage::integer Mesh::FirstLocalIDIter(ElementType etype) const
+	{
+		assert(OneType(etype)); 
+		integer ret = 0, n = ElementNum(etype); 
+		if(n == 5) return 0; 
+		while(ret < static_cast<integer>(links[n].size()) && (links[n][ret] == -1|| Hidden(ComposeHandleNum(n,ret)))) ++ret; 
+		return ret;
+	}
+	Storage::integer Mesh::LastLocalIDIter(ElementType etype) const 
+	{
+		assert(OneType(etype)); 
+		int ret = LastLocalIDNum(ElementNum(etype));
+		return PrevLocalIDIter(etype,ret);
+	}
+
 	Mesh::iteratorStorage Mesh::Begin(ElementType Types)        {return base_iterator<Storage>(Types,this,false);}
 	Mesh::iteratorStorage Mesh::End()                           {return base_iterator<Storage>(this);}
 	Mesh::iteratorElement Mesh::BeginElement(ElementType Types) {return base_iterator<Element>(Types & (NODE | EDGE | FACE | CELL | ESET),this,false);}

Source/Mesh/mesh.cpp

@@ -104,6 +104,9 @@ namespace INMOST
 		errorset = 0;
 		new_element = hide_element = 0;
 
+		memset(hidden_count,0,sizeof(integer)*6);
+		memset(hidden_count_zero,0,sizeof(integer)*6);
+
 		memset(remember,0,sizeof(remember));
 		tag_coords        = CreateTag("PROTECTED_COORD",DATA_REAL, NODE,NONE,dim);
 		tag_high_conn     = CreateTag("PROTECTED_HIGH_CONN",DATA_REFERENCE,ESET|CELL|FACE|EDGE|NODE,NONE);
@@ -1832,7 +1835,7 @@ namespace INMOST
 #ifndef NDEBUG
 		for(int etypenum = 0; etypenum < ElementNum(MESH); ++etypenum)
 		{
-			integer end = LastLocalID(etypenum);
+			integer end = LastLocalIDNum(etypenum);
 			for(integer id = 0; id < end; ++id) 
 				if( isValidElement(etypenum,id) )
 					assert((static_cast<const bulk *>(MGetDenseLink(etypenum,id,MarkersTag()))[n >> MarkerShift] & static_cast<bulk>(n & MarkerMask)) == 0 && "marker was not properly cleared from elements");
@@ -1853,7 +1856,7 @@ namespace INMOST
 #endif
 		for(int etypenum = 0; etypenum < ElementNum(MESH); ++etypenum)
 		{
-			integer end = LastLocalID(etypenum);
+			integer end = LastLocalIDNum(etypenum);
 			for(integer id = 0; id < end; ++id) 
 				if( isValidElement(etypenum,id) )
 					assert((static_cast<const bulk *>(MGetDenseLink(etypenum,id,tag_private_marker[thread]))[n >> MarkerShift] & static_cast<bulk>(n & MarkerMask)) == 0 && "marker was not properly cleared from elements");

Source/Mesh/modify.cpp

@@ -1958,6 +1958,11 @@ public:
 		hide_element = new_element;
 		new_element = temp;
 
+		integer tmp[6];
+		memcpy(tmp,hidden_count,sizeof(integer)*6);
+		memcpy(hidden_count,hidden_count_zero,sizeof(integer)*6);
+		memcpy(hidden_count_zero,tmp,sizeof(integer)*6);
+
 		for(ElementType etype = EDGE; etype <= CELL; etype = etype << 1)
 		{
 			for(integer it = 0; it < LastLocalID(etype); ++it) if( isValidElement(etype,it) )
@@ -2022,8 +2027,11 @@ public:
 		//but may use by retriving lower_bound/higher_bound O(log(n)) operations to narrow performed operations
 		erase->BulkDF(SetComparatorTag()) = ElementSet::HANDLE_COMPARATOR;
 		*/
+
+		//for(integer jt = 0; jt < LastLocalID(ESET); ++jt) if( isValidElementSet(jt) )
 		for(Mesh::iteratorSet it = BeginSet(); it != EndSet(); it++)
 		{
+			//ElementSet it = EsetByLocalID(jt);
 			std::cout << "set name: " << it->GetName() << " size " << it->Size() << " id " << it->LocalID() << std::endl;
 			if( it->HaveParent() && it->GetParent()->Hidden() )
 				it->GetParent()->RemChild(it->self());
@@ -2073,6 +2081,8 @@ public:
 					Destroy(h);
 			}
 		}
+		memset(hidden_count,0,sizeof(integer)*6);
+		memset(hidden_count_zero,0,sizeof(integer)*6);
 		ReleaseMarker(temp_hide_element);
 		ReleaseMarker(new_element);
 		new_element = 0;
@@ -2090,6 +2100,7 @@ public:
 		assert(isValidHandleRange(h));
 		assert(isValidHandle(h));
 		ElementType htype = GetHandleElementType(h);
+		if( Hidden(h) ) Show(h);
 		if( htype & (NODE|EDGE|FACE|CELL) )
 			ElementByHandle(h).Disconnect(true);
 		else if( htype & ESET )
@@ -2109,7 +2120,11 @@ public:
 	{
 		if(HideMarker()) 
 		{
-			SetMarker(h,HideMarker()); 
+			if( !Hidden(h) )
+			{
+				hidden_count[GetHandleElementNum(h)]++;
+				SetMarker(h,HideMarker()); 
+			}
 			return true;
 		} 
 		return false;
@@ -2119,11 +2134,12 @@ public:
 	{
 		if(HideMarker()) 
 		{
-			if( GetMarker(h,HideMarker()) )
+			if( Hidden(h) )
 			{
+				hidden_count[GetHandleElementNum(h)]--;
 				RemMarker(h,HideMarker()); 
-				return true;
 			}
+			return true;
 		} 
 		return false;
 	}
@@ -2135,7 +2151,7 @@ public:
 			if( GetHandleElementType(h) != CELL ) //mark all elements that rely on this that they should be deleted
 			{
 				Element::adj_type & hc = HighConn(h);
-				for(Element::adj_type::size_type it = 0; it < hc.size(); ++it) Delete(hc[it]);
+				for(Element::adj_type::size_type it = 0; it < hc.size(); ++it) if( !Hidden(hc[it]) ) Hide(hc[it]);
 			}
 			return false;
 		}