Commit 61de3bb7 authored by Kirill Terekhov's avatar Kirill Terekhov

Updates to GRDECL/VTU readers

Many updates to GRDECL file format:
- Processing more properties fields, such as pressure, sgas, soil, etc.
- Processing MULTIPLY keyword for properties and transmissibilities
- Introduced field TRANM for transmissibility multipliers set with TRANX/TRANY/TRANZ keywords in MULTIPLY keyword, or with EDITNNC keywords or FAULTS/MULTFLT keyword
- Algorithm that separates faces in place of blocks with zero volume and nonzero ACTNUM. Option ECL_DEGENERATE that replaces algorithm with transmissibility multiplier in TRANM.
-Optimization: bypass algorithm that resolves faces out of set of edges when only one face is expected.
- Processing of NTG keyword.

Fix swap algorithm in dynarray container.

Loading of 2d and mixed 2d/3d grids from VTU files.

Fix Visual Studio warnings in sparse.cpp
parent 1baa532d
......@@ -1018,7 +1018,7 @@ void draw()
glEnable(GL_BLEND);
for(int q = 0; q < inside_face.size() ; q++) inside_face[q].draw();
for (int q = 0; q < inside_face.size(); q++) inside_face[q].draw();
glDisable(GL_BLEND);
......@@ -1271,7 +1271,7 @@ void draw()
{
Element it = mesh->ElementByLocalID(etype,i);
if( !it.isValid() ) continue;
if( it->Boundary() )
if( it->Boundary() || (etype == CELL && Element::GetGeometricDimension(it->GetGeometricType()) == 2))
{
int mmat = color ? (!it->getAsFace()->BackCell().isValid() ? 0 : it->getAsFace()->BackCell()->Integer(colort)) : (!it->getAsFace()->BackCell().isValid() ? 0 : it->getAsFace()->BackCell()->Integer(mat));
if( matfilter == 0 || matfilter-1 == mmat )
......@@ -1285,7 +1285,7 @@ void draw()
{
Element it = mesh->ElementByLocalID(etype,i);
if( !it.isValid() ) continue;
if( it->Boundary() )
if (it->Boundary() || (etype == CELL && Element::GetGeometricDimension(it->GetGeometricType()) == 2))
polygons.push_back(DrawFace(it->self(),0,campos));
}
}
......
......@@ -1812,7 +1812,9 @@ namespace INMOST
void swap(dynarray<element,stacked> & other)
{
size_type k = size(), n = other.size();
if( k <= static_cast<size_type>(stacked) && n <= static_cast<size_type>(stacked) )
bool konstack = (pbegin == stack);
bool nonstack = (other.pbegin == other.stack);
if( konstack && nonstack )
{
char temp[stacked*sizeof(element)];
memcpy(temp,(void *)stack,sizeof(element)*k);
......@@ -1821,7 +1823,7 @@ namespace INMOST
other.pend = other.pbegin+k;
pend = pbegin+n;
}
else if( k <= static_cast<size_type>(stacked) && n > static_cast<size_type>(stacked) )
else if( konstack && !nonstack )
{
memcpy((void *)other.stack,(void *)stack,sizeof(element)*k);
pbegin = other.pbegin;
......@@ -1831,7 +1833,7 @@ namespace INMOST
other.pend = other.stack+k;
other.preserved = other.stack+static_cast<size_type>(stacked);
}
else if( k > static_cast<size_type>(stacked) && n <= static_cast<size_type>(stacked) )
else if( !konstack && nonstack )
{
memcpy((void *)stack,(void *)other.stack,sizeof(element)*n);
other.pbegin = pbegin;
......
......@@ -3015,7 +3015,15 @@ namespace INMOST
/// right now this is "TRUE" by default since results are much better.
/// - "ECL_PROJECT_PERM" - Set "TRUE" to project permeability tensor from grid block coordinates
/// into global coordinates. Otherwise the tensor is considered to be
/// defined on global coordinates.
/// defined on global coordinates. Default "FALSE".
/// - "ECL_DEGENERATE" - In GRDECL format some active grid block may have zero volume, which
/// means there is a fault. Set "TRUE" to introduce a gap between blocks
/// that share degenerate active grid block, set to "TRANM" to introduce
/// zero transmissibility multiplier and keep the grid connected and "FALSE"
/// to simply keep the grid connected. Default: "TRUE".
/// - "ECL_TOPOLOGY" - If "TRUE" checks topology of the grid for errors, this may provide useful
/// warnings if layers of the mesh enter each other and the grid cannot be
/// considered conformal. Default: "FALSE".
///
/// \todo
/// introduce "SET_TAGS_LOAD", "SET_TAGS_SAVE" to explicitly provide set of tags to write
......
This diff is collapsed.
......@@ -1311,12 +1311,13 @@ safe_output:
filled = fscanf(f,"\n");
}
if( grid_is_2d == 2 )
if( grid_is_2d == 2 && ncells )
{
for(i = 0; i < ncells && grid_is_2d == 2; i++)
if( ct[i] > 9 ) grid_is_2d = 0;
//this will create all 2d elements as cells
//grid_is_2d = 1;
if (grid_is_2d == 2) grid_is_2d = 1;
}
if( verbosity > 0 )
......
......@@ -162,6 +162,33 @@ namespace INMOST
newpolyh.push_back(CreateCell(hfaces).first.GetHandle());
}
}
//check grid type
if( grid_is_2d == 2 && ncells) //detect grid type
{
std::stringstream type(v->GetChild("Cells")->GetChildWithAttrib("Name", "types")->GetContents());
int ctype;
bool have_2d = false;
for (int q = 0; q < ncells && grid_is_2d == 2; ++q)
{
type >> ctype;
if( ctype > 9 ) grid_is_2d = 0;
}
if( grid_is_2d == 2 ) grid_is_2d = 1;
}
if (verbosity > 0)
{
switch (grid_is_2d)
{
case 0: std::cout << "Grid has three dimensions" << std::endl; break;
case 1: std::cout << "Grid has two dimensions" << std::endl; break;
case 2: std::cout << "Grid has undetermined dimension" << std::endl; break;
}
}
bool have_faces = false; //some elements go as faces
bool have_edges = false;
bool have_nodes = false;
//read all the cells
{
std::stringstream conn(v->GetChild("Cells")->GetChildWithAttrib("Name", "connectivity")->GetContents());
......@@ -184,7 +211,132 @@ namespace INMOST
RemMarker(newnodes[cconn], unused_marker);
totread++;
}
if (ctype == 10)
if (ctype == 1) //VTK_VERTEX
{
newcells[q] = hnodes.at(0);
have_nodes = true;
}
else if (ctype == 2) //VTK_POLY_VERTEX
{
std::cout << __FILE__ << ":" << __LINE__ << " skipping VTK_POLY_VERTEX" << std::endl;
}
else if (ctype == 3) //VTK_LINE
{
newcells[q] = CreateEdge(hnodes).first.GetHandle();
have_edges = true;
}
else if (ctype == 4)
{
std::cout << __FILE__ << ":" << __LINE__ << " skipping VTK_POLY_LINE" << std::endl;
}
else if (ctype == 5) //VTK_TRIANGLE
{
if (grid_is_2d == 1)
{
ElementArray<Node> e_nodes(this, 1);
ElementArray<Edge> f_edges(this, 2);
ElementArray<Face> c_faces(this);
for (unsigned int k = 0; k < 3; k++)
{
e_nodes.at(0) = hnodes.at(k);
f_edges.at(0) = CreateEdge(e_nodes).first->GetHandle();
e_nodes.at(0) = hnodes.at((k + 1) % 3);
f_edges.at(1) = CreateEdge(e_nodes).first->GetHandle();
c_faces.push_back(CreateFace(f_edges).first);
}
Cell c = CreateCell(c_faces,hnodes).first;
newcells[q] = c->GetHandle();
}
else
{
newcells[q] = CreateFace(hnodes).first->GetHandle();
have_faces = true;
}
break;
}
else if (ctype == 6) //VTK_TRIANGLE_STRIP
{
std::cout << __FILE__ << ":" << __LINE__ << " skipping VTK_TRIANGLE_STRIP" << std::endl;
}
else if (ctype == 7) //VTK_POLYGON
{
if (grid_is_2d == 1)
{
ElementArray<Node> e_nodes(this, 1);
ElementArray<Edge> f_edges(this, 2);
ElementArray<Face> c_faces(this);
for (ElementArray<Node>::size_type k = 0; k < hnodes.size(); k++)
{
e_nodes.at(0) = hnodes.at(k);
f_edges.at(0) = CreateEdge(e_nodes).first->GetHandle();
e_nodes.at(0) = hnodes.at((k + 1) % hnodes.size());
f_edges.at(1) = CreateEdge(e_nodes).first->GetHandle();
c_faces.push_back(CreateFace(f_edges).first);
}
Cell c = CreateCell(c_faces, hnodes).first;
newcells[q] = c->GetHandle();
}
else
{
newcells[q] = CreateFace(hnodes).first->GetHandle();
have_faces = true;
}
break;
}
else if (ctype == 8) //VTK_PIXEL
{
HandleType temp = hnodes.at(2);
hnodes.at(2) = hnodes.at(3);
hnodes.at(3) = temp;
if (grid_is_2d == 1)
{
ElementArray<Node> e_nodes(this, 1);
ElementArray<Edge> f_edges(this, 2);
ElementArray<Face> c_faces(this);
e_nodes.resize(1);
f_edges.resize(2);
for (int k = 0; k < 4; k++)
{
e_nodes.at(0) = hnodes.at(k);
f_edges.at(0) = CreateEdge(e_nodes).first->GetHandle();
e_nodes.at(0) = hnodes.at((k + 1) % 4);
f_edges.at(1) = CreateEdge(e_nodes).first->GetHandle();
c_faces.push_back(CreateFace(f_edges).first);
}
Cell c = CreateCell(c_faces, hnodes).first;
newcells[q] = c->GetHandle();
}
else
{
newcells[q] = CreateFace(hnodes).first->GetHandle();
have_faces = true;
}
}
else if (ctype == 9) //VTK_QUAD
{
if (grid_is_2d == 1)
{
ElementArray<Node> e_nodes(this,1);
ElementArray<Edge> f_edges(this,2);
ElementArray<Face> c_faces(this);
for (int k = 0; k < 4; k++)
{
e_nodes.at(0) = hnodes.at(k);
f_edges.at(0) = CreateEdge(e_nodes).first->GetHandle();
e_nodes.at(0) = hnodes.at((k + 1) % 4);
f_edges.at(1) = CreateEdge(e_nodes).first->GetHandle();
c_faces.push_back(CreateFace(f_edges).first);
}
Cell c = CreateCell(c_faces, hnodes).first;
newcells[q] = c->GetHandle();
}
else
{
newcells[q] = CreateFace(hnodes).first->GetHandle();
have_faces = true;
}
}
else if (ctype == 10) //VTK_TETRA
{
assert(nread == 4);
const integer nodesnum[12] = { 0, 2, 1, 0, 1, 3, 1, 2, 3, 0, 3, 2 };
......@@ -238,6 +390,9 @@ namespace INMOST
};
//5, 6, 7, 8, 9
//0, 1, 2, 3, 4
//9, 8, 7, 6, 5
//0, 1, 2, 3, 4
const integer sizes[7] = { 4, 4, 4, 4, 4, 5, 5 };
newcells[q] = CreateCell(hnodes, nodesnum, sizes, 7).first.GetHandle();
}
......@@ -269,6 +424,22 @@ namespace INMOST
{
std::string dname[2] = { "PointData", "CellData" };
ElementType dtype[2] = { NODE, CELL };
ElementType dsparse[2] = { NONE, NONE };
if (have_faces)
{
dtype[1] |= FACE;
dsparse[1] |= FACE;
}
if (have_edges)
{
dtype[1] |= EDGE;
dsparse[1] |= EDGE;
}
if (have_nodes)
{
dtype[1] |= NODE;
dsparse[1] |= NODE;
}
HandleType * darray[2] = { &newnodes[0], &newcells[0] };
int dsize[2] = { nnodes, ncells };
for (int j = 0; j < 2; ++j)
......@@ -284,7 +455,7 @@ namespace INMOST
int ncomps = 1;
int nca = pd->FindAttrib("NumberOfComponents");
if (nca != pd->NumAttrib()) ncomps = atoi(pd->GetAttrib(nca).value.c_str());
TagRealArray t = CreateTag(pd->GetAttrib("Name"), DATA_REAL, dtype[j], NONE, ncomps);
TagRealArray t = CreateTag(pd->GetAttrib("Name"), DATA_REAL, dtype[j], dsparse[j], ncomps);
std::stringstream inp(pd->GetContents());
for (int l = 0; l < dsize[j]; ++l)
{
......
......@@ -107,7 +107,8 @@ namespace INMOST
}
RowMerger::RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, const std::vector<INMOST_DATA_ENUM_TYPE> & Pre, const std::vector<INMOST_DATA_ENUM_TYPE> & Post, bool Sorted)
: Sorted(Sorted), Nonzeros(0), IntervalBeg(interval_begin), IntervalEnd(interval_end), NonlocalPre(Pre), NonlocalPost(Post), LinkedList(interval_begin-Pre.size(),interval_end+Post.size()+1,Row::make_entry(UNDEF,0.0))
: Sorted(Sorted), Nonzeros(0), IntervalBeg(interval_begin), IntervalEnd(interval_end), NonlocalPre(Pre), NonlocalPost(Post),
LinkedList(static_cast<INMOST_DATA_ENUM_TYPE>(interval_begin - Pre.size()), static_cast<INMOST_DATA_ENUM_TYPE>(interval_end + Post.size() + 1), Row::make_entry(UNDEF, 0.0))
{
//assert(std::is_sorted(Pre.begin(),Pre.end()));
//assert(std::is_sorted(Post.begin(),Post.end()));
......@@ -116,8 +117,8 @@ namespace INMOST
void RowMerger::Resize(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool _Sorted)
{
LinkedList.set_interval_beg(interval_begin-NonlocalPre.size());
LinkedList.set_interval_end(interval_end+1+NonlocalPost.size());
LinkedList.set_interval_beg(static_cast<INMOST_DATA_ENUM_TYPE>(interval_begin - NonlocalPre.size()));
LinkedList.set_interval_end(static_cast<INMOST_DATA_ENUM_TYPE>(interval_end + 1 + NonlocalPost.size()));
IntervalBeg = interval_begin;
IntervalEnd = interval_end;
std::fill(LinkedList.begin(),LinkedList.end(),Row::make_entry(UNDEF,0.0));
......@@ -182,7 +183,7 @@ namespace INMOST
assert(!NonlocalPre.empty()); //there are indices provided
std::vector< INMOST_DATA_ENUM_TYPE >::const_iterator search = std::lower_bound(NonlocalPre.begin(),NonlocalPre.end(),pos);
assert(*search == pos); //is there such index?
return IntervalBeg - NonlocalPre.size() + static_cast<INMOST_DATA_ENUM_TYPE>(search-NonlocalPre.begin());
return static_cast<INMOST_DATA_ENUM_TYPE>(IntervalBeg - NonlocalPre.size() + static_cast<INMOST_DATA_ENUM_TYPE>(search - NonlocalPre.begin()));
}
else if( pos >= IntervalEnd )
{
......@@ -680,7 +681,7 @@ namespace INMOST
MPI_Comm_size(GetCommunicator(),&size);
}
#endif
int * ord = NULL;
INMOST_DATA_ENUM_TYPE * ord = NULL;
if (file_ord != "")
{
std::ifstream input_ord;
......@@ -688,7 +689,7 @@ namespace INMOST
if( input_ord.fail() ) throw -2;
int n;
input_ord >> n;
ord = (int *) malloc(sizeof(int) * n);
ord = (INMOST_DATA_ENUM_TYPE *)malloc(sizeof(INMOST_DATA_ENUM_TYPE)* n);
assert(ord != NULL);
for (int i=0; i<n; i++) input_ord >> ord[i];
int nbl;
......
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