Adding GridTools examples

Various helpful tools for grids.

Examples/CMakeLists.txt

@@ -3,6 +3,7 @@ add_subdirectory(AdaptiveMesh)
 add_subdirectory(DrawGrid)
 add_subdirectory(OldDrawGrid)
 add_subdirectory(GridGen)
+add_subdirectory(GridTools)
 endif(USE_MESH)
 if(USE_SOLVER)
 add_subdirectory(DrawMatrix)

Examples/GridTools/CMakeLists.txt

+project(GridGen)
+
+add_executable(FixFaults fix_faults.cpp)
+add_executable(FixTiny fix_tiny.cpp)
+add_executable(UniteFaces unite_faces.cpp)
+add_executable(Dual dual.cpp)
+add_executable(Slice slice.cpp)
+add_executable(SliceFunc slice_func.cpp)
+
+target_link_libraries(FixFaults inmost)
+if(USE_MPI)
+  message("linking FixFaults with MPI")
+  target_link_libraries(FixFaults ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(FixFaults PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+install(TARGETS FixFaults EXPORT inmost-targets RUNTIME DESTINATION bin)
+
+target_link_libraries(FixTiny inmost)
+if(USE_MPI)
+  message("linking FixTiny with MPI")
+  target_link_libraries(FixTiny ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(FixTiny PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+install(TARGETS FixTiny EXPORT inmost-targets RUNTIME DESTINATION bin)
+
+target_link_libraries(UniteFaces inmost)
+if(USE_MPI)
+  message("linking UniteFaces with MPI")
+  target_link_libraries(UniteFaces ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(UniteFaces PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+install(TARGETS UniteFaces EXPORT inmost-targets RUNTIME DESTINATION bin)
+
+
+target_link_libraries(Dual inmost)
+if(USE_MPI)
+  message("linking Dual with MPI")
+  target_link_libraries(Dual ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(Dual PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+install(TARGETS Dual EXPORT inmost-targets RUNTIME DESTINATION bin)
+
+
+target_link_libraries(Slice inmost)
+if(USE_MPI)
+  message("linking Slice with MPI")
+  target_link_libraries(Slice ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(Slice PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+install(TARGETS Slice EXPORT inmost-targets RUNTIME DESTINATION bin)
+
+
+target_link_libraries(SliceFunc inmost)
+if(USE_MPI)
+  message("linking SliceFunc with MPI")
+  target_link_libraries(SliceFunc ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(SliceFunc PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+install(TARGETS SliceFunc EXPORT inmost-targets RUNTIME DESTINATION bin)
+
+
+
+
+

Examples/GridTools/dual.cpp

+#include <stdio.h>
+#include "inmost.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 A("orig"),B("voro");
+	A.Load(argv[1]);
+	
+	Mesh::GeomParam table;
+	table[ORIENTATION] = FACE;
+	table[CENTROID] = FACE | CELL;
+	A.PrepareGeometricData(table);
+	
+	
+	//create nodes of the voronoi mesh
+	Tag cell2node = A.CreateTag("CELL2NODE",DATA_REMOTE_REFERENCE,CELL|FACE|EDGE|NODE,FACE|EDGE|NODE,1);
+	for(Mesh::iteratorCell it = A.BeginCell(); it != A.EndCell(); ++it)
+	{
+		real cnt[3];
+		it->Centroid(cnt);
+		it->RemoteReference(cell2node) = RemoteHandleType(&B,B.CreateNode(cnt).GetHandle());
+	}
+	for(Mesh::iteratorFace it = A.BeginFace(); it != A.EndFace(); ++it) if( it->Boundary() )
+	{
+		real cnt[3];
+		it->Centroid(cnt);
+		it->RemoteReference(cell2node) = RemoteHandleType(&B,B.CreateNode(cnt).GetHandle());
+	}
+	//add corner nodes
+	int corners = 0;
+	MarkerType corner = A.CreateMarker();
+	for(Mesh::iteratorEdge it = A.BeginEdge(); it != A.EndEdge(); ++it) if( it->Boundary() )
+	{
+		real nrm[2][3], dot, cnt[3];
+		int nbndfaces = 0;
+		ElementArray<Face> faces = it->getFaces();
+		for(ElementArray<Face>::iterator kt = faces.begin(); kt != faces.end(); ++kt) if( kt->Boundary() )
+			kt->UnitNormal(nrm[nbndfaces++]);
+		assert(nbndfaces==2);
+		dot = 0;
+		for(int k = 0; k < 3; ++k) dot += nrm[0][k]*nrm[1][k];
+		if( dot < 0.95 )
+		{
+			corners++;
+			it->SetMarker(corner);
+			it->Centroid(cnt);
+			it->RemoteReference(cell2node) = RemoteHandleType(&B,B.CreateNode(cnt).GetHandle());
+		}
+	}
+	for(Mesh::iteratorNode it = A.BeginNode(); it != A.EndNode(); ++it) if( it->Boundary() )
+	{
+		real cnt[3];
+		int ncorners = it->nbAdjElements(EDGE,corner);
+		if( ncorners > 2 )
+		{
+			corners++;
+			it->SetMarker(corner);
+			it->Centroid(cnt);
+			it->RemoteReference(cell2node) = RemoteHandleType(&B,B.CreateNode(cnt).GetHandle());
+		}
+		else if( ncorners == 2 )
+		{
+			ElementArray<Edge> edges = it->getEdges(corner);
+			Node n1 = edges[0]->getBeg() == it->self() ? edges[0]->getEnd() : edges[0]->getBeg();
+			Node n2 = edges[1]->getBeg() == it->self() ? edges[1]->getEnd() : edges[1]->getBeg();
+			real cnt1[3], cnt2[3], r1[3],r2[3], l, dot;
+			it->Centroid(cnt);
+			n1->Centroid(cnt1);
+			n2->Centroid(cnt2);
+			r1[0] = cnt[0] - cnt1[0];
+			r1[1] = cnt[1] - cnt1[1];
+			r1[2] = cnt[2] - cnt1[2];
+			l = sqrt(r1[0]*r1[0]+r1[1]*r1[1]+r1[2]*r1[2]);
+			r1[0] /= l;
+			r1[1] /= l;
+			r1[2] /= l;
+			r2[0] = cnt2[0] - cnt1[0];
+			r2[1] = cnt2[1] - cnt1[1];
+			r2[2] = cnt2[2] - cnt1[2];
+			l = sqrt(r2[0]*r2[0]+r2[1]*r2[1]+r2[2]*r2[2]);
+			r2[0] /= l;
+			r2[1] /= l;
+			r2[2] /= l;
+			dot = 0;
+			for(int k = 0; k < 3; ++k) dot += r1[k]*r2[k];
+			if( dot < 0.95 )
+			{
+				corners++;
+				it->SetMarker(corner);
+				it->Centroid(cnt);
+				it->RemoteReference(cell2node) = RemoteHandleType(&B,B.CreateNode(cnt).GetHandle());
+			}
+		}
+	}
+	std::cout << "corner nodes: " << corners << std::endl;
+	//create edges of the voronoi mesh
+	Tag face2edge = A.CreateTag("FACE2EDGE",DATA_REMOTE_REFERENCE,FACE|EDGE,EDGE,1);
+	Tag face2edge_corner = A.CreateTag("FACE2EDGE_CORNER",DATA_REMOTE_REFERENCE,EDGE,EDGE,2);
+	Tag face2edge_corner_node = A.CreateTag("FACE2EDGE_CORNER_NODE",DATA_REMOTE_REFERENCE,NODE,NODE);
+	ElementArray<Node> edgenodes(&B,2);
+	for(Mesh::iteratorFace it = A.BeginFace(); it != A.EndFace(); ++it)
+	{
+		Cell c1 = it->BackCell();
+		Cell c2 = it->FrontCell();
+		if( c2.isValid() )
+		{
+			edgenodes[0] = MakeElement(c1.RemoteReference(cell2node)).getAsNode();
+			edgenodes[1] = MakeElement(c2.RemoteReference(cell2node)).getAsNode();
+		}
+		else
+		{
+			edgenodes[0] = MakeElement(c1.RemoteReference(cell2node)).getAsNode();
+			edgenodes[1] = MakeElement(it->RemoteReference(cell2node)).getAsNode();
+		}
+		it->RemoteReference(face2edge) = RemoteHandleType(&B,B.CreateEdge(edgenodes).first.GetHandle());
+	}
+	int process_corners = 0;
+	for(Mesh::iteratorEdge it = A.BeginEdge(); it != A.EndEdge(); ++it) if( it->Boundary() )
+	{
+		ElementArray<Face> faces = it->getFaces();
+		ElementArray<Face> bndfaces(&A);
+		for(ElementArray<Face>::iterator kt = faces.begin(); kt != faces.end(); ++kt) if( kt->Boundary() )
+			bndfaces.push_back(kt->self());
+		assert(bndfaces.size() == 2);
+		
+		if( it->GetMarker(corner) )
+		{
+			process_corners++;
+			edgenodes[0] = MakeElement(bndfaces[0].RemoteReference(cell2node)).getAsNode();
+			edgenodes[1] = MakeElement(it->RemoteReference(cell2node)).getAsNode();
+			
+			it->RemoteReferenceArray(face2edge_corner).at(0) = RemoteHandleType(&B,B.CreateEdge(edgenodes).first.GetHandle());
+			
+			edgenodes[0] = MakeElement(bndfaces[1].RemoteReference(cell2node)).getAsNode();
+			edgenodes[1] = MakeElement(it->RemoteReference(cell2node)).getAsNode();
+			
+			it->RemoteReferenceArray(face2edge_corner).at(1) = RemoteHandleType(&B,B.CreateEdge(edgenodes).first.GetHandle());
+		}
+		else
+		{
+			edgenodes[0] = MakeElement(bndfaces[0].RemoteReference(cell2node)).getAsNode();
+			edgenodes[1] = MakeElement(bndfaces[1].RemoteReference(cell2node)).getAsNode();
+		
+			it->RemoteReference(face2edge) = RemoteHandleType(&B,B.CreateEdge(edgenodes).first.GetHandle());
+		}
+	}
+	//add edges connecting corners
+	for(Mesh::iteratorNode it = A.BeginNode(); it != A.EndNode(); ++it) if( it->Boundary() )
+	{
+		ElementArray<Edge> edges = it->getEdges(corner);
+		if( it->GetMarker(corner) ) //connect each corner node
+		{
+			edgenodes[0] = MakeElement(it->RemoteReference(cell2node)).getAsNode();
+			for(ElementArray<Edge>::iterator kt = edges.begin(); kt != edges.end(); ++kt)
+			{
+				process_corners++;
+				edgenodes[1] = MakeElement(kt->RemoteReference(cell2node)).getAsNode();
+				it->RemoteReferenceArray(face2edge_corner_node).push_back(&B,B.CreateEdge(edgenodes).first.GetHandle());
+			}
+		}
+		else if( edges.size() == 2 ) //connect two corner nodes
+		{
+			process_corners++;
+			edgenodes[0] = MakeElement(edges[0]->RemoteReference(cell2node)).getAsNode();
+			edgenodes[1] = MakeElement(edges[1]->RemoteReference(cell2node)).getAsNode();
+			
+			it->RemoteReferenceArray(face2edge_corner_node).push_back(&B,B.CreateEdge(edgenodes).first.GetHandle());
+		}
+	}
+	std::cout << "corners processed: " << process_corners << std::endl;
+	//create faces of the voronoi mesh
+	Tag edge2face = A.CreateTag("EDGE2FACE",DATA_REMOTE_REFERENCE,EDGE|NODE,NODE,1);
+	ElementArray<Edge> faceedges(&B);
+	process_corners = 0;
+	for(Mesh::iteratorEdge it = A.BeginEdge(); it != A.EndEdge(); ++it)
+	{
+		ElementArray<Face> faces = it->getFaces();
+		ElementArray<Node> nodes = it->getNodes();
+		for(ElementArray<Face>::iterator kt = faces.begin(); kt != faces.end(); ++kt)
+			faceedges.push_back(MakeElement(kt->RemoteReference(face2edge)));
+		
+		if( it->GetMarker(corner) )
+		{
+			process_corners++;
+			faceedges.push_back(it->RemoteReferenceArray(face2edge_corner)[0]);
+			faceedges.push_back(it->RemoteReferenceArray(face2edge_corner)[1]);
+		}
+		else if( it->Boundary() )
+			faceedges.push_back(MakeElement(it->RemoteReference(face2edge)));
+		
+		Face f = B.CreateFace(faceedges).first;
+		f.FixEdgeOrder();
+		it->RemoteReference(edge2face) = RemoteHandleType(&B,B.CreateFace(faceedges).first.GetHandle());
+		faceedges.clear();
+	}
+	//construct boundary faces without corner nodes
+	for(Mesh::iteratorNode it = A.BeginNode(); it != A.EndNode(); ++it) if( it->Boundary() )
+	{
+		ElementArray<Edge> edges = it->getEdges();
+		for(ElementArray<Edge>::iterator kt = edges.begin(); kt != edges.end(); ++kt) if( kt->Boundary() )
+		{
+			if( kt->GetMarker(corner) )
+			{
+				process_corners++;
+				faceedges.push_back(kt->RemoteReferenceArray(face2edge_corner)[0]);
+				faceedges.push_back(kt->RemoteReferenceArray(face2edge_corner)[1]);
+			}
+			else
+				faceedges.push_back(MakeElement(kt->RemoteReference(face2edge)));
+		}
+		Face f = B.CreateFace(faceedges).first;
+		f.FixEdgeOrder();
+		it->RemoteReference(edge2face) = RemoteHandleType(&B,B.CreateFace(faceedges).first.GetHandle());
+		faceedges.clear();
+	}
+	std::cout << "corners processed: " << process_corners << std::endl;
+	//construct cells
+	ElementArray<Face> cellfaces(&B);
+	for(Mesh::iteratorNode it = A.BeginNode(); it != A.EndNode(); ++it)
+	{
+		ElementArray<Edge> edges = it->getEdges();
+		for(ElementArray<Edge>::iterator kt = edges.begin(); kt != edges.end(); ++kt)
+			cellfaces.push_back(MakeElement(kt->RemoteReference(edge2face)));
+		if( it->Boundary() )
+			cellfaces.push_back(MakeElement(it->RemoteReference(edge2face)));
+		B.CreateCell(cellfaces);
+		cellfaces.clear();
+	}
+	
+	//split corner faces by edges connected to corner nodes
+	process_corners= 0;
+	for(Mesh::iteratorNode it = A.BeginNode(); it != A.EndNode(); ++it) if( it->HaveData(face2edge_corner_node) )
+	{
+		process_corners++;
+		Storage::remote_reference_array arr = it->RemoteReferenceArray(face2edge_corner_node);
+		ElementArray<Edge> edges(&B);
+		for(Storage::remote_reference_array::iterator kt = arr.begin(); kt != arr.end(); ++kt)
+			edges.push_back(kt->GetHandle());
+		Face f = MakeElement(it->RemoteReference(edge2face)).getAsFace();
+		Face::SplitFace(f,edges,0);
+	}
+	std::cout << "corners splitted: " << process_corners << std::endl;
+
+	if( A.HaveTag("GRIDNAME") )
+	{
+		Storage::bulk_array nameA = A.self().BulkArray(A.GetTag("GRIDNAME"));
+		Storage::bulk_array nameB = B.self().BulkArray(B.CreateTag("GRIDNAME",DATA_BULK,MESH,NONE));
+		std::string ins = "_dual";
+		nameB.insert(nameB.end(),nameA.begin(),nameA.end());
+		nameB.insert(nameB.end(),ins.c_str(),ins.c_str()+ins.size());
+	}
+	
+	if( argc > 2 )
+	{
+		std::cout << "Save to " << argv[2] << std::endl;
+		B.Save(argv[2]);
+	}
+	else
+	{
+		std::cout << "Save to out.vtk" << std::endl;
+		B.Save("out.vtk");
+	}
+
+	return 0;
+}
\ No newline at end of file

Examples/GridTools/fix_faults.cpp

+#include "inmost.h"
+
+using namespace INMOST;
+
+//todo: non-flat faces
+
+static void GetBox(Element e, float minmax[6])
+{
+	minmax[0] = minmax[2] = minmax[4] = 1.0e20f; //min
+	minmax[1] = minmax[3] = minmax[5] = -1.0e20f; //max
+	ElementArray<Node> nodes = e->getNodes();
+	for (ElementArray<Node>::iterator it = nodes.begin(); it != nodes.end(); it++)
+	{
+		Storage::real_array c = it->Coords();
+		for (int i = 0; i < (int)c.size(); i++)
+		{
+			float v = static_cast<float>(c[i]);
+			if (minmax[1 + 2 * i] < v) minmax[1 + 2 * i] = v; //max
+			if (minmax[0 + 2 * i] > v) minmax[0 + 2 * i] = v; //min
+		}
+	}
+	for (int i = 0; i < 3; ++i)
+	{
+		if (minmax[1 + 2 * i] < minmax[0 + 2 * i])
+		{
+			minmax[1 + 2 * i] = 1.0e-5f;
+			minmax[0 + 2 * i] = -1.0e-5f;
+		}
+		else if (minmax[1 + 2 * i] == minmax[0 + 2 * i])
+		{
+			minmax[1 + 2 * i] += 1.0e-5f;
+			minmax[0 + 2 * i] += -1.0e-5f;
+		}
+	}
+}
+
+static void MergeBox(const float bboxa[6], const float bboxb[6], float bbox[6])
+{
+	for (int k = 0; k < 3; k++)
+	{
+		bbox[0 + 2 * k] = std::min(bboxa[0 + 2 * k], bboxb[0 + 2 * k]);
+		bbox[1 + 2 * k] = std::max(bboxa[1 + 2 * k], bboxb[1 + 2 * k]);
+	}
+}
+
+static bool IntersectBox(const float bboxa[6], const float bboxb[6])
+{
+	for (int k = 0; k < 3; ++k)
+	{
+		float la = bboxa[0 + 2 * k], ra = bboxa[1 + 2 * k];
+		float lb = bboxa[0 + 2 * k], rb = bboxb[1 + 2 * k];
+		if (ra < lb || la > rb) return false;
+	}
+	return true;
+}
+
+static bool MatchPoints(const double v1[3], const double v2[3])
+{
+	double l = 0;
+	for (int k = 0; k < 3; ++k)
+		l += (v1[k] - v2[k])*(v1[k] - v2[k]);
+	return sqrt(l) < 1.0e-7;
+}
+
+//returns distance between line if shortest line is within segments, writes position of distance into last argument
+static double SegmentDistance(const double v1[3], const double v2[3], const double v3[3], const double v4[3], double vout[3])
+{
+	double v13[3], v21[3], v43[3];
+	for (int k = 0; k < 3; ++k)
+	{
+		v13[k] = v1[k] - v3[k];
+		v21[k] = v2[k] - v1[k];
+		v43[k] = v4[k] - v3[k];
+	}
+	double d1321 = 0, d2121 = 0, d4321 = 0, d1343 = 0, d4343 = 0;
+	for (int k = 0; k < 3; ++k)
+	{
+		d1321 += v13[k] * v21[k];
+		d2121 += v21[k] * v21[k];
+		d4321 += v43[k] * v21[k];
+		d1343 += v13[k] * v43[k];
+		d4343 += v43[k] * v43[k];
+	}
+	double mu1 = 0, mu2 = 0;
+	/*
+	double parallel = 0;
+	for (int k = 0; k < 3; ++k)
+	parallel += v21[k] * v43[k];
+	parallel = fabs(parallel);
+	parallel /= sqrt(d2121)*sqrt(d4343);
+	if (fabs(parallel - 1.0) > 1.0e-6)
+	*/
+	if (fabs(d2121*d4343 - d4321*d4321) > 1.0e-6)
+	{
+		mu1 = (d1343*d4321 - d1321*d4343) / (d2121*d4343 - d4321*d4321);
+		mu2 = (d1343 + mu1*d4321) / d4343;
+		if (mu1 > 1) mu1 = 1;
+		if (mu1 < 0) mu1 = 0;
+		if (mu2 > 1) mu2 = 1;
+		if (mu2 < 0) mu2 = 0;
+	}
+	else //parallel lines
+	{
+		//1d problem
+		double la = 0, ra = 0;
+		double lb = 0, rb = 0;
+		for (int k = 0; k < 3; ++k)
+		{
+			la += v1[k] * v21[k];
+			ra += v2[k] * v21[k];
+			lb += v3[k] * v21[k];
+			rb += v4[k] * v21[k];
+		}
+		bool invb = false;
+		if (lb > rb)
+		{
+			std::swap(lb, rb);
+			invb = true;
+		}
+
+		if (ra < lb) // (la,ra) is to the left of (lb,rb), no intersection
+		{
+			mu1 = 1;
+			mu2 = 0;
+		}
+		else if (la > rb) // (la,ra) is to the right of (lb,rb), no intersection
+		{
+			mu1 = 0;
+			mu2 = 1;
+		}
+		else if (lb > la) // (lb,rb) intersects to the right of or contains (la,ra)
+		{
+			mu2 = 0;
+			mu1 = (lb - la) / (ra - la);
+		}
+		else if (rb < ra) // (lb,rb) intersects to the left of or contains (la,ra)
+		{
+			mu2 = 1;
+			mu1 = (rb - la) / (ra - la);
+		}
+		else // (lb,rb) contains (la,ra)
+		{
+			mu1 = 0;
+			mu2 = (la - lb) / (rb - lb);
+		}
+		if (invb) mu2 = 1 - mu2;
+	}
+	double vs1[3], vs2[3], h = 0;
+	for (int k = 0; k < 3; ++k)
+	{
+		vs1[k] = v1[k] + mu1*(v2[k] - v1[k]);
+		vs2[k] = v3[k] + mu2*(v4[k] - v3[k]);
+		vout[k] = (vs1[k] + vs2[k])*0.5;
+		h += (vs2[k] - vs1[k])*(vs2[k] - vs1[k]);
+	}
+	return sqrt(h);
+}
+
+
+class kdtree
+{
+	struct entry
+	{
+		HandleType e;
+		float xyz[3];
+		struct entry & operator =(const struct entry & other)
+		{
+			e = other.e;
+			xyz[0] = other.xyz[0];
+			xyz[1] = other.xyz[1];
+			xyz[2] = other.xyz[2];
+			return *this;
+		}
+	} *set;
+	Mesh * m;
+	INMOST_DATA_ENUM_TYPE size;
+	float bbox[6];
+	kdtree * children;
+	inline static unsigned int flip(const unsigned int * fp)
+	{
+		unsigned int mask = -((int)(*fp >> 31)) | 0x80000000;
+		return *fp ^ mask;
+	}
+	inline static unsigned int _0(unsigned int x)	{ return x & 0x7FF; }
+	inline static unsigned int _1(unsigned int x)	{ return x >> 11 & 0x7FF; }
+	inline static unsigned int _2(unsigned int x)   { return x >> 22; }
+	template<int pos>
+	inline static int cmpElements(const void * a, const void * b)
+	{
+		const entry * ea = ((const entry *)a);
+		const entry * eb = ((const entry *)b);
+		float ad = ea->xyz[pos];
+		float bd = eb->xyz[pos];
+		return (ad > bd) - (ad < bd);
+	}
+	void radix_sort(int dim, struct entry * temp)
+	{
+		unsigned int i;
+		const unsigned int kHist = 2048;
+		unsigned int  b0[kHist * 3];
+		unsigned int *b1 = b0 + kHist;
+		unsigned int *b2 = b1 + kHist;
+		memset(b0, 0, sizeof(unsigned int)*kHist * 3);
+		for (i = 0; i < size; i++)
+		{
+			unsigned int fi = flip((unsigned int *)&set[i].xyz[dim]);
+			++b0[_0(fi)]; ++b1[_1(fi)]; ++b2[_2(fi)];
+		}
+		{
+			unsigned int sum0 = 0, sum1 = 0, sum2 = 0;
+			for (i = 0; i < kHist; i++)
+			{
+				b0[kHist - 1] = b0[i] + sum0; b0[i] = sum0 - 1; sum0 = b0[kHist - 1];
+				b1[kHist - 1] = b1[i] + sum1; b1[i] = sum1 - 1; sum1 = b1[kHist - 1];
+				b2[kHist - 1] = b2[i] + sum2; b2[i] = sum2 - 1; sum2 = b2[kHist - 1];
+			}
+		}
+		for (i = 0; i < size; i++) temp[++b0[_0(flip((unsigned int *)&set[i].xyz[dim]))]] = set[i];
+		for (i = 0; i < size; i++) set[++b1[_1(flip((unsigned int *)&temp[i].xyz[dim]))]] = temp[i];