Commit 9689e3ce authored by Kirill Terekhov's avatar Kirill Terekhov

Fixes and features

Fix dot product of dense matrices, introduce submatrix operator for
dense matrices that allows to manipulate data in part of the matrix.

Do not store automatizator reference in dynamic_variable.

When only PERMX is provided in GRDECL file, consider it as scalar
permeability.
parent f8923774
......@@ -585,8 +585,8 @@ namespace INMOST
assert(Cols() == other.Cols());
assert(Rows() == other.Rows());
typename Promote<Var,typeB>::type ret = 0.0;
for(enumerator i = 0; i < Cols(); ++i)
for(enumerator j = 0; j < Rows(); ++j)
for(enumerator i = 0; i < Rows(); ++i)
for(enumerator j = 0; j < Cols(); ++j)
ret += ((*this)(i,j))*other(i,j);
return ret;
}
......@@ -1243,7 +1243,18 @@ namespace INMOST
/// @param first_col Starting column in the original matrix.
/// @param last_col Last column (excluded) in the original matrix.
/// @return Submatrix of the original matrix.
SubMatrix<Var,storage_type> MakeSubMatrix(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col);
//::INMOST::SubMatrix<Var,storage_type> SubMatrix(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col);
/// Extract submatrix of a matrix for in-place manipulation of elements.
/// Let A = {a_ij}, i in [0,n), j in [0,m) be original matrix.
/// Then the method returns B = {a_ij}, i in [ibeg,iend),
/// and j in [jbeg,jend).
/// @param first_row Starting row in the original matrix.
/// @param last_row Last row (excluded) in the original matrix.
/// @param first_col Starting column in the original matrix.
/// @param last_col Last column (excluded) in the original matrix.
/// @return Submatrix of the original matrix.
::INMOST::SubMatrix<Var,storage_type> operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col);
};
/// This class allows for in-place operations on submatrix of the matrix elements.
template<typename Var, typename Storage>
......@@ -1327,9 +1338,9 @@ namespace INMOST
/// not affect elements of the submatrix or original matrix
/// used to create submatrix.
/// @return Matrix with same entries as submatrix.
Matrix<Var> MakeMatrix()
Matrix<Var> Matrix()
{
Matrix<Var> ret(Rows(),Cols());
::INMOST::Matrix<Var> ret(Rows(),Cols());
for(enumerator i = 0; i < Rows(); ++i)
for(enumerator j = 0; j < Cols(); ++j)
ret(i,j) = (*this)(i,j);
......@@ -1550,11 +1561,12 @@ namespace INMOST
enumerator * order = new enumerator [m];
std::pair<Matrix<typename Promote<Var,typeB>::type>,bool>
ret = std::make_pair(Matrix<typename Promote<Var,typeB>::type>(m,l),true);
Var max, temp;
typename Promote<Var,typeB>::type tempb;
for(enumerator i = 0; i < m; ++i) order[i] = i;
for(enumerator i = 0; i < m; i++)
{
enumerator maxk = i, maxq = i, temp2;
Var max, temp;
max = fabs(AtA(maxk,maxq));
//Find best pivot
for(enumerator k = i; k < m; k++) // over rows
......@@ -1581,9 +1593,9 @@ namespace INMOST
//exchange rhs
for(enumerator q = 0; q < l; q++) // over columns of B
{
temp = AtB(maxk,q);
tempb = AtB(maxk,q);
AtB(maxk,q) = AtB(i,q);
AtB(i,q) = temp;
AtB(i,q) = tempb;
}
}
//Exchange columns
......@@ -1730,10 +1742,16 @@ namespace INMOST
return ret;
}
//template<typename Var,typename storage_type>
//SubMatrix<Var,storage_type> Matrix<Var,storage_type>::SubMatrix(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col)
//{
// return ::INMOST::SubMatrix<Var,storage_type>(*this,first_row,last_row,first_col,last_col);
//}
template<typename Var,typename storage_type>
SubMatrix<Var,storage_type> Matrix<Var,storage_type>::MakeSubMatrix(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col)
SubMatrix<Var,storage_type> Matrix<Var,storage_type>::operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col)
{
return SubMatrix<Var,storage_type>(*this,first_row,last_row,first_col,last_col);
return ::INMOST::SubMatrix<Var,storage_type>(*this,first_row,last_row,first_col,last_col);
}
/// shortcut for matrix of real values.
typedef Matrix<INMOST_DATA_REAL_TYPE> rMatrix;
......
......@@ -253,24 +253,23 @@ namespace INMOST
class dynamic_variable : public shell_dynamic_variable<var_expression,dynamic_variable >
{
private:
Automatizator & aut;
Tag index_tag, value_tag;
MarkerType mask;
INMOST_DATA_ENUM_TYPE comp;
public:
dynamic_variable(Automatizator & paut, INMOST_DATA_ENUM_TYPE pregval, INMOST_DATA_ENUM_TYPE pcomp = 0) : aut(paut), comp(pcomp)
dynamic_variable() : index_tag(),value_tag(),mask(0),comp(ENUMUNDEF) {}
dynamic_variable(Automatizator & paut, INMOST_DATA_ENUM_TYPE pregval, INMOST_DATA_ENUM_TYPE pcomp = 0) : comp(pcomp)
{
if( pregval != ENUMUNDEF )
{
mask = aut.GetMask(pregval);
value_tag = aut.GetValueTag(pregval);
index_tag = aut.GetIndexTag(pregval);
mask = paut.GetMask(pregval);
value_tag = paut.GetValueTag(pregval);
index_tag = paut.GetIndexTag(pregval);
}
}
dynamic_variable(const dynamic_variable & other) : aut(other.aut), index_tag(other.index_tag), value_tag(other.value_tag), mask(other.mask), comp(other.comp) {}
dynamic_variable(const dynamic_variable & other) : index_tag(other.index_tag), value_tag(other.value_tag), mask(other.mask), comp(other.comp) {}
dynamic_variable & operator =(const dynamic_variable & other)
{
aut = other.aut;
index_tag = other.index_tag;
value_tag = other.value_tag;
mask = other.mask;
......@@ -351,6 +350,7 @@ namespace INMOST
Tag variable_tag;
INMOST_DATA_ENUM_TYPE comp;
public:
stored_variable() : variable_tag(), comp(ENUMUNDEF) {}
stored_variable(Tag t, INMOST_DATA_ENUM_TYPE pcomp = 0) : variable_tag(t), comp(pcomp)
{
assert(t.GetDataType() == DATA_REAL || t.GetDataType() == DATA_VARIABLE);
......
......@@ -87,7 +87,9 @@
#define SEG_START 0
#define SEG_END 1
#define HAVE_PERM_X 1
#define HAVE_PERM_Y 2
#define HAVE_PERM_Z 4
namespace INMOST
{
......@@ -768,11 +770,12 @@ namespace INMOST
void Mesh::LoadECL(std::string File)
{
char have_perm = 0;
std::cout << std::scientific;
bool perform_splitting = false;
bool curvilinear_edges = true;
bool triangulated_edges = false; //this is not working, yet
bool project_perm = true;
bool project_perm = false;
for(INMOST_DATA_ENUM_TYPE k = 0; k < file_options.size(); ++k)
{
if( file_options[k].first == "ECL_SPLIT_GLUED" )
......@@ -943,6 +946,7 @@ namespace INMOST
downread = totread = dims[0]*dims[1]*dims[2];
argtype = ECL_VAR_REAL;
offset = state = ECL_PERMX;
have_perm |= HAVE_PERM_X;
}
else if( !ECLSTRCMP(p,"PERMY") )
{
......@@ -954,6 +958,7 @@ namespace INMOST
argtype = ECL_VAR_REAL;
offset = ECL_PERMX;
state = ECL_PERMY;
have_perm |= HAVE_PERM_Y;
}
else if( !ECLSTRCMP(p,"PERMZ") )
{
......@@ -965,6 +970,7 @@ namespace INMOST
argtype = ECL_VAR_REAL;
offset = ECL_PERMX;
state = ECL_PERMZ;
have_perm |= HAVE_PERM_Z;
}
else if( !ECLSTRCMP(p,"PORO") )
{
......@@ -1163,6 +1169,58 @@ ecl_exit_loop:
{
std::cout << __FILE__ << ":" << __LINE__ << " radial grids not supported yet " << std::endl;
}
if( !have_perm ) perm.clear();
if( !(HAVE_PERM_X & have_perm) )
{
if( HAVE_PERM_Y & have_perm )
{
for(int i = 0; i < dims[0]; ++i)
for(int j = 0; j < dims[1]; ++j)
for(int k = 0; k < dims[2]; ++k)
perm[3*ECL_IJK_DATA(i,j,k)+0] = perm[3*ECL_IJK_DATA(i,j,k)+1];
}
if( HAVE_PERM_Z & have_perm )
{
for(int i = 0; i < dims[0]; ++i)
for(int j = 0; j < dims[1]; ++j)
for(int k = 0; k < dims[2]; ++k)
perm[3*ECL_IJK_DATA(i,j,k)+0] = perm[3*ECL_IJK_DATA(i,j,k)+2];
}
}
if( !(HAVE_PERM_Y & have_perm) )
{
if( HAVE_PERM_X & have_perm )
{
for(int i = 0; i < dims[0]; ++i)
for(int j = 0; j < dims[1]; ++j)
for(int k = 0; k < dims[2]; ++k)
perm[3*ECL_IJK_DATA(i,j,k)+1] = perm[3*ECL_IJK_DATA(i,j,k)+0];
}
if( HAVE_PERM_Z & have_perm )
{
for(int i = 0; i < dims[0]; ++i)
for(int j = 0; j < dims[1]; ++j)
for(int k = 0; k < dims[2]; ++k)
perm[3*ECL_IJK_DATA(i,j,k)+1] = perm[3*ECL_IJK_DATA(i,j,k)+2];
}
}
if( !(HAVE_PERM_Z & have_perm) )
{
if( HAVE_PERM_X & have_perm )
{
for(int i = 0; i < dims[0]; ++i)
for(int j = 0; j < dims[1]; ++j)
for(int k = 0; k < dims[2]; ++k)
perm[3*ECL_IJK_DATA(i,j,k)+2] = perm[3*ECL_IJK_DATA(i,j,k)+0];
}
if( HAVE_PERM_Y & have_perm )
{
for(int i = 0; i < dims[0]; ++i)
for(int j = 0; j < dims[1]; ++j)
for(int k = 0; k < dims[2]; ++k)
perm[3*ECL_IJK_DATA(i,j,k)+2] = perm[3*ECL_IJK_DATA(i,j,k)+1];
}
}
int beg_dims[2], end_dims[2];
beg_dims[0] = 0;
beg_dims[1] = 0;
......
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