Commit 27edd03c authored by Kirill Terekhov's avatar Kirill Terekhov

Feature and fixes

Added option “ECL_SPLIT_GLUED” to triangulate faces for blocks that
degenerate on three pillars.

Corrected that Face::SplitFace will not connect new nodes to existing
cells.

PMF file format now restores tag on mesh that indicates topology error.
parent 7d7e2fbb
......@@ -3956,6 +3956,18 @@ void draw_screen()
bclipper->clip_plane(p,n);
bclipper->gen_clip(clip_boundary,n);
clipboxupdate = false;
for(int k = 0; k < (int)orphans.size(); ++k)
{
if( orphans[k]->GetElementType() == FACE )
clip_boundary.push_back(DrawFace(orphans[k]->getAsFace()));
else if( orphans[k]->GetElementType() == CELL )
{
ElementArray<Face> faces = orphans[k].getFaces();
for(int q = 0; q < faces.size(); ++q)
clip_boundary.push_back(DrawFace(faces[q]));
}
}
if( current_picker != NULL ) {delete current_picker; current_picker = NULL;}
if( !clip_boundary.empty() ) current_picker = new picker(clip_boundary);
......@@ -4107,13 +4119,15 @@ void draw_screen()
}
glEnd();
if( disp_e.isValid() )
DrawElement(disp_e,color_t(1,1,0),color_t(1,0,0),color_t(0,0,1));
if( draw_orphan )
for(int k = 0; k < (int)orphans.size(); ++k)
DrawElement(orphans[k],color_t(1,0,0),color_t(0,1,0),color_t(0,0,1));
DrawElement(orphans[k],color_t(1,0,1),color_t(0,1,1),color_t(0,0,1));
if( disp_e.isValid() )
DrawElement(disp_e,color_t(1,1,0),color_t(1,0,0),color_t(0,0,1));
double top = 0.96;
if( picked != -1 )
{
......@@ -4694,17 +4708,33 @@ int main(int argc, char ** argv)
}
{
std::map<ElementType,int> num_orphans;
std::map<ElementType,int> num_orphans, num_topo;
/*
for(Mesh::iteratorElement it = mesh->BeginElement(FACE|EDGE|NODE); it != mesh->EndElement(); ++it)
if( it->nbAdjElements(CELL) == 0 )
{
orphans.push_back(it->self());
num_orphans[it->GetElementType()]++;
}
*/
printf("number of orphan elements: %d\n",orphans.size());
for(std::map<ElementType,int>::iterator it = num_orphans.begin(); it != num_orphans.end(); ++it)
printf("%s %d\n",ElementTypeName(it->first),it->second);
int was = orphans.size();
if(mesh->TopologyErrorTag().isValid())
{
for(Mesh::iteratorElement it = mesh->BeginElement(FACE|EDGE|NODE); it != mesh->EndElement(); ++it)
if( it->HaveData(mesh->TopologyErrorTag()) )
{
orphans.push_back(it->self());
num_topo[it->GetElementType()]++;
}
printf("number of elements with topology error: %d\n",orphans.size()-was);
for(std::map<ElementType,int>::iterator it = num_topo.begin(); it != num_topo.end(); ++it)
printf("%s %d\n",ElementTypeName(it->first),it->second);
for(int k = was; k < orphans.size(); ++k)
std::cout << ElementTypeName(orphans[k]->GetElementType()) << ":" << orphans[k]->LocalID() << std::endl;
}
}
quatinit();
......
......@@ -2964,6 +2964,7 @@ namespace INMOST
/// Current availible file options:
/// - "VTK_GRID_DIMS" - set "2" for two-dimensional vtk grids, "3" for three-dimensional vtk grids
/// - "VERBOSITY" - set "2" for progress messages, "1" for reports, "0" for silence
/// - "ECL_SPLIT_GLUED" - set "TRUE" to triangulate faces of the blocks that degenerate on three pillars.
///
/// \todo
/// introduce "SET_TAGS_LOAD", "SET_TAGS_SAVE" to explicitly provide set of tags to write
......
......@@ -518,6 +518,17 @@ namespace INMOST
void Mesh::LoadECL(std::string File)
{
std::cout << std::scientific;
bool perform_splitting = false;
for(INMOST_DATA_ENUM_TYPE k = 0; k < file_options.size(); ++k)
{
if( file_options[k].first == "ECL_SPLIT_GLUED" )
{
if( file_options[k].second == "TRUE" )
perform_splitting = true;
else
perform_splitting = false;
}
}
FILE * f = fopen(File.c_str(),"r");
if( f == NULL )
{
......@@ -1884,18 +1895,20 @@ ecl_exit_loop:
//DeleteTag(block_number, NODE);
//now construct top and bottom interfaces
printf("started tops/bottoms/cells\n");
//Tag split_face = CreateTag("SPLIT_FACE",DATA_REFERENCE,FACE,NONE,1); //points to edge that should be used to split face
Tag split_face;
if( perform_splitting )
split_face = CreateTag("SPLIT_FACE",DATA_REFERENCE,FACE,NONE,1); //points to edge that should be used to split face
#if defined(USE_OMP)
#pragma omp parallel
#endif
{
ElementArray<Node> edge_nodes(this,2);
for(int q = 0; q < 2; ++q) //uneven cell construction for OpenMP
for(int uneven = 0; uneven < 2; ++uneven) //uneven cell construction for OpenMP
{
#if defined(USE_OMP)
#pragma omp for
#endif
for(int i = q; i < dims[0]; i+=2)
for(int i = uneven; i < dims[0]; i+=2)
{
//printf("top/bottom/cells %6.2f%%\r", ((Storage::real)i)/((Storage::real)dims[0]-1)*100);
//fflush(stdout);
......@@ -1952,54 +1965,90 @@ ecl_exit_loop:
{
Face f = CreateFace(face_edges).first;
if( TopologyErrorTag().isValid() && f->HaveData(TopologyErrorTag()) ) std::cout << __FILE__ <<":"<<__LINE__ << " topology error on face " << f->GetHandle() << std::endl;
/*
int split = 0;
if( block_nodes[cur*8+0] == block_nodes[cur*8+4] &&
block_nodes[cur*8+3] == block_nodes[cur*8+7] )
split = 1; //SE-NW diagonal
else
if( block_nodes[cur*8+1] == block_nodes[cur*8+5] &&
block_nodes[cur*8+2] == block_nodes[cur*8+6] )
split = 2; //SW-NE diagonal
HandleType diag = f->Reference(split_face);
//split up quad faces into triagnles if they degenerate
if( split )
{
if( diag == InvalidHandle() )
{
if( split == 1 ) //SE-NW diagonal
{
edge_nodes.at(0) = block_nodes[cur*8+1+q*4];
edge_nodes.at(1) = block_nodes[cur*8+2+q*4];
f->Reference(split_face) = CreateEdge(edge_nodes).first.GetHandle();
}
else //SW-NE diagonal
{
edge_nodes.at(0) = block_nodes[cur*8+0+q*4];
edge_nodes.at(1) = block_nodes[cur*8+3+q*4];
f->Reference(split_face) = CreateEdge(edge_nodes).first.GetHandle();
}
}
else //the face was already split
{
int was_split = 1; //was SE-NW
Edge e(this,diag);
if( e->getBeg()->GetHandle() == block_nodes[cur*8+0+q*4] &&
e->getEnd()->GetHandle() == block_nodes[cur*8+3+q*4])
was_split = 2;
if( split != was_split )
if( perform_splitting )
{
std::cout << "Got different split direction" << std::endl;
std::cout << "was " << e->getBeg()->GetHandle() << " <-> " << e->getEnd()->GetHandle() << (was_split == 1 ? "SE-NW" : "SW-NE") << std::endl;
std::cout << (q? "top" : "bottom") << " nodes: " << std::endl;
for(int l = 0; l < 4; ++l)
std::cout << block_nodes[cur*8+l+q*4];
std::cout << "want direction " << (split == 1? "SE-NW" : "SW-NE") << std::endl;
}
}
int split = 0;
bool senw = block_nodes[cur*8+0] == block_nodes[cur*8+4] && block_nodes[cur*8+3] == block_nodes[cur*8+7];
bool swne = block_nodes[cur*8+1] == block_nodes[cur*8+5] && block_nodes[cur*8+2] == block_nodes[cur*8+6];
if( senw && swne )
split = 0;
else if( senw )
split = 1; //SE-NW diagonal
else if( swne )
split = 2; //SW-NE diagonal
HandleType diag = f->Reference(split_face);
//split up quad faces into triagnles if they degenerate
if( split )
{
if( diag == InvalidHandle() )
{
if( split == 1 ) //SE-NW diagonal
{
edge_nodes.at(0) = block_nodes[cur*8+1+q*4];
edge_nodes.at(1) = block_nodes[cur*8+2+q*4];
f->Reference(split_face) = CreateEdge(edge_nodes).first.GetHandle();
}
else if( split == 2 ) //SW-NE diagonal
{
edge_nodes.at(0) = block_nodes[cur*8+0+q*4];
edge_nodes.at(1) = block_nodes[cur*8+3+q*4];
f->Reference(split_face) = CreateEdge(edge_nodes).first.GetHandle();
}
}
else //the face was already split
{
int was_split = 0;
Edge e(this,diag);
if( e->getBeg()->GetHandle() == block_nodes[cur*8+1+q*4] &&
e->getEnd()->GetHandle() == block_nodes[cur*8+2+q*4])
was_split = 1; //was SE-NW
else if( e->getBeg()->GetHandle() == block_nodes[cur*8+0+q*4] &&
e->getEnd()->GetHandle() == block_nodes[cur*8+3+q*4])
was_split = 2; //was SW-NE
if( !was_split )
{
std::cout << "Cannot detect how the face was priviously split" << std::endl;
std::cout << "edge: " << e->getBeg()->GetHandle() << " <-> " << e->getEnd()->GetHandle() << std::endl;
std::cout << "SE-NW: " << block_nodes[cur*8+0+q*4] << " <-> " << block_nodes[cur*8+3+q*4] << std::endl;
std::cout << "SW-NE: " << block_nodes[cur*8+1+q*4] << " <-> " << block_nodes[cur*8+2+q*4] << std::endl;
}
if( split != was_split )
{
//replace with node
e->Delete();
Storage::real xyz[3] = {0,0,0};
for(int l = 0; l < 4; ++l)
{
Storage::real_array c = Node(this,block_nodes[cur*8+l+q*4]).Coords();
xyz[0] += c[0]*0.25;
xyz[1] += c[1]*0.25;
xyz[2] += c[2]*0.25;
}
f->Reference(split_face) = CreateNode(xyz)->GetHandle();
/*
std::cout << "Got different split direction" << std::endl;
std::cout << "was " << e->getBeg()->GetHandle() << " <-> " << e->getEnd()->GetHandle() << " " << (was_split == 1 ? "SE-NW" : "SW-NE") << std::endl;
std::cout << (q? "top" : "bottom") << " nodes: " << std::endl;
for(int l = 0; l < 4; ++l)
std::cout << l+q*4 << ":" << block_nodes[cur*8+l+q*4] << " ";
std::cout << std::endl;
std::cout << (1-q? "top" : "bottom") << " nodes: " << std::endl;
for(int l = 0; l < 4; ++l)
std::cout << l+(1-q)*4 << ":" << block_nodes[cur*8+l+(1-q)*4] << " ";
std::cout << std::endl;
std::cout << "SE-NW: " << block_nodes[cur*8+0+q*4] << " <-> " << block_nodes[cur*8+3+q*4] << std::endl;
std::cout << "SW-NE: " << block_nodes[cur*8+1+q*4] << " <-> " << block_nodes[cur*8+2+q*4] << std::endl;
std::cout << "want direction " << (split == 1? "SE-NW" : "SW-NE") << std::endl;
*/
}
}
}
}
*/
f->Integer(face_origin) = 2;
......@@ -2102,23 +2151,40 @@ ecl_exit_loop:
} //q
}
printf("finished tops/bottoms/cells\n");
/*
printf("started splitting faces of degenerate cells\n");
if( perform_splitting )
{
printf("started splitting faces of degenerate cells\n");
#if defined(USE_OMP)
#pragma omp parallel for
#endif
for(integer it = 0; it < FaceLastLocalID(); ++it) if( isValidFace(it) )
{
Face f = FaceByLocalID(it);
HandleType diag = f->Reference(split_face);
if( diag != InvalidHandle() )
for(integer it = 0; it < FaceLastLocalID(); ++it) if( isValidFace(it) )
{
ElementArray<Edge> split_edge(this,1,diag);
Face::SplitFace(f,split_edge,0);
Face f = FaceByLocalID(it);
HandleType diag = f->Reference(split_face);
if( diag != InvalidHandle() )
{
if( GetHandleElementType(diag) == NODE )
{
ElementArray<Edge> split_edges(this);
ElementArray<Node> edge_nodes(this,2);
ElementArray<Node> face_nodes = f->getNodes();
edge_nodes.at(0) = diag;
for(int k = 0; k < (int)face_nodes.size(); ++k)
{
edge_nodes[1] = face_nodes[k];
split_edges.push_back(CreateEdge(edge_nodes).first);
}
Face::SplitFace(f,split_edges,0);
}
else
{
if( diag != InvalidHandle() )
Face::SplitFace(f,ElementArray<Edge>(this,1,diag),0);
}
}
}
printf("finished splitting faces of degenerate cells\n");
}
printf("finished splitting faces of degenerate cells\n");
*/
//printf("\n");
//cleanup data
//DeleteTag(split_face);
......
......@@ -1413,6 +1413,8 @@ namespace INMOST
for(ElementType etype = NODE; etype <= CELL; etype = etype << 1)
if( tag_global_id.isDefined(etype) ) have_global_id |= etype;
}
if( HaveTag("TOPOLOGY_ERROR_TAG") )
tag_topologyerror = GetTag("TOPOLOGY_ERROR_TAG");
#if defined(USE_MPI)
if( m_state == Mesh::Parallel )
{
......
......@@ -1119,15 +1119,25 @@ namespace INMOST
for(dynarray<HandleType,2>::size_type it = 0; it < cells.size(); ++it)
{
adj_type & hc = m->HighConn(cells[it]);
adj_type & hc = m->HighConn(cells[it]); //cell nodes
hc.clear(); //have to recompute cell nodes
adj_type & lc = m->LowConn(cells[it]);
adj_type & lc = m->LowConn(cells[it]); //cell faces
lc.insert(lc.end(),ret.begin(),ret.end());
m->ComputeGeometricType(cells[it]);
ElementArray<Node> nn(m);
m->RestoreCellNodes(cells[it],nn);
hc.insert(hc.begin(),nn.begin(),nn.end());
m->RecomputeGeometricData(cells[it]);
//have to add cell to nodes that do not yet have connection to the cell
for(ElementArray<Node>::iterator jt = nn.begin(); jt != nn.end(); ++jt)
{
adj_type & nlc = m->LowConn(*jt); //node cells
bool have_cell = false;
for(adj_type::iterator kt = nlc.begin(); kt != nlc.end() && !have_cell; ++kt)
if( *kt == cells[it] ) have_cell = true;
if( !have_cell )
nlc.push_back(cells[it]);
}
}
return ret;
......
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