Commit 2c7f5178 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Fixed ResolveShared

parent 0a11ce1e
......@@ -294,7 +294,7 @@ namespace INMOST
AdaptiveMesh::AdaptiveMesh() : Mesh()
{
//create a tag that stores maximal refinement level of each element
level = CreateTag("REFINEMENT_LEVEL",DATA_INTEGER,CELL|FACE|EDGE|NODE|ESET,NONE,1);
level = CreateTag("REFINEMENT_LEVEL",DATA_INTEGER,CELL|FACE|EDGE|ESET,NONE,1);
tag_status = CreateTag("TAG_STATUS",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
tag_an = CreateTag("TAG_AN",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1);
ref_tag = CreateTag("REF",DATA_REFERENCE,CELL|FACE|EDGE|NODE,NONE);
......@@ -376,6 +376,7 @@ namespace INMOST
scheduled = Integrate(scheduled);
if( scheduled ) schedule_counter++;
}
ExchangeData(indicator,CELL | FACE | EDGE,0);
//6.Refine
BeginModification();
while(schedule_counter)
......@@ -396,7 +397,7 @@ namespace INMOST
//create middle node
Node n = CreateNode(xyz);
//set increased level for new node
level[n] = level[e.getBeg()] = level[e.getEnd()] = level[e]+1;
//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);
......@@ -427,7 +428,7 @@ namespace INMOST
//create middle node
Node n = CreateNode(xyz);
//set increased level for the new node
level[n] = level[f]+1;
//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);
......@@ -471,7 +472,7 @@ namespace INMOST
//create middle node
Node n = CreateNode(xyz);
//set increased level for the new node
level[n] = level[c]+1;
//level[n] = level[c]+1;
//retrive all edges of current face to mark them
ElementArray<Edge> cell_edges = c.getEdges();
//mark all edges so that we can retive them later
......@@ -594,7 +595,7 @@ namespace INMOST
//11. Restore parallel connectivity, global ids
ResolveShared(true);
//if (call_counter == 0)
ResolveModification();
// ResolveModification(false,1);
//12. Let the user update their data
//todo: call back function for New() cells
//13. Delete old elements of the mesh
......@@ -602,7 +603,8 @@ namespace INMOST
//14. Done
//cout << rank << ": Before end " << endl;
EndModification();
ExchangeData(hanging_nodes,CELL | FACE,0);
CheckCentroids();
//ExchangeData(hanging_nodes,CELL | FACE,0);
//cout << rank << ": After end " << endl;
//reorder element's data to free up space
ReorderEmpty(CELL|FACE|EDGE|NODE);
......@@ -648,11 +650,13 @@ namespace INMOST
indicator[e] = INT_MAX;
for(ElementArray<Element>::size_type kt = 0; kt < adj.size(); ++kt)
if( level[e] == level[adj[kt]]) indicator[e] = std::min(indicator[e],indicator[adj[kt]]);
assert(indicator[e] != INT_MAX);
//if (indicator[e] == INT_MAX) cout << ro() << rank << ": " << ElementTypeName(e.GetElementType()) << e.GlobalID() << endl;
//assert(indicator[e] != INT_MAX);
}
}
//2.Communicate indicator on faces and edges
ReduceData(indicator,FACE|EDGE,0,ReduceMin);
ExchangeData(indicator,FACE|EDGE,0);
//3.If there is adjacent finer edge that are not marked for coarsening
// then this cell should not be coarsened
unscheduled = scheduled = 0;
......@@ -734,6 +738,7 @@ namespace INMOST
}
//c) data reduction to get maximum over mesh partition
ReduceData(coarse_indicator,EDGE,0,ReduceMax);
ExchangeData(coarse_indicator,EDGE,0);
//d) look from cells if any edge is coarsened earlier
for(Storage::integer it = 0; it < CellLastLocalID(); ++it) if( isValidCell(it) )
{
......@@ -752,6 +757,7 @@ namespace INMOST
}
}
}
ExchangeData(indicator,CELL|FACE|EDGE,0);
//5.Go back to 1 until no new elements scheduled
scheduled = Integrate(scheduled);
unscheduled = Integrate(unscheduled);
......@@ -895,6 +901,20 @@ namespace INMOST
assert(visited);
}
}
/*
for(Storage::integer it = 0; it < NodeLastLocalID(); ++it) if( isValidNode(it) )
{
Node e = NodeByLocalID(it);
if( !e.Hidden() )
{
int my_level = -1;
ElementArray<Edge> edges = e.getEdges();
for(ElementArray<Edge>::iterator kt = edges.begin(); kt != edges.end(); ++kt)
my_level = std::max(my_level,level[*kt]);
level[e] = my_level;
}
}
*/
//jump to later schedule
schedule_counter--;
}
......@@ -902,13 +922,14 @@ namespace INMOST
DeleteTag(indicator,FACE|EDGE);
//todo:
ResolveShared(true);
ResolveModification();
//ResolveModification(false,1);
//todo:
//let the user update their data
ApplyModification();
//done
EndModification();
ExchangeData(hanging_nodes,CELL | FACE,0);
CheckCentroids();
//ExchangeData(hanging_nodes,CELL | FACE,0);
//cleanup null links to hanging nodes
for(ElementType etype = FACE; etype <= CELL; etype = NextElementType(etype))
{
......
......@@ -57,6 +57,7 @@ int main(int argc, char ** argv)
q += (xyz[d]-cnt[d])*(xyz[d]-cnt[d]);
q = sqrt(q);
if( q < r*(k+1) && q > r*k)
// if( q < 0.4 && q > 0.3)
//if( q < 0.02 )
{
indicator[it->self()] = 1;
......@@ -91,6 +92,7 @@ int main(int argc, char ** argv)
}
while(numref);
refcnt = 0;
//if (0)
do
{
numref = 0;
......@@ -111,7 +113,7 @@ int main(int argc, char ** argv)
numref = m.Integrate(numref);
if( numref )
{
std::cout << "k " << k << " crscnt " << refcnt << " " << r*k << " < r < " << r*(k+1) << " numcrs " << numref << " cells " << m.NumberOfCells() << std::endl;
std::cout << ": k " << k << " crscnt " << refcnt << " " << r*k << " < r < " << r*(k+1) << " numcrs " << numref << " cells " << m.NumberOfCells() << std::endl;
/*
{
std::stringstream file;
......@@ -119,6 +121,8 @@ int main(int argc, char ** argv)
m.Save(file.str());
}
*/
if( !m.Coarse(indicator) ) break;
{
......@@ -127,10 +131,25 @@ int main(int argc, char ** argv)
m.Save(file.str());
}
std::cout << " Output " << k << std::endl;
m.UpdateStatus();
TagInteger isface = m.CreateTag("isface",DATA_INTEGER,CELL|FACE,NONE,1);
for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) isface[*it] = 1;
std::stringstream file;
//m.SetFileOption("VTK_OUTPUT_FACES","1");
file << "step_" << k << "_" << m.GetProcessorRank() << ".xml";
//file << "step_" << k << "_" << refcnt << ".pvtk";
m.Save(file.str());
if (k == 6 && refcnt==0) exit(-1);
//cleanup indicator
for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) indicator[it->self()] = 0;
}
{
}
//return -1;
refcnt++;
......@@ -138,7 +157,13 @@ int main(int argc, char ** argv)
while(numref);
{
std::cout << " Output " << k << std::endl;
m.UpdateStatus();
TagInteger isface = m.CreateTag("isface",DATA_INTEGER,CELL|FACE,NONE,1);
for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) isface[*it] = 1;
std::stringstream file;
//m.SetFileOption("VTK_OUTPUT_FACES","1");
//file << "step_" << k << "_" << m.GetProcessorRank() << ".xml";
file << "step_" << k << ".pvtk";
m.Save(file.str());
}
......
......@@ -2113,6 +2113,12 @@ int main(int argc, char ** argv)
if( argc > 2 ) mesh->SetFileOption("VTK_GRID_DIMS",argv[2]);
mesh->Load(argv[1]);
}
TagInteger global_id = mesh->GetTag("GLOBAL_ID");
for(Mesh::iteratorElement e = mesh->BeginElement(FACE); e != mesh->EndElement(); ++e)
{
if( global_id[*e] == 201 || global_id[*e] == 203 || global_id[*e] == 208 ) std::cout << e->LocalID() << " " << global_id[*e] << std::endl;
}
/*
catch(...)
{
......
......@@ -54,7 +54,7 @@
// output xml files for debugging of parallel algorithms
// search for style.xsl within examples for comfortable
// view of generated xml files
//#define USE_PARALLEL_WRITE_TIME
#define USE_PARALLEL_WRITE_TIME
// this will revert Mesh::PrepareReceiveInner to always
// use MPI point to point functionality disregarding problem type
......
......@@ -3163,11 +3163,12 @@ namespace INMOST
/// newly created elements, provide them valid global identificators, resolve owners of
/// the elements potentially optimized using information from BridgeTag and LayersTag
/// May use ResolveShared function as basis but instead the whole mesh run the same algorithm for subset.
void ResolveModification(int metric=0); //resolve parallel state of newly created elements, restore ghost layers; not implemented, resuse ResolveShared code
void ResolveModification(bool only_new, int metric=0); //resolve parallel state of newly created elements, restore ghost layers; not implemented, resuse ResolveShared code
void EndModification (); //delete hidden elements
enumerator getNext (const HandleType * arr, enumerator size, enumerator k, MarkerType marker) const;
enumerator Count (const HandleType * arr, enumerator size, MarkerType marker) const;
void FindMinDijkstra (Cell c, double& dist, int& owner);
void CheckCentroids();
//implemented in mesh.cpp
private:
Tag tag_topologyerror;
......
......@@ -1658,7 +1658,8 @@ namespace INMOST
assert(oc.GetHandle() != *it);
if (tag_mark[oc] == 1) continue;
if (GetMarker(oc.GetHandle(),NewMarker()) == false)
//if (GetMarker(oc.GetHandle(),NewMarker()) == false)
if (oc.GetStatus() == Element::Owned)
{
dist = tag_dist[*it] + 1;
owner = oc.IntegerDF(tag_owner);
......@@ -1686,15 +1687,23 @@ namespace INMOST
}
void Mesh::ResolveModification(int metric)
void Mesh::ResolveModification(bool only_new, int metric)
{
int rank = GetProcessorRank();
int size = GetProcessorsNumber();
if (size == 1) return;
Tag tag = CreateTag("TEMP_DISTANSE",DATA_REAL,CELL,CELL,2);
TagInteger tag_d = CreateTag("DIJ",DATA_INTEGER,CELL,NONE,1);
TagInteger tag_ow = CreateTag("DIJ_OW",DATA_INTEGER,CELL,NONE,1);
TagReal tag_dis = CreateTag("DIJ_DIS",DATA_REAL,CELL,NONE,1);
stringstream ss;
for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++) if (GetMarker(*it,NewMarker()))
for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++)
//if (GetMarker(*it,NewMarker()))
if (!only_new || it->nbAdjElements(NODE | EDGE | FACE | CELL, NewMarker()) != 0)
{
if (it->GetStatus() == Element::Owned) continue;
double min_distance = 0;
......@@ -1718,7 +1727,10 @@ namespace INMOST
}
else
{
tag_d[*it] = 1;
FindMinDijkstra(it->getAsCell(), min_distance, owner);
tag_ow[*it] = owner;
tag_dis[*it] = min_distance;
}
it->RealArray(tag)[0] = owner;
it->RealArray(tag)[1] = min_distance;
......@@ -1727,7 +1739,8 @@ namespace INMOST
ReduceData(tag, CELL, 0, OperationMinDistance);
ExchangeData(tag, CELL, 0);
for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++) if (GetMarker(*it,NewMarker()) && (it->GetStatus() == Element::Ghost || it->GetStatus() == Element::Shared) )
for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++) //if (GetMarker(*it,NewMarker()) && (it->GetStatus() == Element::Ghost || it->GetStatus() == Element::Shared) )
if ( !only_new || it->nbAdjElements(NODE | EDGE | FACE | CELL, NewMarker()) != 0)
{
if (it->GetStatus() == Element::Owned) continue;
int new_owner = (int)it->RealArray(tag)[0];
......
......@@ -750,14 +750,21 @@ namespace INMOST
void determine_my_procs_low(Mesh * m, HandleType h, dynarray<Storage::integer,64> & result, dynarray<Storage::integer,64> & intersection)
{
MarkerType hm = m->HideMarker();
Element::adj_type const & subelements = m->LowConn(h);
Element::adj_type::const_iterator i = subelements.begin();
Storage::integer_array p = m->IntegerArrayDV(*i,m->ProcessorsTag());
result.clear();
while(i != subelements.end())
{
if( hm && m->GetMarker(*i,hm) ) {i++; continue;}
Storage::integer_array p = m->IntegerArrayDV(*i,m->ProcessorsTag());
result.insert(result.end(),p.begin(),p.end());
i++;
break;
}
while(i != subelements.end())
{
if( hm && m->GetMarker(*i,hm) ) {i++; continue;}
Storage::integer_array q = m->IntegerArrayDV(*i,m->ProcessorsTag());
intersection.resize(std::max(static_cast<unsigned>(result.size()),static_cast<unsigned>(q.size())));
dynarray<Storage::integer,64>::iterator qt = std::set_intersection(result.begin(),result.end(),q.begin(),q.end(),intersection.begin());
......@@ -769,12 +776,14 @@ namespace INMOST
void determine_my_procs_high(Mesh * m, HandleType h, const Tag & procs, dynarray<Storage::integer,64> & result, dynarray<Storage::integer,64> & intersection)
{
MarkerType hm = m->HideMarker();
Element::adj_type const & overelements = m->HighConn(h);
if( overelements.empty() ) return;
Element::adj_type::const_iterator i = overelements.begin();
result.clear();
while(i != overelements.end())
{
if( hm && m->GetMarker(*i,hm) ) {i++; continue;}
Storage::integer_array q = m->IntegerArrayDV(*i,procs);
intersection.resize(result.size()+q.size());
dynarray<Storage::integer,64>::iterator qt = std::set_union(result.begin(),result.end(),q.begin(),q.end(),intersection.begin());
......@@ -1047,7 +1056,7 @@ namespace INMOST
}
else
{
/*
if (only_new)
{
for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++) if (GetMarker(*it,NewMarker()))
......@@ -1059,6 +1068,7 @@ namespace INMOST
}
}
}
*/
double time = Timer();
Storage::real epsilon = GetEpsilon();
......@@ -1081,6 +1091,7 @@ namespace INMOST
for(integer nit = 0; nit < NodeLastLocalID(); ++nit) if( isValidNode(nit) )
{
Node it = NodeByLocalID(nit);
if (only_new && GetMarker(it->GetHandle(),NewMarker()) == false) continue;
Storage::integer_array arr = it->IntegerArrayDV(tag_processors);
arr.resize(1);
arr[0] = mpirank;
......@@ -1100,6 +1111,7 @@ namespace INMOST
for(iteratorNode n = BeginNode(); n != EndNode(); n++)
{
if (only_new && GetMarker(*n,NewMarker()) == false) continue;
real_array c = n->Coords();
for(real_array::size_type k = 0; k < procs.size(); k++)
if( point_in_bbox(c.data(),bboxs.data()+procs[k]*dim*2,dim,GetEpsilon()) )
......@@ -1316,7 +1328,12 @@ namespace INMOST
REPORT_VAL("time",time);
MarkerType new_lc;
if( only_new ) new_lc = CreateMarker();
TagInteger lc_mrk = CreateTag("LC_MARKER",DATA_INTEGER,CELL|FACE|EDGE,NONE,1);
MarkerType hm = HideMarker();
for (iteratorElement it = BeginElement(CELL|FACE|EDGE); it != EndElement(); it++) lc_mrk[*it] = 0;
for(ElementType current_mask = EDGE; current_mask <= CELL; current_mask = NextElementType(current_mask))
{
......@@ -1330,6 +1347,32 @@ namespace INMOST
Element::Status estat;
time = Timer();
if(only_new)
{
for(integer eit = 0; eit < LastLocalID(current_mask); ++eit) if( isValidElement(current_mask,eit) )
{
Element it = ElementByLocalID(current_mask,eit);
if( it->GetMarker(hm) ) continue;
if (GetMarker(it.GetHandle(),NewMarker()) == true)
{
SetMarker(ComposeHandle(current_mask,eit),new_lc);
lc_mrk[ComposeHandle(current_mask,eit)] = 1;
continue;
}
Element::adj_type low_conn = LowConn(ComposeHandle(current_mask,eit));
for (Element::adj_type::iterator it = low_conn.begin(); it != low_conn.end(); it++)
if (GetMarker(*it,NewMarker()))
{
SetMarker(ComposeHandle(current_mask,eit),new_lc);
lc_mrk[ComposeHandle(current_mask,eit)] = 1;
break;
}
}
}
stringstream ss;
ss << "file_" << GetProcessorRank() << ".txt";
ofstream ofs(ss.str().c_str());
//Determine what processors potentially share the element
#if defined(USE_OMP)
#pragma omp parallel
......@@ -1344,10 +1387,37 @@ namespace INMOST
if( isValidElement(current_mask,eit) )
{
Element it = ElementByLocalID(current_mask,eit);
if (only_new && GetMarker(it.GetHandle(),NewMarker()) == false) continue;
if( it->GetMarker(hm) ) continue;
if (only_new && GetMarker(it.GetHandle(),new_lc) == false) continue;
determine_my_procs_low(this,it->GetHandle(), result, intersection);
Storage::integer_array p = it->IntegerArrayDV(tag_processors);
//if (current_mask == CELL)
//if (GetProcessorRank() == 0 && it->LocalID() == 59 ||GetProcessorRank() == 1 && it->LocalID() == 49)
{
Element::adj_type const & subelements = LowConn(it->GetHandle());
Element::adj_type::const_iterator i = subelements.begin();
double cnt[3];
it->Centroid(cnt);
ofs << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << " " << cnt[0] << " " << cnt[1] << " " << cnt[2] << std::endl;
for (int k = 0; k < result.size(); k++) ofs << result[k] << " ";
ofs << std::endl;
while( i != subelements.end() )
{
Storage::integer_array p = IntegerArrayDV(*i,ProcessorsTag());
ofs << "| face:" << GetHandleID(*i) << " " << Element::StatusName(Bulk(*i,SharedTag())) << " ";
for (int k = 0; k < p.size(); k++) ofs << p[k] << " ";
ofs << std::endl;
i++;
}
//out << ro() << GetProcessorRank() << "!!!!! " << ss.str() << endl;
}
if( result.empty() )
{
p.clear();
......@@ -1414,18 +1484,23 @@ namespace INMOST
message_send.push_back(0);
for(Mesh::iteratorElement it = BeginElement(current_mask); it != EndElement(); it++)
{
if (only_new && GetMarker(*it,NewMarker()) == false) continue;
if (only_new && GetMarker(*it,new_lc) == false) continue;
Storage::integer_array pr = it->IntegerArrayDV(tag_processors);
if( std::binary_search(pr.begin(),pr.end(),*p) )
{
Element::adj_type & sub = LowConn(*it);
if( sub.size() == 0 ) throw Impossible;
message_send.push_back(static_cast<int>(sub.size()));
int message_size_pos = message_send.size();
message_send.push_back(0);
REPORT_VAL("number of connections",sub.size());
for(Element::adj_type::iterator kt = sub.begin(); kt != sub.end(); kt++)
REPORT_STR("element " << ElementTypeName(current_mask) << ":" << it->LocalID());
for(Element::adj_type::iterator kt = sub.begin(); kt != sub.end(); kt++) if( !hm || !GetMarker(*kt,hm) )
{
message_send.push_back(GlobalID(*kt));
REPORT_STR("global id " << GlobalID(*kt) << " local id " << GetHandleID(*kt) << " " << Element::StatusName(Element(this,*kt)->GetStatus()));
message_send[message_size_pos]++;
double cnt[3];
ElementByLocalID(PrevElementType(current_mask),GetHandleID(*kt))->Centroid(cnt);
REPORT_STR("global id " << GlobalID(*kt) << " local id " << GetHandleID(*kt) << " " << Element::StatusName(Element(this,*kt)->GetStatus()) << " cnt " << cnt[0] << " " << cnt[1] << " " << cnt[2]);
}
message_send[1]++;
elements[m].push_back(*it);
......@@ -1512,6 +1587,7 @@ namespace INMOST
}
int find_local_id = mapping[find].second;
REPORT_STR("global id " << global_id << " local id " << find_local_id << " " << Element::StatusName(ElementByLocalID(PrevElementType(current_mask), find_local_id)->GetStatus()));
if( GetMarker(ComposeHandle(PrevElementType(current_mask), find_local_id),hm) ) std::cout << "Found hidden element" << std::endl;
sub_elements.push_back(ComposeHandle(PrevElementType(current_mask), find_local_id));
}
......@@ -1568,7 +1644,7 @@ namespace INMOST
//Now mark all the processors status
for(Mesh::iteratorElement it = BeginElement(current_mask); it != EndElement(); it++)
{
if (only_new && (GetMarker(*it,NewMarker()) == false)) continue;
if (only_new && (GetMarker(*it,new_lc) == false)) continue;
Storage::integer_array pr = it->IntegerArrayDV(tag_processors);
if( pr.empty() )
......@@ -1614,6 +1690,8 @@ namespace INMOST
REPORT_STR("Set parallel info");
REPORT_VAL("time",time);
}
ReleaseMarker(new_lc,EDGE | FACE|CELL);
}
RemoveGeometricData(table);
}
......@@ -4841,6 +4919,28 @@ namespace INMOST
#endif
}
void Mesh::CheckCentroids()
{
TagRealArray checker = CreateTag("CHECK_CENTROIDS",DATA_REAL,CELL|FACE|EDGE|NODE,CELL|FACE|EDGE|NODE,3);
for(Mesh::iteratorElement it = BeginElement(CELL|FACE|EDGE|NODE); it != EndElement(); ++it) if( it->GetStatus() & (Element::Ghost | Element::Shared) )
it->Centroid(checker[*it].data());
ExchangeData(checker,CELL|FACE|EDGE|NODE,0);
for(Mesh::iteratorElement it = BeginElement(CELL|FACE|EDGE|NODE); it != EndElement(); ++it) if( it->GetStatus() & (Element::Ghost | Element::Shared) )
{
double cnt[3];
it->Centroid(cnt);
bool problem = false;
for(int k = 0; k < 3; ++k)
if( fabs(cnt[k] - checker[*it][k]) > 1.0e-9 )
problem = true;
if( problem )
{
std::cout << "bad centroid on " << GetProcessorRank() << " " << ElementTypeName(it->GetElementType()) << ":" << it->LocalID() << " " << cnt[0] << " " << cnt[1] << " " << cnt[2] << " != " << checker[*it][0] << " " << checker[*it][1] << " " << checker[*it][2] << std::endl;
}
}
}
void Mesh::SortParallelStorage(parallel_storage & ghost, parallel_storage & shared, ElementType mask)
{
ENTER_FUNC();
......
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