Commit c75eb406 authored by Kirill Terekhov's avatar Kirill Terekhov

add more timers for parallel adaptation code

parent be438831
......@@ -476,9 +476,10 @@ namespace INMOST
//for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
// if( it->getNodes().size() != 2 ) {REPORT_STR("edge " << it->LocalID() << " has " << it->getNodes().size() << " nodes ");}
//EXIT_BLOCK();
ENTER_BLOCK();
m->CheckSetLinks(__FILE__,__LINE__);
CheckParentSet(__FILE__,__LINE__);
EXIT_BLOCK();
/*
ENTER_BLOCK();
m->ResolveSets();
......@@ -599,318 +600,411 @@ namespace INMOST
Storage::real xyz[3] = {0,0,0};
//7.split all edges of the current schedule
ENTER_BLOCK();
for(Storage::integer it = 0; it < m->EdgeLastLocalID(); ++it) if( m->isValidEdge(it) )
{
Edge e = m->EdgeByLocalID(it);
if( !e.Hidden() && indicator[e] == schedule_counter )
int new_nodes = 0, splits = 0;
double t1, t2, tadj = 0, tcreate = 0, thanging = 0, tdata = 0, tsplit = 0;
for(Storage::integer it = 0; it < m->EdgeLastLocalID(); ++it) if( m->isValidEdge(it) )
{
//ENTER_BLOCK();
//for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
// if( it->getNodes().size() != 2 ) {REPORT_STR("edge " << it->LocalID() << " has " << it->getNodes().size() << " nodes ");}
//EXIT_BLOCK();
//REPORT_STR("split edge " << e.LocalID() << " nodes " << e.getBeg().LocalID() << "," << e.getEnd().LocalID() << " level " << level[e] << " lc size " << m->LowConn(e.GetHandle()).size() );
//ElementArray<Node> nodes = e.getNodes();
//for(int q = 0; q < nodes.size(); ++q) REPORT_STR("node " << nodes[q].GetHandle() << " " << nodes[q].LocalID() << (nodes[q].Hidden()?" hidden " : " good ") );
//remember adjacent faces that should get information about new hanging node
ElementArray<Face> edge_faces = e.getFaces();
//location on the center of the edge
for(Storage::integer d = 0; d < m->GetDimensions(); ++d)
xyz[d] = (e.getBeg().Coords()[d]+e.getEnd().Coords()[d])*0.5;
//todo: request transformation of node location according to geometrical model
//create middle node
Node n = m->CreateNode(xyz);
//set increased level for new node
level[n] = level[e.getBeg()] = level[e.getEnd()] = level[e]+1;
//for each face provide link to a new hanging node
for(ElementArray<Face>::size_type kt = 0; kt < edge_faces.size(); ++kt)
hanging_nodes[edge_faces[kt]].push_back(n);
//CheckClosure(__FILE__,__LINE__);
//split the edge by the middle node
ElementArray<Edge> new_edges = Edge::SplitEdge(e,ElementArray<Node>(m,1,n.GetHandle()),0);
for(ElementArray<Face>::size_type kt = 0; kt < edge_faces.size(); ++kt) assert(edge_faces[kt].Closure());
//set increased level for new edges
level[new_edges[0]] = level[new_edges[1]] = level[e]+1;
Edge e = m->EdgeByLocalID(it);
if( !e.Hidden() && indicator[e] == schedule_counter )
{
//ENTER_BLOCK();
//for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
// if( it->getNodes().size() != 2 ) {REPORT_STR("edge " << it->LocalID() << " has " << it->getNodes().size() << " nodes ");}
//EXIT_BLOCK();
//REPORT_STR("split edge " << e.LocalID() << " nodes " << e.getBeg().LocalID() << "," << e.getEnd().LocalID() << " level " << level[e] << " lc size " << m->LowConn(e.GetHandle()).size() );
//ElementArray<Node> nodes = e.getNodes();
//for(int q = 0; q < nodes.size(); ++q) REPORT_STR("node " << nodes[q].GetHandle() << " " << nodes[q].LocalID() << (nodes[q].Hidden()?" hidden " : " good ") );
//remember adjacent faces that should get information about new hanging node
t1 = Timer();
ElementArray<Face> edge_faces = e.getFaces();
t2 = Timer(), tadj += t2 - t1, t1 = t2;
//location on the center of the edge
for(Storage::integer d = 0; d < m->GetDimensions(); ++d)
xyz[d] = (e.getBeg().Coords()[d]+e.getEnd().Coords()[d])*0.5;
//todo: request transformation of node location according to geometrical model
//create middle node
Node n = m->CreateNode(xyz);
new_nodes++;
//set increased level for new node
level[n] = level[e.getBeg()] = level[e.getEnd()] = level[e]+1;
t2 = Timer(), tcreate += t2 - t1, t1 = t2;
//for each face provide link to a new hanging node
for(ElementArray<Face>::size_type kt = 0; kt < edge_faces.size(); ++kt)
hanging_nodes[edge_faces[kt]].push_back(n);
t2 = Timer(), thanging += t2 - t1, t1 = t2;
//CheckClosure(__FILE__,__LINE__);
//split the edge by the middle node
ElementArray<Edge> new_edges = Edge::SplitEdge(e,ElementArray<Node>(m,1,n.GetHandle()),0);
splits++;
for(ElementArray<Face>::size_type kt = 0; kt < edge_faces.size(); ++kt) assert(edge_faces[kt].Closure());
t2 = Timer(), tsplit += t2 - t1, t1 = t2;
//set increased level for new edges
level[new_edges[0]] = level[new_edges[1]] = level[e]+1;
#if defined(USE_AUTODIFF) && defined(USE_SOLVER)
if( model ) model->EdgeRefinement(e,new_edges);
if( model ) model->EdgeRefinement(e,new_edges);
#endif
//for(int q = 0; q < 2; ++q)
//{
// REPORT_STR("new edges["<<q<<"]" << new_edges[q].LocalID() << " nodes " << new_edges[q].getBeg().LocalID() << "," << new_edges[q].getEnd().LocalID() << " level " << level[new_edges[q]]);
//}
//CheckClosure(__FILE__,__LINE__);
//if( !Element::CheckConnectivity(m) ) std::cout << __FILE__ << ":" << __LINE__ << " broken connectivity" << std::endl;
t2 = Timer(), tdata += t2 - t1, t1 = t2;
//for(int q = 0; q < 2; ++q)
//{
// REPORT_STR("new edges["<<q<<"]" << new_edges[q].LocalID() << " nodes " << new_edges[q].getBeg().LocalID() << "," << new_edges[q].getEnd().LocalID() << " level " << level[new_edges[q]]);
//}
//CheckClosure(__FILE__,__LINE__);
//if( !Element::CheckConnectivity(m) ) std::cout << __FILE__ << ":" << __LINE__ << " broken connectivity" << std::endl;
}
}
REPORT_VAL("adjacencies", tadj);
REPORT_VAL("create (new nodes)", tcreate);
REPORT_VAL("hanging", thanging);
REPORT_VAL("data", tdata);
REPORT_VAL("split", tsplit);
REPORT_VAL("new nodes", new_nodes);
REPORT_VAL("splits", splits);
}
EXIT_BLOCK();
ENTER_BLOCK();
assert(Element::CheckConnectivity(m));
CheckClosure(__FILE__,__LINE__);
EXIT_BLOCK();
//8.split all faces of the current schedule, using hanging nodes on edges
ENTER_BLOCK();
for(Storage::integer it = 0; it < m->FaceLastLocalID(); ++it) if( m->isValidFace(it) )
{
Face f = m->FaceByLocalID(it);
if( !f.Hidden() && indicator[f] == schedule_counter )
int new_nodes = 0, new_edges = 0, splits = 0;
double t1, t2, tadj = 0, tcreate = 0, thanging = 0, tdata = 0, tsplit = 0;
for(Storage::integer it = 0; it < m->FaceLastLocalID(); ++it) if( m->isValidFace(it) )
{
#if !defined(NDEBUG)
ElementArray<Edge> face_edges = f.getEdges();
for(ElementArray<Edge>::iterator jt = face_edges.begin(); jt != face_edges.end(); ++jt)
Face f = m->FaceByLocalID(it);
if( !f.Hidden() && indicator[f] == schedule_counter )
{
if( level[*jt] != level[f]+1 )
#if !defined(NDEBUG)
ElementArray<Edge> face_edges = f.getEdges();
for(ElementArray<Edge>::iterator jt = face_edges.begin(); jt != face_edges.end(); ++jt)
{
std::cout << m->GetProcessorRank() << " face " << f.LocalID();
std::cout << " " << Element::StatusName(f.GetStatus()) << " owner " << f.Integer(m->OwnerTag());
std::cout << " lvl " << level[f] << " ind " << indicator[f];
std::cout << " edge " << jt->LocalID();
std::cout << " " << Element::StatusName(jt->GetStatus()) << " owner " << jt->Integer(m->OwnerTag());
std::cout << " lvl " << level[*jt] << " ind " << indicator[*jt];
std::cout << std::endl;
ElementArray<Cell> adj_cells;
adj_cells = f.getCells();
std::cout << "face cells ";
for(ElementArray<Cell>::iterator kt = adj_cells.begin(); kt != adj_cells.end(); ++kt)
{
std::cout << m->GetProcessorRank() << " face " << kt->LocalID();
std::cout << " " << Element::StatusName(kt->GetStatus()) << " owner " << kt->Integer(m->OwnerTag());
std::cout << " lvl " << level[*kt] << " ind " << indicator[*kt];
std::cout << std::endl;
}
std::cout << "edge cells ";
adj_cells = jt->getCells();
for(ElementArray<Cell>::iterator kt = adj_cells.begin(); kt != adj_cells.end(); ++kt)
if( level[*jt] != level[f]+1 )
{
std::cout << m->GetProcessorRank() << " face " << kt->LocalID();
std::cout << " " << Element::StatusName(kt->GetStatus()) << " owner " << kt->Integer(m->OwnerTag());
std::cout << " lvl " << level[*kt] << " ind " << indicator[*kt];
std::cout << m->GetProcessorRank() << " face " << f.LocalID();
std::cout << " " << Element::StatusName(f.GetStatus()) << " owner " << f.Integer(m->OwnerTag());
std::cout << " lvl " << level[f] << " ind " << indicator[f];
std::cout << " edge " << jt->LocalID();
std::cout << " " << Element::StatusName(jt->GetStatus()) << " owner " << jt->Integer(m->OwnerTag());
std::cout << " lvl " << level[*jt] << " ind " << indicator[*jt];
std::cout << std::endl;
ElementArray<Cell> adj_cells;
adj_cells = f.getCells();
std::cout << "face cells ";
for(ElementArray<Cell>::iterator kt = adj_cells.begin(); kt != adj_cells.end(); ++kt)
{
std::cout << m->GetProcessorRank() << " face " << kt->LocalID();
std::cout << " " << Element::StatusName(kt->GetStatus()) << " owner " << kt->Integer(m->OwnerTag());
std::cout << " lvl " << level[*kt] << " ind " << indicator[*kt];
std::cout << std::endl;
}
std::cout << "edge cells ";
adj_cells = jt->getCells();
for(ElementArray<Cell>::iterator kt = adj_cells.begin(); kt != adj_cells.end(); ++kt)
{
std::cout << m->GetProcessorRank() << " face " << kt->LocalID();
std::cout << " " << Element::StatusName(kt->GetStatus()) << " owner " << kt->Integer(m->OwnerTag());
std::cout << " lvl " << level[*kt] << " ind " << indicator[*kt];
std::cout << std::endl;
}
}
assert(level[*jt] == level[f]+1);
}
assert(level[*jt] == level[f]+1);
}
#endif //NDEBUG
//connect face center to hanging nodes of the face
Storage::reference_array face_hanging_nodes = hanging_nodes[f];
//remember adjacent cells that should get information about new hanging node
//and new hanging edges
ElementArray<Cell> face_cells = f.getCells();
//create node at face center
//f->Centroid(xyz);
for(int d = 0; d < 3; ++d) xyz[d] = 0.0;
for(Storage::reference_array::size_type kt = 0; kt < face_hanging_nodes.size(); ++kt)
for(int d = 0; d < 3; ++d) xyz[d] += face_hanging_nodes[kt].getAsNode().Coords()[d];
for(int d = 0; d < 3; ++d) xyz[d] /= (Storage::real)face_hanging_nodes.size();
//todo: request transformation of node location according to geometrical model
//create middle node
Node n = m->CreateNode(xyz);
//set increased level for the new node
level[n] = level[f]+1;
//for each cell provide link to new hanging node
for(ElementArray<Face>::size_type kt = 0; kt < face_cells.size(); ++kt)
hanging_nodes[face_cells[kt]].push_back(n);
ElementArray<Node> edge_nodes(m,2); //to create new edges
ElementArray<Edge> hanging_edges(m,face_hanging_nodes.size());
edge_nodes[0] = n;
for(Storage::reference_array::size_type kt = 0; kt < face_hanging_nodes.size(); ++kt)
{
edge_nodes[1] = face_hanging_nodes[kt].getAsNode();
hanging_edges[kt] = m->CreateEdge(edge_nodes).first;
//set increased level for new edges
level[hanging_edges[kt]] = level[f]+1;
}
//split the face by these edges
ElementArray<Face> new_faces = Face::SplitFace(f,hanging_edges,0);
//set increased level to new faces
for(ElementArray<Face>::size_type kt = 0; kt < new_faces.size(); ++kt)
level[new_faces[kt]] = level[f]+1;
t1 = Timer();
//remember adjacent cells that should get information about new hanging node
//and new hanging edges
ElementArray<Cell> face_cells = f.getCells();
t2 = Timer(), tadj += t2 - t1, t1 = t2;
//create node at face center
//f->Centroid(xyz);
for(int d = 0; d < 3; ++d) xyz[d] = 0.0;
//connect face center to hanging nodes of the face
Storage::reference_array face_hanging_nodes = hanging_nodes[f];
for(Storage::reference_array::size_type kt = 0; kt < face_hanging_nodes.size(); ++kt)
for(int d = 0; d < 3; ++d) xyz[d] += face_hanging_nodes[kt].getAsNode().Coords()[d];
for(int d = 0; d < 3; ++d) xyz[d] /= (Storage::real)face_hanging_nodes.size();
//todo: request transformation of node location according to geometrical model
//create middle node
Node n = m->CreateNode(xyz);
new_nodes++;
//set increased level for the new node
level[n] = level[f]+1;
t2 = Timer(), tcreate += t2 - t1, t1 = t2;
//for each cell provide link to new hanging node
for(ElementArray<Face>::size_type kt = 0; kt < face_cells.size(); ++kt)
hanging_nodes[face_cells[kt]].push_back(n);
ElementArray<Node> edge_nodes(m,2); //to create new edges
ElementArray<Edge> hanging_edges(m,face_hanging_nodes.size());
edge_nodes[0] = n;
for(Storage::reference_array::size_type kt = 0; kt < face_hanging_nodes.size(); ++kt)
{
edge_nodes[1] = face_hanging_nodes[kt].getAsNode();
hanging_edges[kt] = m->CreateEdge(edge_nodes).first;
new_edges++;
//set increased level for new edges
level[hanging_edges[kt]] = level[f]+1;
}
t2 = Timer(), thanging += t2 - t1, t1 = t2;
//split the face by these edges
ElementArray<Face> new_faces = Face::SplitFace(f,hanging_edges,0);
splits++;
t2 = Timer(), tsplit += t2 - t1, t1 = t2;
//set increased level to new faces
for(ElementArray<Face>::size_type kt = 0; kt < new_faces.size(); ++kt)
level[new_faces[kt]] = level[f]+1;
#if defined(USE_AUTODIFF) && defined(USE_SOLVER)
if( model ) model->FaceRefinement(f,new_faces);
if( model ) model->FaceRefinement(f,new_faces);
#endif
t2 = Timer(), tdata += t2 - t1, t1 = t2;
}
}
REPORT_VAL("adjacencies", tadj);
REPORT_VAL("create (new nodes)", tcreate);
REPORT_VAL("hanging (new edges)", thanging);
REPORT_VAL("data", tdata);
REPORT_VAL("split", tsplit);
REPORT_VAL("new nodes", new_nodes);
REPORT_VAL("new edges", new_edges);
REPORT_VAL("splits", splits);
}
EXIT_BLOCK();
ENTER_BLOCK();
assert(Element::CheckConnectivity(m));
CheckClosure(__FILE__,__LINE__);
EXIT_BLOCK();
TagReferenceArray internal_face_edges;
MarkerType mark_cell_edges, mark_hanging_nodes;
ENTER_BLOCK();
//this tag helps recreate internal face
TagReferenceArray internal_face_edges = m->CreateTag("INTERNAL_FACE_EDGES",DATA_REFERENCE,NODE,NODE,4);
internal_face_edges = m->CreateTag("INTERNAL_FACE_EDGES",DATA_REFERENCE,NODE,NODE,4);
//this marker helps detect edges of current cell only
MarkerType mark_cell_edges = m->CreateMarker();
mark_cell_edges = m->CreateMarker();
//this marker helps detect nodes hanging on edges of unrefined cell
MarkerType mark_hanging_nodes = m->CreateMarker();
mark_hanging_nodes = m->CreateMarker();
EXIT_BLOCK();
//9.split all cells of the current schedule
ENTER_BLOCK();
for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
{
Cell c = m->CellByLocalID(it);
if( !c.Hidden() && indicator[c] == schedule_counter )
int new_nodes = 0, new_edges = 0, new_faces = 0, new_sets = 0, splits = 0;
double t1, t2, tadj = 0, tcreate = 0, thanging1 = 0, thanging2 = 0, tset = 0, tdata = 0, tsplit = 0, tsort = 0, tpconn = 0;
for(Storage::integer it = 0; it < m->CellLastLocalID(); ++it) if( m->isValidCell(it) )
{
Storage::reference_array cell_hanging_nodes = hanging_nodes[c]; //nodes to be connected
//create node at cell center
for(int d = 0; d < 3; ++d) xyz[d] = 0.0;
for(Storage::reference_array::size_type kt = 0; kt < cell_hanging_nodes.size(); ++kt)
for(int d = 0; d < 3; ++d) xyz[d] += cell_hanging_nodes[kt].getAsNode().Coords()[d];
for(int d = 0; d < 3; ++d) xyz[d] /= (Storage::real)cell_hanging_nodes.size();
//c->Centroid(xyz);
//todo: request transformation of node location according to geometrical model
//create middle node
Node n = m->CreateNode(xyz);
//set increased level for the new node
level[n] = level[c]+1;
//retrive all edges of current face to mark them
ElementArray<Edge> cell_edges = c.getEdges();
Cell c = m->CellByLocalID(it);
if( !c.Hidden() && indicator[c] == schedule_counter )
{
t1 = Timer();
Storage::reference_array cell_hanging_nodes = hanging_nodes[c]; //nodes to be connected
//create node at cell center
for(int d = 0; d < 3; ++d) xyz[d] = 0.0;
for(Storage::reference_array::size_type kt = 0; kt < cell_hanging_nodes.size(); ++kt)
for(int d = 0; d < 3; ++d) xyz[d] += cell_hanging_nodes[kt].getAsNode().Coords()[d];
for(int d = 0; d < 3; ++d) xyz[d] /= (Storage::real)cell_hanging_nodes.size();
//c->Centroid(xyz);
//todo: request transformation of node location according to geometrical model
//create middle node
Node n = m->CreateNode(xyz);
new_nodes++;
//set increased level for the new node
level[n] = level[c]+1;
t2 = Timer(), tcreate += t2 - t1, t1 = t2;
//retrive all edges of current face to mark them
ElementArray<Edge> cell_edges = c.getEdges();
#if !defined(NDEBUG)
for(ElementArray<Edge>::iterator jt = cell_edges.begin(); jt != cell_edges.end(); ++jt) assert(level[*jt] == level[c]+1);
ElementArray<Face> cell_faces = c.getFaces();
for(ElementArray<Face>::iterator jt = cell_faces.begin(); jt != cell_faces.end(); ++jt) assert(level[*jt] == level[c]+1);
for(ElementArray<Edge>::iterator jt = cell_edges.begin(); jt != cell_edges.end(); ++jt) assert(level[*jt] == level[c]+1);
ElementArray<Face> cell_faces = c.getFaces();
for(ElementArray<Face>::iterator jt = cell_faces.begin(); jt != cell_faces.end(); ++jt) assert(level[*jt] == level[c]+1);
#endif //NDEBUG
//mark all edges so that we can retive them later
cell_edges.SetMarker(mark_cell_edges);
//connect face center to centers of faces by edges
ElementArray<Node> edge_nodes(m,2);
ElementArray<Edge> edges_to_faces(m,cell_hanging_nodes.size());
edge_nodes[0] = n;
for(Storage::reference_array::size_type kt = 0; kt < cell_hanging_nodes.size(); ++kt)
{
assert(cell_hanging_nodes[kt].isValid());
//todo: unmark hanging node on edge if no more cells depend on it
edge_nodes[1] = cell_hanging_nodes[kt].getAsNode();
edges_to_faces[kt] = m->CreateEdge(edge_nodes).first;
//set increased level for new edges
level[edges_to_faces[kt]] = level[c]+1;
//for each node other then the hanging node of the face
//(this is hanging node on the edge)
//we record a pair of edges to reconstruct internal faces
ElementArray<Edge> hanging_edges = cell_hanging_nodes[kt].getEdges(mark_cell_edges,0);
for(ElementArray<Edge>::size_type lt = 0; lt < hanging_edges.size(); ++lt)
t2 = Timer(), tadj += t2 - t1, t1 = t2;
//mark all edges so that we can retive them later
cell_edges.SetMarker(mark_cell_edges);
//connect face center to centers of faces by edges
ElementArray<Node> edge_nodes(m,2);
ElementArray<Edge> edges_to_faces(m,cell_hanging_nodes.size());
edge_nodes[0] = n;
for(Storage::reference_array::size_type kt = 0; kt < cell_hanging_nodes.size(); ++kt)
{
//get hanging node on the edge
assert(hanging_edges[lt].getBeg() == cell_hanging_nodes[kt] || hanging_edges[lt].getEnd() == cell_hanging_nodes[kt]);
Node v = hanging_edges[lt].getBeg() == cell_hanging_nodes[kt]? hanging_edges[lt].getEnd() : hanging_edges[lt].getBeg();
//mark so that we can collect all of them
v.SetMarker(mark_hanging_nodes);
//fill the edges
Storage::reference_array face_edges = internal_face_edges[v];
//fill first two in forward order
//this way we make a closed loop
assert(face_edges[0] == InvalidElement() || face_edges[2] == InvalidElement());
if( face_edges[0] == InvalidElement() )
assert(cell_hanging_nodes[kt].isValid());
//todo: unmark hanging node on edge if no more cells depend on it
edge_nodes[1] = cell_hanging_nodes[kt].getAsNode();
edges_to_faces[kt] = m->CreateEdge(edge_nodes).first;
new_edges++;
//set increased level for new edges
level[edges_to_faces[kt]] = level[c]+1;
//for each node other then the hanging node of the face
//(this is hanging node on the edge)
//we record a pair of edges to reconstruct internal faces
ElementArray<Edge> hanging_edges = cell_hanging_nodes[kt].getEdges(mark_cell_edges,0);
for(ElementArray<Edge>::size_type lt = 0; lt < hanging_edges.size(); ++lt)
{
face_edges[0] = edges_to_faces[kt];
face_edges[1] = hanging_edges[lt];
//get hanging node on the edge
assert(hanging_edges[lt].getBeg() == cell_hanging_nodes[kt] || hanging_edges[lt].getEnd() == cell_hanging_nodes[kt]);
Node v = hanging_edges[lt].getBeg() == cell_hanging_nodes[kt]? hanging_edges[lt].getEnd() : hanging_edges[lt].getBeg();
//mark so that we can collect all of them
v.SetMarker(mark_hanging_nodes);
//fill the edges
Storage::reference_array face_edges = internal_face_edges[v];
//fill first two in forward order
//this way we make a closed loop
assert(face_edges[0] == InvalidElement() || face_edges[2] == InvalidElement());
if( face_edges[0] == InvalidElement() )
{
face_edges[0] = edges_to_faces[kt];
face_edges[1] = hanging_edges[lt];
}
else //last two in reverse
{
assert(face_edges[2] ==InvalidElement());
face_edges[2] = hanging_edges[lt];
face_edges[3] = edges_to_faces[kt];
}
}
else //last two in reverse
}
//remove marker from cell edges
cell_edges.RemMarker(mark_cell_edges);
t2 = Timer(), thanging1 += t2 - t1, t1 = t2;
//now we have to create internal faces
ElementArray<Node> edge_hanging_nodes = c.getNodes(mark_hanging_nodes,0);
ElementArray<Face> internal_faces(m,edge_hanging_nodes.size());
//unmark hanging nodes on edges
edge_hanging_nodes.RemMarker(mark_hanging_nodes);
for(ElementArray<Node>::size_type kt = 0; kt < edge_hanging_nodes.size(); ++kt)
{
//create a face based on collected edges
Storage::reference_array face_edges = internal_face_edges[edge_hanging_nodes[kt]];
assert(face_edges[0].isValid());
assert(face_edges[1].isValid());
assert(face_edges[2].isValid());
assert(face_edges[3].isValid());
internal_faces[kt] = m->CreateFace(ElementArray<Edge>(m,face_edges.begin(),face_edges.end())).first;
new_faces++;
//set increased level
level[internal_faces[kt]] = level[c]+1;
//clean up structure, so that other cells can use it
edge_hanging_nodes[kt].DelData(internal_face_edges);
}
t2 = Timer(), thanging2 += t2 - t1, t1 = t2;
//if( c.GlobalID() == 228 )
//{
// double cnt[3];
// c.Centroid(cnt);
// std::cout << "Split CELL:" << c.LocalID() << " " << c.GlobalID() << " " << Element::StatusName(c.GetStatus()) << " " << cnt[0] << " " << cnt[1] << " " << cnt[2] << std::endl;
//
//}
//split the cell
//retrive parent set
ElementSet parent(m,parent_set[c]);
//create set corresponding to old coarse cell
//Storage::real cnt[3];
//c.Centroid(cnt);
std::stringstream set_name;
//set_name << parent.GetName() << "_C" << c.GlobalID(); //rand may be unsafe
if( parent == root )
set_name << "AM_R" << set_id[c];
else
set_name << parent.GetName() << "C" << set_id[c];
//set_name << base64_encode_((unsigned char *)cnt,3*sizeof(double)/sizeof(unsigned char));
#if !defined(NDEBUG)
ElementSet check_set = m->GetSet(set_name.str());
if( check_set.isValid() )
{
std::cout << rank << " set " << set_name.str() << " for cell " << c.GlobalID() << " " << Element::StatusName(c.GetStatus()) << " already exists" << std::endl;
if( check_set->HaveParent() )
std::cout << rank << " parent is " << check_set->GetParent()->GetName() << " cell parent is " << parent.GetName() << std::endl;
std::cout << rank << " Elements of " << check_set.GetName() << ": ";
for(ElementSet::iterator it = check_set.Begin(); it != check_set.End(); ++it)
std::cout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << "," << it->GlobalID() << "," << Element::StatusName(c.GetStatus()) << "," << level[*it] << "," << indicator[*it] << " ";
std::cout << std::endl;
std::cout << rank << " Elements of " << parent.GetName() << ": ";
for(ElementSet::iterator it = parent.Begin(); it != parent.End(); ++it)
std::cout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << "," << it->GlobalID() << "," << Element::StatusName(c.GetStatus()) << "," << level[*it] << "," << indicator[*it] << " ";
std::cout << std::endl;
if( parent.HaveChild() )
{
assert(face_edges[2] ==InvalidElement());
face_edges[2] = hanging_edges[lt];
face_edges[3] = edges_to_faces[kt];
std::cout << rank << " Children of " << parent.GetName() << ": ";
for(ElementSet jt = parent.GetChild(); jt.isValid(); jt = jt.GetSibling() )
std::cout << jt.GetName() << " size " << jt.Size() << " ";
std::cout << std::endl;
}
exit(-1);
}
}
//remove marker from cell edges
cell_edges.RemMarker(mark_cell_edges);
//now we have to create internal faces
ElementArray<Node> edge_hanging_nodes = c.getNodes(mark_hanging_nodes,0);
ElementArray<Face> internal_faces(m,edge_hanging_nodes.size());
//unmark hanging nodes on edges
edge_hanging_nodes.RemMarker(mark_hanging_nodes);
for(ElementArray<Node>::size_type kt = 0; kt < edge_hanging_nodes.size(); ++kt)
{
//create a face based on collected edges
Storage::reference_array face_edges = internal_face_edges[edge_hanging_nodes[kt]];
assert(face_edges[0].isValid());
assert(face_edges[1].isValid());
assert(face_edges[2].isValid());
assert(face_edges[3].isValid());
internal_faces[kt] = m->CreateFace(ElementArray<Edge>(m,face_edges.begin(),face_edges.end())).first;
//set increased level
level[internal_faces[kt]] = level[c]+1;
//clean up structure, so that other cells can use it
edge_hanging_nodes[kt].DelData(internal_face_edges);
}
//if( c.GlobalID() == 228 )
//{
// double cnt[3];
// c.Centroid(cnt);
// std::cout << "Split CELL:" << c.LocalID() << " " << c.GlobalID() << " " << Element::StatusName(c.GetStatus()) << " " << cnt[0] << " " << cnt[1] << " " << cnt[2] << std::endl;
//
//}
//split the cell
//retrive parent set
ElementSet parent(m,parent_set[c]);
//create set corresponding to old coarse cell
//Storage::real cnt[3];
//c.Centroid(cnt);
std::stringstream set_name;
//set_name << parent.GetName() << "_C" << c.GlobalID(); //rand may be unsafe
if( parent == root )
set_name << "AM_R" << set_id[c];
else
set_name << parent.GetName() << "C" << set_id[c];
//set_name << base64_encode_((unsigned char *)cnt,3*sizeof(double)/sizeof(unsigned char));
ElementSet check_set = m->GetSet(set_name.str());
if( check_set.isValid() )
{
std::cout << rank << " set " << set_name.str() << " for cell " << c.GlobalID() << " " << Element::StatusName(c.GetStatus()) << " already exists" << std::endl;
if( check_set->HaveParent() )
std::cout << rank << " parent is " << check_set->GetParent()->GetName() << " cell parent is " << parent.GetName() << std::endl;
std::cout << rank << " Elements of " << check_set.GetName() << ": ";
for(ElementSet::iterator it = check_set.Begin(); it != check_set.End(); ++it)
std::cout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << "," << it->GlobalID() << "," << Element::StatusName(c.GetStatus()) << "," << level[*it] << "," << indicator[*it] << " ";
std::cout << std::endl;
std::cout << rank << " Elements of " << parent.GetName() << ": ";
for(ElementSet::iterator it = parent.Begin(); it != parent.End(); ++it)
std::cout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << "," << it->GlobalID() << "," << Element::StatusName(c.GetStatus()) << "," << level[*it] << "," << indicator[*it] << " ";
std::cout << std::endl;
if( parent.HaveChild() )
#endif
ElementSet cell_set = m->CreateSetUnique(set_name.str()).first;
new_sets++;
//cell_set->SetExchange(ElementSet::SYNC_ELEMENTS_ALL);
level[cell_set] = level[c]+1;
set_id[cell_set] = set_id[c];
t2 = Timer(), tset += t2 - t1, t1 = t2;
ElementArray<Cell> new_cells = Cell::SplitCell(c,internal_faces,0);
splits++;
t2 = Timer(), tsplit += t2 - t1, t1 = t2;
std::sort(new_cells.begin(),new_cells.end(),Mesh::CentroidComparator(m));
t2 = Timer(), tsort += t2 - t1, t1 = t2;
//set up increased level for the new cells
for(ElementArray<Cell>::size_type kt = 0; kt < new_cells.size(); ++kt)
{
std::cout << rank << " Children of " << parent.GetName() << ": ";
for(ElementSet jt = parent.GetChild(); jt.isValid(); jt = jt.GetSibling() )
std::cout << jt.GetName() << " size " << jt.Size() << " ";
set_id[new_cells[kt]] = (int)kt;
level[new_cells[kt]] = level[c]+1;
cell_set.PutElement(new_cells[kt]);
parent_set[new_cells[kt]] = cell_set.GetHandle();
}
/*
if( check_set.isValid() )
{
std::cout << rank << " Elements: ";
for(ElementSet::iterator it = check_set.Begin(); it != check_set.End(); ++it)
std::cout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << "," << it->GlobalID() << "," << Element::StatusName(c.GetStatus()) << "," << level[*it] << " ";
std::cout << std::endl;
}
exit(-1);
}
ElementSet cell_set = m->CreateSetUnique(set_name.str()).first;
//cell_set->SetExchange(ElementSet::SYNC_ELEMENTS_ALL);
level[cell_set] = level[c]+1;
set_id[cell_set] = set_id[c];
ElementArray<Cell> new_cells = Cell::SplitCell(c,internal_faces,0);
std::sort(new_cells.begin(),new_cells.end(),Mesh::CentroidComparator(m));
//set up increased level for the new cells
for(ElementArray<Cell>::size_type kt = 0; kt < new_cells.size(); ++kt)
{
set_id[new_cells[kt]] = (int)kt;
level[new_cells[kt]] = level[c]+1;
cell_set.PutElement(new_cells[kt]);
parent_set[new_cells[kt]] = cell_set.GetHandle();
}
/*
if( check_set.isValid() )
{
std::cout << rank << " Elements: ";
for(ElementSet::iterator it = check_set.Begin(); it != check_set.End(); ++it)
std::cout << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << "," << it->GlobalID() << "," << Element::StatusName(c.GetStatus()) << "," << level[*it] << " ";
std::cout << std::endl;
}
*/
//if( !cell_set->HaveParent() )
parent.AddChild(cell_set);
*/
//if( !cell_set->HaveParent() )
parent.AddChild(cell_set);
t2 = Timer(), tpconn += t2 - t1, t1 = t2;
#if defined(USE_AUTODIFF) && defined(USE_SOLVER)
if( model ) model->CellRefinement(c,new_cells);
if( model ) model->CellRefinement(c,new_cells);
#endif
//else assert(cell_set->GetParent() == parent);
//increment number of refined cells
ret++;
t2 = Timer(), tdata += t2 - t1, t1 = t2;
//else assert(cell_set->GetParent() == parent);
//increment number of refined cells
ret++;
}
}
REPORT_VAL("adjacencies", tadj);
REPORT_VAL("create (new nodes)", tcreate);
REPORT_VAL("hanging1 (new edges)", thanging1);
REPORT_VAL("hanging2 (new faces)", thanging2);
REPORT_VAL("create set", tset);