Commit d9ce2066 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Fixes and features

Added operations for matrix with hessian calculation

Fix a bug in submatrix calculation

Fix a bug in reading data of dense type with variable size from files
in xml format

Fix an incorrect warning when reading variables with variations from
files in xml format
parent f6feda49
...@@ -797,7 +797,7 @@ public: ...@@ -797,7 +797,7 @@ public:
colors.push_back(color_t(0.45,0,0.55)); colors.push_back(color_t(0.45,0,0.55));
colors.push_back(color_t(0,0,0)); colors.push_back(color_t(0,0,0));
samples = 2048; samples = 4096;
float * pixel_array = new float[(samples+2)*4]; float * pixel_array = new float[(samples+2)*4];
...@@ -808,8 +808,8 @@ public: ...@@ -808,8 +808,8 @@ public:
float t = 1.0f*q/static_cast<float>(samples+1); float t = 1.0f*q/static_cast<float>(samples+1);
color_t c = pick_color(t); color_t c = pick_color(t);
//countour lines //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); // c = pick_color(1-t) + color_t(0,2*t*(1-t),0);
pixel_array[(q)*4+0] = c.r(); pixel_array[(q)*4+0] = c.r();
pixel_array[(q)*4+1] = c.g(); pixel_array[(q)*4+1] = c.g();
...@@ -4093,11 +4093,13 @@ double display_elem_info(Element e, double top, double left, double interval) ...@@ -4093,11 +4093,13 @@ double display_elem_info(Element e, double top, double left, double interval)
char str[4096]; char str[4096];
char temp[4096]; char temp[4096];
str[0] = '\0'; str[0] = '\0';
int dsize;
switch(t->GetDataType()) switch(t->GetDataType())
{ {
case DATA_INTEGER: case DATA_INTEGER:
{ {
Storage::integer_array arr = e->IntegerArray(*t); Storage::integer_array arr = e->IntegerArray(*t);
dsize = arr.size();
for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++) for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
{ {
sprintf(temp,"%s %d",str,arr[k]); sprintf(temp,"%s %d",str,arr[k]);
...@@ -4108,6 +4110,7 @@ double display_elem_info(Element e, double top, double left, double interval) ...@@ -4108,6 +4110,7 @@ double display_elem_info(Element e, double top, double left, double interval)
case DATA_REAL: case DATA_REAL:
{ {
Storage::real_array arr = e->RealArray(*t); Storage::real_array arr = e->RealArray(*t);
dsize = arr.size();
for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++) for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
{ {
sprintf(temp,"%s %lf",str,arr[k]); sprintf(temp,"%s %lf",str,arr[k]);
...@@ -4118,6 +4121,7 @@ double display_elem_info(Element e, double top, double left, double interval) ...@@ -4118,6 +4121,7 @@ double display_elem_info(Element e, double top, double left, double interval)
case DATA_BULK: case DATA_BULK:
{ {
Storage::bulk_array arr = e->BulkArray(*t); Storage::bulk_array arr = e->BulkArray(*t);
dsize = arr.size();
for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++) for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
{ {
sprintf(temp,"%s %d",str,arr[k]); sprintf(temp,"%s %d",str,arr[k]);
...@@ -4128,6 +4132,7 @@ double display_elem_info(Element e, double top, double left, double interval) ...@@ -4128,6 +4132,7 @@ double display_elem_info(Element e, double top, double left, double interval)
case DATA_REFERENCE: case DATA_REFERENCE:
{ {
Storage::reference_array arr = e->ReferenceArray(*t); Storage::reference_array arr = e->ReferenceArray(*t);
dsize = arr.size();
for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++) for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
{ {
if(arr.at(k) == InvalidHandle()) sprintf(temp,"%s NULL",str); 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) ...@@ -4139,6 +4144,7 @@ double display_elem_info(Element e, double top, double left, double interval)
case DATA_REMOTE_REFERENCE: case DATA_REMOTE_REFERENCE:
{ {
Storage::remote_reference_array arr = e->RemoteReferenceArray(*t); Storage::remote_reference_array arr = e->RemoteReferenceArray(*t);
dsize = arr.size();
for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++) 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); 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) ...@@ -4151,6 +4157,7 @@ double display_elem_info(Element e, double top, double left, double interval)
case DATA_VARIABLE: case DATA_VARIABLE:
{ {
Storage::var_array arr = e->VariableArray(*t); Storage::var_array arr = e->VariableArray(*t);
dsize = arr.size();
for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++) for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
{ {
std::stringstream stream; std::stringstream stream;
...@@ -4167,7 +4174,7 @@ double display_elem_info(Element e, double top, double left, double interval) ...@@ -4167,7 +4174,7 @@ double display_elem_info(Element e, double top, double left, double interval)
} }
#endif #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); strcpy(str,temp);
top -= interval; top -= interval;
glRasterPos2f(left,top); glRasterPos2f(left,top);
...@@ -5036,7 +5043,7 @@ int main(int argc, char ** argv) ...@@ -5036,7 +5043,7 @@ int main(int argc, char ** argv)
//if( false ) //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()); printf("preparing octree around mesh, was sets %d\n",mesh->NumberOfSets());
Octree octsearch = Octree(mesh->CreateSet("octsearch").first); Octree octsearch = Octree(mesh->CreateSet("octsearch").first);
......
...@@ -30,6 +30,11 @@ namespace INMOST ...@@ -30,6 +30,11 @@ namespace INMOST
template<> struct Promote<INMOST_DATA_REAL_TYPE, variable> {typedef variable type;}; 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, INMOST_DATA_REAL_TYPE> {typedef variable type;};
template<> struct Promote<variable, variable> {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 #else
__INLINE INMOST_DATA_REAL_TYPE get_value(INMOST_DATA_REAL_TYPE x) {return x;} __INLINE INMOST_DATA_REAL_TYPE get_value(INMOST_DATA_REAL_TYPE x) {return x;}
#endif #endif
...@@ -868,16 +873,30 @@ namespace INMOST ...@@ -868,16 +873,30 @@ namespace INMOST
break; break;
case 21: //symmetric elasticity tensor (note - diagonal first, then off-diagonal rows) case 21: //symmetric elasticity tensor (note - diagonal first, then off-diagonal rows)
{ {
int shift[5] = {6,11,15,18,20}; Kc(0,0) = K[0]; //c11
for(int i = 0; i < 6; ++i) Kc(0,1) = Kc(1,0) = K[1]; //c12
{ Kc(0,2) = Kc(2,0) = K[2]; //c13
Kc(i,i) = K[i]; //diagonal term first Kc(0,3) = Kc(3,0) = K[3]; //c14
for(int j = i+1; j < 6; ++j) Kc(0,4) = Kc(4,0) = K[4]; //c15
Kc(i,j) = Kc(j,i) = K[shift[i]+j-i-1]; //off diagonal parts 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,3) = K[19]; //c56
Kc(5,5) = K[20]; //c66
break; break;
} }
case 36: //full permeability tensor case 36: //full elasticity tensor
for(int i = 0; i < 6; ++i) for(int i = 0; i < 6; ++i)
for(int j = 0; j < 6; ++j) for(int j = 0; j < 6; ++j)
Kc(i,j) = K[6*i+j]; Kc(i,j) = K[6*i+j];
...@@ -1082,7 +1101,7 @@ namespace INMOST ...@@ -1082,7 +1101,7 @@ namespace INMOST
for(enumerator i = ibeg; i < iend; ++i) for(enumerator i = ibeg; i < iend; ++i)
{ {
for(enumerator j = jbeg; j < jend; ++j) 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; return ret;
} }
...@@ -1091,6 +1110,7 @@ namespace INMOST ...@@ -1091,6 +1110,7 @@ namespace INMOST
typedef Matrix<INMOST_DATA_REAL_TYPE> rMatrix; //shortcut for real matrix typedef Matrix<INMOST_DATA_REAL_TYPE> rMatrix; //shortcut for real matrix
#if defined(USE_AUTODIFF) #if defined(USE_AUTODIFF)
typedef Matrix<variable> vMatrix; //shortcut for matrix with variations typedef Matrix<variable> vMatrix; //shortcut for matrix with variations
typedef Matrix<hessian_variable> hMatrix; //shortcut for matrix with second variations
#endif #endif
} }
......
...@@ -334,7 +334,11 @@ namespace INMOST ...@@ -334,7 +334,11 @@ namespace INMOST
{ {
return 1 + v[0].first; return 1 + v[0].first;
} }
void Print() const
{
std::cout << value << std::endl;
entries.Print();
}
friend class multivar_expression_reference; friend class multivar_expression_reference;
}; };
......
...@@ -726,7 +726,7 @@ namespace INMOST ...@@ -726,7 +726,7 @@ namespace INMOST
Storage::real_array data = RealArray(*it,t); Storage::real_array data = RealArray(*it,t);
std::vector<Storage::real> Vector; int Repeat; std::vector<Storage::real> Vector; int Repeat;
reader.ParseReal(val,Vector,Repeat,set_size); reader.ParseReal(val,Vector,Repeat,set_size);
if( t.GetSize() == ENUMUNDEF ) data.resize((enumerator)Vector.size()*Repeat); if( t.GetSize() == ENUMUNDEF ) data.resize((enumerator)Vector.size()*Repeat);
else if( t.GetSize() != Vector.size()*Repeat && sparse_read ) else if( t.GetSize() != Vector.size()*Repeat && sparse_read )
{ {
reader.Report("Cannot write record of size %d into tag data of size %d",Vector.size()*Repeat,t.GetSize()); reader.Report("Cannot write record of size %d into tag data of size %d",Vector.size()*Repeat,t.GetSize());
...@@ -745,6 +745,11 @@ namespace INMOST ...@@ -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; break;
case DATA_INTEGER: case DATA_INTEGER:
...@@ -770,6 +775,11 @@ namespace INMOST ...@@ -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; break;
case DATA_BULK: case DATA_BULK:
...@@ -795,6 +805,11 @@ namespace INMOST ...@@ -795,6 +805,11 @@ namespace INMOST
} }
} }
} }
if( !sparse_read && t.GetSize() == ENUMUNDEF )
{
++it;
if( ((int)(it-set_elems)) < (int)set_size ) data = BulkArray(*it,t);
}
} }
break; break;
case DATA_REFERENCE: case DATA_REFERENCE:
...@@ -820,6 +835,11 @@ namespace INMOST ...@@ -820,6 +835,11 @@ namespace INMOST
} }
} }
} }
if( !sparse_read && t.GetSize() == ENUMUNDEF )
{
++it;
if( ((int)(it-set_elems)) < (int)set_size ) data = ReferenceArray(*it,t);
}
} }
break; break;
case DATA_REMOTE_REFERENCE: case DATA_REMOTE_REFERENCE:
...@@ -848,6 +868,11 @@ namespace INMOST ...@@ -848,6 +868,11 @@ namespace INMOST
} }
} }
} }
if( !sparse_read && t.GetSize() == ENUMUNDEF )
{
++it;
if( ((int)(it-set_elems)) < (int)set_size ) data = RemoteReferenceArray(*it,t);
}
} }
else else
{ {
...@@ -874,6 +899,11 @@ namespace INMOST ...@@ -874,6 +899,11 @@ namespace INMOST
} }
} }
} }
if( !sparse_read && t.GetSize() == ENUMUNDEF )
{
++it;
if( ((int)(it-set_elems)) < (int)set_size ) data = RemoteReferenceArray(*it,t);
}
} }
} }
break; break;
...@@ -901,6 +931,11 @@ namespace INMOST ...@@ -901,6 +931,11 @@ namespace INMOST
} }
} }
} }
if( !sparse_read && t.GetSize() == ENUMUNDEF )
{
++it;
if( ((int)(it-set_elems)) < (int)set_size ) data = VariableArray(*it,t);
}
} }
break; break;
#endif #endif
......
...@@ -1181,7 +1181,7 @@ namespace INMOST ...@@ -1181,7 +1181,7 @@ namespace INMOST
{ {
INMOST::Sparse::Row entries; INMOST::Sparse::Row entries;
INMOST::Storage::real val; INMOST::Storage::real val;
if( !(_str[0] == '(' && _str[1] == ')' ) ) Report("Expected scopes for variable %s",_str); if( !(_str[0] == '(' && _str[strlen(_str)-1] == ')' ) ) Report("Expected scopes for variable %s",_str);
std::string str(_str+1,strlen(_str)-2); std::string str(_str+1,strlen(_str)-2);
std::vector<std::string> decomposed; std::vector<std::string> decomposed;
ParseCommaSeparated(str,decomposed,';'); ParseCommaSeparated(str,decomposed,';');
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment