Commit 35fde864 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Some fixes in GridTools

parent cc18b29b
......@@ -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-7;
return sqrt(l) < 1.0e-6;
}
//returns distance between line if shortest line is within segments, writes position of distance into last argument
......@@ -316,31 +316,60 @@ class kdtree
kdtree() : set(NULL), size(0), children(NULL) {}
void clear_children() { if (children) { children[0].clear_children(); children[1].clear_children(); free(children); } }
public:
kdtree(Mesh * m) : m(m), children(NULL)
kdtree(Mesh * m, MarkerType mrk = 0) : m(m), children(NULL)
{
double tt;
size = m->NumberOfFaces();
assert(size > 1);
set = new entry[size];
INMOST_DATA_ENUM_TYPE k = 0;
tt = Timer();
printf("Prepearing face set.\n");
m->ReorderEmpty(FACE);
if( mrk == 0 )
{
size = m->NumberOfFaces();
assert(size > 1);
set = new entry[size];
INMOST_DATA_ENUM_TYPE k = 0;
tt = Timer();
printf("Prepearing face set.\n");
m->ReorderEmpty(FACE);
#if defined(USE_OMP)
#pragma omp parallel for
#endif
for (int i = 0; i < m->FaceLastLocalID(); ++i) if (m->isValidFace(i))
for (int i = 0; i < m->FaceLastLocalID(); ++i) if (m->isValidFace(i))
{
Face it = m->FaceByLocalID(i);
Storage::real cnt[3];
it->Centroid(cnt);
int k = it->DataLocalID();
set[k].e = it->GetHandle();
set[k].xyz[0] = static_cast<float>(cnt[0]);
set[k].xyz[1] = static_cast<float>(cnt[1]);
set[k].xyz[2] = static_cast<float>(cnt[2]);
}
printf("Done. Time %lg\n", Timer() - tt);
}
else
{
Face it = m->FaceByLocalID(i);
Storage::real cnt[3];
it->Centroid(cnt);
int k = it->DataLocalID();
set[k].e = it->GetHandle();
set[k].xyz[0] = static_cast<float>(cnt[0]);
set[k].xyz[1] = static_cast<float>(cnt[1]);
set[k].xyz[2] = static_cast<float>(cnt[2]);
size = 0;
for (int i = 0; i < m->FaceLastLocalID(); ++i)
if( m->isValidFace(i) && m->FaceByLocalID(i).GetMarker(mrk) ) size++;
assert(size > 1);
set = new entry[size];
INMOST_DATA_ENUM_TYPE k = 0;
tt = Timer();
printf("Prepearing face set.\n");
for (int i = 0; i < m->FaceLastLocalID(); ++i) if (m->isValidFace(i))
{
Face it = m->FaceByLocalID(i);
if( it->GetMarker(mrk) )
{
Storage::real cnt[3];
it->Centroid(cnt);
set[k].e = it->GetHandle();
set[k].xyz[0] = static_cast<float>(cnt[0]);
set[k].xyz[1] = static_cast<float>(cnt[1]);
set[k].xyz[2] = static_cast<float>(cnt[2]);
k++;
}
}
printf("Done. Time %lg\n", Timer() - tt);
}
printf("Done. Time %lg\n", Timer() - tt);
int done = 0, total = size;
printf("Building KD-tree.\n");
tt = Timer();
......@@ -448,7 +477,7 @@ struct coords
void FixFaults::FixMeshFaults()
void FixFaults::FixMeshFaults(MarkerType mrk)
{
std::cout << "Start:" << std::endl;
......@@ -468,7 +497,7 @@ void FixFaults::FixMeshFaults()
std::cout << "Total face normals fixed: " << fixed << std::endl;
kdtree t(&m);
kdtree t(&m,mrk);
TagReferenceArray new_points = m.CreateTag("EDGEPOINTS", DATA_REFERENCE, EDGE, NONE);
TagReferenceArray new_edges = m.CreateTag("FACEEDGES", DATA_REFERENCE, FACE, NONE);
......@@ -484,9 +513,13 @@ void FixFaults::FixMeshFaults()
for (int it = 0; it < m.FaceLastLocalID(); ++it) if (m.isValidFace(it))
{
Face f = m.FaceByLocalID(it);
if( !f.GetMarker(mrk) ) continue;
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;
if (!ifaces.empty())
{
marked++;
......@@ -873,6 +906,7 @@ void FixFaults::FixMeshFaults()
//start splitting
m.BeginModification();
//split all edges
int nedges = 0;
for (Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it) if (!it->New())
{
Storage::reference_array nodes = new_points[it->self()];
......@@ -892,20 +926,25 @@ void FixFaults::FixMeshFaults()
ElementArray<Edge> new_edges = Edge::SplitEdge(it->self(), split_nodes, 0);
for (int k = 0; k < (int)new_edges.size(); ++k)
mark[new_edges[k]] = 1;
nedges++;
}
}
std::cout << "split edges " << nedges << std::endl;
//split all faces
int nfaces = 0;
for (Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) if (!it->New())
{
Storage::reference_array edges = new_edges[it->self()];
if (!edges.empty())
{
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;
for (int k = 0; k < (int)new_faces.size(); ++k)
mark[new_faces[k]] = 1;
nfaces++;
}
}
std::cout << "split faces " << nedges << std::endl;
m.EndModification();
m.DeleteTag(new_points);
......
......@@ -14,7 +14,7 @@ public:
FixFaults & operator =(FixFaults const & b) {m = b.m; return *this;}
~FixFaults() {}
void FixMeshFaults();
void FixMeshFaults(MarkerType mrk = 0);
};
#endif //_FIX_FAULTS_H
......@@ -75,8 +75,17 @@ void Fracture::FaceCenter(Face f, INMOST_DATA_REAL_TYPE cnt[3]) const
void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
{
fracture_aperture = aperture;
std::cout << "Cells: " << m->NumberOfCells() << std::endl;
std::cout << "Faces: " << m->NumberOfFaces() << std::endl;
std::cout << "Edges: " << m->NumberOfEdges() << std::endl;
std::cout << "Nodes: " << m->NumberOfNodes() << std::endl;
m->BeginModification();
fracture_marker = m->CreateMarker();
std::cout << "create marker fracture_marker " << fracture_marker << std::endl;
for(Mesh::iteratorFace f = m->BeginFace(); f != m->EndFace(); ++f)
if( f->GetMarker(fracture_marker) ) std::cout << "Face " << f->LocalID() << " already marked fracture " << fracture_marker << std::endl;
for(Mesh::iteratorCell f = m->BeginCell(); f != m->EndCell(); ++f)
if( f->GetMarker(fracture_marker) ) std::cout << "Cell " << f->LocalID() << " already marked fracture " << fracture_marker << std::endl;
for(Mesh::iteratorFace f = m->BeginFace(); f != m->EndFace(); ++f)
if( f->HaveData(aperture) ) f->SetMarker(fracture_marker);
m->self()->Integer(m->CreateTag("FRACTURE_MARKER",DATA_INTEGER,MESH,NONE,1)) = fracture_marker;
......@@ -92,31 +101,41 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
Tag origcoords = m->CreateTag("ORIGINAL_FRACTURE_COORDS",DATA_REAL,NODE,NODE,3);
//std::cout << "Initial number of nodes: " << m->NumberOfNodes() << std::endl;
std::cout << "mult is " << mult << std::endl;
//Unmark any boundary faces
int unmrk = 0, mrk = 0;
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it)
{
if( it->GetMarker(isFracture()) ) mrk++;
if( !it->FrontCell().isValid() && it->GetMarker(isFracture()) )
{
it->RemMarker(isFracture());
unmrk++;
}
}
std::cout << "fracture marked " << mrk << " boundary unmarked " << unmrk << std::endl;
std::vector<Tag> transfer_node_real_tags;
//std::vector<Tag> transfer_node_integer_tags;
std::vector<Tag> transfer_node_integer_tags;
for(Mesh::iteratorTag t = m->BeginTag(); t != m->EndTag(); ++t)
if( t->isDefined(NODE) && *t != m->CoordsTag() )
{
if( t->GetDataType() == DATA_REAL )
transfer_node_real_tags.push_back(*t);
//else if( t->GetDataType() == DATA_INTEGER )
// transfer_node_integer_tags.push_back(*t);
else if( t->GetDataType() == DATA_INTEGER )
transfer_node_integer_tags.push_back(*t);
}
//std::cout << "Transfer real tags:" << std::endl;
//for(int q = 0; q < (int)transfer_node_real_tags.size(); ++q)
// std::cout << transfer_node_real_tags[q].GetTagName() << std::endl;
//std::cout << "Transfer integer tags:" << std::endl;
//for(int q = 0; q < (int)transfer_node_integer_tags.size(); ++q)
// std::cout << transfer_node_integer_tags[q].GetTagName() << std::endl;
std::cout << "Transfer real tags:" << std::endl;
for(int q = 0; q < (int)transfer_node_real_tags.size(); ++q)
std::cout << transfer_node_real_tags[q].GetTagName() << std::endl;
std::cout << "Transfer integer tags:" << std::endl;
for(int q = 0; q < (int)transfer_node_integer_tags.size(); ++q)
std::cout << transfer_node_integer_tags[q].GetTagName() << std::endl;
//For each node create several new ones that will be used in open fracture
std::cout << "Nodes in mesh " << m->NumberOfNodes() << std::endl;
int nfrac_nodes =0, new_nodes=0;
for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it)
{
ElementArray<Face> n_faces = it->getFaces(); //retrive all faces joining at the node
......@@ -146,6 +165,7 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
nums.resize(std::unique(nums.begin(),nums.end()) - nums.begin());
if( nums.size() > 1 ) //at least two distinctive nodes
{
nfrac_nodes++;
for(int k = 0; k < nums.size(); k++)
{
Node image;
......@@ -169,6 +189,7 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
for( Storage::real_array::size_type q = 0; q < coords.size(); ++q)
xyz[q] = coords[q]*mult + (1-mult)*xyz[q];
image = m->CreateNode(xyz);
new_nodes++;
Storage::real_array save_coords = image->RealArray(origcoords);
save_coords[0] = coords[0];
save_coords[1] = coords[1];
......@@ -182,7 +203,7 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
for(int qq = 0; qq < source.size(); ++qq)
target[qq] = source[qq];
}
/*
for(int q = 0; q < (int)transfer_node_integer_tags.size(); ++q)
if( it->HaveData(transfer_node_integer_tags[q]) )
{
......@@ -192,7 +213,7 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
for(int qq = 0; qq < source.size(); ++qq)
target[qq] = source[qq];
}
*/
for(int q = 0; q < n_cells.size(); q++)
{
if( n_cells[q].IntegerDF(indexes) == nums[k] )
......@@ -209,6 +230,11 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
}
}
std::cout << "Nodes in mesh " << m->NumberOfNodes() << " frac nodes " << nfrac_nodes << " new nodes " << new_nodes << std::endl;
std::cout << "Cells: " << m->NumberOfCells() << std::endl;
std::cout << "Faces: " << m->NumberOfFaces() << std::endl;
std::cout << "Edges: " << m->NumberOfEdges() << std::endl;
std::cout << "Nodes: " << m->NumberOfNodes() << std::endl;
m->DeleteTag(indexes);
//this tag is used to transfer position of fracture edge center onto fracture face center
......@@ -234,16 +260,24 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
{
fracture_edge_length = m->CreateTag("LENGTH_FRACTURE",DATA_REAL,FACE,FACE,1);
//Mark edges that will be incorporated into fracture cells
int nedges = 0;
for(Mesh::iteratorEdge it = m->BeginEdge(); it != m->EndEdge(); ++it)
{
if( it->nbAdjElements(FACE,isFracture()) > 2) it->SetMarker(isFracture());
if( it->nbAdjElements(FACE,isFracture()) > 2)
{
it->SetMarker(isFracture());
nedges++;
}
}
std::cout << "fracture edges with multiple joints " << nedges << std::endl;
//create new matrix-matrix faces, adjacent to fractures
//for all non-fracture faces that have any fracture edge create new face
int nfaces = 0;
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( !it->GetMarker(isFracture()) && it->nbAdjElements(NODE,isFracture()) )
{
ElementArray<Edge> edges = it->getEdges(); //go through edges of current face
......@@ -280,6 +314,7 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
else newedges.push_back(*jt);
}
std::pair<Face,bool> f = m->CreateFace(newedges);
nfaces++;
if( !f.second ) std::cout << __FILE__ << ":" << __LINE__ << " Face already exists!!! source " << *it << " new " << f.first.GetHandle() << std::endl;
else
{
......@@ -307,6 +342,7 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
if(it->BackCell().isValid() ) it->BackCell()->ReferenceArray(connface).push_back(f.first);
if(it->FrontCell().isValid() ) it->FrontCell()->ReferenceArray(connface).push_back(f.first);
}
std::cout << "new fracture faces " << nfaces << std::endl;
if( !fracture_aperture.isDefined(CELL) ) // this is to extend volume factor
m->CreateTag(fracture_aperture.GetTagName(),fracture_aperture.GetDataType(),CELL,CELL,fracture_aperture.GetSize());
......@@ -345,7 +381,7 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
//now create fracture-fracture control volumes, add fracture-matrix faces to matrix cells
int kid = 0;
int kid = 0, new_edges = 0, new_faces = 0, new_cells = 0;
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->GetMarker(isFracture()) && it->nbAdjElements(NODE,isFracture()) )
{
//I need all elements adjacent to neighbouring cells of current face
......@@ -377,11 +413,13 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
enodes[0] = nodesb[k];
enodes[1] = nodesb[(k+1)%N];
edgesb[k] = m->CreateEdge(enodes).first;
new_edges++;
}
else edgesb.data()[k] = edges.data()[k];
}
//This is matrix-fracture face
std::pair<Face,bool> facesb = m->CreateFace(edgesb);
new_faces++;
if(!facesb.second ) std::cout << __FILE__ << ":" << __LINE__ << " Face already exists!!! " << *it << " " << facesb.first->GetHandle() << std::endl;
else
{
......@@ -432,11 +470,13 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
enodes[0] = nodesf[k];
enodes[1] = nodesf[(k+1)%N];
edgesf[k] = m->CreateEdge(enodes).first;
new_edges++;
}
else edgesf[k] = edges[k];
}
//This is matrix-fracture face
std::pair<Face,bool> facesf = m->CreateFace(edgesf);
new_faces++;
if( !facesf.second ) std::cout << __FILE__ << ":" << __LINE__ << " Face already exists!!! " << *it << " " << facesf.first->GetHandle() << std::endl;
else
{
......@@ -480,12 +520,13 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
enodes2[1] = nodes[(k+1)%N];
f_edges.push_back(edgesb[k]);
if( enodes1[0] != enodes1[1] ) f_edges.push_back(m->CreateEdge(enodes1).first);
if( enodes1[0] != enodes1[1] ) {f_edges.push_back(m->CreateEdge(enodes1).first); new_edges++;}
f_edges.push_back(edges[k]);
if( enodes2[0] != enodes2[1] ) f_edges.push_back(m->CreateEdge(enodes2).first);
if( enodes2[0] != enodes2[1] ) {f_edges.push_back(m->CreateEdge(enodes2).first); new_edges++;}
assert(f_edges.size() > 2 );
Face f1 = m->CreateFace(f_edges).first;
new_faces++;
f1->Real(fracture_edge_length) = edges[k]->Length();
if( centroid_tag.isValid() )
{
......@@ -506,15 +547,16 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
enodes4[0] = nodesf[(k+1)%N];
enodes4[1] = nodes[(k+1)%N];
if( enodes3[0] != enodes3[1] ) f_edges.push_back(m->CreateEdge(enodes3).first);
if( enodes3[0] != enodes3[1] ) {f_edges.push_back(m->CreateEdge(enodes3).first); new_edges++;}
f_edges.push_back(edgesf[k]);
if( enodes4[0] != enodes4[1] ) f_edges.push_back(m->CreateEdge(enodes4).first);
if( enodes4[0] != enodes4[1] ) {f_edges.push_back(m->CreateEdge(enodes4).first); new_edges++;}
f_edges.push_back(edges[k]);
assert(f_edges.size() > 2 );
Face f2 = m->CreateFace(f_edges).first;
new_faces++;
if( centroid_tag.isValid() )
{
edges[k]->Centroid(ecnt);
......@@ -539,16 +581,19 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
enodes[0] = nodesb[k];
enodes[1] = nodes[k];
f_edges.push_back(m->CreateEdge(enodes).first);
new_edges++;
enodes[0] = nodes[k];
enodes[1] = nodesf[k];
f_edges.push_back(m->CreateEdge(enodes).first);
new_edges++;
}
else if( nodes[k].GetMarker(isFracture()) )
{
enodes[0] = nodesb[k];
enodes[1] = nodesf[k];
f_edges.push_back(m->CreateEdge(enodes).first);
new_edges++;
}
f_edges.push_back(edgesf[k]);
......@@ -557,16 +602,19 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
enodes[0] = nodesf[(k+1)%N];
enodes[1] = nodes[(k+1)%N];
f_edges.push_back(m->CreateEdge(enodes).first);
new_edges++;
enodes[0] = nodes[(k+1)%N];
enodes[1] = nodesb[(k+1)%N];
f_edges.push_back(m->CreateEdge(enodes).first);
new_edges++;
}
else if( nodes[(k+1)%N].GetMarker(isFracture()) )
{
enodes[0] = nodesf[(k+1)%N];
enodes[1] = nodesb[(k+1)%N];
f_edges.push_back(m->CreateEdge(enodes).first);
new_edges++;
}
f_edges.push_back(edgesb[k]);
......@@ -574,6 +622,7 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
assert(f_edges.size() > 2 );
Face f = m->CreateFace(f_edges).first;
new_faces++;
if( centroid_tag.isValid() )
{
edges[k]->Centroid(ecnt);
......@@ -608,14 +657,15 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
enodes2[0] = nodesb[(k+1)%N];
enodes2[1] = nodes[(k+1)%N];
if( enodes1[0] != enodes1[1] ) f_edges.push_back(m->CreateEdge(enodes1).first);
if( enodes1[0] != enodes1[1] ) {f_edges.push_back(m->CreateEdge(enodes1).first); new_edges++;}
f_edges.push_back(edgesb.data()[k]);
if( enodes2[0] != enodes2[1] ) f_edges.push_back(m->CreateEdge(enodes2).first);
if( enodes2[0] != enodes2[1] ) {f_edges.push_back(m->CreateEdge(enodes2).first); new_edges++;}
f_edges.push_back(edges.data()[k]);
assert(f_edges.size() > 2);
Face f = m->CreateFace(f_edges).first;
new_faces++;
fracfaces.push_back(f);
f->Real(fracture_edge_length) = edges[k]->Length();
f_edges.clear();
......@@ -639,6 +689,7 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
fout.close();
*/
Cell fcell = m->CreateCell(fracfaces).first;
new_cells++;
fcell->Real(fracture_volume) = it->Area()*it->Real(fracture_aperture);
//move permiability and volume factor to cell
/*
......@@ -670,7 +721,13 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
}
kid++;
}
std::cout << "kid " << kid << " new edges " << new_edges << " new faces " << new_faces << " new_cells " << new_cells << std::endl;
std::cout << "Cells: " << m->NumberOfCells() << std::endl;
std::cout << "Faces: " << m->NumberOfFaces() << std::endl;
std::cout << "Edges: " << m->NumberOfEdges() << std::endl;
std::cout << "Nodes: " << m->NumberOfNodes() << std::endl;
//now reconnect matrix cells to new faces
int reconnect = 0;
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) if( !it->GetMarker(isFracture()) && it->nbAdjElements(NODE,isFracture()) )
{
ElementArray<Face> change(m);
......@@ -682,12 +739,18 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
}
if( !change.empty() )
{
reconnect++;
Storage::reference_array newfaces = it->ReferenceArray(connface);
assert(change.size() == newfaces.size());
it->Disconnect(change.data(),(Storage::enumerator)change.size());
it->Connect(newfaces.data(),(Storage::enumerator)newfaces.size());
}
}
std::cout << "reconnect " << reconnect << std::endl;
std::cout << "Cells: " << m->NumberOfCells() << std::endl;
std::cout << "Faces: " << m->NumberOfFaces() << std::endl;
std::cout << "Edges: " << m->NumberOfEdges() << std::endl;
std::cout << "Nodes: " << m->NumberOfNodes() << std::endl;
}
else //number of dimensions is 2
{
......@@ -1101,7 +1164,10 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
it->SetMarker(multiple_fracture_joints);
}
//precalculate area for fractures
std::cout << "Cells: " << m->NumberOfCells() << std::endl;
std::cout << "Faces: " << m->NumberOfFaces() << std::endl;
std::cout << "Edges: " << m->NumberOfEdges() << std::endl;
std::cout << "Nodes: " << m->NumberOfNodes() << std::endl;
fracture_area = m->CreateTag("AREA_FRACTURE",DATA_REAL,FACE,FACE,1);
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->GetMarker(isFracture()) )
{
......@@ -1133,5 +1199,9 @@ void Fracture::Open(Tag aperture, bool fill_fracture, double gap_multiplier)
}
m->ApplyModification();
m->EndModification();
std::cout << "Cells: " << m->NumberOfCells() << std::endl;
std::cout << "Faces: " << m->NumberOfFaces() << std::endl;
std::cout << "Edges: " << m->NumberOfEdges() << std::endl;
std::cout << "Nodes: " << m->NumberOfNodes() << std::endl;
//m->Save("fracture_test.vtk");
}
......@@ -122,6 +122,7 @@ void Slice::SliceMesh(Mesh & m, bool remove_material_zero)
MarkerType original = m.CreateMarker();
std::cout << "marker original " << original << std::endl;
#if defined(USE_OMP)
#pragma omp parallel for
#endif
......@@ -135,7 +136,9 @@ void Slice::SliceMesh(Mesh & m, bool remove_material_zero)
for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it) it->SetMarker(original);
MarkerType slice = m.CreateMarker();
std::cout << "marker slice " << slice << std::endl;
MarkerType mrk = m.CreateMarker();
std::cout << "marker mrk " << mrk << std::endl;
int nslice = 0, nmark = 0;
for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it) if( !it->GetMarker(mrk) )
{
......@@ -290,6 +293,7 @@ void Slice::SliceMesh(Mesh & m, bool remove_material_zero)
nslice = 0;
MarkerType unique = m.CreateMarker();
std::cout << "marker unique " << unique << std::endl;
for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) if( !it->GetMarker(mrk) )
{
ElementArray<Node> nodes = it->getNodes(slice); //those nodes should be ordered so that each pair forms an edge
......@@ -523,6 +527,7 @@ void Slice::SliceMesh(Mesh & m, bool remove_material_zero)
//else std::cout << "Only one adjacent slice node, face " << it->LocalID() << std::endl;
}
m.ReleaseMarker(unique);
std::cout << "release marker unique " << unique << std::endl;
nmark = 0;
......@@ -571,8 +576,11 @@ void Slice::SliceMesh(Mesh & m, bool remove_material_zero)
nslice = 0;
TagInteger indx = m.CreateTag("TEMP_INDX",DATA_INTEGER,NODE|EDGE,NONE,1);
MarkerType visited = m.CreateMarker();
std::cout << "marker visited " << visited << std::endl;
MarkerType cmrk = m.CreateMarker();
std::cout << "marker cmrk " << cmrk << std::endl;
MarkerType isolate = m.CreateMarker();
std::cout << "marker isolate " << isolate << std::endl;
for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( !it->GetMarker(mrk) )
{
ElementArray<Edge> edges = it->getEdges(slice);
......@@ -1159,8 +1167,11 @@ void Slice::SliceMesh(Mesh & m, bool remove_material_zero)
}
}
m.ReleaseMarker(isolate);
std::cout << "release marker isolate " << isolate << std::endl;
m.ReleaseMarker(cmrk);
std::cout << "release marker cmrk " << cmrk << std::endl;
m.ReleaseMarker(visited);
std::cout << "release marker visited " << visited << std::endl;
m.DeleteTag(indx);
std::cout << "sliced cells: " << nslice << std::endl;
......@@ -1228,8 +1239,11 @@ void Slice::SliceMesh(Mesh & m, bool remove_material_zero)
}
m.ReleaseMarker(slice,NODE|EDGE);
std::cout << "release marker slice " << slice << std::endl;
m.ReleaseMarker(original,NODE|EDGE);
m.ReleaseMarker(mrk,EDGE|FACE|NODE);
std::cout << "release marker original " << original << std::endl;
m.ReleaseMarker(mrk,EDGE|FACE|NODE|CELL);
std::cout << "release marker mrk " << mrk << std::endl;
// m.DeleteTag(material);
// m.DeleteTag(level);
......
......@@ -8,14 +8,11 @@ class SliceFault : public Slice
{
double v[3]; // vector for projection in z direction
std::vector< std::pair<double,double> > curvxy; //a curve that defines a fault
double func(double x, double y, double z)
{