Commit 68ec2014 authored by Kirill Terekhov's avatar Kirill Terekhov

Updates

Improvements for eclipse grid reader.

Fixed a bug in ElementArray::iterator.
parent 9e65b7d1
......@@ -1223,7 +1223,10 @@ std::vector<face2gl> clip_boundary;
void draw_faces_nc(std::vector<face2gl> & set, int highlight = -1)
{
if( drawedges == 2 || drawedges == 3 ) return;
if( drawedges == 2 || drawedges == 3 ) return;
glEnable(GL_POLYGON_OFFSET_FILL);
glPolygonOffset(1.0, 1.0);
glColor4f(0,1,0,0.1);
glBegin(GL_TRIANGLES);
for(INMOST_DATA_ENUM_TYPE q = 0; q < set.size() ; q++) set[q].draw();
......@@ -1235,12 +1238,17 @@ void draw_faces_nc(std::vector<face2gl> & set, int highlight = -1)
set[highlight].draw();
glEnd();
}
glDisable(GL_POLYGON_OFFSET_FILL);
}
void draw_faces(std::vector<face2gl> & set, int highlight = -1)
{
if( drawedges == 2 || drawedges == 3) return;
glEnable(GL_POLYGON_OFFSET_FILL);
glPolygonOffset(1.0, 1.0);
if( visualization_tag.isValid() ) CommonColorBar->BindTexture();
glBegin(GL_TRIANGLES);
for(INMOST_DATA_ENUM_TYPE q = 0; q < set.size() ; q++) set[q].draw_colour();
glEnd();
......@@ -1252,16 +1260,25 @@ void draw_faces(std::vector<face2gl> & set, int highlight = -1)
set[highlight].draw();
glEnd();
}
glDisable(GL_POLYGON_OFFSET_FILL);
}
void draw_faces_alpha(std::vector<face2gl> & set, double alpha)
{
if( drawedges == 2 || drawedges==3) return;
if( visualization_tag.isValid() ) CommonColorBar->BindTexture();
glEnable(GL_POLYGON_OFFSET_FILL);
glPolygonOffset(1.0, 1.0);
glBegin(GL_TRIANGLES);
for(INMOST_DATA_ENUM_TYPE q = 0; q < set.size() ; q++) set[q].draw_colour_alpha(alpha);
glEnd();
if( visualization_tag.isValid() ) CommonColorBar->UnbindTexture();
glDisable(GL_POLYGON_OFFSET_FILL);
}
void draw_edges(std::vector<face2gl> & set, int highlight = -1)
......@@ -1283,31 +1300,53 @@ void draw_edges(std::vector<face2gl> & set, int highlight = -1)
void draw_faces_interactive_nc(std::vector<face2gl> & set)
{
if( drawedges == 2 || drawedges==3 ) return;
if( drawedges == 2 || drawedges==3 ) return;
glColor4f(0,1,0,0.1);
glEnable(GL_POLYGON_OFFSET_FILL);
glPolygonOffset(1.0, 1.0);
glBegin(GL_TRIANGLES);
for(INMOST_DATA_ENUM_TYPE q = 0; q < set.size() ; q++) if( set[q].get_flag() ) set[q].draw();
glEnd();
glDisable(GL_POLYGON_OFFSET_FILL);
}
void draw_faces_interactive(std::vector<face2gl> & set)
{
if( drawedges == 2 || drawedges==3 ) return;
if( visualization_tag.isValid() ) CommonColorBar->BindTexture();
glEnable(GL_POLYGON_OFFSET_FILL);
glPolygonOffset(1.0, 1.0);
glBegin(GL_TRIANGLES);
for(INMOST_DATA_ENUM_TYPE q = 0; q < set.size() ; q++) if( set[q].get_flag() ) set[q].draw_colour();
glEnd();
if( visualization_tag.isValid() ) CommonColorBar->UnbindTexture();
glDisable(GL_POLYGON_OFFSET_FILL);
}
void draw_faces_interactive_alpha(std::vector<face2gl> & set, double alpha)
{
if( drawedges == 2 || drawedges==3 ) return;
if( visualization_tag.isValid() ) CommonColorBar->BindTexture();
glEnable(GL_POLYGON_OFFSET_FILL);
glPolygonOffset(1.0, 1.0);
glBegin(GL_TRIANGLES);
for(INMOST_DATA_ENUM_TYPE q = 0; q < set.size() ; q++) if( set[q].get_flag() ) set[q].draw_colour_alpha(alpha);
glEnd();
if( visualization_tag.isValid() ) CommonColorBar->UnbindTexture();
glDisable(GL_POLYGON_OFFSET_FILL);
}
void draw_edges_interactive(std::vector<face2gl> & set)
......@@ -2329,6 +2368,9 @@ public:
void draw_clip(INMOST_DATA_ENUM_TYPE pace, Storage::real n[3])
{
if( visualization_tag.isValid() ) CommonColorBar->BindTexture();
glEnable(GL_POLYGON_OFFSET_FILL);
glPolygonOffset(1.0, 1.0);
glBegin(GL_TRIANGLES);
for(INMOST_DATA_ENUM_TYPE k = 0; k < cells.size(); k+=pace) if( cells[k]->GetMarker(marker))
{
......@@ -2405,6 +2447,8 @@ public:
}
}
glEnd();
glDisable(GL_POLYGON_OFFSET_FILL);
if( visualization_tag.isValid() ) CommonColorBar->UnbindTexture();
}
void draw_clip_edges(INMOST_DATA_ENUM_TYPE pace, Storage::real n[3])
......@@ -2608,6 +2652,11 @@ public:
if( state == CLIP_FACE_INSIDE )
{
ElementArray<Node> nodes = Face(mm,faces[k])->getNodes();
glEnable(GL_POLYGON_OFFSET_FILL);
glPolygonOffset(1.0, 1.0);
glBegin(GL_POLYGON);
if( visualization_tag.isValid() )
{
......@@ -2624,6 +2673,8 @@ public:
glVertex3dv(&nodes[q].Coords()[0]);
}
glEnd();
glDisable(GL_POLYGON_OFFSET_FILL);
}
}
mm->IntegerDF(faces[k],clip_state) = state;
......@@ -2899,7 +2950,7 @@ void set_matrix3d()
else
{
const double pi = 3.1415926535897932384626433832795;
const double znear = 1;
const double znear = zoom*2*std::max(std::max( sright-sleft, stop-sbottom ), sfar-snear)*0.1;
const double zfar = 1000000.0;
const double fH = znear * tan( 60.0 / 360.0 * pi);
double fW = fH * aspect;
......@@ -3775,6 +3826,7 @@ double display_elem_info(Element e, double top, double left, double interval)
void draw_screen()
{
//glDepthMask(GL_TRUE);
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glLoadIdentity();
set_matrix3d();
......
......@@ -75,9 +75,9 @@ namespace INMOST
_iterator(etype * i):e(i){}
_iterator(const _iterator & other){e = other.e;}
~_iterator() {};
_iterator operator -(size_t n) { return _iterator(e-n); }
_iterator operator -(size_t n) const { return _iterator(e-n); }
_iterator & operator -=(size_t n) { e-=n; return *this; }
_iterator operator +(size_t n) { return _iterator(e+n); }
_iterator operator +(size_t n) const { return _iterator(e+n); }
_iterator & operator +=(size_t n) { e+=n; return *this; }
_iterator & operator ++(){ ++e; return *this;}
_iterator operator ++(int){ return _iterator(e++); }
......@@ -94,7 +94,7 @@ namespace INMOST
bool operator <=(const _iterator & other) const { return e <= other.e;}
bool operator >=(const _iterator & other) const { return e >= other.e;}
operator void *() {return static_cast<void *> (e);}
operator const void *() {return const_cast<const void *> (e);}
operator const void *() const {return const_cast<const void *> (e);}
};
typedef _iterator<element> iterator;
typedef _iterator<const element> const_iterator;
......@@ -113,9 +113,9 @@ namespace INMOST
_reverse_iterator(etype * i):e(i){}
_reverse_iterator(const _reverse_iterator & other){e = other.e;}
~_reverse_iterator() {};
_reverse_iterator operator -(size_t n) { return _reverse_iterator(e+n); }
_reverse_iterator operator -(size_t n) const { return _reverse_iterator(e+n); }
_reverse_iterator & operator -=(size_t n) { e+=n; return *this; }
_reverse_iterator operator +(size_t n) {return _reverse_iterator(e-n); }
_reverse_iterator operator +(size_t n) const {return _reverse_iterator(e-n); }
_reverse_iterator & operator +=(size_t n) { e-=n; return *this; }
_reverse_iterator & operator ++(){ --e; return *this;}
_reverse_iterator operator ++(int){ return _reverse_iterator(e--); }
......@@ -132,6 +132,7 @@ namespace INMOST
bool operator <=(const _reverse_iterator & other) const { return e <= other.e;}
bool operator >=(const _reverse_iterator & other) const { return e >= other.e;}
operator void *() {return static_cast<void *> (e);}
operator const void *() const {return static_cast<const void *> (e);}
};
typedef _reverse_iterator<element> reverse_iterator;
typedef _reverse_iterator<const element> const_reverse_iterator;
......@@ -536,9 +537,9 @@ namespace INMOST
_iterator(dtype * i):e(i){}
_iterator(const _iterator & other){e = other.e;}
~_iterator() {};
_iterator operator -(size_t n) { return _iterator(e-n); }
_iterator operator -(size_t n) const { return _iterator(e-n); }
_iterator & operator -=(size_t n) { e-=n; return *this; }
_iterator operator +(size_t n) { return _iterator(e+n); }
_iterator operator +(size_t n) const { return _iterator(e+n); }
_iterator & operator +=(size_t n) { e+=n; return *this; }
_iterator & operator ++(){ ++e; return *this;}
_iterator operator ++(int){ return _iterator(e++); }
......@@ -573,9 +574,9 @@ namespace INMOST
_reverse_iterator(dtype * i):e(i){}
_reverse_iterator(const _reverse_iterator & other){e = other.e;}
~_reverse_iterator() {};
_reverse_iterator operator -(size_t n) { return _reverse_iterator(e+n); }
_reverse_iterator operator -(size_t n) const { return _reverse_iterator(e+n); }
_reverse_iterator & operator -=(size_t n) { e+=n; return *this; }
_reverse_iterator operator +(size_t n) {return _reverse_iterator(e-n); }
_reverse_iterator operator +(size_t n) const {return _reverse_iterator(e-n); }
_reverse_iterator & operator +=(size_t n) { e-=n; return *this; }
_reverse_iterator & operator ++(){ --e; return *this;}
_reverse_iterator operator ++(int){ return _reverse_iterator(e--); }
......@@ -1547,9 +1548,9 @@ namespace INMOST
_iterator(dtype * i):e(i){}
_iterator(const _iterator & other){e = other.e;}
~_iterator() {};
_iterator operator -(size_type n) { return _iterator(e-n); }
_iterator operator -(size_type n) const { return _iterator(e-n); }
_iterator & operator -=(size_type n) { e-=n; return *this; }
_iterator operator +(size_type n) { return _iterator(e+n); }
_iterator operator +(size_type n) const { return _iterator(e+n); }
_iterator & operator +=(size_type n) { e+=n; return *this; }
_iterator & operator ++(){ ++e; return *this;}
_iterator operator ++(int){ return _iterator(e++); }
......@@ -1559,14 +1560,14 @@ namespace INMOST
dtype & operator *() { return *e; }
dtype * operator ->() { return e; }
_iterator & operator =(_iterator const & other) { e = other.e; return *this; }
bool operator ==(const _iterator & other) { return e == other.e;}
bool operator !=(const _iterator & other) { return e != other.e;}
bool operator <(const _iterator & other) { return e < other.e;}
bool operator >(const _iterator & other) { return e > other.e;}
bool operator <=(const _iterator & other) { return e <= other.e;}
bool operator >=(const _iterator & other) { return e >= other.e;}
bool operator ==(const _iterator & other) const { return e == other.e;}
bool operator !=(const _iterator & other) const { return e != other.e;}
bool operator <(const _iterator & other) const { return e < other.e;}
bool operator >(const _iterator & other) const { return e > other.e;}
bool operator <=(const _iterator & other) const { return e <= other.e;}
bool operator >=(const _iterator & other) const { return e >= other.e;}
operator void *() {return static_cast<void *> (e);}
void * cast_to_void(){return static_cast<void *>(e);}
operator const void *() const {return static_cast<const void *> (e);}
};
typedef _iterator<element> iterator;
typedef _iterator<const element> const_iterator;
......@@ -1585,9 +1586,9 @@ namespace INMOST
_reverse_iterator(dtype * i):e(i){}
_reverse_iterator(const _reverse_iterator & other){e = other.e;}
~_reverse_iterator() {};
_reverse_iterator operator -(size_type n) { return _reverse_iterator(e+n); }
_reverse_iterator operator -(size_type n) const { return _reverse_iterator(e+n); }
_reverse_iterator & operator -=(size_type n) { e+=n; return *this; }
_reverse_iterator operator +(size_type n) {return _reverse_iterator(e-n); }
_reverse_iterator operator +(size_type n) const {return _reverse_iterator(e-n); }
_reverse_iterator & operator +=(size_type n) { e-=n; return *this; }
_reverse_iterator & operator ++(){ --e; return *this;}
_reverse_iterator operator ++(int){ return _reverse_iterator(e--); }
......@@ -1597,14 +1598,14 @@ namespace INMOST
dtype & operator *() { return *e; }
dtype * operator ->() { return e; }
_reverse_iterator & operator =(_reverse_iterator const & other) { e = other.e; return *this;}
bool operator ==(const _reverse_iterator & other) { return e == other.e;}
bool operator !=(const _reverse_iterator & other) { return e != other.e;}
bool operator <(const _reverse_iterator & other) { return e < other.e;}
bool operator >(const _reverse_iterator & other) { return e > other.e;}
bool operator <=(const _reverse_iterator & other) { return e <= other.e;}
bool operator >=(const _reverse_iterator & other) { return e >= other.e;}
bool operator ==(const _reverse_iterator & other) const { return e == other.e;}
bool operator !=(const _reverse_iterator & other) const { return e != other.e;}
bool operator <(const _reverse_iterator & other) const { return e < other.e;}
bool operator >(const _reverse_iterator & other) const { return e > other.e;}
bool operator <=(const _reverse_iterator & other) const { return e <= other.e;}
bool operator >=(const _reverse_iterator & other) const { return e >= other.e;}
operator void *() {return static_cast<void *> (e);}
void * cast_to_void(){return static_cast<void *>(e);}
operator const void *() const {return static_cast<const void *> (e);}
};
typedef _reverse_iterator<element> reverse_iterator;
typedef _reverse_iterator<const element> const_reverse_iterator;
......@@ -2037,12 +2038,12 @@ namespace INMOST
dtype & operator *() { return *it; }
dtype * operator ->() { return &*it; }
_iterator & operator =(_iterator const & other) {it = other.it; cur_list = other.cur_list; lists = other.lists; return *this; }
bool operator ==(const _iterator & other) { return it == other.it;}
bool operator !=(const _iterator & other) { return it != other.it;}
bool operator <(const _iterator & other) { return it < other.it;}
bool operator >(const _iterator & other) { return it > other.it;}
bool operator <=(const _iterator & other) { return it <= other.it;}
bool operator >=(const _iterator & other) { return it >= other.it;}
bool operator ==(const _iterator & other) const { return it == other.it;}
bool operator !=(const _iterator & other) const { return it != other.it;}
bool operator <(const _iterator & other) const { return it < other.it;}
bool operator >(const _iterator & other) const { return it > other.it;}
bool operator <=(const _iterator & other) const { return it <= other.it;}
bool operator >=(const _iterator & other) const { return it >= other.it;}
};
typedef _iterator<inner_type> iterator;
typedef _iterator<const inner_type> const_iterator;
......@@ -2189,9 +2190,9 @@ namespace INMOST
iterator(chunk_array<element,block_bits> * c, size_type pos):link(c),pos(pos){}
iterator(const iterator & other){link = other.link; pos = other.pos;}
~iterator() {};
iterator operator -(size_type n) { return iterator(link,pos-n); }
iterator operator -(size_type n) const { return iterator(link,pos-n); }
iterator & operator -=(size_type n) { pos-=n; return *this; }
iterator operator +(size_type n) { return iterator(link,pos+n); }
iterator operator +(size_type n) const { return iterator(link,pos+n); }
iterator & operator +=(size_type n) { pos+=n; return *this; }
iterator & operator ++(){ ++pos; return *this;}
iterator operator ++(int){ return iterator(link,pos++); }
......@@ -2201,12 +2202,12 @@ namespace INMOST
element & operator *() { return link->at(pos); }
element * operator ->() { return &link->at(pos); }
iterator & operator =(iterator const & other) { link = other.link; pos = other.pos; return *this; }
bool operator ==(const iterator & other) { assert(link == other.link); return pos == other.pos;}
bool operator !=(const iterator & other) { assert(link == other.link); return pos != other.pos;}
bool operator <(const iterator & other) { assert(link == other.link); return pos < other.pos;}
bool operator >(const iterator & other) { assert(link == other.link); return pos > other.pos;}
bool operator <=(const iterator & other) { assert(link == other.link); return pos <= other.pos;}
bool operator >=(const iterator & other) { assert(link == other.link); return pos >= other.pos;}
bool operator ==(const iterator & other) const { assert(link == other.link); return pos == other.pos;}
bool operator !=(const iterator & other) const { assert(link == other.link); return pos != other.pos;}
bool operator <(const iterator & other) const { assert(link == other.link); return pos < other.pos;}
bool operator >(const iterator & other) const { assert(link == other.link); return pos > other.pos;}
bool operator <=(const iterator & other) const { assert(link == other.link); return pos <= other.pos;}
bool operator >=(const iterator & other) const { assert(link == other.link); return pos >= other.pos;}
};
class const_iterator
......@@ -2224,9 +2225,9 @@ namespace INMOST
const_iterator(const chunk_array<element,block_bits> * c, size_type pos):link(c),pos(pos){}
const_iterator(const const_iterator & other) :link(other.link) {pos = other.pos;}
~const_iterator() {};
const_iterator operator -(size_type n) { return const_iterator(link,pos-n); }
const_iterator operator -(size_type n) const { return const_iterator(link,pos-n); }
const_iterator & operator -=(size_type n) { pos-=n; return *this; }
const_iterator operator +(size_type n) { return const_iterator(link,pos+n); }
const_iterator operator +(size_type n) const { return const_iterator(link,pos+n); }
const_iterator & operator +=(size_type n) { pos+=n; return *this; }
const_iterator & operator ++(){ ++pos; return *this;}
const_iterator operator ++(int){ return const_iterator(link,pos++); }
......@@ -2236,12 +2237,12 @@ namespace INMOST
const element & operator *() { return link->at(pos); }
const element * operator ->() { return &link->at(pos); }
const_iterator & operator =(const_iterator const & other) { link = other.link; pos = other.pos; return *this; }
bool operator ==(const const_iterator & other) { assert(link == other.link); return pos == other.pos;}
bool operator !=(const const_iterator & other) { assert(link == other.link); return pos != other.pos;}
bool operator <(const const_iterator & other) { assert(link == other.link); return pos < other.pos;}
bool operator >(const const_iterator & other) { assert(link == other.link); return pos > other.pos;}
bool operator <=(const const_iterator & other) { assert(link == other.link); return pos <= other.pos;}
bool operator >=(const const_iterator & other) { assert(link == other.link); return pos >= other.pos;}
bool operator ==(const const_iterator & other) const { assert(link == other.link); return pos == other.pos;}
bool operator !=(const const_iterator & other) const { assert(link == other.link); return pos != other.pos;}
bool operator <(const const_iterator & other) const { assert(link == other.link); return pos < other.pos;}
bool operator >(const const_iterator & other) const { assert(link == other.link); return pos > other.pos;}
bool operator <=(const const_iterator & other) const { assert(link == other.link); return pos <= other.pos;}
bool operator >=(const const_iterator & other) const { assert(link == other.link); return pos >= other.pos;}
};
void inner_resize(size_type new_size)
......
......@@ -163,8 +163,8 @@ namespace INMOST
iterator(Mesh * m, const cont_t::iterator & other ) : cont_t::iterator(other) , m_link(m){assert(m_link);}
iterator(const iterator & other ) : cont_t::iterator(other), m_link(other.m_link) {assert(m_link);}
ptrdiff_t operator -(const iterator & other) const {return static_cast<const cont_t::iterator>(*this)-static_cast<const cont_t::iterator>(other);}
iterator operator +(size_t n){iterator ret(*this); ret.cont_t::iterator::operator+(n); return ret;}
iterator operator -(size_t n){iterator ret(*this); ret.cont_t::iterator::operator-(n); return ret;}
iterator operator +(size_t n) const{return iterator(m_link,cont_t::iterator::operator +(n));}
iterator operator -(size_t n) const{return iterator(m_link,cont_t::iterator::operator -(n));}
iterator & operator ++() {cont_t::iterator::operator++(); return *this;}
iterator operator ++(int) {iterator ret(*this); cont_t::iterator::operator++(); return ret;}
iterator & operator --() {cont_t::iterator::operator--(); return *this;}
......@@ -180,8 +180,8 @@ namespace INMOST
reverse_iterator(Mesh * m, const cont_t::reverse_iterator & other ) : cont_t::reverse_iterator(other), m_link(m) {assert(m_link);}
reverse_iterator(const reverse_iterator & other ) : cont_t::reverse_iterator(other), m_link(other.m_link) {assert(m_link);}
ptrdiff_t operator -(const reverse_iterator & other) const {return static_cast<const cont_t::reverse_iterator>(other)-static_cast<const cont_t::reverse_iterator>(*this);}
reverse_iterator operator +(size_t n){reverse_iterator ret(*this); ret.cont_t::reverse_iterator::operator+(n); return ret;}
reverse_iterator operator -(size_t n){reverse_iterator ret(*this); ret.cont_t::reverse_iterator::operator-(n); return ret;}
reverse_iterator operator +(size_t n) const {return reverse_iterator(m_link,cont_t::reverse_iterator::operator+(n));}
reverse_iterator operator -(size_t n) const {return reverse_iterator(m_link,cont_t::reverse_iterator::operator-(n));}
reverse_iterator & operator ++() {cont_t::reverse_iterator::operator++(); return *this;}
reverse_iterator operator ++(int) {reverse_iterator ret(*this); cont_t::reverse_iterator::operator++(); return ret;}
reverse_iterator & operator --() {cont_t::reverse_iterator::operator--(); return *this;}
......
......@@ -17,9 +17,18 @@
// 7) add top and bottom interface
//todo
// 1. remove edges on pillars that do not get any associated block
// 1. (ok) remove edges on pillars that do not get any associated block
// 2. (ok) do not create cells with actnum = 0
// 2. fix fast intersect() algorithm, check against intersect_naive()
// 3. FixEdgeOrder cannot fix order of edges for certain interfaces - investigate
// 3. CheckEdgeOrder reports bad order of edges for certain interfaces - investigate
// 4. when encounter cells with one face with actnum = 1 should disconnect adjacent cells
// 5. (ok) populate keywords data, poro, perm, actnum, satnum, etc...
// 5.1 add saturation, pressure
// 6. local grid refinement
// 7. omp-parallel loading
// 8. mpi-parallel loading
// 9. do not convert arrays xyz and zcorn
//10. read wells into sets
//eclipse states
#define ECLSTRCMP(x,y) strncmp(x,y,8)
......@@ -42,6 +51,7 @@
#define ECL_COORDS 14
#define ECL_ZCORN 15
#define ECL_ACTNUM 16
#define ECL_SATNUM 17
#define ECL_GTYPE_NONE 0
#define ECL_GTYPE_TOPS 1
......@@ -86,18 +96,18 @@ namespace INMOST
Point(double _x, double _y) : x(_x), y(_y) {}
bool operator <(const Point & b) const
{
if (y < b.y - 1.0e-9) return true;
else if (y > b.y + 1.0e-9) return false;
else if (x < b.x - 1.0e-9) return true;
if (y < b.y - 1.0e-5) return true;
else if (y > b.y + 1.0e-5) return false;
else if (x < b.x - 1.0e-5) return true;
else return false;
}
bool operator ==(const Point & b) const
{
return fabs(y - b.y) < 1.0e-9 && fabs(x - b.x) < 1.0e-9;
return fabs(y - b.y) < 1.0e-5 && fabs(x - b.x) < 1.0e-5;
}
bool operator !=(const Point & b) const
{
return fabs(y - b.y) > 1.0e-9 || fabs(x - b.x) > 1.0e-9;
return fabs(y - b.y) > 1.0e-5 || fabs(x - b.x) > 1.0e-5;
}
};
//Comparator for map of nodes
......@@ -108,29 +118,41 @@ namespace INMOST
{
for(int k = 0; k < 3; ++k)
{
if( a[k] < b[k] - 1.0e-6)
if( a[k] < b[k] - 1.0e-5)
return true;
else if( a[k] > b[k] + 1.0e-6 )
else if( a[k] > b[k] + 1.0e-5 )
return false;
}
return false;
}
};
//Comparator for events in line sweep algorithm
class event_less
{
public:
bool operator()(const std::pair<double, int> & a, const std::pair<double, int> & b) const
{
if (a.first < b.first - 1.0e-6)
if (a.first < b.first - 1.0e-5)
return true;
else if (a.first > b.first + 1.0e-6)
else if (a.first > b.first + 1.0e-5)
return false;
else if (a.second < b.second)
return true;
return false;
}
};
//Comparator for depth of nodes in pillar
class pillar_less
{
public:
bool operator()(Storage::real a, Storage::real b) const
{
if( a < b - 1.0e-5)
return true;
return false;
}
};
class index_comparator
{
Storage::integer_array & data;
......@@ -180,7 +202,7 @@ namespace INMOST
__INLINE void normalize(Storage::real v[3])
__INLINE Storage::real normalize(Storage::real v[3])
{
Storage::real l = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
if( l )
......@@ -189,6 +211,7 @@ namespace INMOST
v[1] /= l;
v[2] /= l;
}
return l;
}
void compute_conormal(const Edge & a, const Edge & b, Storage::real vn[3])
......@@ -276,10 +299,12 @@ namespace INMOST
Point project(Storage::real pos[3], const Storage::real p0beg[3], const Storage::real p0end[3], const Storage::real p1beg[3], const Storage::real p1end[3])
{
Storage::real alpha, beta, c;
Storage::real alpha, beta, c,l, zmax, zmin;
Storage::real v[3], v0[3], v1[3];
alpha = (pos[2]-p0end[2])/(p0beg[2]-p0end[2]);
beta = (pos[2]-p1end[2])/(p1beg[2]-p1end[2]);
zmax = std::max(std::max(p0beg[2],p0end[2]),std::max(p1beg[2],p1end[2]));
zmin = std::min(std::min(p0beg[2],p0end[2]),std::min(p1beg[2],p1end[2]));
//get slice of pillar at z position of beginning of the segment
v0[0] = p0end[0] + alpha*(p0beg[0]-p0end[0]);
v0[1] = p0end[1] + alpha*(p0beg[1]-p0end[1]);
......@@ -289,10 +314,10 @@ namespace INMOST
v1[2] = pos[2];
//get vector connecting pillars
make_vec(v1,v0,v);
normalize(v);
l=normalize(v);
//project coordinates
c = dot_prod(pos,v);
return Point(pos[2],c);
c = ((pos[0]-v0[0])*v[0]+(pos[1]-v0[1])*v[1])/l;//dot_prod(pos,v);
return Point((pos[2]-zmin)/(zmax-zmin),c);
}
//intersect a pair of segments
......@@ -323,27 +348,27 @@ namespace INMOST
if (fabs(paend.x - pabeg.x) > 1.0e-9)
{
t1 = (pfind.x - pabeg.x) / (paend.x - pabeg.x);
if (t1 < 1.0e-9 || t1 > 1.0 - 1.0e-9) { if (print) std::cout << "out of bound: " << t1 << std::endl; return std::make_pair(false,InvalidNode()); }
if (t1 < 1.0e-4 || t1 > 1.0 - 1.0e-4) { if (print) std::cout << "out of bound: " << t1 << std::endl; return std::make_pair(false,InvalidNode()); }
}
if (fabs(paend.y - pabeg.y) > 1.0e-9)
{
t1 = (pfind.y - pabeg.y) / (paend.y - pabeg.y);
if (t1 < 1.0e-9 || t1 > 1.0 - 1.0e-9) { if (print) std::cout << "out of bound: " << t1 << std::endl; return std::make_pair(false,InvalidNode()); }
if (t1 < 1.0e-4 || t1 > 1.0 - 1.0e-4) { if (print) std::cout << "out of bound: " << t1 << std::endl; return std::make_pair(false,InvalidNode()); }
}
if (fabs(pbend.x - pbbeg.x) > 1.0e-9)
{
t2 = (pfind.x - pbbeg.x) / (pbend.x - pbbeg.x);
if (t2 < 1.0e-9 || t2 > 1.0 - 1.0e-9) { if (print) std::cout << "out of bound: " << t2 << std::endl; return std::make_pair(false,InvalidNode()); }
if (t2 < 1.0e-4 || t2 > 1.0 - 1.0e-4) { if (print) std::cout << "out of bound: " << t2 << std::endl; return std::make_pair(false,InvalidNode()); }
}
if (fabs(pbend.y - pbbeg.y) > 1.0e-9)
{
t2 = (pfind.y - pbbeg.y) / (pbend.y - pbbeg.y);
if (t2 < 1.0e-9 || t2 > 1.0 - 1.0e-9) { if (print) std::cout << "out of bound: " << t2 << std::endl; return std::make_pair(false,InvalidNode()); }
if (t2 < 1.0e-4 || t2 > 1.0 - 1.0e-4) { if (print) std::cout << "out of bound: " << t2 << std::endl; return std::make_pair(false,InvalidNode()); }
}
//restore third coordinate
for(int k = 0; k < 3; ++k)
find[k] = 0.5*((1-t1)*abeg[k]+t1*aend[k] + (1-t2)*bbeg[k]+t2*bend[k]);
if (print) std::cout << "intersection accepted (" << find[0] << "," << find[1] << "," << find[2] << ")" << std::endl;
if (print) std::cout << "intersection accepted (" << find[0] << "," << find[1] << "," << find[2] << ") t1 " << t1 << " t2 " << t2 << std::endl;
Node I;
std::map<position,Node,position_less>::iterator search = intersections.find(find);
//check whether intersection already exists
......@@ -703,33 +728,39 @@ namespace INMOST
void block_number_intersection(ElementArray<Element> & adj, Tag block, std::vector<int> & out)
{
Storage::integer_array be = adj[0]->IntegerArray(block);
std::vector<int> inter(be.begin(),be.end()), tmp(inter.size());
for(ElementArray<Edge>::size_type k = 1; k < adj.size(); ++k)
out.clear();
if( !adj.empty() )
{
be = adj[k]->IntegerArray(block);
tmp.resize(std::set_intersection(inter.begin(),inter.end(),be.begin(),be.end(),tmp.begin())-tmp.begin());
inter.swap(tmp);
Storage::integer_array be = adj[0]->IntegerArray(block);
std::vector<int> inter(be.begin(),be.end()), tmp(inter.size());
for(ElementArray<Edge>::size_type k = 1; k < adj.size(); ++k)
{
be = adj[k]->IntegerArray(block);
tmp.resize(std::set_intersection(inter.begin(),inter.end(),be.begin(),be.end(),tmp.begin())-tmp.begin());
inter.swap(tmp);
}
out.insert(out.end(),inter.begin(),inter.end());
}
out.clear();
out.insert(out.end(),inter.begin(),inter.end());
//Storage::integer_array bn = n->IntegerArray(block);
//bn.replace(bn.begin(),bn.end(),inter.begin(),inter.end());
}
void block_number_union(Element n, ElementArray<Element> & adj, Tag block)
{
Storage::integer_array be = adj[0]->IntegerArray(block);
std::vector<int> uni(be.begin(),be.end()), tmp(uni.size());
for(ElementArray<Edge>::size_type k = 1; k < adj.size(); ++k)
if( !adj.empty() )
{
be = adj[k]->IntegerArray(block);
tmp.resize(uni.size()+be.size());
tmp.resize(std::set_union(uni.begin(),uni.end(),be.begin(),be.end(),tmp.begin())-tmp.begin());
uni.swap(tmp);
Storage::integer_array be = adj[0]->IntegerArray(block);
std::vector<int> uni(be.begin(),be.end()), tmp(uni.size());
for(ElementArray<Edge>::size_type k = 1; k < adj.size(); ++k)
{
be = adj[k]->IntegerArray(block);
tmp.resize(uni.size()+be.size());
tmp.resize(std::set_union(uni.begin(),uni.end(),be.begin(),be.end(),tmp.begin())-tmp.begin());
uni.swap(tmp);
}
Storage::integer_array bn = n->IntegerArray(block);
bn.replace(bn.begin(),bn.end(),uni.begin(),uni.end());
}
Storage::integer_array bn = n->IntegerArray(block);
bn.replace(bn.begin(),bn.end(),uni.begin(),uni.end());
}
......@@ -763,7 +794,7 @@ namespace INMOST
Storage::integer dims[3], mapaxis[6] = {0,1,0,0,1,0};
Storage::real inrad = 0;
std::vector<Storage::real> xyz,perm,poro, tops,zcorn;
std::vector<Storage::integer> actnum;
std::vector<Storage::integer> actnum, satnum;
while(!fs.empty())
{
while(fgets(readline,2048,fs.back().first.first) != NULL)
......@@ -923,6 +954,16 @@ namespace INMOST
argtype = ECL_VAR_INT;
offset = state = ECL_ACTNUM;
}
else if( !ECLSTRCMP(p,"SATNUM") )