Add Examples/GridTools/CubeTransform performing bilinear transformation of...

Add Examples/GridTools/CubeTransform performing bilinear transformation of cube in space; add Examples/TestSuite/test_ns_turek for Naver-Stokes Schafer-Turek tests; correct Examples/GridTools/Sector

Examples/ADMFD/diffusion.cpp

@@ -430,6 +430,8 @@ int main(int argc,char ** argv)
             m->Save("out.vtk");
         else
             m->Save("out.pvtk");
+            
+        m->Save("solution.pmf");
 
         delete m; //clean up the mesh
     }

Examples/GridTools/CMakeLists.txt

@@ -16,6 +16,7 @@ add_executable(Sew sew.cpp)
 add_executable(Scale scale.cpp)
 add_executable(Move move.cpp)
 add_executable(Convert convert.cpp)
+add_executable(CubeTransform cube_transform.cpp)
 add_executable(SameMeshDifference difference_same.cpp)
 add_executable(MatchSameMeshDifference difference_same_match.cpp)
 add_executable(MeshDifference difference_map.cpp)
@@ -264,6 +265,17 @@ endif(USE_MPI)
 install(TARGETS Convert EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
 
 
+target_link_libraries(CubeTransform inmost)
+if(USE_MPI)
+  message("linking CubeTransform with MPI")
+  target_link_libraries(CubeTransform ${MPI_LIBRARIES})
+  if(MPI_LINK_FLAGS)
+    set_target_properties(CubeTransform PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif()
+endif(USE_MPI)
+install(TARGETS CubeTransform EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
+
+
 target_link_libraries(SameMeshDifference inmost)
 if(USE_MPI)
   message("linking SameMeshDifference with MPI")

Examples/GridTools/cube_transform.cpp

+#include "inmost.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+
+
+std::string mesh_name = "_TRANSFORMED";
+
+using namespace INMOST;
+
+
+
+
+int main(int argc, char *argv[]) 
+{
+	
+	if( argc < 3 )
+	{
+		printf("Usage: %s mapping_file mesh_file [output=out.pmf] [dims=3]\n",argv[0]);
+		printf("mapping file contains 8 (for 3D) 4 (for 2D) 2 (for 1D) coordinates mapping each corner of the cube:\n");
+		for(int k = 0; k < 8; ++k)
+			printf("x%d y%d z%d\n",k,k,k);
+		printf("      (6)----------(7)\n");
+		printf("      /|           /|\n");
+		printf("     / |          / |\n");
+		printf("    /  |         /  |\n");
+		printf("  (4)----------(5)  |\n");
+		printf("   |   |        |   |\n");
+		printf("   |   |        |   |\n");
+		printf("   |  (2)-------|--(3)\n");
+		printf("   |  /         |  /\n");
+		printf("   | /          | /\n");
+		printf("   |/           |/\n");
+		printf("z (0)----------(1)\n");
+		printf("^ y\n");
+		printf("|/\n");
+		printf("*-->x\n");
+		return -1;
+	}
+	std::string mapfile(argv[1]);
+	std::string infile(argv[2]);
+	std::string outfile = "out.pmf";
+	int dims = 3;
+	if( argc > 3 ) outfile = std::string(argv[3]);
+	if( argc > 4 ) dims = atoi(argv[4]);
+	  
+	std::fstream map(mapfile.c_str(),std::ios::in);
+	double xyz[8][3] =  {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};
+	for(int k = 0; k < (1<<dims); ++k)
+		for(int d = 0; d < dims; ++d)
+			map >> xyz[k][d];
+	map.close();
+	
+	for(int k = (1<<dims); k < 8; k+=(1<<dims))
+	{
+		for(int q = 0; q < (1<<dims); ++q)
+			for(int d = 0; d < dims; ++d)
+				xyz[k+q][d] = xyz[k-(1<<dims)+q][d];
+	}
+	
+	for(int k = 0; k < 8; ++k)
+		printf("%g %g %g\n",xyz[k][0],xyz[k][1],xyz[k][2]);
+	
+	Mesh mesh; 
+	mesh.Load(infile);
+	
+	double cmin[3], cmax[3];
+	for(int d = 0; d < 3; ++d)
+	{
+		cmin[d]=1.0e+20;
+		cmax[d]=-1.0e+20;
+	}
+	
+	for(Mesh::iteratorNode n = mesh.BeginNode(); n != mesh.EndNode(); ++n)
+	{
+		for(int d = 0; d < 3; ++d)
+		{
+			if( cmin[d] > n->Coords()[d] ) cmin[d] = n->Coords()[d];
+			if( cmax[d] < n->Coords()[d] ) cmax[d] = n->Coords()[d];
+		}
+	}
+	printf("before\n");
+	for(int d = 0; d < 3; ++d)
+		printf("%d %g:%g\n",d,cmin[d],cmax[d]);
+		
+	  
+	for(Mesh::iteratorNode n = mesh.BeginNode(); n != mesh.EndNode(); ++n)
+	{
+		double c[3] = {0,0,0};
+		for(int d = 0; d < 3; ++d)
+			c[d] = (n->Coords()[d]-cmin[d])/(cmax[d]-cmin[d]);
+		for(int d = 0; d < dims; ++d)
+			n->Coords()[d] = ((xyz[0][d]*(1-c[0]) + xyz[1][d]*c[0])*(1-c[1])+
+			                  (xyz[2][d]*(1-c[0]) + xyz[3][d]*c[0])*c[1])*(1-c[2])+
+			                 ((xyz[4][d]*(1-c[0]) + xyz[5][d]*c[0])*(1-c[1])+
+			                  (xyz[6][d]*(1-c[0]) + xyz[7][d]*c[0])*c[1])*c[2];
+	} 
+	
+	for(int d = 0; d < 3; ++d)
+	{
+		cmin[d]=1.0e+20;
+		cmax[d]=-1.0e+20;
+	}
+	
+	for(Mesh::iteratorNode n = mesh.BeginNode(); n != mesh.EndNode(); ++n)
+	{
+		for(int d = 0; d < 3; ++d)
+		{
+			if( cmin[d] > n->Coords()[d] ) cmin[d] = n->Coords()[d];
+			if( cmax[d] < n->Coords()[d] ) cmax[d] = n->Coords()[d];
+		}
+	}
+	
+	printf("after\n");
+	for(int d = 0; d < 3; ++d)
+		printf("%d %g:%g\n",d,cmin[d],cmax[d]);
+		
+	
+	Storage::bulk_array name = mesh->self()->BulkArray(mesh.CreateTag("GRIDNAME",DATA_BULK,MESH,NONE));
+	name.insert(name.end(),mesh_name.begin(),mesh_name.end());
+	
+	  
+	printf("I'm ready!\n");
+	mesh.Save(outfile);
+	  
+	
+	return 0;
+}

Examples/GridTools/sector.cpp

@@ -74,12 +74,15 @@ int main(int argc, char ** argv)
 		else if( stype == 1 )
 		{
 			mx += 1 + inner;
-			my += x*sin(alpha)*sqrt(2);
+			//~ my += x*sin(alpha)*sqrt(2);
+			my = 3*y-1 - (1-x)*(2*y-1);
 		}
 		else if( stype == 2 )
 		{
 			mx += 1;
-			my += x*sin(alpha)*sqrt(2);
+			//~ my += x*sin(alpha)*sqrt(2);
+			//~ my += x*(2*y-1);
+			my = 3*y-1 - (1-x)*(2*y-1);
 		}
 		else if( stype == 3 )
 		{

Examples/TestSuite/CMakeLists.txt

@@ -31,7 +31,7 @@ add_executable(test_ns_linear ns_linear.cpp)
 add_executable(test_ns_cavity ns_cavity.cpp)
 add_executable(test_ns_pousielle ns_pousielle.cpp)
 add_executable(test_ns_backward_step ns_backward_step.cpp)
-
+add_executable(test_ns_turek ns_turek.cpp)
 
 
 target_link_libraries(test_ns_linear inmost)
@@ -62,7 +62,17 @@ if(USE_MPI)
   endif() 
 endif(USE_MPI)
 install(TARGETS test_ns_pousielle EXPORT projects-targets RUNTIME DESTINATION bin/TestSuite)
-
+
+
+target_link_libraries(test_ns_turek inmost)
+if(USE_MPI)
+  message("linking test_ns_cavity with MPI")
+  target_link_libraries(test_ns_turek ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(test_ns_turek PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+install(TARGETS test_ns_turek EXPORT projects-targets RUNTIME DESTINATION bin/TestSuite)
 
 target_link_libraries(test_ns_backward_step inmost)
 if(USE_MPI)

Examples/TestSuite/ns_turek.cpp

+#include "inmost.h"
+
+using namespace INMOST;
+
+// use ADMFD or ADVDIFF to solve for phi
+
+//setup BC
+//setup on mesh:
+//
+// FORCE - 3 entries
+//
+// BOUNDARY_CONDITION_VELOCITY - 7 entries
+// r = (a5,a6,a7)
+// n.(a1*u + a2*t) = n.r
+// (I - nn.)(a3*u + a4*t) = (I-nn.)r
+//
+// BOUNDARY_CONDITION_PRESSURE - 1 entry
+//
+// REFERENCE_VELOCITY - 3 entries
+//
+// REFERENCE_PRESSURE - 1 entry
+//
+// BOUNDARY_CONDTION - 3 entries (for calculation of phi and measures of drag and lift)
+
+
+typedef Storage::real real;
+const double eps = 1.0e-5;
+
+
+int main(int argc, char ** argv)
+{
+	
+	if( argc < 2 )
+	{
+		std::cout << "Usage: " << argv[0] << " mesh [mesh_out=grid_out.pmf] [Umax=2.25]" << std::endl;
+		return 0;
+	}
+	
+	Mesh * m = new Mesh;
+	m->SetFileOption("VERBOSITY","2");
+	try
+	{
+		m->Load(argv[1]);
+	}
+	catch(...)
+	{
+		std::cout << "Cannot load the mesh " << argv[1] << std::endl;
+		return -1;
+	}
+	
+	std::string fout = "grid_out.pmf";
+	double Umax = 2.25;
+	if( argc > 2 ) fout = std::string(argv[2]);
+	if( argc > 3 ) Umax = atof(argv[3]);
+	
+	
+	
+	
+	double cmax[3] = {-1.0e20,-1.0e20,-1.0e20}, cmin[3] = {1.0e20,1.0e20,1.0e20};
+	for(Mesh::iteratorNode n = m->BeginNode(); n != m->EndNode(); ++n)
+	{
+		Storage::real_array c = n->Coords();
+		for(int k = 0; k < 3; ++k)
+		{
+			if( cmax[k] < c[k] ) cmax[k] = c[k];
+			if( cmin[k] > c[k] ) cmin[k] = c[k];
+		}
+	}
+	for(int d = 0; d < 3; ++d)
+		std::cout << d << " " << cmin[d] << ":" << cmax[d] << std::endl;
+			
+	
+	
+	//TagRealArray force = m->CreateTag("FORCE",DATA_REAL,CELL,NONE,3);
+	TagRealArray  bc      = m->CreateTag("BOUNDARY_CONDITION_VELOCITY",DATA_REAL,FACE,FACE,7);
+	TagRealArray  bcphi   = m->CreateTag("BOUNDARY_CONDITION",DATA_REAL,FACE,FACE,3);
+	TagReal       bcp     = m->CreateTag("BOUNDARY_CONDITION_PRESSURE",DATA_REAL,FACE,FACE,1);
+	TagRealArray  uvw     = m->CreateTag("UVW",DATA_REAL,CELL,NONE,3);
+	TagReal       p       = m->CreateTag("P",DATA_REAL,CELL,NONE,1);
+	
+	
+	
+	//this should not be needed?
+	for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->Boundary() )
+		it->FixNormalOrientation();
+	
+	
+	
+	
+	for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->Boundary() )
+	{
+		double  n[3], c[3];
+		it->UnitNormal(n);
+		it->Barycenter(c);
+		bcphi[*it][0] = 1;
+		bcphi[*it][1] = 0;
+		bcphi[*it][2] = 0;
+		if( fabs(n[0]-1) < 1.0e-3 && c[0] > 2.5-eps) // outflow
+		{
+			bcp[*it] = 0;
+		}
+		else if(  fabs(n[0]+1) < 1.0e-3 && c[0] < 0.0+eps) //inflow
+		{
+			bc[*it][0] = 1;
+			bc[*it][1] = 0;
+			bc[*it][2] = 1;
+			bc[*it][3] = 0;
+			bc[*it][4] = 16*Umax*c[1]*c[2]*(0.41-c[1])*(0.41-c[2])/pow(0.41,4);
+			bc[*it][5] = 0;
+			bc[*it][6] = 0;
+		}
+		else //no-slip walls
+		{
+			bc[*it][0] = 1;
+			bc[*it][1] = 0;
+			bc[*it][2] = 1;
+			bc[*it][3] = 0;
+			bc[*it][4] = 0;
+			bc[*it][5] = 0;
+			bc[*it][6] = 0;
+			if( fabs(n[2]) < eps && c[0] > 0.45-eps && c[0] < 0.55+eps && c[1] > 0.15-eps && c[1] < 0.25+eps )
+				bcphi[*it][2] = 1;
+		}
+	}
+	
+	std::cout << "Saving output to " << fout << std::endl;
+	
+	m->Save(fout);
+	
+	return 0;
+}