Commit f52d6054 authored by Kirill Terekhov's avatar Kirill Terekhov

Updates and fixes

Examples/AdaptiveMesh do not inherit Mesh anymore - easier to work with
the method

Examples/Fracture got the option to not to fill the space

SearchKDTree class to search for elements by coordinates in mesh (some
rearrangements and improvements may follow)

Code to support mesh adaptation in models

Fix boundary condition in TestSuite/cavity
parent 9b984278
Pipeline #170 failed with stages
in 10 minutes and 9 seconds
This diff is collapsed.
......@@ -4,8 +4,10 @@
namespace INMOST
{
class AdaptiveMesh : public Mesh
class AdaptiveMesh
{
Mesh * m;
Model * model;
ElementSet root; //< Root set that links all the other sets for coarsements
TagInteger level; //< Refinement level of the cell
TagReference parent_set; //<Link to the set that contains an element.
......@@ -15,7 +17,7 @@ namespace INMOST
void PrepareSet();
public:
Storage::integer GetLevel(const Storage & e) {return level[e];}
AdaptiveMesh();
AdaptiveMesh(Mesh & m);
~AdaptiveMesh();
/// Indicator must be 1 on cells to be refined
/// and 0 on all other cells
......@@ -23,6 +25,7 @@ namespace INMOST
bool Coarse(TagInteger & indicator);
/// Delete all data related to mesh refinement-coarsement.
void ClearData();
void SetModel(Model * mm) {model = mm;}
};
}
......
......@@ -8,9 +8,10 @@ int main(int argc, char ** argv)
if( argc > 1 )
{
AdaptiveMesh m;
Mesh m;
m.SetCommunicator(INMOST_MPI_COMM_WORLD);
m.Load(argv[1]);
AdaptiveMesh am(m);
//m.SetTopologyCheck(NEED_TEST_CLOSURE);
//m.SetTopologyCheck(PROHIBIT_MULTILINE);
//m.SetTopologyCheck(PROHIBIT_MULTIPOLYGON);
......@@ -56,7 +57,7 @@ int main(int argc, char ** argv)
do
{
numref = 0;
for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( m.GetLevel(it->self()) < max_levels )
for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( am.GetLevel(it->self()) < max_levels )
{
it->Barycenter(xyz);
//refine a circle
......@@ -82,7 +83,7 @@ int main(int argc, char ** argv)
m.Save(file.str());
}
*/
if( !m.Refine(indicator) ) break;
if( !am.Refine(indicator) ) break;
{
std::stringstream file;
......@@ -100,7 +101,7 @@ int main(int argc, char ** argv)
do
{
numref = 0;
for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( m.GetLevel(it->self()) > 0 )
for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( am.GetLevel(it->self()) > 0 )
{
it->Barycenter(xyz);
//refine a circle
......@@ -125,7 +126,7 @@ int main(int argc, char ** argv)
m.Save(file.str());
}
*/
if( !m.Coarse(indicator) ) break;
if( !am.Coarse(indicator) ) break;
{
std::stringstream file;
......
......@@ -8,9 +8,10 @@ int main(int argc, char ** argv)
if( argc > 1 )
{
AdaptiveMesh m;
m.Load(argv[1]);
TagInteger indicator = m.CreateTag("INDICATOR",DATA_INTEGER,CELL,NONE,1);
Mesh mm;
mm.Load(argv[1]);
TagInteger indicator = mm.CreateTag("INDICATOR",DATA_INTEGER,CELL,NONE,1);
AdaptiveMesh m(mm);
int max_levels = 2;
if( argc > 2 ) max_levels = atoi(argv[2]);
......@@ -19,7 +20,7 @@ int main(int argc, char ** argv)
do
{
numref = 0;
for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
for(Mesh::iteratorCell it = mm.BeginCell(); it != mm.EndCell(); ++it)
if( m.GetLevel(it->self()) < max_levels )
{
double x[3];
......@@ -33,13 +34,13 @@ int main(int argc, char ** argv)
if( numref )
{
if( !m.Refine(indicator) ) break;
for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) indicator[it->self()] = 0;
for(Mesh::iteratorCell it = mm.BeginCell(); it != mm.EndCell(); ++it) indicator[it->self()] = 0;
}
}
while(numref);
std::string file = "out.vtk";
if( argc > 3 ) file = std::string(argv[3]);
m.Save(file);
mm.Save(file);
}
else std::cout << "Usage: " << argv[0] << " mesh_file [max_levels=2] [mesh_out=out.vtk]" << std::endl;
}
\ No newline at end of file
}
......@@ -8,9 +8,10 @@ int main(int argc, char ** argv)
if( argc > 1 )
{
AdaptiveMesh m;
Mesh m;
m.Load(argv[1]);
TagInteger indicator = m.CreateTag("INDICATOR",DATA_INTEGER,CELL,NONE,1);
AdaptiveMesh am(m);
int max_levels = 2;
if( argc > 2 ) max_levels = atoi(argv[2]);
......@@ -20,7 +21,7 @@ int main(int argc, char ** argv)
{
numref = 0;
for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
if( m.GetLevel(it->self()) < max_levels )
if( am.GetLevel(it->self()) < max_levels )
{
double x[3];
it->Centroid(x);
......@@ -32,7 +33,7 @@ int main(int argc, char ** argv)
}
if( numref )
{
if( !m.Refine(indicator) ) break;
if( !am.Refine(indicator) ) break;
for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) indicator[it->self()] = 0;
}
}
......@@ -42,4 +43,4 @@ int main(int argc, char ** argv)
m.Save(file);
}
else std::cout << "Usage: " << argv[0] << " mesh_file [max_levels=2] [mesh_out=out.vtk]" << std::endl;
}
\ No newline at end of file
}
......@@ -29,9 +29,22 @@ add_executable(shestakov_mesh shestakov_grid.cpp)
add_executable(sinusoidal_mesh sinusoidal_grid.cpp)
add_executable(split_faces split_faces.cpp)
add_executable(glue_faces glue_faces.cpp)
add_executable(check_collapse check_collapse.cpp)
add_executable(check_collapse check_collapse.cpp)
add_executable(test_fracture test_fracture.cpp)
add_library(FractureLib fracture.cpp fracture.h)
target_link_libraries(test_fracture inmost)
target_link_libraries(test_fracture FractureLib)
if(USE_MPI)
message("linking test_fracture with MPI")
target_link_libraries(test_fracture ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(test_fracture PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS test_fracture EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(FixTinyCollapse inmost)
if(USE_MPI)
......
......@@ -72,11 +72,13 @@ void Fracture::FaceCenter(Face f, INMOST_DATA_REAL_TYPE cnt[3]) const
else f->Centroid(cnt);
}
void Fracture::Open(Tag aperture)
void Fracture::Open(Tag aperture, bool fill_fracture)
{
fracture_aperture = aperture;
m->BeginModification();
fracture_marker = m->CreateMarker();
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;
//break up all the faces with fractures into volumes?
fracture_volume = m->CreateTag("VOLUME_FRACTURE",DATA_REAL,CELL,CELL,1);
......@@ -213,7 +215,7 @@ void Fracture::Open(Tag aperture)
Tag centroid_tag;
if( m->HaveTag("GEOM_UTIL_CENTROID") )
centroid_tag = m->GetTag("GEOM_UTIL_CENTROID");
if( !centroid_tag.isDefined(FACE) )
if( centroid_tag.isValid() && !centroid_tag.isDefined(FACE) )
centroid_tag = Tag();
////////
std::vector<Tag> transfer_face_real_tags;
......@@ -620,7 +622,7 @@ void Fracture::Open(Tag aperture)
}
}
}
if( fill_fracture )
{
assert(fracfaces.size() >= 4); //tetrahedra or more
/*
......@@ -986,6 +988,7 @@ void Fracture::Open(Tag aperture)
}
}
if( fill_fracture )
{
assert(fracfaces.size() >= 3); //triangle or more
Cell fcell = m->CreateCell(fracfaces).first;
......
......@@ -29,7 +29,7 @@ public:
INMOST_DATA_REAL_TYPE Volume(Cell c) const;
INMOST_DATA_REAL_TYPE Area(Face f) const;
void FaceCenter(Face f, INMOST_DATA_REAL_TYPE cnt[3]) const;
void Open(Tag aperture);
void Open(Tag aperture, bool fill_fracture);
};
#endif //__FRACTURE_H
\ No newline at end of file
#endif //__FRACTURE_H
#include "inmost.h"
#include "fracture.h"
using namespace INMOST;
// cell-> node, edge, face
// face-> node, edge
// edge-> node
// based on eigenvalues in bounding ellipse
int main(int argc, char ** argv)
{
if (argc < 2)
{
std::cout << "Usage: " << argv[0] << " mesh [fill_fracture=true] [mesh_out=grid.vtk]" << std::endl;
return -1;
}
std::string grid_out = "grid.vtk";
bool fill_fracture = true;
if( argc > 2 ) fill_fracture = atoi(argv[2]) ? true : false;
if (argc > 3) grid_out = std::string(argv[3]);
Mesh m;
m.Load(argv[1]);
//Mesh::GeomParam table;
//table[MEASURE] = CELL;
//m.PrepareGeometricData(table);
TagReal aperture = m.CreateTag("aperture",DATA_REAL,FACE,FACE,1);
TagBulk sliced;
if( m.HaveTag("SLICED") )
sliced = m.GetTag("SLICED");
else
{
std::cout << "No tag named SLICED on the grid" << std::endl;
return -1;
}
for(int k = 0; k < m.FaceLastLocalID(); ++k) if( m.isValidFace(k) )
{
Face f = m.FaceByLocalID(k);
if( f.HaveData(sliced) )
aperture[f] = 0.001;
}
std::cout << "Start:" << std::endl;
std::cout << "Cells: " << m.NumberOfCells() << std::endl;
std::cout << "Faces: " << m.NumberOfFaces() << std::endl;
std::cout << "Edges: " << m.NumberOfEdges() << std::endl;
double tt = Timer();
Fracture f(m);
f.Open(aperture,fill_fracture);
std::cout << "Time to open fracture:" << Timer() - tt << std::endl;
std::cout << "Cells: " << m.NumberOfCells() << std::endl;
std::cout << "Faces: " << m.NumberOfFaces() << std::endl;
std::cout << "Edges: " << m.NumberOfEdges() << std::endl;
m.Save(grid_out);
return 0;
}
......@@ -67,8 +67,13 @@ int main(int argc, char ** argv)
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->Boundary() )
{
it->UnitNormal(n.data());
if( is2D && (fabs(n(2,0)-1) < 1.0-3 || fabs(n(2,0)+1) < 1.0e-3 ) ) //roller bc
//std::cout << "face " << it->LocalID() << std::endl;
//std::cout << std::scientific;
//std::cout << "n " << (n(2,0)-1) << " " << (fabs(n(2,0)-1) < 1.0e-3) << " " << (n(2,0)+1) << " " << (fabs(n(2,0)+1) < 1.0e-3) << std::endl;
//n.Transpose().Print();
if( is2D && (fabs(n(2,0)-1) < 1.0e-3 || fabs(n(2,0)+1) < 1.0e-3) ) //roller bc
{
//std::cout << "roller bc" << std::endl;
bc[*it][0] = 1;
bc[*it][1] = 0;
bc[*it][2] = 0;
......@@ -85,12 +90,14 @@ int main(int argc, char ** argv)
bc[*it][3] = 0;
if( fabs(n(1,0)-1) < 1.0e-3 ) //top
{
//std::cout << "dirichlet bc top" << std::endl;
bc[*it][4] = 1;
bc[*it][5] = 0;
bc[*it][6] = 0;
}
else
{
//std::cout << "dirichlet bc" << std::endl;
bc[*it][4] = 0;
bc[*it][5] = 0;
bc[*it][6] = 0;
......
......@@ -272,5 +272,20 @@ namespace INMOST
#endif
return ret;
}
void Model::Adaptation(Mesh & m) const
{
array<HandleType> old_cells;
for(int k = 0; k < m.CellLastLocalID(); ++k) if( m.isValidCell(k) )
{
Cell c = m.CellByLocalID(k);
if( c.Hidden() ) old_cells.push_back(c.GetHandle());
}
//std::cout << "old cells: " << old_cells.size() << std::endl;
SearchKDTree tree(&m,old_cells.data(),old_cells.size());
for(std::vector< std::pair<std::string, AbstractSubModel *> >::const_iterator it = SubModels.begin();
it != SubModels.end(); ++it)
it->second->Adaptation(m,tree);
}
}
#endif
......@@ -3388,6 +3388,259 @@ namespace INMOST
/// Find mesh by name.
static Mesh * GetMesh(std::string name);
};
class SearchKDTree
{
public:
inline static bool cell_point(const Cell & c, const Storage::real p[3])
{
return c.Inside(p);
}
template<typename bbox_type>
inline static int bbox_point(const Storage::real p[3], const bbox_type bbox[6])
{
for(int i = 0; i < 3; i++)
{
if( p[i] < bbox[i*2] || p[i] > bbox[i*2+1] )
return 0;
}
return 1;
}
private:
struct entry
{
HandleType e;
float xyz[3];
struct entry & operator =(const struct entry & other)
{
e = other.e;
xyz[0] = other.xyz[0];
xyz[1] = other.xyz[1];
xyz[2] = other.xyz[2];
return *this;
}
} * set;
Mesh * m;
INMOST_DATA_ENUM_TYPE size;
float bbox[6];
SearchKDTree * children;
static int cmpElements0(const void * a,const void * b)
{
const entry * ea = ((const entry *)a);
const entry * eb = ((const entry *)b);
float ad = ea->xyz[0];
float bd = eb->xyz[0];
return (ad > bd) - (ad < bd);
}
static int cmpElements1(const void * a,const void * b)
{
const entry * ea = ((const entry *)a);
const entry * eb = ((const entry *)b);
float ad = ea->xyz[1];
float bd = eb->xyz[1];
return (ad > bd) - (ad < bd);
}
static int cmpElements2(const void * a,const void * b)
{
const entry * ea = ((const entry *)a);
const entry * eb = ((const entry *)b);
float ad = ea->xyz[2];
float bd = eb->xyz[2];
return (ad > bd) - (ad < bd);
}
inline static unsigned int flip(const unsigned int * fp)
{
unsigned int mask = -((int)(*fp >> 31)) | 0x80000000;
return *fp ^ mask;
}
#define _m0(x) (x & 0x7FF)
#define _m1(x) (x >> 11 & 0x7FF)
#define _m2(x) (x >> 22 )
void radix_sort(int dim, struct entry * temp)
{
unsigned int i;
const unsigned int kHist = 2048;
unsigned int b0[kHist * 3];
unsigned int *b1 = b0 + kHist;
unsigned int *b2 = b1 + kHist;
memset(b0,0,sizeof(unsigned int)*kHist*3);
for (i = 0; i < size; i++)
{
unsigned int fi = flip((unsigned int *)&set[i].xyz[dim]);
++b0[_m0(fi)]; ++b1[_m1(fi)]; ++b2[_m2(fi)];
}
{
unsigned int sum0 = 0, sum1 = 0, sum2 = 0;
for (i = 0; i < kHist; i++)
{
b0[kHist-1] = b0[i] + sum0; b0[i] = sum0 - 1; sum0 = b0[kHist-1];
b1[kHist-1] = b1[i] + sum1; b1[i] = sum1 - 1; sum1 = b1[kHist-1];
b2[kHist-1] = b2[i] + sum2; b2[i] = sum2 - 1; sum2 = b2[kHist-1];
}
}
for (i = 0; i < size; i++) temp[++b0[_m0(flip((unsigned int *)&set[i].xyz[dim]))]] = set[i];
for (i = 0; i < size; i++) set[++b1[_m1(flip((unsigned int *)&temp[i].xyz[dim]))]] = temp[i];
for (i = 0; i < size; i++) temp[++b2[_m2(flip((unsigned int *)&set[i].xyz[dim]))]] = set[i];
for (i = 0; i < size; i++) set[i] = temp[i];
}
void kdtree_build(int dim, int & done, int total, struct entry * temp)
{
if( size > 1 )
{
if( size > 128 ) radix_sort(dim,temp); else
switch(dim)
{
case 0: qsort(set,size,sizeof(entry),cmpElements0);break;
case 1: qsort(set,size,sizeof(entry),cmpElements1);break;
case 2: qsort(set,size,sizeof(entry),cmpElements2);break;
}
children = static_cast<SearchKDTree *>(malloc(sizeof(SearchKDTree)*2));//new kdtree[2];
children[0].children = NULL;
children[0].set = set;
children[0].size = size/2;
children[0].m = m;
children[1].children = NULL;
children[1].set = set+size/2;
children[1].size = size - size/2;
children[1].m = m;
children[0].kdtree_build((dim+1)%3,done,total,temp);
children[1].kdtree_build((dim+1)%3,done,total,temp);
for(int k = 0; k < 3; k++)
{
bbox[0+2*k] = std::min<float>(children[0].bbox[0+2*k],children[1].bbox[0+2*k]);
bbox[1+2*k] = std::max<float>(children[0].bbox[1+2*k],children[1].bbox[1+2*k]);
}
}
else
{
assert(size == 1);
{
ElementArray<Node> nodes = Element(m,set[0].e)->getNodes();
bbox[0] = bbox[2] = bbox[4] = 1.0e20f;
bbox[1] = bbox[3] = bbox[5] = -1.0e20f;
for(ElementArray<Node>::size_type k = 0; k < nodes.size(); ++k)
{
Storage::real_array coords = nodes[k].Coords();
for(Storage::real_array::size_type q = 0; q < coords.size(); q++)
{
bbox[q*2+0] = std::min<float>(bbox[q*2+0],(float)coords[q]);
bbox[q*2+1] = std::max<float>(bbox[q*2+1],(float)coords[q]);
}
}
for(int k = m->GetDimensions(); k < 3; ++k)
{
bbox[k*2+0] = -1.0e20f;
bbox[k*2+1] = 1.0e20f;
}
}
}
}
SearchKDTree() : set(NULL), size(0), children(NULL) {}
Cell SubSearchCell(const Storage::real p[3])
{
Cell ret = InvalidCell();
if( size == 1 )
{
if( cell_point(Cell(m,set[0].e),p) )
ret = Cell(m,set[0].e);
}
else
{
assert(size > 1);
if( bbox_point(p,bbox) )
{
ret = children[0].SubSearchCell(p);
if( !ret.isValid() )
ret = children[1].SubSearchCell(p);
}
}
return ret;
}
void clear_children() { if( children ) {children[0].clear_children(); children[1].clear_children(); free(children);}}
public:
SearchKDTree(Mesh * m) : m(m), children(NULL)
{
size = m->NumberOfCells();
if( size )
{
set = new entry[size];
INMOST_DATA_ENUM_TYPE k = 0;
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
{
set[k].e = it->GetHandle();
Storage::real cnt[3] = {0,0,0};
m->GetGeometricData(set[k].e,CENTROID,cnt);
set[k].xyz[0] = (float)cnt[0];
set[k].xyz[1] = (float)cnt[1];
set[k].xyz[2] = (float)cnt[2];
++k;
}
int done = 0, total = size;
struct entry * temp = new entry[size];
kdtree_build(0,done,total,temp);
delete [] temp;
if( size > 1 )
{
for(int k = 0; k < 3; k++)
{
bbox[0+2*k] = std::min<float>(children[0].bbox[0+2*k],children[1].bbox[0+2*k]);
bbox[1+2*k] = std::max<float>(children[0].bbox[1+2*k],children[1].bbox[1+2*k]);
}
}
}
}
SearchKDTree(Mesh * m, HandleType * _set, unsigned set_size) : m(m), children(NULL)
{
size = set_size;
if( size )
{
set = new entry[size];
for(unsigned k = 0; k < set_size; ++k)
{
set[k].e = _set[k];
Storage::real cnt[3] = {0,0,0};
m->GetGeometricData(set[k].e,CENTROID,cnt);
set[k].xyz[0] = (float)cnt[0];
set[k].xyz[1] = (float)cnt[1];
set[k].xyz[2] = (float)cnt[2];
}
int done = 0, total = size;
struct entry * temp = new entry[size];
kdtree_build(0,done,total,temp);
delete [] temp;
if( size > 1 )
{
for(int k = 0; k < 3; k++)
{
bbox[0+2*k] = std::min<float>(children[0].bbox[0+2*k],children[1].bbox[0+2*k]);
bbox[1+2*k] = std::max<float>(children[0].bbox[1+2*k],children[1].bbox[1+2*k]);
}
}
}
else set = NULL;
}
~SearchKDTree()
{
if( size )
{
delete [] set;
clear_children();
}
}
Cell SearchCell(const Storage::real * point)
{
return SubSearchCell(point);
}
};
//////////////////////////////////////////////////////////////////////
......
......@@ -38,9 +38,13 @@ namespace INMOST
/// Roll back to previous step.
virtual bool RestoreTimeStep() = 0;
/// Calculate multiplier for update for this model. Can simply return 1.
virtual double UpdateMultiplier(const Sparse::Vector & sol) const {(void)sol; return 1;}
virtual double UpdateMultiplier(const Sparse::Vector & sol) const {(void)sol; return 1;}
/// Calculate time step for this model. Can simply return dt.
virtual double AdjustTimeStep(double dt) const {return dt;}
/// Adapt the data of the model after the mesh refinement/coarsement.
/// No algorithm by default
/// If this submodel depends on provided adpated mesh, it should update it's data
virtual void Adaptation(Mesh & m, SearchKDTree & search_old_cells) const {};
};
/// A class to organize a model.
......@@ -140,6 +144,9 @@ namespace INMOST
double UpdateMultiplier(const Sparse::Vector & sol) const;
/// Calculate optimal time step for submodels.
double AdjustTimeStep(double dt) const;
/// Adapt the data of the model after the mesh refinement/coarsement.
/// Those model that use the adapted mesh should update their data
virtual void Adaptation(Mesh & m) const;
};
}
......
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