Merge branch 'master' into andre_branch

Conflicts:
	Source/Mesh/modify.cpp

.gitignore

@@ -5,10 +5,14 @@ SyncToy*
 
 *.iml
 
-.idea/*.xml
+.idea/
 
 Source/Solver/solver_fcbiilu2/fcbiilu2.cpp
 
 Source/Solver/solver_k3biilu2/k3d.cpp
 
 Source/Solver/solver_k3biilu2/k3d.h
+
+.vscode/*
+
+.DS_Store

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/ADFVDiscr/main.cpp

+#include "inmost.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
 #include <string.h>
 
-#include "inmost.h"
 using namespace INMOST;
 
 #ifndef M_PI

Examples/ADMFD/main.cpp

+#include "inmost.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
 #include <string.h>
 
-#include "inmost.h"
 using namespace INMOST;
 
 #ifndef M_PI

Examples/AdaptiveMesh/amesh.cpp

 #include "amesh.h"
+#include <iomanip>
+
+using namespace std;
 
 /// todo:
 /// 1. coarsment
@@ -31,6 +34,233 @@ namespace INMOST
 		(void) size;
 		element->Integer(tag) = std::min(element->Integer(tag),*((const INMOST_DATA_INTEGER_TYPE *)data));
 	}
+
+
+    void AdaptiveMesh::PrintSetLocal(string offset, ElementSet it, stringstream& ss)
+    {
+        stringstream ss1;
+        ss1 << offset << rank << ": Set : " << setw(5) << it.GetName() << " ";
+        for (int i = ss1.str().length(); i < 23; i++) ss1 << " ";
+        ss << ss1.str();
+        ss << setw(6);
+        if (it.GetStatus() == Element::Shared) ss << "shared";
+        else if (it.GetStatus() == Element::Ghost) ss << "ghost";
+        else if (it.GetStatus() == Element::Owned) ss << "owned";
+        else ss << "none";
+
+        ss << " tag_owner (" << it.IntegerDF(OwnerTag()) << ")";
+
+        //ss << "   level (" << level[it.self()] << ")  ";
+        ss << " tag_processors (";
+        stringstream ss2;
+        Storage::integer_array arr = it.IntegerArrayDV(tag_processors);
+        for (int i = 0; i < arr.size(); i++)
+            ss2 << arr[i] << " ";
+        ss << setw(5) << ss2.str() <<")";
+        
+
+        ElementSet::iterator p = it.Begin();
+        ss << "     | Refs: ";
+        int first = 0;
+        while(p != it.End())
+        {
+        //    if (first++ == 0) ss << endl << offset << "   ";
+            string type = "unknw";
+            if (p->GetElementType() == CELL) type = "cell";
+            if (p->GetElementType() == FACE) type = "face";
+            if (p->GetElementType() == EDGE) type = "edge";
+            if (p->GetElementType() == NODE) type = "node";
+            ss << type << "-" << p->GlobalID() << " ";
+            p++;
+        }
+        ss << endl;
+
+        for(ElementSet child = it.GetChild(); child.isValid(); child = child.GetSibling())
+        {
+            PrintSetLocal(offset + "   ",child,ss);
+        }
+    }
+
+
+    void AdaptiveMesh::PrintSet()
+    {
+        stringstream ss;
+        for(Mesh::iteratorSet it = BeginSet(); it != EndSet(); ++it) 
+        {
+            if (it->HaveParent()) continue;
+            PrintSetLocal("",ElementSet(this,*it),ss);
+        }
+        cout << ss.str() << endl;
+    }
+
+    void AdaptiveMesh::PrintMesh(ostream& os, int cell, int face, int edge, int node)
+    {
+        if (cell + face + edge + node == 0) return;
+        stringstream ss;
+        ss << "================= " << rank << " =====================" << endl;
+        if (cell)
+        {
+            ss << "Cells: " << NumberOfCells() <<  endl;
+            for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); ++it) 
+            {
+                ss << rank << ": " << it->GlobalID() << " - " ;            
+                if (it->GetStatus() == Element::Shared) ss << "shared";
+                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
+                else ss << "none";
+
+                Storage::reference_array refs = ref_tag[it->self()];
+                if (refs.size() > 0) ss << ". Ref: ";
+                for(Storage::reference_array::size_type i = 0; i < refs.size(); ++i)
+                {
+                    string type = "unknw";
+                    if (refs[i].GetElementType() == CELL) type = "cell";
+                    if (refs[i].GetElementType() == FACE) type = "face";
+                    if (refs[i].GetElementType() == EDGE) type = "edge";
+                    if (refs[i].GetElementType() == NODE) type = "node";
+                    ss << "(" << type << "," << refs[i]->GlobalID() << ") ";
+                }
+
+                ss << endl;
+            }
+        }
+
+        if (face)
+        {
+            ss << "Faces: " << NumberOfFaces() <<  endl;
+            for(Mesh::iteratorFace it = BeginFace(); it != EndFace(); ++it) 
+            {
+                ss << rank << ": " << it->GlobalID() << " - " ;            
+                if (it->GetStatus() == Element::Shared) ss << "shared";
+                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
+                else ss << "none";
+
+                Storage::reference_array refs = ref_tag[it->self()];
+                if (refs.size() > 0) ss << ". Ref: ";
+                for(Storage::reference_array::size_type i = 0; i < refs.size(); ++i)
+                {
+                    string type = "unknw";
+                    if (refs[i].GetElementType() == CELL) type = "cell";
+                    if (refs[i].GetElementType() == FACE) type = "face";
+                    if (refs[i].GetElementType() == EDGE) type = "edge";
+                    if (refs[i].GetElementType() == NODE) type = "node";
+                    ss << "(" << type << "," << refs[i]->GlobalID() << ") ";
+                }
+
+                ss << endl;
+            }
+        }
+
+        if (edge)
+        {
+            ss << "Edges: " << NumberOfEdges() <<  endl;
+            for(Mesh::iteratorEdge it = BeginEdge(); it != EndEdge(); ++it) 
+            {
+                ss << rank << ": " << it->GlobalID() << " - " ;            
+                if (it->GetStatus() == Element::Shared) ss << "shared";
+                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
+                else ss << "none";
+
+                Storage::reference_array refs = ref_tag[it->self()];
+                if (refs.size() > 0) ss << ". Ref: ";
+                for(Storage::reference_array::size_type i = 0; i < refs.size(); ++i)
+                {
+                    string type = "unknw";
+                    if (refs[i].GetElementType() == CELL) type = "cell";
+                    if (refs[i].GetElementType() == FACE) type = "face";
+                    if (refs[i].GetElementType() == EDGE) type = "edge";
+                    if (refs[i].GetElementType() == NODE) type = "node";
+                    ss << "(" << type << "," << refs[i]->GlobalID() << ") ";
+                }
+
+                ss << endl;
+            }
+        }
+
+        if (node)
+        {
+            ss << "Nodes:" << endl;
+            for(Mesh::iteratorNode it = BeginNode(); it != EndNode(); ++it) 
+            {
+                ss << rank << ": " << it->GlobalID() << " - " ;            
+                if (it->GetStatus() == Element::Shared) ss << "shared";
+                else if (it->GetStatus() == Element::Ghost) ss << "ghost";
+                else ss << "none";
+
+                {
+                    Storage::reference_array refs = ref_tag[it->self()];
+                    if (refs.size() > 0) ss << ". Ref: ";
+                    for(Storage::reference_array::size_type i = 0; i < refs.size(); ++i)
+                    {
+                        string type = "unknw";
+                        if (refs[i].GetElementType() == CELL) type = "cell";
+                    if (refs[i].GetElementType() == FACE) type = "face";
+                    if (refs[i].GetElementType() == EDGE) type = "edge";
+                        if (refs[i].GetElementType() == NODE) type = "node";
+                        ss << "(" << type << "," << refs[i]->GlobalID() << ") ";
+                        
+                    }
+                }
+                ss << endl;
+            }
+        }
+
+        ss << "=========================================" << endl;
+        os << ss.str() << endl;
+    }
+
+    void AdaptiveMesh::UpdateStatus()
+    {
+
+        for(ElementType mask = CELL; mask >= NODE; mask = PrevElementType(mask))
+        {
+            for(iteratorElement it = BeginElement(mask); it != EndElement(); it++)
+            {
+                int stat = 0;
+                if (it->GetStatus() == Element::Shared) stat = 1;
+                else if (it->GetStatus() == Element::Ghost)  stat = 2;
+
+                tag_status[it->self()] = stat;
+            }
+        }
+
+
+		for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); ++it) 
+        {
+                }
+    }
+    
+    void AdaptiveMesh::PrintSet(ElementSet set, std::string offset)
+    {
+        std::cout << offset << "Set: " << std::endl;
+
+		ElementSet::iterator it = set.Begin();
+		while(it != set.End())
+        {
+		    ElementSet parent(this,parent_set[*it]);
+            std::cout << offset << it->GlobalID() << " - " << level[*it] << " : ";
+		    
+            ElementSet::iterator p = parent.Begin();
+    		while(p != parent.End())
+            {
+                std::cout << p->GlobalID() << " ";
+                p++;
+            }
+            std::cout << std::endl;
+            it++;
+        }
+
+		for(ElementSet child = set.GetChild(); child.isValid(); child = child.GetSibling())
+        {
+            PrintSet(child,offset + "  ");
+        }
+
+    }
+
+    void AdaptiveMesh::Test()
+    {
+        std::cout << rank << ": ================" << endl;
+        PrintSet(root,"");
+    }
 	
 	
 	void AdaptiveMesh::ClearData()
@@ -57,7 +287,7 @@ namespace INMOST
 					parent_set[it->self()] = root.GetHandle();
 				}
 			}
-		}
+        }
 		if( !HaveGlobalID(CELL) ) AssignGlobalID(CELL); //for unique set names
 	}
 	
@@ -65,10 +295,15 @@ namespace INMOST
 	{
 		//create a tag that stores maximal refinement level of each element
 		level = CreateTag("REFINEMENT_LEVEL",DATA_INTEGER,CELL|FACE|EDGE|NODE|ESET,NONE,1);
+		tag_status = CreateTag("TAG_STATUS",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
+		tag_an = CreateTag("TAG_AN",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
+		ref_tag = CreateTag("REF",DATA_REFERENCE,CELL|FACE|EDGE|NODE,NONE);
 		//create a tag that stores links to all the hanging nodes of the cell
 		hanging_nodes = CreateTag("HANGING_NODES",DATA_REFERENCE,CELL|FACE,NONE);
 		//create a tag that stores links to sets
 		parent_set = CreateTag("PARENT_SET",DATA_REFERENCE,CELL,NONE,1);
+	    size = GetProcessorsNumber();
+    	rank = GetProcessorRank();
 	}
 	
 	AdaptiveMesh::~AdaptiveMesh()
@@ -79,6 +314,7 @@ namespace INMOST
 	
 	bool AdaptiveMesh::Refine(TagInteger & indicator)
 	{
+        cout << ro() << rank << ": IN REFINE" << endl;
 		static int call_counter = 0;
 		int ret = 0; //return number of refined cells
 		//initialize tree structure
@@ -326,31 +562,56 @@ namespace INMOST
 			ReleaseMarker(mark_hanging_nodes);
 			ReleaseMarker(mark_cell_edges);
 			DeleteTag(internal_face_edges);
+			//ExchangeData(hanging_nodes,CELL | FACE,0);
 			//10.jump to later schedule, and go to 7.
 			schedule_counter--;
 		}
+
 		//free created tag
 		DeleteTag(indicator,FACE|EDGE);
+
+
+        stringstream ss;
+        ss << ro() << rank << ": during refine: ";
+        for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); ++it) 
+        {
+            int marked = (it->GetMarker(NewMarker())) ? 1 : 0;
+            int st = 0;
+            if (it->GetStatus() == Element::Shared) st = 1;
+            else if (it->GetStatus() == Element::Ghost) st = 2;
+            //tag_an[it->self()] = it->IntegerDF(OwnerTag());
+            tag_an[it->self()] = marked;
+            ss << "(" << it->GlobalID() << "," << level[it->self()] << "," << marked << ") ";
+
+        }
+        for(Mesh::iteratorNode it = BeginNode(); it != EndNode(); ++it) 
+        {
+            int marked = (it->GetMarker(NewMarker())) ? 1 : 0;
+            tag_an[it->self()] = marked;
+        }
+        //cout << ss.str() << endl;
+
+
 		//11. Restore parallel connectivity, global ids
-		//ResolveModification();
+        ResolveShared(true);
+        //if (call_counter == 0)
+		ResolveModification();
 		//12. Let the user update their data
 		//todo: call back function for New() cells
 		//13. Delete old elements of the mesh
 		ApplyModification();
 		//14. Done
+        //cout << rank << ": Before end " << endl;
 		EndModification();
+        //cout << rank << ": After end " << endl;
 		//reorder element's data to free up space
 		ReorderEmpty(CELL|FACE|EDGE|NODE);
 		//return number of refined cells
 		call_counter++;
+        cout << ro() << rank << ": END REFINE " << (ret != 0) << endl;
 		return ret != 0;
 	}
 	
-	
-	
-	
-	
-	
 	bool AdaptiveMesh::Coarse(TagInteger & indicator)
 	{
 		static int call_counter = 0;
@@ -635,12 +896,14 @@ namespace INMOST
 				}
 			}
 			//jump to later schedule
+			//ExchangeData(hanging_nodes,CELL | FACE,0);
 			schedule_counter--;
 		}
 		//free created tag
 		DeleteTag(indicator,FACE|EDGE);
 		//todo:
-		//ResolveModification();
+    	ResolveShared(true);
+		ResolveModification();
 		//todo:
 		//let the user update their data
 		ApplyModification();
@@ -666,4 +929,4 @@ namespace INMOST
 		call_counter++;
 		return ret != 0;
 	}
-}
\ No newline at end of file
+}

Examples/AdaptiveMesh/amesh.h

@@ -7,13 +7,19 @@ 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
+		TagInteger tag_status; 
+		TagInteger tag_an; 
 		TagReference parent_set; //<Link to the set that contains an element.
 		TagReferenceArray hanging_nodes; //< Link to current hanging nodes of the cell.
+        int rank;
+        int size;
 		/// Prepare sets for coarsements.
 		/// Do not do this in constructor, since mesh may contain no cells.
 		void PrepareSet();
+        void PrintSetLocal(std::string offset, ElementSet it, std::stringstream& ss);
 	public:
+		TagInteger level; //< Refinement level of the cell
+		TagReferenceArray ref_tag; //<Link to the set that contains an element.
 		Storage::integer GetLevel(const Storage & e) {return level[e];}
 		AdaptiveMesh();
 		~AdaptiveMesh();
@@ -23,6 +29,12 @@ namespace INMOST
 		bool Coarse(TagInteger & indicator);
 		/// Delete all data related to mesh refinement-coarsement.
 		void ClearData();
+        void PrintSet(ElementSet set, std::string offset);
+        void Test();
+        void PrintMesh(std::ostream& os, int cell = 0, int face = 0, int edge = 0, int node = 0);
+        void PrintSet();
+        void UpdateStatus();
+        void test_sets();
 	};
 }
 

Examples/AdaptiveMesh/main.cpp

@@ -9,6 +9,7 @@ int main(int argc, char ** argv)
 	if( argc > 1 )
 	{
 		AdaptiveMesh m;
+		m.SetCommunicator(INMOST_MPI_COMM_WORLD);
 		m.Load(argv[1]);
 		//m.SetTopologyCheck(NEED_TEST_CLOSURE);
 		//m.SetTopologyCheck(PROHIBIT_MULTILINE);
@@ -29,6 +30,8 @@ int main(int argc, char ** argv)
 				if( it->Coords()[d] < cmin[d] ) cmin[d] = it->Coords()[d];
 			}
 		}
+		m.AggregateMax(cmax,3);
+		m.AggregateMin(cmin,3);
 		r = 1;
 		for(int d = 0; d < 3; ++d)
 		{
@@ -47,18 +50,20 @@ int main(int argc, char ** argv)
 				numref = 0;
 				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( m.GetLevel(it->self()) < max_levels )
 				{
-					it->Centroid(xyz);
+					it->Barycenter(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)
+					//if( q < 0.02 )
 					{
 						indicator[it->self()] = 1;
 						numref++;
 					}
 				}
+				numref = m.Integrate(numref);
 				if( numref )
 				{
 					std::cout << "k " << k << " refcnt " << refcnt << " " <<  r*k << " < r < " << r*(k+1) << " numref " << numref << " cells " << m.NumberOfCells() << std::endl;
@@ -89,7 +94,7 @@ int main(int argc, char ** argv)
 				numref = 0;
 				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( m.GetLevel(it->self()) > 0 )
 				{
-					it->Centroid(xyz);
+					it->Barycenter(xyz);
 					//refine a circle
 					q = 0;
 					for(int d = 0; d < 3; ++d)
@@ -101,6 +106,7 @@ int main(int argc, char ** argv)
 						numref++;
 					}
 				}
+				numref = m.Integrate(numref);
 				if( numref )
 				{
 					std::cout << "k " << k << " crscnt " << refcnt << " " << r*k << " < r < " << r*(k+1) << " numcrs " << numref << " cells " << m.NumberOfCells() <<  std::endl;
@@ -131,11 +137,19 @@ int main(int argc, char ** argv)
 			
 			{
 				std::stringstream file;
-				file << "step_" << k << ".vtk";
+				file << "step_" << k << ".pvtk";
 				m.Save(file.str());
 			}
 			
 			{
+				TagInteger tag_owner = m.CreateTag("OWN",DATA_INTEGER,CELL,NONE,1);
+				TagInteger tag_owner0 = m.GetTag("OWNER_PROCESSOR");
+				TagInteger tag_stat = m.CreateTag("STAT",DATA_INTEGER,CELL,NONE,1);
+				for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
+				{
+					tag_owner[*it] = tag_owner0[*it];
+					tag_stat[*it] = it->GetStatus();
+				}
 				std::stringstream file;
 				file << "step_" << k << ".pmf";
 				m.Save(file.str());
@@ -143,4 +157,6 @@ int main(int argc, char ** argv)
 		}
 	}
 	else std::cout << "Usage: " << argv[0] << " mesh_file [max_levels=2]" << std::endl;
-}
\ No newline at end of file
+	
+	Mesh::Finalize();
+}

Examples/DrawMatrix/main.cpp

@@ -17,8 +17,9 @@ Sparse::Matrix * m = NULL;
 
 int zoom = 200;
 int block_size = 0;
+int draw_color = 0;
 
-//~ Storage::real min = 1e20,max = -1e20;
+double amin = 1e20,amax = -1e20;
 
 void printtext(const char * fmt, ...)
 {
@@ -379,6 +380,18 @@ void keyboard(unsigned char key, int x, int y)
 
 		}
 	}
+	if( key == 'd' )
+	{
+		draw_color = (draw_color+1)%5;
+		switch(draw_color)
+		{
+			case 0: std::cout << "color diagonal" << std::endl; break;
+			case 1: std::cout << "color by min:max" << std::endl; break;
+			case 2: std::cout << "color by min:max in log scale" << std::endl; break;
+			case 3: std::cout << "color by min:max per row" << std::endl; break;
+			case 4: std::cout << "color by min:max per row in log scale" << std::endl; break;
+		}
+	}
 	if (key == 'c') ord->clear();
 	if (key == 'r')
 	{
@@ -415,53 +428,128 @@ void draw()
 	glVertex2i(m->Size(), m->Size() - ord->GetMatrixPartSize());
 	glEnd();
 	*/
-
-	glColor3f(0,0,0);
-	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() )
-			//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();
-
+	
 
 	/*
-	glColor3f(0.0, 1.0, 0);
-	glBegin(GL_QUADS);
-	for (Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it)
+	 glColor3f(0.0, 1.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]) > row_sum[it-m->Begin()]) //DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(ind));
+	 DrawEntry((it - m->Begin()), m->Size() - ind);