Commit 04aaeaca authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Sync changes

Restored OctreeCutcell example;
Added FixEdgeOrder, CheckEdgeOrder algorithms to check of consistency of
2d elements;
Element::Connect now uses FixEdgeOrder for faces;
Automatic grid dimension detection for vtk;
Fixed bug that sets may contain invalid links to elements even after
ApplyModification;
Corrected condition estimator in ddpqiluc;
Switched off pseudoinverse in bicgstab(l) since svd algorithm is bad.
parent da9fed1e
......@@ -8,12 +8,281 @@ namespace INMOST
bool Cell::CheckEdgeOrder() const
{
assert(GetHandleElementType(GetHandle())==CELL);
Mesh * m = GetMeshLink();
if( Element::GetGeometricDimension(m->GetGeometricType(GetHandle())) == 2 )
{
if( !m->HideMarker() )
{
HandleType last, first;
adj_type const & lc = m->LowConn(GetHandle());
if( lc.size() < 3 ) return false;
HandleType q = lc[0]; //edge 0
adj_type const & qlc = m->LowConn(q);
first = qlc[0]; //node 0
last = qlc[1]; //node 1
HandleType r = lc[1]; //edge 1
adj_type & rlc = m->LowConn(r);
if( first == rlc[0] || first == rlc[1] )
{
HandleType temp = first;
first = last;
last = temp;
}
else if ( !(last == rlc[0] || last == rlc[1]) )
return false;
adj_type::size_type it = 1, iend = lc.size()-1;
while(it < iend) //loop over edges
{
adj_type const & ilc = m->LowConn(lc[it]);
if( last == ilc[0] ) last = ilc[1];
else if( last == ilc[1] ) last = ilc[0];
else return false;
++it;
}
}
else
{
MarkerType hm = m->HideMarker();
HandleType first, last;
enumerator i = ENUMUNDEF, k = ENUMUNDEF, k1 = ENUMUNDEF, k2;
adj_type const & lc = m->LowConn(GetHandle());
if( m->Count(lc.data(),lc.size(),hm) < 3 ) return false;
i = m->getNext(lc.data(),static_cast<enumerator>(lc.size()),i,hm);
HandleType q = lc[i]; //edge 0
adj_type const & qlc = m->LowConn(q);
k = m->getNext(qlc.data(),static_cast<enumerator>(qlc.size()),k,hm);
first = qlc[k]; //node 0
k = m->getNext(qlc.data(),static_cast<enumerator>(qlc.size()),k,hm);
last = qlc[k]; //node 1
i = m->getNext(lc.data(),static_cast<enumerator>(lc.size()),i,hm);
HandleType r = lc[i]; //edge 1
adj_type const & rlc = m->LowConn(r);
k1 = m->getNext(rlc.data(),static_cast<enumerator>(rlc.size()),k1,hm);
k2 = m->getNext(rlc.data(),static_cast<enumerator>(rlc.size()),k1,hm);
if( first == rlc[k1] || first == rlc[k2] )
{
HandleType temp = first;
first = last;
last = temp;
}
else if ( !(last == rlc[k1] || last == rlc[k2]) )
return false;
adj_type::size_type it = 1, iend = lc.size()-1;
while(it != iend) if( !m->GetMarker(lc[it],hm) ) //loop over edges
{
adj_type const & ilc = m->LowConn(lc[it]);
k1 = ENUMUNDEF;
k1 = m->getNext(ilc.data(),static_cast<enumerator>(ilc.size()),k1,hm);
k2 = m->getNext(ilc.data(),static_cast<enumerator>(ilc.size()),k1,hm);
if( last == ilc[k1] ) last = ilc[k2];
else if( last == ilc[k2] ) last = ilc[k1];
else return false;
++it;
}
}
}
return true;
}
bool Cell::FixEdgeOrder() const
{
assert(GetHandleElementType(GetHandle())==CELL);
Mesh * m = GetMeshLink();
if( Element::GetGeometricDimension(m->GetGeometricType(GetHandle())) == 2 )
{
if( !m->HideMarker() )
{
HandleType last, first;
adj_type & lc = m->LowConn(GetHandle());
if( lc.size() < 3 ) return false;
HandleType q = lc[0]; //edge 0
adj_type const & qlc = m->LowConn(q);
if( qlc.size() != 2 ) return false;
first = qlc[0]; //node 0
last = qlc[1]; //node 1
HandleType r = lc[1]; //edge 1
adj_type & rlc = m->LowConn(r);
if( rlc.size() != 2 ) return false;
if( first == rlc[0] || first == rlc[1] )
{
HandleType temp = first;
first = last;
last = temp;
}
else if ( !(last == rlc[0] || last == rlc[1]) )
{
adj_type::size_type jt = 2, jend = lc.size();
while(jt < jend)
{
adj_type const & ilc = m->LowConn(lc[jt]);
if( ilc.size() != 2 ) return false;
if( first == ilc[0] || first == ilc[1] )
{
HandleType temp = lc[1];
lc[1] = lc[jt];
lc[jt] = lc[1];
temp = first;
first = last;
last = temp;
break;
}
else if( last == ilc[0] || last == ilc[1] )
{
HandleType temp = lc[1];
lc[1] = lc[jt];
lc[jt] = lc[1];
break;
}
++jt;
}
if( jt == jend ) return false; //no matching edge
}
adj_type::size_type it = 1, iend = lc.size()-1;
while(it < iend) //loop over edges
{
adj_type const & ilc = m->LowConn(lc[it]);
if( ilc.size() != 2 ) return false;
if( last == ilc[0] )
{
last = ilc[1];
++it;
}
else if( last == ilc[1] )
{
last = ilc[0];
++it;
}
else //search for the connected edge and swap with current
{
adj_type::size_type jt = it+1, jend = lc.size();
while(jt < jend)
{
adj_type const & ilc = m->LowConn(lc[jt]);
if( ilc.size() != 2 ) return false;
if( last == ilc[0] || last == ilc[1] )
{
HandleType temp = lc[it];
lc[it] = lc[jt];
lc[jt] = temp;
break;
}
}
if( jt == jend ) return false; //no matching edge
}
}
//check that the loop is closed
adj_type const & ilc = m->LowConn(lc[iend]);
if( ilc.size() != 2 ) return false;
if( !( ilc[0] == last && ilc[1] == first || ilc[0] == first && ilc[1] == last ) )
return false;
}
else
{
HandleType first, last;
enumerator i = ENUMUNDEF, k = ENUMUNDEF, k1 = ENUMUNDEF, k2;
MarkerType hm = m->HideMarker();
adj_type & lc = m->LowConn(GetHandle());
if( m->Count(lc.data(),lc.size(),hm) < 3 ) return false;
i = m->getNext(lc.data(),static_cast<enumerator>(lc.size()),i,hm);
HandleType q = lc[i]; //edge 0
adj_type const & qlc = m->LowConn(q);
if(m->Count(qlc.data(),qlc.size(),hm)!=2) return false;
k = m->getNext(qlc.data(),static_cast<enumerator>(qlc.size()),k,hm);
first = qlc[k]; //node 0
k = m->getNext(qlc.data(),static_cast<enumerator>(qlc.size()),k,hm);
last = qlc[k]; //node 1
i = m->getNext(lc.data(),static_cast<enumerator>(lc.size()),i,hm);
HandleType r = lc[i]; //edge 1
adj_type const & rlc = m->LowConn(r);
if(m->Count(rlc.data(),rlc.size(),hm)!=2) return false;
k1 = m->getNext(rlc.data(),static_cast<enumerator>(rlc.size()),k1,hm);
k2 = m->getNext(rlc.data(),static_cast<enumerator>(rlc.size()),k1,hm);
if( first == rlc[k1] || first == rlc[k2] )
{
HandleType temp = first;
first = last;
last = temp;
}
else if( !(last == rlc[k1] || last == rlc[k2]) )
{
adj_type::size_type jt = i+1, jend = lc.size();
while(jt < jend) if( !m->GetMarker(lc[jt],hm) ) //loop over edges
{
adj_type const & ilc = m->LowConn(lc[jt]);
if( m->Count(ilc.data(),ilc.size(),hm) != 2 ) return false;
k1 = ENUMUNDEF;
k1 = m->getNext(ilc.data(),static_cast<enumerator>(ilc.size()),k1,hm);
k2 = m->getNext(ilc.data(),static_cast<enumerator>(ilc.size()),k1,hm);
if( first == ilc[k1] || first == ilc[k2] )
{
HandleType temp = lc[i];
lc[i] = lc[jt];
lc[jt] = lc[i];
temp = first;
first = last;
last = temp;
break;
}
else if( last == ilc[k1] || last == ilc[k2] )
{
HandleType temp = lc[i];
lc[i] = lc[jt];
lc[jt] = lc[i];
break;
}
++jt;
}
if( jt == jend ) return false; //no matching edge
}
adj_type::size_type it = i, iend = lc.size()-1;
while(it != iend) if( !m->GetMarker(lc[it],hm) ) //loop over edges
{
adj_type const & ilc = m->LowConn(lc[it]);
if(m->Count(ilc.data(),ilc.size(),hm)!=2) return false;
k1 = ENUMUNDEF;
k1 = m->getNext(ilc.data(),static_cast<enumerator>(ilc.size()),k1,hm);
k2 = m->getNext(ilc.data(),static_cast<enumerator>(ilc.size()),k1,hm);
if( last == ilc[k1] )
{
last = ilc[k2];
++it;
}
else if( last == ilc[k2] )
{
last = ilc[k1];
++it;
}
else//search for the connected edge and swap with current
{
adj_type::size_type jt = it+1, jend = lc.size();
while(jt < jend) if( !m->GetMarker(lc[jt],hm) ) //loop over edges
{
adj_type const & ilc = m->LowConn(lc[jt]);
if(m->Count(ilc.data(),ilc.size(),hm)!=2) return false;
k1 = ENUMUNDEF;
k1 = m->getNext(ilc.data(),static_cast<enumerator>(ilc.size()),k1,hm);
k2 = m->getNext(ilc.data(),static_cast<enumerator>(ilc.size()),k1,hm);
if( last == ilc[k1] || last == ilc[k2] )
{
HandleType temp = lc[it];
lc[it] = lc[jt];
lc[jt] = temp;
break;
}
}
if( jt == jend ) return false; //no matching edge
}
}
//check that the loop is closed
adj_type const & ilc = m->LowConn(lc[iend]);
if(m->Count(ilc.data(),ilc.size(),hm)!=2) return false;
k1 = ENUMUNDEF;
k1 = m->getNext(ilc.data(),static_cast<enumerator>(ilc.size()),k1,hm);
k2 = m->getNext(ilc.data(),static_cast<enumerator>(ilc.size()),k1,hm);
if( !( ilc[k1] == last && ilc[k2] == first || ilc[k1] == first && ilc[k2] == last ) )
return false;
}
}
return false;
}
......
......@@ -4,5 +4,5 @@ add_subdirectory(DrawMatrix)
add_subdirectory(MatSolve)
add_subdirectory(GridGen)
add_subdirectory(FVDiscr)
#add_subdirectory(OctreeCutcell)
add_subdirectory(OctreeCutcell)
add_subdirectory(Solver)
......@@ -80,9 +80,9 @@ int main(int argc, char ** argv)
s.SetParameterEnum("rescale_iterations",8);
s.SetParameterEnum("adapt_ddpq_tolerance",0);
s.SetParameterReal("drop_tolerance",1.0e-4);
s.SetParameterReal("reuse_tolerance",1.0e-8);
s.SetParameterReal("ddpq_tolerance",0.8);
s.SetParameterReal("drop_tolerance",1.0e-2);
s.SetParameterReal("reuse_tolerance",1.0e-3);
s.SetParameterReal("ddpq_tolerance",0.4);
s.SetParameterEnum("condition_estimation",1);
......
......@@ -17,12 +17,10 @@ set(SOURCE main.cpp
add_executable(OctreeCutcell ${SOURCE})
target_link_libraries(OctreeCutcell inmost)
if(USE_MPI)
message("linking with mpi")
message("linking OctreeCutcell with mpi")
target_link_libraries(OctreeCutcell ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
message("adding link flags")
......@@ -30,19 +28,23 @@ if(USE_MPI)
endif()
endif(USE_MPI)
find_package(OpenGL)
find_package(GLUT)
if(OPENGL_FOUND)
if(GLUT_FOUND)
include_directories(${OPENGL_INCLUDE_DIR})
include_directories(${GLUT_INCLUDE_DIR})
add_definitions(-D__GRAPHICS__)
target_link_libraries(OctreeCutcell ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES})
else(GLUT_FOUND)
message("GLUT not found")
endif(GLUT_FOUND)
if(GLUT_FOUND)
message("linking OctreeCutcell with GLUT and OpenGL")
include_directories(${OPENGL_INCLUDE_DIR})
include_directories(${GLUT_INCLUDE_DIR})
add_definitions(-D__GRAPHICS__)
target_link_libraries(OctreeCutcell ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES})
else(GLUT_FOUND)
message("GLUT not found")
endif(GLUT_FOUND)
else(OPENGL_FOUND)
message("OpenGL not found")
message("OpenGL not found")
endif(OPENGL_FOUND)
install(TARGETS OctreeCutcell EXPORT inmost-targets RUNTIME DESTINATION bin)
......@@ -137,12 +137,12 @@ void cell_init_data(struct grid * g, int cell)
for(int i = 0; i < g->cells[cell].mr->size(); i++)
{
Storage::real vol = (*g->cells[cell].mr)[i]->Volume();
Storage::integer mat = (*g->cells[cell].mr)[i]->Integer(g->cell_material);
Storage::real vol = Cell(g->mesh,(*g->cells[cell].mr)[i])->Volume();
Storage::integer mat = Cell(g->mesh,(*g->cells[cell].mr)[i])->Integer(g->cell_material);
vol_and_data_by_mat & current_mat = (*g->cells[cell].data)[mat];
for(int j = 0; j < ntags; j++)
{
Storage::real val = (*g->cells[cell].mr)[i]->Real(tags[j]);
Storage::real val = Cell(g->mesh,(*g->cells[cell].mr)[i])->Real(tags[j]);
current_mat.second[tags[j]] += vol*val;
//~ if( isnan(current_mat.second[tags[j]]) ) throw -1;
}
......@@ -279,7 +279,7 @@ void fill_K(Storage::real * center, Storage::real * Kvec, Storage::real * K)
K[i*3+j] += qmat[i*3+k]*qrmat[k*3+j];
}
void cell_to_INMOST(struct grid * g, int cell, Cell * r)
void cell_to_INMOST(struct grid * g, int cell, Cell r)
{
Storage::integer mat = r->Integer(g->cell_material);
Storage::real center[3];
......@@ -344,7 +344,7 @@ void init_mesh(struct grid * g)
}
std::map<Tag,Storage::real> cell_small_unite(dynarray<Cell *,32> & unite)
std::map<Tag,Storage::real> cell_small_unite(ElementArray<Cell> & unite)
{
std::map<Tag,Storage::real> ret;
Storage::real vol = 0;
......@@ -439,13 +439,13 @@ public:
face2gl DrawFace(Element * f, int mmat, double campos[3])
face2gl DrawFace(Element f, int mmat, double campos[3])
{
double cnt[3];
face2gl ret;
f->Centroid(cnt);
ret.set_center(cnt,campos);
adjacent<Node> nodes = f->getNodes();
ElementArray<Node> nodes = f->getNodes();
if( f->nbAdjElements(CELL) == 0 )
ret.set_color(1,0,0,0.25);
......@@ -459,7 +459,7 @@ face2gl DrawFace(Element * f, int mmat, double campos[3])
ret.set_color(r,1-r,0.5,0.25);
}
for( adjacent<Node>::iterator kt = nodes.begin(); kt != nodes.end(); kt++)
for( ElementArray<Node>::iterator kt = nodes.begin(); kt != nodes.end(); kt++)
ret.add_vert(&(kt->Coords()[0]));
......@@ -511,8 +511,8 @@ void draw()
{
if( it->nbAdjElements(CELL) <= 1 )
{
int mmat = (it->BackCell() == NULL ? -1 : it->BackCell()->Integer(thegrid.cell_material));
polygons.push_back(DrawFace(&*it,mmat,campos));
int mmat = (!it->BackCell().isValid() ? -1 : it->BackCell()->Integer(thegrid.cell_material));
polygons.push_back(DrawFace(it->self(),mmat,campos));
}
}
std::sort(polygons.begin(),polygons.end());
......@@ -591,12 +591,12 @@ void draw2()
//if( m == show_region )
{
int i = 0;
adjacent<Face> faces = it->getFaces();
for(adjacent<Face>::iterator f = faces.begin(); f != faces.end(); f++) if( f->Boundary() )
ElementArray<Face> faces = it->getFaces();
for(ElementArray<Face>::iterator f = faces.begin(); f != faces.end(); f++) if( f->Boundary() )
{
adjacent<Node> nodes = f->getNodes();
ElementArray<Node> nodes = f->getNodes();
glBegin(GL_POLYGON);
for(adjacent<Node>::iterator n = nodes.begin(); n != nodes.end(); n++)
for(ElementArray<Node>::iterator n = nodes.begin(); n != nodes.end(); n++)
glVertex3dv(&n->RealArray(thegrid.mesh->CoordsTag())[0]);
glEnd();
}
......@@ -614,7 +614,7 @@ void draw2()
glBegin(GL_LINES);
for(Mesh::iteratorEdge f = thegrid.mesh->BeginEdge(); f != thegrid.mesh->EndEdge(); f++)
{
adjacent<Node> nodes = f->getNodes();
ElementArray<Node> nodes = f->getNodes();
glVertex3dv(&nodes[0].RealArray(thegrid.mesh->CoordsTag())[0]);
glVertex3dv(&nodes[1].RealArray(thegrid.mesh->CoordsTag())[0]);
}
......@@ -627,10 +627,10 @@ void draw2()
{
for(Mesh::iteratorFace f = thegrid.mesh->BeginFace(); f != thegrid.mesh->EndFace(); f++)
{
adjacent<Node> nodes = f->getNodes();
ElementArray<Node> nodes = f->getNodes();
glColor3f(0.0f,0.5f,0.0f);
glBegin(GL_POLYGON);
for(adjacent<Node>::iterator n = nodes.begin(); n != nodes.end(); n++)
for(ElementArray<Node>::iterator n = nodes.begin(); n != nodes.end(); n++)
glVertex3dv(&n->RealArray(thegrid.mesh->CoordsTag())[0]);
glEnd();
}
......@@ -812,7 +812,7 @@ void idle(void)
void vert_to_INMOST(struct grid * g, int vert, Node * v)
void vert_to_INMOST(struct grid * g, int vert, Node v)
{
(void) g; (void) vert;
Storage::real_array c = v->Coords();
......
......@@ -23,8 +23,8 @@ void default_vert_init_data(struct grid * g, int vert) {(void) g; (void) vert;};
void default_cell_init_data(struct grid * g, int cell) {(void) g; (void) cell;};
void default_vert_destroy_data(struct grid * g, int vert) {(void) g; (void) vert;};
void default_cell_destroy_data(struct grid * g, int cell) {(void) g; (void) cell;};
void default_vert_to_INMOST(struct grid * g, int vert, Node * v) {(void) g; (void) vert; (void) v;};
void default_cell_to_INMOST(struct grid * g, int cell, Cell * r) {(void) g; (void) cell; (void) r;};
void default_vert_to_INMOST(struct grid * g, int vert, Node v) {(void) g; (void) vert; (void) v;};
void default_cell_to_INMOST(struct grid * g, int cell, Cell r) {(void) g; (void) cell; (void) r;};
void default_init_mesh(struct grid * g) {(void ) g;} ;
......@@ -156,17 +156,16 @@ int cellAround(struct grid * g, int m, int side, int neighbours[1<<(DIM-1)])
void vertDestroyINMOST(struct grid *g, int m)
{
//printf("%s\n",__FUNCTION__);
if( g->verts[m].mv != NULL )
if( g->verts[m].mv.isValid() && !g->verts[m].mv.Hidden() )
{
delete g->verts[m].mv;
g->verts[m].mv = NULL;
g->verts[m].mv->Delete();
}
}
void vertCreateINMOST(struct grid * g, int m)
{
Storage::real xyz[3];
if( g->verts[m].mv != NULL ) return;
if( g->verts[m].mv.isValid() ) return;
vertGetCoord(g,m,xyz);
g->transformation(xyz);
g->verts[m].mv = g->mesh->CreateNode(xyz);
......@@ -183,7 +182,8 @@ void cellDestroyINMOST(struct grid * g, int m)
std::vector<Storage::integer> other_del;
for(int k = 0; k < g->cells[m].mr->size(); k++)
{
Storage::integer_array p = (*g->cells[m].mr)[k]->IntegerArray(g->parent);
Cell ck = Cell(g->mesh,(*g->cells[m].mr)[k]);
Storage::integer_array p = ck->IntegerArray(g->parent);
for(int j = 0; j < p.size(); j++) if( p[j] != m )
{
other_del.push_back(p[j]); //this cell is united, should delete it's other parents
......@@ -194,24 +194,24 @@ void cellDestroyINMOST(struct grid * g, int m)
break;
}
}
adjacent<Face> faces = (*g->cells[m].mr)[k]->getFaces();
adjacent<Edge> edges = (*g->cells[m].mr)[k]->getEdges();
adjacent<Node> nodes = (*g->cells[m].mr)[k]->getNodes();
delete (*g->cells[m].mr)[k];
for(adjacent<Face>::iterator f = faces.begin(); f != faces.end(); f++)
ElementArray<Face> faces = ck->getFaces();
ElementArray<Edge> edges = ck->getEdges();
ElementArray<Node> nodes = ck->getNodes();
ck->Delete();
for(ElementArray<Face>::iterator f = faces.begin(); f != faces.end(); f++)
{
if( f->nbAdjElements(CELL) == 0 )
delete &*f;
f->Delete();
}
for(adjacent<Edge>::iterator e = edges.begin(); e != edges.end(); e++)
for(ElementArray<Edge>::iterator e = edges.begin(); e != edges.end(); e++)
{
if( e->nbAdjElements(FACE) == 0 )
delete &*e;
e->Delete();
}
for(adjacent<Node>::iterator e = nodes.begin(); e != nodes.end(); e++)
for(ElementArray<Node>::iterator e = nodes.begin(); e != nodes.end(); e++)
{
if( !e->GetMarker(g->octree_node) && e->nbAdjElements(EDGE) == 0)
delete &*e;
e->Delete();
}
}
g->cells[m].mr->clear();
......@@ -290,7 +290,7 @@ int vertGetMiddle(struct grid * g, int m, int side, int edge)
return -1;
}
void cellGetFaceVerts(struct grid * g, int m, int side, int * nverts, Node * verts[54], bool * mid,int * faces, int reverse)
void cellGetFaceVerts(struct grid * g, int m, int side, int * nverts, Node verts[54], bool * mid,int * faces, int reverse)
{
int i,middle;
const int nvf[6][4] = {{0,4,6,2},{1,3,7,5},{0,1,5,4},{2,6,7,3},{0,2,3,1},{4,5,7,6}};
......@@ -319,32 +319,32 @@ void cellGetFaceVerts(struct grid * g, int m, int side, int * nverts, Node * ver
std::vector<Edge *> traverse_edges_sub(Edge * start, Edge * current, MarkerType edgeset, MarkerType visited_bridge, MarkerType visited_edge)
std::vector<Edge> traverse_edges_sub(Edge start, Edge current, MarkerType edgeset, MarkerType visited_bridge, MarkerType visited_edge)
{
//~ if( current == start ) return std::vector<Edge *> (1,start);
std::vector< std::vector<Edge *> > paths;
adjacent<Node> n = current->getNodes();
std::vector< std::vector<Edge> > paths;
ElementArray<Node> n = current->getNodes();
for(int j = 0; j < n.size(); j++)
{
if( !n[j].GetMarker(visited_bridge) )
{
n[j].SetMarker(visited_bridge);
adjacent<Edge> e = n[j].getEdges();
ElementArray<Edge> e = n[j].getEdges();
for(int i = 0; i < e.size(); i++)
{
if( &e[i] == start )
if( e[i] == start )
{
n[j].RemMarker(visited_bridge);
return std::vector<Edge *> (1,start);
return std::vector<Edge> (1,start);
}
if( e[i].GetMarker(edgeset) && !e[i].GetMarker(visited_edge) )
{
e[i].SetMarker(visited_edge);
std::vector<Edge *> ret