Merge remote-tracking branch 'INMOST-DEV/master'

Examples/OldDrawGrid/main.cpp

@@ -797,7 +797,7 @@ public:
 		colors.push_back(color_t(0.45,0,0.55));
 		colors.push_back(color_t(0,0,0));
 
-		samples = 2048;
+		samples = 4096;
 
 		float * pixel_array = new float[(samples+2)*4];
 
@@ -808,7 +808,7 @@ public:
 			float t = 1.0f*q/static_cast<float>(samples+1);
 			color_t c = pick_color(t);
 			//countour lines
-			//if( ((q+1) % 64 == 0 || (q+1) % 64 == 63) && (q+1) < samples ) 
+			//if( ((q+1) % 128 == 0 || (q+1) % 128 == 127) && (q+1) < samples )
 			//	c = pick_color(1-t) + color_t(0,2*t*(1-t),0);
 
 			pixel_array[(q)*4+0] = c.r();
@@ -4093,11 +4093,13 @@ double display_elem_info(Element e, double top, double left, double interval)
 			char str[4096];
 			char temp[4096];
 			str[0] = '\0';
+			int dsize;
 			switch(t->GetDataType())
 			{
 				case DATA_INTEGER:
 				{
 					Storage::integer_array arr = e->IntegerArray(*t);
+					dsize = arr.size();
 					for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
 					{
 						sprintf(temp,"%s %d",str,arr[k]);
@@ -4108,6 +4110,7 @@ double display_elem_info(Element e, double top, double left, double interval)
 				case DATA_REAL:
 				{
 					Storage::real_array arr = e->RealArray(*t);
+					dsize = arr.size();
 					for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
 					{
 						sprintf(temp,"%s %lf",str,arr[k]);
@@ -4118,6 +4121,7 @@ double display_elem_info(Element e, double top, double left, double interval)
 				case DATA_BULK:
 				{
 					Storage::bulk_array arr = e->BulkArray(*t);
+					dsize = arr.size();
 					for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
 					{
 						sprintf(temp,"%s %d",str,arr[k]);
@@ -4128,6 +4132,7 @@ double display_elem_info(Element e, double top, double left, double interval)
 				case DATA_REFERENCE:
 				{
 					Storage::reference_array arr = e->ReferenceArray(*t);
+					dsize = arr.size();
 					for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
 					{
 						if(arr.at(k) == InvalidHandle()) sprintf(temp,"%s NULL",str);
@@ -4139,6 +4144,7 @@ double display_elem_info(Element e, double top, double left, double interval)
 				case DATA_REMOTE_REFERENCE:
 				{
 					Storage::remote_reference_array arr = e->RemoteReferenceArray(*t);
+					dsize = arr.size();
 					for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
 					{
 						if(arr.at(k).first == NULL || arr.at(k).second == InvalidHandle()) sprintf(temp,"%s NULL",str);
@@ -4151,6 +4157,7 @@ double display_elem_info(Element e, double top, double left, double interval)
 				case DATA_VARIABLE:
 				{
 					Storage::var_array arr = e->VariableArray(*t);
+					dsize = arr.size();
 					for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
 					{
 						std::stringstream stream;
@@ -4167,7 +4174,7 @@ double display_elem_info(Element e, double top, double left, double interval)
 				}
 #endif
 			}
-			sprintf(temp,"%s %s %s",t->GetTagName().c_str(),DataTypeName(t->GetDataType()),str);
+			sprintf(temp,"%s[%d] %s %s",t->GetTagName().c_str(),dsize,DataTypeName(t->GetDataType()),str);
 			strcpy(str,temp);
 			top -= interval;
 			glRasterPos2f(left,top);
@@ -5036,7 +5043,7 @@ int main(int argc, char ** argv)
 
 
   //if( false )
-  //if( mesh->HaveTag("VELOCITY") && mesh->GetTag("VELOCITY").isDefined(CELL) )
+  if( mesh->HaveTag("VELOCITY") && mesh->GetTag("VELOCITY").isDefined(CELL) )
   {
     printf("preparing octree around mesh, was sets %d\n",mesh->NumberOfSets());
     Octree octsearch = Octree(mesh->CreateSet("octsearch").first);

Source/Headers/inmost_data.h

@@ -409,6 +409,7 @@ namespace INMOST
 			Element operator[] (size_type n) const;
 			class iterator : public shell<remote_reference>::iterator 
 			{
+			public:
 				iterator() :shell<remote_reference>::iterator() {}
 				iterator(const shell<remote_reference>::iterator & other) : shell<remote_reference>::iterator(other) {}
 				iterator(const iterator & other) : shell<remote_reference>::iterator(other) {}
@@ -416,6 +417,7 @@ namespace INMOST
 			};
 			class const_iterator : public shell<remote_reference>::const_iterator 
 			{
+			public:
 				const_iterator() :shell<remote_reference>::const_iterator() {}
 				const_iterator(const shell<remote_reference>::const_iterator & other) : shell<remote_reference>::const_iterator(other) {}
 				const_iterator(const const_iterator & other) : shell<remote_reference>::const_iterator(other) {}
@@ -423,6 +425,7 @@ namespace INMOST
 			};
 			class reverse_iterator : public shell<remote_reference>::reverse_iterator 
 			{
+			public:
 				reverse_iterator() :shell<remote_reference>::reverse_iterator() {}
 				reverse_iterator(const shell<remote_reference>::reverse_iterator & other) : shell<remote_reference>::reverse_iterator(other) {}
 				reverse_iterator(const reverse_iterator & other) : shell<remote_reference>::reverse_iterator(other) {}
@@ -430,6 +433,7 @@ namespace INMOST
 			};
 			class const_reverse_iterator : public shell<remote_reference>::const_reverse_iterator 
 			{
+			public:
 				const_reverse_iterator() :shell<remote_reference>::const_reverse_iterator() {}
 				const_reverse_iterator(const shell<remote_reference>::const_reverse_iterator & other) : shell<remote_reference>::const_reverse_iterator(other) {}
 				const_reverse_iterator(const const_reverse_iterator & other) : shell<remote_reference>::const_reverse_iterator(other) {}

Source/Headers/inmost_dense.h

@@ -30,6 +30,11 @@ namespace INMOST
 	template<> struct Promote<INMOST_DATA_REAL_TYPE, variable>  {typedef variable type;};
 	template<> struct Promote<variable, INMOST_DATA_REAL_TYPE>  {typedef variable type;};
 	template<> struct Promote<variable, variable> {typedef variable type;};
+	template<> struct Promote<INMOST_DATA_REAL_TYPE, hessian_variable>  {typedef hessian_variable type;};
+	template<> struct Promote<hessian_variable, INMOST_DATA_REAL_TYPE>  {typedef hessian_variable type;};
+	template<> struct Promote<variable, hessian_variable>  {typedef hessian_variable type;};
+	template<> struct Promote<hessian_variable, variable>  {typedef hessian_variable type;};
+	template<> struct Promote<hessian_variable, hessian_variable> {typedef hessian_variable type;};
 #else
 	__INLINE INMOST_DATA_REAL_TYPE get_value(INMOST_DATA_REAL_TYPE x) {return x;}
 #endif
@@ -868,16 +873,30 @@ namespace INMOST
 						break;
 					case 21: //symmetric elasticity tensor (note - diagonal first, then off-diagonal rows)
 					{
-						int shift[5] = {6,11,15,18,20};
-						for(int i = 0; i < 6; ++i)
-						{
-							Kc(i,i) = K[i]; //diagonal term first
-							for(int j = i+1; j < 6; ++j)
-								Kc(i,j) = Kc(j,i) = K[shift[i]+j-i-1]; //off diagonal parts
-						}
+						Kc(0,0) = K[0]; //c11
+						Kc(0,1) = Kc(1,0) = K[1]; //c12
+						Kc(0,2) = Kc(2,0) = K[2]; //c13
+						Kc(0,3) = Kc(3,0) = K[3]; //c14
+						Kc(0,4) = Kc(4,0) = K[4]; //c15
+						Kc(0,5) = Kc(5,0) = K[5]; //c16
+						Kc(1,1) = K[6]; //c22
+						Kc(1,2) = Kc(2,1) = K[7]; //c23
+						Kc(1,3) = Kc(3,1) = K[8]; //c24
+						Kc(1,4) = Kc(4,1) = K[9]; //c25
+						Kc(1,5) = Kc(5,1) = K[10]; //c26
+						Kc(2,2) = K[11]; //c33
+						Kc(2,3) = Kc(3,2) = K[12]; //c34
+						Kc(2,4) = Kc(4,2) = K[13]; //c35
+						Kc(2,5) = Kc(5,2) = K[14]; //c36
+						Kc(3,3) = K[15]; //c44
+						Kc(3,4) = Kc(4,3) = K[16]; //c45
+						Kc(3,5) = Kc(5,3) = K[17]; //c46
+						Kc(4,4) = K[18]; //c55
+						Kc(4,5) = Kc(5,4) = K[19]; //c56
+						Kc(5,5) = K[20]; //c66
 						break;
 					}
-					case 36: //full permeability tensor
+					case 36: //full elasticity tensor
 						for(int i = 0; i < 6; ++i)
 							for(int j = 0; j < 6; ++j)
 								Kc(i,j) = K[6*i+j];
@@ -1082,7 +1101,7 @@ namespace INMOST
 			for(enumerator i = ibeg; i < iend; ++i)
 			{
 				for(enumerator j = jbeg; j < jend; ++j)
-					ret(i-ibeg,j-ibeg) = (*this)(i,j);
+					ret(i-ibeg,j-jbeg) = (*this)(i,j);
 			}
 			return ret;
 		}
@@ -1091,6 +1110,7 @@ namespace INMOST
 	typedef Matrix<INMOST_DATA_REAL_TYPE> rMatrix; //shortcut for real matrix
 #if defined(USE_AUTODIFF)
 	typedef Matrix<variable> vMatrix; //shortcut for matrix with variations
+	typedef Matrix<hessian_variable> hMatrix; //shortcut for matrix with second variations
 #endif
 	
 }

Source/Headers/inmost_expression.h

@@ -334,7 +334,11 @@ namespace INMOST
 		{
 			return 1 + v[0].first;
 		}
-		
+		void Print() const
+		{
+			std::cout << value << std::endl;
+			entries.Print();
+		}
 		friend class multivar_expression_reference;
 	};
 	
@@ -401,12 +405,12 @@ namespace INMOST
 		__INLINE multivar_expression GetVariable(INMOST_DATA_ENUM_TYPE index)
 		{
 			multivar_expression ret(0);
-			for(int k = 0; k < entries.Size(); ++k)
+			for(int k = 0; k < (int)entries.Size(); ++k)
 				if( entries.GetIndex(k) == index )
 				{
 					ret.SetValue(entries.GetValue(k));
 					Sparse::Row & r = ret.GetRow();
-					for(int q = 0; q < hessian_entries.Size(); ++q)
+					for(int q = 0; q < (int)hessian_entries.Size(); ++q)
 					{
 						Sparse::HessianRow::index & i = hessian_entries.GetIndex(q);
 						if( i.first == index )
@@ -456,20 +460,24 @@ namespace INMOST
 		__INLINE hessian_multivar_expression & operator +=(basic_expression const & expr)
 		{
 			value += expr.GetValue();
-			Sparse::Row tmp(entries);
-			Sparse::HessianRow htmp(hessian_entries);
-			expr.GetHessian(1.0,tmp,1.0,htmp);
+			Sparse::Row tmpr, tmp;
+			Sparse::HessianRow htmpr, htmp;
+			expr.GetHessian(1.0,tmpr,1.0,htmpr);
+			Sparse::Row::MergeSortedRows(1.0,entries,1.0,tmpr,tmp);
 			entries.Swap(tmp);
+			Sparse::HessianRow::MergeSortedRows(1.0,hessian_entries,1.0,htmpr,htmp);
 			hessian_entries.Swap(htmp);
 			return *this;
 		}
 		__INLINE hessian_multivar_expression & operator -=(basic_expression const & expr)
 		{
 			value -= expr.GetValue();
-			Sparse::Row tmp(entries);
-			Sparse::HessianRow htmp(hessian_entries);
-			expr.GetHessian(-1.0,tmp,-1.0,htmp);
+			Sparse::Row tmpr, tmp;
+			Sparse::HessianRow htmpr, htmp;
+			expr.GetHessian(1.0,tmpr,1.0,htmpr);
+			Sparse::Row::MergeSortedRows(1.0,entries,-1.0,tmpr,tmp);
 			entries.Swap(tmp);
+			Sparse::HessianRow::MergeSortedRows(1.0,hessian_entries,-1.0,htmpr,htmp);
 			hessian_entries.Swap(htmp);
 			return *this;
 		}
@@ -1041,7 +1049,13 @@ namespace INMOST
     }
     __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
     {
-        arg.GetHessian(multJ*dmult,J,-multH*value,H);
+		Sparse::HessianRow htmp;
+        arg.GetHessian(1,J,1,htmp);
+		assert(J.isSorted());
+		assert(htmp.isSorted());
+		Sparse::HessianRow::MergeJacobianHessian(-value*multH,J,J,dmult*multH,htmp,H);
+		for(Sparse::Row::iterator it = J.Begin(); it != J.End(); ++it) it->second*=dmult*multJ;
+		assert(H.isSorted());
     }
 	};
 
@@ -1069,7 +1083,11 @@ namespace INMOST
     }
     __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
     {
-        arg.GetHessian(multJ*dmult,J,-multH*value,H);
+        //arg.GetHessian(multJ*dmult,J,-multH*value,H);
+		Sparse::HessianRow htmp;
+        arg.GetHessian(1,J,1,htmp);
+		Sparse::HessianRow::MergeJacobianHessian(-value*multH,J,J,dmult*multH,htmp,H);
+		for(Sparse::Row::iterator it = J.Begin(); it != J.End(); ++it) it->second*=dmult*multJ;
     }
 	};
 
@@ -1095,7 +1113,13 @@ namespace INMOST
     }
     __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
     {
-        arg.GetHessian(0.5*multJ/value,J,-0.25*multH/::pow(value,3),H);
+		//general formula:
+		// (F(G))'' = F'(G) G'' + F''(G) G'.G'
+		Sparse::HessianRow htmp;
+		arg.GetHessian(1,J,1,htmp);
+		Sparse::HessianRow::MergeJacobianHessian(-0.25/::pow(value,3.0)*multH,J,J,0.5/value*multH,htmp,H);
+		for(Sparse::Row::iterator it = J.Begin(); it != J.End(); ++it) it->second *= 0.5/value*multJ;
+        //arg.GetHessian(0.5*multJ/value,J,-0.25*multH/::pow(value,3),H);
     }
 	};
 
@@ -1253,18 +1277,25 @@ namespace INMOST
     }
     __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
     {
+		// (F*G)'' = (F'G+G'F)' = (F''G + F'G' + G''F + G'F') = (F''G + G''F + 2F'G')
         Sparse::Row JL, JR; //temporary jacobian rows from left and right expressions
         Sparse::HessianRow HL, HR; //temporary hessian rows form left and right expressions
-        left.GetHessian(multJ,JL,multH,HL); //retrive jacobian row and hessian matrix of the left expression
-        right.GetHessian(multJ,JR,multH,HR); //retrive jacobian row and hessian matrix of the right expression
+        left.GetHessian(1,JL,1,HL); //retrive jacobian row and hessian matrix of the left expression
+		assert(JL.isSorted());
+		assert(HL.isSorted());
+        right.GetHessian(1,JR,1,HR); //retrive jacobian row and hessian matrix of the right expression
+		assert(JR.isSorted());
+		assert(HR.isSorted());
         //assume rows are sorted (this is to be ensured by corresponding GetHessian functions)
         //preallocate J to JL.Size+JR.Size
         //perform merging of two sorted arrays
         //resize to correct size
-        Sparse::Row::MergeSortedRows(right.GetValue(),JL,left.GetValue(),JR,J);
+        Sparse::Row::MergeSortedRows(right.GetValue()*multJ,JL,left.GetValue()*multJ,JR,J);
+		assert(J.isSorted());
         //preallocate H to HL.Size+HR.Size+JL.Size*JR.Size
         //merge sorted
-        Sparse::HessianRow::MergeJacobianHessian(2.0,JL,JR,right.GetValue(),HL,left.GetValue(),HR,H);
+        Sparse::HessianRow::MergeJacobianHessian(2.0*multH,JL,JR,right.GetValue()*multH,HL,left.GetValue()*multH,HR,H);
+		assert(H.isSorted());
     }
   };
 
@@ -1340,10 +1371,16 @@ namespace INMOST
     {
         Sparse::Row JL, JR; //temporary jacobian rows from left and right expressions
         Sparse::HessianRow HL, HR; //temporary hessian rows form left and right expressions
-        left.GetHessian(multJ,JL,multH,HL); //retrive jacobian row and hessian matrix of the left expression
-        right.GetHessian(multJ,JR,multH,HR); //retrive jacobian row and hessian matrix of the right expression
-        Sparse::Row::MergeSortedRows(1.0,JL,1.0,JR,J);
-        Sparse::HessianRow::MergeSortedRows(1.0,HL,1.0,HR,H);
+        left.GetHessian(1,JL,1,HL); //retrive jacobian row and hessian matrix of the left expression
+		assert(JL.isSorted());
+		assert(HL.isSorted());
+        right.GetHessian(1,JR,1,HR); //retrive jacobian row and hessian matrix of the right expression
+		assert(JR.isSorted());
+		assert(HR.isSorted());
+        Sparse::Row::MergeSortedRows(multJ,JL,multJ,JR,J);
+		assert(J.isSorted());
+        Sparse::HessianRow::MergeSortedRows(multH,HL,multH,HR,H);
+		assert(H.isSorted());
     }
   };
 
@@ -1373,12 +1410,13 @@ namespace INMOST
     }
     __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
     {
+		//(F-G)'' = (F'-G')' = (F''-G'')
         Sparse::Row JL, JR; //temporary jacobian rows from left and right expressions
         Sparse::HessianRow HL, HR; //temporary hessian rows form left and right expressions
-        left.GetHessian(multJ,JL,multH,HL); //retrive jacobian row and hessian matrix of the left expression
-        right.GetHessian(multJ,JR,multH,HR); //retrive jacobian row and hessian matrix of the right expression
-        Sparse::Row::MergeSortedRows(1.0,JL,-1.0,JR,J);
-        Sparse::HessianRow::MergeSortedRows(1.0,HL,-1.0,HR,H);
+        left.GetHessian(1,JL,1,HL); //retrive jacobian row and hessian matrix of the left expression
+        right.GetHessian(1,JR,1,HR); //retrive jacobian row and hessian matrix of the right expression
+        Sparse::Row::MergeSortedRows(multJ,JL,-multJ,JR,J);
+        Sparse::HessianRow::MergeSortedRows(multH,HL,-multH,HR,H);
     }
   };
 

Source/Headers/inmost_sparse.h

@@ -319,9 +319,9 @@ namespace INMOST
             /// that allow for the modification of individual entries.
             /// @param size New size of the row.
             void                    Resize(INMOST_DATA_ENUM_TYPE size) {data.resize(size);}
-            void                    Print()
+            void                    Print() const
             {
-                for(iterator it = Begin(); it != End(); ++it) std::cout << "(" << it->first.first << "," << it->first.second << "," << it->second << ") ";
+                for(const_iterator it = Begin(); it != End(); ++it) std::cout << "(" << it->first.first << "," << it->first.second << "," << it->second << ") ";
                 std::cout << std::endl;
             }
             bool                    isSorted() const;
@@ -506,17 +506,10 @@ namespace INMOST
 		/// This class may be used to sum multiple sparse rows.
         /// \warning
         /// In parallel column indices of the matrix may span wider then
-        /// local row indices, to prevent any problem you are currently
-        /// advised to set total size of the matrix as interval of the
-        /// RowMerger. In future this may change, see todo 2 below.
-        /// \todo
-        /// 1. (testing!) Add iterators over entries.
-        /// 2. Implement multiple intervals for distributed computation,
-        ///    then in parallel the user may specify additional range of indexes
-        ///    for elements that lay on the borders between each pair of processors.
-        ///    Or even better implement mapping that will remap nonlocal entries to the
-        ///    end of current linked list when added and put them back in correct places
-        ///    when retrived. May use algorithm from class OrderInfo.
+        /// local row indices, to prevent any problem you can safely
+        /// set total size of the matrix as interval of the RowMerger.
+        /// For efficiency you should use RowMerger::SetNonlocal function
+		/// to provide information about non-local elements.
         class RowMerger
         {
         public:

Source/IO/mesh_xml_file.cpp

@@ -745,6 +745,11 @@ namespace INMOST
                         }
                       }
                     }
+					if( !sparse_read && t.GetSize() == ENUMUNDEF )
+					{
+						++it;
+						if( ((int)(it-set_elems)) < (int)set_size ) data = RealArray(*it,t);
+					}
                   }
                   break;
                 case DATA_INTEGER:
@@ -770,6 +775,11 @@ namespace INMOST
                         }
                       }
                     }
+					if( !sparse_read && t.GetSize() == ENUMUNDEF )
+					{
+						++it;
+						if( ((int)(it-set_elems)) < (int)set_size ) data = IntegerArray(*it,t);
+					}
                   }
                   break;
                 case DATA_BULK:
@@ -795,6 +805,11 @@ namespace INMOST
                         }
                       }