Skip to content

GitLab

Commit 73f0c811 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

fix Face::SplitFace to properly connect already existing faces to cells; updates to grid tools

parent 35fde864
Hide whitespace changes
Inline Side-by-side
Examples/GridTools/fix_faults.cpp
View file @ 73f0c811
......@@ -23,13 +23,13 @@ static void GetBox(Element e, float minmax[6])
{
if (minmax[1 + 2 * i] < minmax[0 + 2 * i])
{
minmax[1 + 2 * i] = 1.0e-5f;
minmax[0 + 2 * i] = -1.0e-5f;
minmax[1 + 2 * i] = 1.0e-9f;
minmax[0 + 2 * i] = -1.0e-9f;
}
else if (minmax[1 + 2 * i] == minmax[0 + 2 * i])
{
minmax[1 + 2 * i] += 1.0e-5f;
minmax[0 + 2 * i] += -1.0e-5f;
minmax[1 + 2 * i] += 1.0e-9f;
minmax[0 + 2 * i] += -1.0e-9f;
}
}
}
......@@ -59,7 +59,7 @@ static bool MatchPoints(const double v1[3], const double v2[3])
double l = 0;
for (int k = 0; k < 3; ++k)
l += (v1[k] - v2[k])*(v1[k] - v2[k]);
return sqrt(l) < 1.0e-6;
return sqrt(l) < 1.0e-9;
}
//returns distance between line if shortest line is within segments, writes position of distance into last argument
......@@ -90,7 +90,7 @@ static double SegmentDistance(const double v1[3], const double v2[3], const doub
parallel /= sqrt(d2121)*sqrt(d4343);
if (fabs(parallel - 1.0) > 1.0e-6)
*/
if (fabs(d2121*d4343 - d4321*d4321) > 1.0e-6)
if (fabs(d2121*d4343 - d4321*d4321) > 1.0e-9)
{
mu1 = (d1343*d4321 - d1321*d4343) / (d2121*d4343 - d4321*d4321);
mu2 = (d1343 + mu1*d4321) / d4343;
......@@ -278,7 +278,7 @@ class kdtree
Storage::real fnrm[3], ffnrm[3];
f.UnitNormal(fnrm);
ff.UnitNormal(ffnrm);
if (fabs(fnrm[0] * ffnrm[0] + fnrm[1] * ffnrm[1] + fnrm[2] * ffnrm[2]) >= 1.0 - 1.0e-3)
if (fabs(fnrm[0] * ffnrm[0] + fnrm[1] * ffnrm[1] + fnrm[2] * ffnrm[2]) >= 0.85)
{
//todo: check that faces intersect
int fid = f.LocalID(), ffid = ff.LocalID();
......@@ -303,7 +303,7 @@ class kdtree
}
l1 = sqrt(l1);
l2 = sqrt(l2);
if (h < (l1 + l2)*1.0e-13)
if (h < (l1 + l2)*1.0e-9)
intersect = true;
}
}
......@@ -496,6 +496,9 @@ void FixFaults::FixMeshFaults(MarkerType mrk)
std::cout << "Time to fix normals: " << Timer() - tt << std::endl;
std::cout << "Total face normals fixed: " << fixed << std::endl;
if( !Element::CheckConnectivity(&m) )
std::cout << "Connectivity is broken on entry to " << __FUNCTION__ << std::endl;
kdtree t(&m,mrk);
......@@ -517,9 +520,9 @@ void FixFaults::FixMeshFaults(MarkerType mrk)
int fid = f.LocalID();
ElementArray<Face> ifaces(&m);
t.intersect_nonadj_face(f, ifaces);
std::cout << "face " << f.LocalID() << " finds: ";
for(int kk = 0; kk < ifaces.size(); ++kk) std::cout << ifaces[kk].LocalID() << " ";
std::cout << std::endl;
//std::cout << "face " << f.LocalID() << " finds: ";
//for(int kk = 0; kk < ifaces.size(); ++kk) std::cout << ifaces[kk].LocalID() << " ";
//std::cout << std::endl;
if (!ifaces.empty())
{
marked++;
......@@ -584,12 +587,12 @@ void FixFaults::FixMeshFaults(MarkerType mrk)
parallel /= sqrt(d2121)*sqrt(d4343);
if( fabs(parallel - 1.0) > 1.0e-6 )
*/
if (fabs(d2121*d4343 - d4321*d4321) > 1.0e-6)
if (fabs(d2121*d4343 - d4321*d4321) > 1.0e-9)
{
mu1 = (d1343*d4321 - d1321*d4343) / (d2121*d4343 - d4321*d4321);
mu2 = (d1343 + mu1*d4321) / d4343;
//here only one or zero intersections, but we have to introduce an edge
if (mu1 >= 0 - 1.0e-13 && mu1 <= 1 + 1.0e-13 && mu2 >= 0 - 1.0e-13 && mu2 <= 1 + 1.0e-13)
if (mu1 >= 0 - 1.0e-9 && mu1 <= 1 + 1.0e-9 && mu2 >= 0 - 1.0e-9 && mu2 <= 1 + 1.0e-9)
{
//find distance and point
double vs1[3], vs2[3], h = 0, vv[3];
......@@ -600,7 +603,7 @@ void FixFaults::FixMeshFaults(MarkerType mrk)
vv[k] = (vs1[k] + vs2[k])*0.5;
h += (vs2[k] - vs1[k])*(vs2[k] - vs1[k]);
}
if (h < (l1 + l2)*1.0e-13)
if (h < (l1 + l2)*1.0e-9)
{
Node matchnode = InvalidNode();
if (MatchPoints(vv, v1))
......@@ -738,7 +741,7 @@ void FixFaults::FixMeshFaults(MarkerType mrk)
double vs2 = v3[k] + mu2*(v4[k] - v3[k]);
h += (vs2 - vs1)*(vs2 - vs1);
}
if (h < (l1 + l2)*1.0e-13)
if (h < (l1 + l2)*1.0e-9)
{
Storage::reference_array nodes = new_points[it->self()];
for (int k = 0; k < nt; ++k)
......@@ -900,6 +903,9 @@ void FixFaults::FixMeshFaults(MarkerType mrk)
Storage::real_array c = it->second->Coords();
std::cout << it->second->LocalID() << " " << it->first.v[0] << "," << it->first.v[1] << "," << it->first.v[2] << " " << c[0] << "," << c[1] << "," << c[2] << std::endl;
}
if( !Element::CheckConnectivity(&m) )
std::cout << "Connectivity is broken before modification in " << __FUNCTION__ << std::endl;
tt = Timer();
TagInteger mark = m.CreateTag("SPLITTED", DATA_INTEGER, FACE | EDGE, NONE, 1);
......@@ -937,16 +943,42 @@ void FixFaults::FixMeshFaults(MarkerType mrk)
Storage::reference_array edges = new_edges[it->self()];
if (!edges.empty())
{
ElementArray<Cell> adj_cells = it->getCells();
ElementArray<Face> new_faces = Face::SplitFace(it->self(), ElementArray<Edge>(&m, edges.begin(), edges.end()), 0);
if( new_faces.size() == 1 || new_faces.size() == 0 ) std::cout << "split " << it->LocalID() << " gets " << (new_faces.empty()?-1:new_faces[0].LocalID()) << " edges " << edges.size() << std::endl;
//std::cout << "split " << it->LocalID();
//std::cout << " adj cells: ";
//for (int k = 0; k < (int)adj_cells.size(); ++k) std::cout << adj_cells[k].LocalID() << " ";
//std::cout << " new faces: ";
for (int k = 0; k < (int)new_faces.size(); ++k)
{
mark[new_faces[k]] = 1;
//std::cout << new_faces[k].LocalID() << " (";
//ElementArray<Cell> adj_cells2 = new_faces[k].getCells();
//for(int q = 0; q < (int)adj_cells2.size(); ++q)
// std::cout << adj_cells2[q].LocalID() << " ";
//std::cout << ") ";
}
//std::cout << std::endl;
nfaces++;
}
}
std::cout << "split faces " << nedges << std::endl;
std::cout << "split faces " << nfaces << std::endl;
m.ApplyModification();
m.EndModification();
/*
for (Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) if( !it->CheckElementConnectivity() )
{
std::cout << "connectivity problem for face " << it->LocalID() << std::endl;
it->PrintElementConnectivity();
}
for (Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( !it->CheckElementConnectivity() )
{
std::cout << "connectivity problem for cell " << it->LocalID() << std::endl;
it->PrintElementConnectivity();
}
*/
m.DeleteTag(new_points);
m.DeleteTag(new_edges);
......@@ -957,4 +989,13 @@ void FixFaults::FixMeshFaults(MarkerType mrk)
std::cout << "Cells: " << m.NumberOfCells() << std::endl;
std::cout << "Faces: " << m.NumberOfFaces() << std::endl;
std::cout << "Edges: " << m.NumberOfEdges() << std::endl;
m.ReorderEmpty(CELL|FACE|EDGE|NODE);
if( !Element::CheckConnectivity(&m) )
{
std::cout << "Connectivity is broken on end of " << __FUNCTION__ << std::endl;
m.Save("problem.pmf");
m.Save("problem.xml");
exit(-1);
}
}
Examples/GridTools/fracture.cpp
View file @ 73f0c811
......@@ -79,6 +79,10 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
std::cout << "Faces: " << m->NumberOfFaces() << std::endl;
std::cout << "Edges: " << m->NumberOfEdges() << std::endl;
std::cout << "Nodes: " << m->NumberOfNodes() << std::endl;
if( !Element::CheckConnectivity(m) )
std::cout << "Connectivity is broken on entry to " << __FUNCTION__ << std::endl;
m->BeginModification();
fracture_marker = m->CreateMarker();
std::cout << "create marker fracture_marker " << fracture_marker << std::endl;
......@@ -1199,6 +1203,9 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
}
m->ApplyModification();
m->EndModification();
if( !Element::CheckConnectivity(m) )
std::cout << "Connectivity is broken on exit from " << __FUNCTION__ << std::endl;
std::cout << "Cells: " << m->NumberOfCells() << std::endl;
std::cout << "Faces: " << m->NumberOfFaces() << std::endl;
std::cout << "Edges: " << m->NumberOfEdges() << std::endl;
......
Examples/GridTools/test_make_faults.cpp
View file @ 73f0c811
......@@ -16,13 +16,13 @@ public:
v[2] = 1;
double l = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
v[0] /= l; v[1] /= l; v[2] /= l;
double dx = (cmax[0]-cmin[0])/25.0;
double dy = (cmax[1]-cmin[1])/25.0;
double dx = (cmax[0]-cmin[0])/50.0;
double dy = (cmax[1]-cmin[1])/50.0;
double dl = sqrt(dx*dx+dy*dy);
double t = 2.0*(rand()*1.0/RAND_MAX)-1.0;
double a = M_PI*0.5 + t*M_PI*0.2;
double a = M_PI*0.5 + t*M_PI*0.15;
double y = cmin[1];
double x = (0.1 + (rand()*1.0/RAND_MAX)*0.8)*(cmax[0]-cmin[0]) + cmin[0];
double x = (0.2 + (rand()*1.0/RAND_MAX)*0.6)*(cmax[0]-cmin[0]) + cmin[0];
std::cout << " starting point " << x << ", " << y << std::endl;
curvxy.push_back(std::make_pair(x,y));
bool ended = false;
......@@ -36,7 +36,7 @@ public:
if( y <= cmin[1] ) {y = cmin[1]; ended = true;}
curvxy.push_back(std::make_pair(x,y));
t = 2.0*(rand()*1.0/RAND_MAX)-1.0;
a = a + M_PI*0.1*t;
a = a + M_PI*0.05*t;
}
std::cout << " ending point " << x << ", " << y << std::endl;
std::cout << " generated points " << curvxy.size() << std::endl;
......@@ -77,13 +77,13 @@ int main(int argc, char ** argv)
{
if (argc < 2)
{
std::cout << "Usage: " << argv[0] << " mesh [mesh_out=grid.vtk]" << std::endl;
std::cout << "Usage: " << argv[0] << " mesh [mesh_out=grid.vtk] [nfaults=1]" << std::endl;
return -1;
}
std::string grid_out = "grid.vtk";
if (argc > 2) grid_out = std::string(argv[2]);
int nfaults = argc > 3 ? atoi(argv[3]) : 1;
{
Mesh m;
......@@ -104,57 +104,67 @@ int main(int argc, char ** argv)
if( cmax[d] < n->Coords()[d] ) cmax[d] = n->Coords()[d];
}
}
SliceFault sf(cmin,cmax);
sf.SliceMesh(m,false);
//create fake aperture tag
TagReal aperture = m.CreateTag("APERTURE",DATA_REAL,FACE,FACE,1);
TagBulk sliced = m.GetTag("SLICED"); //this was created by Slice
TagInteger mat = m.GetTag("MATERIAL");
int nmrk = 0;
for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) if( it->HaveData(sliced) ) {aperture[*it] = 1; nmrk++;}
std::cout << "marked for fracture " << nmrk << std::endl;
Fracture fr(m);
//fr.Open(aperture,false,0.5);
fr.Open(aperture,false,1.0);
//shift nodes
double v[3];
sf.Getv(v);
double t = 0.1+0.9*rand()*1.0/RAND_MAX;
double sc = t*(cmax[2]-cmin[2])*0.4;
std::cout << "shift nodes by " << sc << std::endl;
for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
for(int kk = 0; kk < nfaults; ++kk)
{
ElementArray<Cell> adj_cells = n->getCells();
int mats[3] = {0,0,0};
for(ElementArray<Cell>::iterator kt = adj_cells.begin(); kt != adj_cells.end(); ++kt)
mats[mat[*kt]]++;
if( mats[1] == 0 )
SliceFault sf(cmin,cmax);
sf.SliceMesh(m,false);
//create fake aperture tag
TagReal aperture = m.CreateTag("APERTURE",DATA_REAL,FACE,FACE,1);
TagBulk sliced = m.GetTag("SLICED"); //this was created by Slice
TagInteger mat = m.GetTag("MATERIAL");
int nmrk = 0;
for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) if( it->HaveData(sliced) ) {aperture[*it] = 1; nmrk++;}
std::cout << "marked for fracture " << nmrk << std::endl;
Fracture fr(m);
//fr.Open(aperture,false,0.5);
fr.Open(aperture,false,1.0);
//shift nodes
double v[3];
sf.Getv(v);
double t = 0.1+0.9*rand()*1.0/RAND_MAX;
double sc = t*(cmax[2]-cmin[2])*0.05;
std::cout << "shift nodes by " << sc << std::endl;
for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n)
{
n->Coords()[0] += v[0]*sc;
n->Coords()[1] += v[1]*sc;
n->Coords()[2] += v[2]*sc;
ElementArray<Cell> adj_cells = n->getCells();
int mats[3] = {0,0,0};
for(ElementArray<Cell>::iterator kt = adj_cells.begin(); kt != adj_cells.end(); ++kt)
mats[mat[*kt]]++;
if( mats[1] == 0 )
{
n->Coords()[0] += v[0]*sc;
n->Coords()[1] += v[1]*sc;
n->Coords()[2] += v[2]*sc;
}
}
FixFaults fix(m);
MarkerType fault = m.CreateMarker();
std::cout << "fault marker: " << fault << std::endl;
for(Mesh::iteratorFace f = m.BeginFace(); f != m.EndFace(); ++f)
if( mat[*f] == 2 ) f->SetMarker(fault);
fix.FixMeshFaults(fault);
std::cout << "release fault marker: " << fault << std::endl;
m.ReleaseMarker(fault,FACE);
//delete to clear data
m.DeleteTag(aperture);
m.DeleteTag(sliced);
m.DeleteTag(mat);
}
FixFaults fix(m);
MarkerType fault = m.CreateMarker();
for(Mesh::iteratorFace f = m.BeginFace(); f != m.EndFace(); ++f)
if( mat[*f] == 2 ) f->SetMarker(fault);
fix.FixMeshFaults(fault);
m.ReleaseMarker(fault,FACE);
m.Save(grid_out);
}
return 0;
......
Source/Mesh/modify.cpp
View file @ 73f0c811
......@@ -1464,9 +1464,31 @@ namespace INMOST
{
if( loop.size() > 2 )
{
ret.push_back(m->CreateFace(loop).first);
std::pair<Face,bool> ff = m->CreateFace(loop);
ret.push_back(ff.first); //collect new faces
adj_type & hc = m->HighConn(ret.back()->GetHandle());
hc.replace(hc.begin(),hc.end(),cells.begin(),cells.end());
if( ff.second ) //connect brand new face to cells
{
//std::cout << "hello1 from " << ff.first.LocalID() << std::endl;
hc.replace(hc.begin(),hc.end(),cells.begin(),cells.end());
}
else //face already existed before and may be connected to some cells
{
//std::cout << "hello2 from " << ff.first.LocalID() << std::endl;
//std::cout << "cells: ";
//for(int k = 0; k < (int)cells.size(); ++k) std::cout << GetHandleID(cells[k]) << " ";
//std::cout << std::endl;
//std::cout << "conns: ";
//for(int k = 0; k < (int)hc.size(); ++k) std::cout << GetHandleID(hc[k]) << " ";
//std::cout << std::endl;
for(int k = 0; k < (int)cells.size(); ++k)
{
bool add = true;
for(int j = 0; j < (int)hc.size() && add; ++j)
if( hc[j] == cells[k] ) add = false;
if( add ) hc.push_back(cells[k]);
} //FIXME: what if more then 2 connections? have to fire topology exception, unless there are hidden cells
}
}
}
else break;
......
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