Many fixes and updates for model composition

.gitignore

@@ -12,3 +12,5 @@ Source/Solver/solver_fcbiilu2/fcbiilu2.cpp
 Source/Solver/solver_k3biilu2/k3d.cpp
 
 Source/Solver/solver_k3biilu2/k3d.h
+
+.vscode/*

CMakeLists.txt

@@ -265,6 +265,9 @@ endif(COMPILE_TESTS)
 
 set(INMOST_INSTALL_HEADERS Source/Headers/inmost.h
                            Source/Headers/inmost_autodiff.h
+                           Source/Headers/inmost_residual.h
+                           Source/Headers/inmost_model.h
+                           Source/Headers/inmost_operator.h
                            Source/Headers/inmost_common.h
                            Source/Headers/inmost_data.h
                            Source/Headers/inmost_dense.h
@@ -276,6 +279,7 @@ set(INMOST_INSTALL_HEADERS Source/Headers/inmost.h
                            Source/Headers/inmost_sparse.h
                            Source/Headers/inmost_xml.h
                            Source/Headers/inmost_variable.h
+                           Source/Headers/inmost_block_variable.h
                            Source/Headers/container.hpp)
 
 
@@ -295,6 +299,9 @@ set_property(TARGET inmost PROPERTY PUBLIC_HEADER
 	"${PROJECT_BINARY_DIR}/inmost_options.h"
 	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost.h"
 	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_autodiff.h"
+	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_residual.h"
+	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_operator.h"
+	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_model.h"
 	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_common.h"
 	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_data.h"
 	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_dense.h"
@@ -305,6 +312,7 @@ set_property(TARGET inmost PROPERTY PUBLIC_HEADER
 	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_solver.h"
 	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_sparse.h"
 	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_variable.h"
+	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_block_variable.h"
 	"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_xml.h"
 	"${PROJECT_SOURCE_DIR}/Source/Headers/container.hpp")
 

Examples/DrawMatrix/main.cpp

@@ -420,7 +420,7 @@ void draw()
 	glBegin(GL_QUADS);
 	for(Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
 		for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
-			if( jt->first != it - m->Begin() )
+			if( jt->first != it - m->Begin() && jt->second)
 			//DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(jt->first));//, sqrt((jt->second-min)/(max-min)));
 			DrawEntry(jt->first, m->Size() - (it - m->Begin()));//, sqrt((jt->second-min)/(max-min)));
 	glEnd();
@@ -451,17 +451,17 @@ void draw()
 	}
 	glEnd();
 
-	zoom += 1;
+	//zoom += 1;
 	glColor3f(1.0, 0, 0);
 	glBegin(GL_QUADS);
 	for (Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
 	{
 		int ind = it - m->Begin();
-		if (fabs((*it)[ind]) < 1e-9) //DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(ind));
+		if (fabs((*it)[ind]) < 1e-13) //DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(ind));
 			DrawEntry((it - m->Begin()), m->Size() - ind);
 	}
 	glEnd();
-	zoom -= 1;
+	//zoom -= 1;
 
 	if (CommonInput != NULL)
 	{

Examples/GridTools/CMakeLists.txt

@@ -9,6 +9,10 @@ add_executable(SliceFunc slice_func.cpp)
 add_executable(SplitNonplanar split_nonplanar.cpp)
 add_executable(CollapseDegenerate collapse_degenerate.cpp)
 add_executable(Bnd2Stl bnd2stl.cpp)
+add_executable(Sector sector.cpp)
+add_executable(Sew sew.cpp)
+add_executable(Scale scale.cpp)
+add_executable(Move move.cpp)
 
 target_link_libraries(FixFaults inmost)
 if(USE_MPI)
@@ -106,4 +110,44 @@ if(USE_MPI)
 endif(USE_MPI)
 install(TARGETS Bnd2Stl EXPORT inmost-targets RUNTIME DESTINATION bin)
 
+target_link_libraries(Sector inmost)
+if(USE_MPI)
+  message("linking Sector with MPI")
+  target_link_libraries(Sector ${MPI_LIBRARIES})
+  if(MPI_LINK_FLAGS)
+    set_target_properties(Sector PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif()
+endif(USE_MPI)
+install(TARGETS Sector EXPORT inmost-targets RUNTIME DESTINATION bin)
+
+target_link_libraries(Sew inmost)
+if(USE_MPI)
+  message("linking Sew with MPI")
+  target_link_libraries(Sew ${MPI_LIBRARIES})
+  if(MPI_LINK_FLAGS)
+    set_target_properties(Sew PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif()
+endif(USE_MPI)
+install(TARGETS Sew EXPORT inmost-targets RUNTIME DESTINATION bin)
+
+target_link_libraries(Scale inmost)
+if(USE_MPI)
+  message("linking Scale with MPI")
+  target_link_libraries(Scale ${MPI_LIBRARIES})
+  if(MPI_LINK_FLAGS)
+    set_target_properties(Scale PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif()
+endif(USE_MPI)
+install(TARGETS Scale EXPORT inmost-targets RUNTIME DESTINATION bin)
+
+
+target_link_libraries(Move inmost)
+if(USE_MPI)
+  message("linking Move with MPI")
+  target_link_libraries(Move ${MPI_LIBRARIES})
+  if(MPI_LINK_FLAGS)
+    set_target_properties(Move PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif()
+endif(USE_MPI)
+install(TARGETS Move EXPORT inmost-targets RUNTIME DESTINATION bin)
 

Examples/GridTools/move.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 [move_x=1] [move_y=1] [move_z=1] [output_mesh]\n",argv[0]);
+		return -1;
+	}
+
+	Mesh m;
+	m.Load(argv[1]);
+	double mx = argc > 2 ? atof(argv[2]) : 1;
+	double my = argc > 3 ? atof(argv[3]) : 1;
+	double mz = argc > 4 ? atof(argv[4]) : 1;
+	std::string save_file = argc > 5 ? std::string(argv[5]) : "out.vtk";
+	
+	double xmax = -1.0e20, xmin = 1.0e20;
+	double ymax = -1.0e20, ymin = 1.0e20;
+	double zmax = -1.0e20, zmin = 1.0e20;
+	
+	for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
+	{
+		double x = n->Coords()[0], y = n->Coords()[1], z = n->Coords()[2];
+		if( x > xmax ) xmax = x;
+		if( x < xmin ) xmin = x;
+		if( y > ymax ) ymax = y;
+		if( y < ymin ) ymin = y;
+		if( z > zmax ) zmax = z;
+		if( z < zmin ) zmin = z;
+	}
+	
+	std::cout << "x: " << xmin << ":" << xmax << std::endl;
+	std::cout << "y: " << ymin << ":" << ymax << std::endl;
+	std::cout << "z: " << zmin << ":" << zmax << std::endl;
+	
+	for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
+	{
+		n->Coords()[0] += mx;
+		n->Coords()[1] += my;
+		n->Coords()[2] += mz;
+	}
+	
+	xmax = -1.0e20, xmin = 1.0e20;
+	ymax = -1.0e20, ymin = 1.0e20;
+	zmax = -1.0e20, zmin = 1.0e20;
+	
+	for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
+	{
+		double x = n->Coords()[0], y = n->Coords()[1], z = n->Coords()[2];
+		if( x > xmax ) xmax = x;
+		if( x < xmin ) xmin = x;
+		if( y > ymax ) ymax = y;
+		if( y < ymin ) ymin = y;
+		if( z > zmax ) zmax = z;
+		if( z < zmin ) zmin = z;
+	}
+	
+	std::cout << "x: " << xmin << ":" << xmax << std::endl;
+	std::cout << "y: " << ymin << ":" << ymax << std::endl;
+	std::cout << "z: " << zmin << ":" << zmax << std::endl;
+	
+	std::cout << "Save to " << save_file << std::endl;
+	m.Save(save_file);
+	
+	return 0;
+}

Examples/GridTools/scale.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 [scale_x=1] [scale_y=1] [scale_z=1] [output_mesh]\n",argv[0]);
+		return -1;
+	}
+
+	Mesh m;
+	m.Load(argv[1]);
+	double sx = argc > 2 ? atof(argv[2]) : 1;
+	double sy = argc > 3 ? atof(argv[3]) : 1;
+	double sz = argc > 4 ? atof(argv[4]) : 1;
+	std::string save_file = argc > 5 ? std::string(argv[5]) : "out.vtk";
+	
+	double xmax = -1.0e20, xmin = 1.0e20;
+	double ymax = -1.0e20, ymin = 1.0e20;
+	double zmax = -1.0e20, zmin = 1.0e20;
+	
+	for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
+	{
+		double x = n->Coords()[0], y = n->Coords()[1], z = n->Coords()[2];
+		if( x > xmax ) xmax = x;
+		if( x < xmin ) xmin = x;
+		if( y > ymax ) ymax = y;
+		if( y < ymin ) ymin = y;
+		if( z > zmax ) zmax = z;
+		if( z < zmin ) zmin = z;
+	}
+	
+	std::cout << "x: " << xmin << ":" << xmax << std::endl;
+	std::cout << "y: " << ymin << ":" << ymax << std::endl;
+	std::cout << "z: " << zmin << ":" << zmax << std::endl;
+	
+	for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
+	{
+		n->Coords()[0] *= sx;
+		n->Coords()[1] *= sy;
+		n->Coords()[2] *= sz;
+	}
+	
+	xmax = -1.0e20, xmin = 1.0e20;
+	ymax = -1.0e20, ymin = 1.0e20;
+	zmax = -1.0e20, zmin = 1.0e20;
+	
+	for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
+	{
+		double x = n->Coords()[0], y = n->Coords()[1], z = n->Coords()[2];
+		if( x > xmax ) xmax = x;
+		if( x < xmin ) xmin = x;
+		if( y > ymax ) ymax = y;
+		if( y < ymin ) ymin = y;
+		if( z > zmax ) zmax = z;
+		if( z < zmin ) zmin = z;
+	}
+	
+	std::cout << "x: " << xmin << ":" << xmax << std::endl;
+	std::cout << "y: " << ymin << ":" << ymax << std::endl;
+	std::cout << "z: " << zmin << ":" << zmax << std::endl;
+	
+	std::cout << "Save to " << save_file << std::endl;
+	m.Save(save_file);
+	
+	return 0;
+}

Examples/GridTools/sector.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 [rotation_angle=0 degrees] [output_mesh]\n",argv[0]);
+		return -1;
+	}
+
+	Mesh m;
+	m.Load(argv[1]);
+	
+	const double pi = 3.1415926536;
+	double theta = 0, ct, st;
+	if( argc > 2 ) theta = atof(argv[2])/180.0*pi;
+	ct = cos(theta);
+	st = sin(theta);
+	
+	double xmin = 1.0e20, xmax = -1.0e20;
+	double ymin = 1.0e20, ymax = -1.0e20;
+	
+	for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
+	{
+		double x = n->Coords()[0];
+		double y = n->Coords()[1];
+		if( x > xmax ) xmax = x;
+		if( x < xmin ) xmin = x;
+		if( y > ymax ) ymax = y;
+		if( y < ymin ) ymin = y;
+	}
+	
+	std::cout << "x: " << xmin << ":" << xmax << std::endl;
+	std::cout << "y: " << ymin << ":" << ymax << std::endl;
+	
+	for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
+	{
+		double x = n->Coords()[0], mx = x;
+		double y = n->Coords()[1], my = y;
+		double alpha = pi/4*(2*y-1);
+		mx += 1 - x*(1-cos(alpha)) + x*cos(2*alpha)*0.2;
+		my += x*sin(alpha);
+		n->Coords()[0] = (mx-0.5)*ct - (my-0.5)*st + 0.5;
+		n->Coords()[1] = (mx-0.5)*st + (my-0.5)*ct + 0.5;
+	}
+	
+	xmin = 1.0e20, xmax = -1.0e20;
+	ymin = 1.0e20, ymax = -1.0e20;
+	
+	for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
+	{
+		double x = n->Coords()[0];
+		double y = n->Coords()[1];
+		if( x > xmax ) xmax = x;
+		if( x < xmin ) xmin = x;
+		if( y > ymax ) ymax = y;
+		if( y < ymin ) ymin = y;
+	}
+	
+	std::cout << "x: " << xmin << ":" << xmax << std::endl;
+	std::cout << "y: " << ymin << ":" << ymax << std::endl;
+	
+	if( argc > 3 )
+	{
+		std::cout << "Save to " << argv[3] << std::endl;
+		m.Save(argv[3]);
+	}
+	else
+	{
+		std::cout << "Save to out.vtk" << std::endl;
+		m.Save("out.vtk");
+	}
+
+	return 0;
+}

Examples/GridTools/sew.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_mesh1 [input_mesh2] ... [output_mesh]\n",argv[0]);
+		return -1;
+	}
+
+	Mesh m;
+	for(int k = 1 ; k < argc; ++k)
+		m.Load(argv[k]);
+	
+	double xmax = -1.0e20, xmin = 1.0e20;
+	double ymax = -1.0e20, ymin = 1.0e20;
+	double zmax = -1.0e20, zmin = 1.0e20;
+	
+	for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
+	{
+		double x = n->Coords()[0], y = n->Coords()[1], z = n->Coords()[2];
+		if( x > xmax ) xmax = x;
+		if( x < xmin ) xmin = x;
+		if( y > ymax ) ymax = y;
+		if( y < ymin ) ymin = y;
+		if( z > zmax ) zmax = z;
+		if( z < zmin ) zmin = z;
+	}
+	
+	std::cout << "x: " << xmin << ":" << xmax << std::endl;
+	std::cout << "y: " << ymin << ":" << ymax << std::endl;
+	std::cout << "z: " << zmin << ":" << zmax << std::endl;
+	
+	std::cout << "Save to out.pmf" << std::endl;
+	m.Save("out.pmf");
+	
+	return 0;
+}

Examples/GridTools/slice_func.cpp

@@ -34,6 +34,10 @@ double func(double x, double y, double z, int n)
 		if (y < Ly) return y - Ly;
 		return fmin( fmin(x-Lx,Rx-x), fmin(y-Ly, Ry-y) );
 	}
+	else if( n == 6 )
+	{
+		return std::min(std::min(sqrt((x-0.48)*(x-0.48) + (z-1)*(z-1))-0.25,sqrt((y-0.53)*(y-0.53) + (z-3)*(z-3))-0.24),sqrt((x-0.6)*(x-0.6) + (z-5)*(z-5))-0.3);
+	}
 	return 1;
 }
 
@@ -1057,7 +1061,7 @@ int main(int argc, char ** argv)
 			nvv++;
 		}
 		vv /= nvv;
-		if( v < vv*0.15 )
+		if( v < vv*0.01 )
 		{
 			if( false )
 			{
@@ -1117,7 +1121,7 @@ int main(int argc, char ** argv)
 	}
 	
 	for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
-		if( material[*it] == 0 )
+		if( material[*it] == 0 || it->Integer(collapse) )
 			it->Delete();
 	
 	for(ElementType etype = FACE; etype >= NODE; etype = PrevElementType(etype) )
@@ -1130,4 +1134,4 @@ int main(int argc, char ** argv)
 	m.ReleaseMarker(mrk,EDGE|FACE|NODE);
 	m.Save(grid_out);
 	return 0;
-}
+}
\ No newline at end of file

Examples/OldDrawGrid/CMakeLists.txt

@@ -29,7 +29,9 @@ set(SOURCE main.cpp
 		   picker.h
 		   picker.cpp
 		   clipper.h
-		   clipper.cpp)
+		   clipper.cpp
+		   vector.h
+		   vector.cpp)
 find_package(OpenGL)
 find_package(GLUT)
 

Examples/OldDrawGrid/clipboard.cpp

@@ -73,14 +73,14 @@ bool setTextToPasteboard(std::string str)
     PasteboardCreate( kPasteboardClipboard, &pasteboard );
     
     err = PasteboardClear( pasteboard );
-    require_noerr( err, PasteboardClear_FAILED );
+    //require_noerr( err, PasteboardClear_FAILED );
     
     CFDataRef data;
     
     data =   CFDataCreate(kCFAllocatorDefault, (UInt8*)byteArrayIndex, strlen(byteArrayIndex)+1);
     
     err = PasteboardPutItemFlavor( pasteboard, (PasteboardItemID)1, kUTTypeUTF8PlainText, data, 0);
-    require_noerr( err, PasteboardPutItemFlavor_FAILED );
+    //require_noerr( err, PasteboardPutItemFlavor_FAILED );
     
 PasteboardPutItemFlavor_FAILED:
 PasteboardClear_FAILED:
@@ -101,14 +101,14 @@ std::string getTextFromPasteboard()
     err = badPasteboardSyncErr;
     
     err = PasteboardGetItemCount( inPasteboard, &itemCount );
-    require_noerr( err, CantGetPasteboardItemCount );
+    //require_noerr( err, CantGetPasteboardItemCount );
     
     for( int itemIndex = 1; itemIndex <= itemCount; itemIndex++ ) 
 	{
         PasteboardItemID itemID;
         CFDataRef flavorData;
         err = PasteboardGetItemIdentifier( inPasteboard, itemIndex, &itemID );
-        require_noerr( err, CantGetPasteboardItemIdentifier );
+        //require_noerr( err, CantGetPasteboardItemIdentifier );
         err = PasteboardCopyItemFlavorData( inPasteboard, itemID, CFSTR("public.utf8-plain-text"), &flavorData );
         if(err==noErr)data = (char*)CFDataGetBytePtr(flavorData);
         if( data!=NULL && err==noErr )

Examples/OldDrawGrid/main.cpp

@@ -23,6 +23,7 @@
 #include "tga.h"
 #include "screenshot.h"
 #include "volumetric.h"
+#include "vector.h"
 #include "input.h"
 #include "picker.h"
 #include "clipper.h"
@@ -206,6 +207,7 @@ std::vector<face2gl> clip_boundary;
 
 
 volumetric * CommonVolumetricView;
+Vectors * CommonVectors = NULL;
 
 
 
@@ -1294,6 +1296,9 @@ void draw_screen()
 			if( bndupdate ) CommonVolumetricView->camera(campos,interactive);
 			CommonVolumetricView->draw(interactive);
 		}
+		
+		if( CommonVectors )
+			CommonVectors->Draw(interactive);
 
 
 		for(int k = 0; k < streamlines.size(); ++k)
@@ -1504,6 +1509,11 @@ void draw_screen()
 							correct_input = true;
 							streamlines.clear();
 							glutPostRedisplay();
+							if( CommonVectors )
+							{
+								delete CommonVectors;
+								CommonVectors = NULL;
+							}
 
 						}
 					}
@@ -1515,79 +1525,90 @@ void draw_screen()
 						visualization_type = NONE;
 						for (size_t q = 0; q < stype.size(); ++q)
 							stype[q] = tolower(stype[q]);
-						if (stype == "node") visualization_type = NODE;
-						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 (stype == "scale")
+						{
+							double scale = atof(visualization_prompt + k + 1);
+							std::cout << "input scale: " << scale << std::endl;
+							if( CommonVectors )
+								CommonVectors->SetScale(scale);
+							correct_input = true;
+							glutPostRedisplay();
+						}
 						else
 						{
-							disp_e = InvalidElement();
-							bool is_number = true;
-							for (l = k + 1; l < slen; ++l)
-							if (!isdigit(visualization_prompt[l]))
-								is_number = false;
-							if (is_number)
+							if (stype == "node") visualization_type = NODE;
+							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);
+							else
 							{
-								comp = atoi(visualization_prompt + k + 1);
-								visualization_prompt[k] = ':';
+								disp_e = InvalidElement();
+								bool is_number = true;
+								for (l = k + 1; l < slen; ++l)
+								if (!isdigit(visualization_prompt[l]))
+									is_number = false;
+								if (is_number)
+								{
+									comp = atoi(visualization_prompt + k + 1);
+									visualization_prompt[k] = ':';
 
 
-								correct_input = true;
+									correct_input = true;
 
-								if (mesh->isValidElement(visualization_type, comp))
-								{
-									printf("Display data for %s:%d\n", typen, comp);
-									disp_e = mesh->ElementByLocalID(visualization_type, comp);
+									if (mesh->isValidElement(visualization_type, comp))
+									{
+										printf("Display data for %s:%d\n", typen, comp);
+										disp_e =