Skip to content

GitLab

Commit d2320a57 authored by Dmitry Bagaev's avatar Dmitry Bagaev
Browse files

Merge remote-tracking branch 'INMOST-DEV/master' 

# Conflicts:
#	Examples/ADMFD/main.cpp
#	Examples/GridGen/main.cpp
#	Tests/solver_test000/CMakeLists.txt
#	Tests/solver_test002/CMakeLists.txt
parents c7a00014 f6feda49
Expand all Hide whitespace changes
Inline Side-by-side
Examples/ADFVDiscr/main.cpp
View file @ d2320a57
......@@ -146,10 +146,10 @@ int main(int argc,char ** argv)
//~ BARRIER
//~ if( m->GetProcessorRank() == 0 ) std::cout << "Prepare geometric data: " << Timer()-ttt << std::endl;
{
Automatizator aut(m);
Automatizator aut;
Automatizator::MakeCurrent(&aut);
INMOST_DATA_ENUM_TYPE iphi = aut.RegisterDynamicTag(phi,CELL);
aut.EnumerateDynamicTags();
INMOST_DATA_ENUM_TYPE iphi = aut.RegisterTag(phi,CELL);
aut.EnumerateTags();
// Set the indeces intervals for the matrix and vectors
R.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
......@@ -175,7 +175,7 @@ int main(int argc,char ** argv)
Cell r2 = face->FrontCell();
if( ((!r1->isValid() || (s1 = r1->GetStatus()) == Element::Ghost)?0:1) +
((!r2->isValid() || (s2 = r2->GetStatus()) == Element::Ghost)?0:1) == 0) continue;
Storage::integer i1 = aut.GetDynamicIndex(r1,iphi), i2;
Storage::integer i1 = aut.GetIndex(r1,iphi), i2;
Storage::real f_nrm[3], r1_cnt[3], r2_cnt[3], f_cnt[3], d1, d2, D, v[3], T;
Storage::real f_area = face->Area(); // Get the face area
face->UnitNormal(f_nrm); // Get the face normal
......@@ -198,7 +198,7 @@ int main(int argc,char ** argv)
}
else
{
i2 = aut.GetDynamicIndex(r2,iphi);
i2 = aut.GetIndex(r2,iphi);
r2->Centroid(r2_cnt);
D = dot_prod(f_nrm,f_cnt);
d1 = fabs(dot_prod(r1_cnt,f_nrm) - D);
......@@ -276,7 +276,7 @@ int main(int argc,char ** argv)
{
Storage::real old = cell->Real(phi);
Storage::real exact = cell->Mean(func, 0); // Compute the mean value of the function over the cell
Storage::real res = Update[aut.GetDynamicIndex(cell->self(),iphi)];
Storage::real res = Update[aut.GetIndex(cell->self(),iphi)];
Storage::real sol = old-res;
Storage::real err = fabs (sol - exact);
if (err > local_err_C) local_err_C = err;
......
Examples/ADMFD/main.cpp
View file @ d2320a57
......@@ -547,11 +547,10 @@ int main(int argc,char ** argv)
ttt = Timer();
{ //Main loop for problem solution
Automatizator aut(m); // declare class to help manage unknowns
Automatizator aut; // declare class to help manage unknowns
Automatizator::MakeCurrent(&aut);
dynamic_variable P(aut,aut.RegisterDynamicTag(tag_P,CELL|FACE)); //register pressure as primary unknown
aut.EnumerateDynamicTags(); //enumerate all primary variables
dynamic_variable P(aut,aut.RegisterTag(tag_P,CELL|FACE)); //register pressure as primary unknown
aut.EnumerateTags(); //enumerate all primary variables
std::cout << "Enumeration done, size " << aut.GetLastIndex() - aut.GetFirstIndex() << std::endl;
Residual R("",aut.GetFirstIndex(),aut.GetLastIndex());
......
Examples/ADMFD/main_nln.cpp 0 → 100644
View file @ d2320a57
This diff is collapsed.
Examples/GridGen/main.cpp
View file @ d2320a57
......@@ -25,11 +25,11 @@ using namespace INMOST;
(1)*-------*(3) */
void CreateCubeElement(Mesh *m, ElementArray<Node> verts)
{
// Define six cube faces assuming verts are numerated in the way presented above
const INMOST_DATA_INTEGER_TYPE face_nodes[24] = {0,4,6,2, 1,3,7,5, 0,1,5,4, 2,6,7,3, 0,2,3,1, 4,5,7,6};
const INMOST_DATA_INTEGER_TYPE num_nodes[6] = {4, 4, 4, 4, 4, 4};
// Define six cube faces assuming verts are numerated in the way presented above
const INMOST_DATA_INTEGER_TYPE face_nodes[24] = {0,4,6,2, 1,3,7,5, 0,1,5,4, 2,6,7,3, 0,2,3,1, 4,5,7,6};
const INMOST_DATA_INTEGER_TYPE num_nodes[6] = {4, 4, 4, 4, 4, 4};
m->CreateCell(verts,face_nodes,num_nodes,6); // Create the cubic cell in the mesh
m->CreateCell(verts,face_nodes,num_nodes,6); // Create the cubic cell in the mesh
}
/* (4)*-------*(6)
......@@ -47,14 +47,14 @@ void CreateCubeElement(Mesh *m, ElementArray<Node> verts)
(1)*-------*(3) */
void CreateNEPrismElements(Mesh *m, ElementArray<Node> verts)
{
// Define prism faces assuming verts are numerated in the way presented above
const INMOST_DATA_INTEGER_TYPE nw_face_nodes[18] = {0,4,6,2, 2,6,5,1, 1,5,4,0, 0,2,1, 4,5,6};
const INMOST_DATA_INTEGER_TYPE nw_num_nodes[5] = {4, 4, 4, 3, 3};
const INMOST_DATA_INTEGER_TYPE se_face_nodes[18] = {1,5,6,2, 1,3,7,5, 7,3,2,6, 1,2,3, 6,5,7};
const INMOST_DATA_INTEGER_TYPE se_num_nodes[5] = {4, 4, 4, 3, 3};
m->CreateCell(verts,nw_face_nodes,nw_num_nodes,5); // Create north-west prismatic cell
m->CreateCell(verts,se_face_nodes,se_num_nodes,5); // Create south-east prismatic cell
// Define prism faces assuming verts are numerated in the way presented above
const INMOST_DATA_INTEGER_TYPE nw_face_nodes[18] = {0,4,6,2, 2,6,5,1, 1,5,4,0, 0,2,1, 4,5,6};
const INMOST_DATA_INTEGER_TYPE nw_num_nodes[5] = {4, 4, 4, 3, 3};
const INMOST_DATA_INTEGER_TYPE se_face_nodes[18] = {1,5,6,2, 1,3,7,5, 7,3,2,6, 1,2,3, 6,5,7};
const INMOST_DATA_INTEGER_TYPE se_num_nodes[5] = {4, 4, 4, 3, 3};
m->CreateCell(verts,nw_face_nodes,nw_num_nodes,5); // Create north-west prismatic cell
m->CreateCell(verts,se_face_nodes,se_num_nodes,5); // Create south-east prismatic cell
}
/* (4)*-------*(6)
......@@ -72,14 +72,14 @@ void CreateNEPrismElements(Mesh *m, ElementArray<Node> verts)
(1)*-------*(3) */
void CreateNWPrismElements(Mesh *m, ElementArray<Node> verts)
{
// Define prism faces assuming verts are numerated in the way presented above
const INMOST_DATA_INTEGER_TYPE ne_face_nodes[18] = {0,4,6,2, 0,3,7,4, 2,6,7,3, 0,2,3, 4,7,6};
const INMOST_DATA_INTEGER_TYPE ne_num_nodes[5] = {4, 4, 4, 3, 3};
const INMOST_DATA_INTEGER_TYPE sw_face_nodes[18] = {0,4,7,3, 1,3,7,5, 0,1,5,4, 0,3,1, 4,5,7};
const INMOST_DATA_INTEGER_TYPE sw_num_nodes[5] = {4, 4, 4, 3, 3};
m->CreateCell(verts,ne_face_nodes,ne_num_nodes,5); // Create north-east prismatic cell
m->CreateCell(verts,sw_face_nodes,sw_num_nodes,5); // Create south-west prismatic cell
// Define prism faces assuming verts are numerated in the way presented above
const INMOST_DATA_INTEGER_TYPE ne_face_nodes[18] = {0,4,6,2, 0,3,7,4, 2,6,7,3, 0,2,3, 4,7,6};
const INMOST_DATA_INTEGER_TYPE ne_num_nodes[5] = {4, 4, 4, 3, 3};
const INMOST_DATA_INTEGER_TYPE sw_face_nodes[18] = {0,4,7,3, 1,3,7,5, 0,1,5,4, 0,3,1, 4,5,7};
const INMOST_DATA_INTEGER_TYPE sw_num_nodes[5] = {4, 4, 4, 3, 3};
m->CreateCell(verts,ne_face_nodes,ne_num_nodes,5); // Create north-east prismatic cell
m->CreateCell(verts,sw_face_nodes,sw_num_nodes,5); // Create south-west prismatic cell
}
......@@ -98,47 +98,48 @@ void CreateNWPrismElements(Mesh *m, ElementArray<Node> verts)
(1)*-------*(3) */
void CreateNWTetElements(Mesh *m, ElementArray<Node> verts)
{
// Define prism faces assuming verts are numerated in the way presented above
const INMOST_DATA_INTEGER_TYPE ne_face_nodes1[12] = {0,1,5, 5,1,3, 1,0,3, 3,0,5};
const INMOST_DATA_INTEGER_TYPE ne_num_nodes1[4] = {3,3,3,3};
// Define prism faces assuming verts are numerated in the way presented above
const INMOST_DATA_INTEGER_TYPE ne_face_nodes1[12] = {0,1,5, 5,1,3, 1,0,3, 3,0,5};
const INMOST_DATA_INTEGER_TYPE ne_num_nodes1[4] = {3,3,3,3};
const INMOST_DATA_INTEGER_TYPE ne_face_nodes2[12] = {0,3,5, 0,7,3, 5,3,7, 0,5,7};
const INMOST_DATA_INTEGER_TYPE ne_num_nodes2[4] = {3,3,3,3};
const INMOST_DATA_INTEGER_TYPE ne_face_nodes2[12] = {0,3,5, 0,7,3, 5,3,7, 0,5,7};
const INMOST_DATA_INTEGER_TYPE ne_num_nodes2[4] = {3,3,3,3};
const INMOST_DATA_INTEGER_TYPE ne_face_nodes3[12] = {0,7,5, 4,5,7, 0,5,4, 0,4,7};
const INMOST_DATA_INTEGER_TYPE ne_num_nodes3[4] = {3,3,3,3};
const INMOST_DATA_INTEGER_TYPE ne_face_nodes3[12] = {0,7,5, 4,5,7, 0,5,4, 0,4,7};
const INMOST_DATA_INTEGER_TYPE ne_num_nodes3[4] = {3,3,3,3};
const INMOST_DATA_INTEGER_TYPE sw_face_nodes1[12] = {0,3,7, 2,7,3, 0,7,2, 0,2,3};
const INMOST_DATA_INTEGER_TYPE sw_num_nodes1[4] = {3,3,3,3};
const INMOST_DATA_INTEGER_TYPE sw_face_nodes1[12] = {0,3,7, 2,7,3, 0,7,2, 0,2,3};
const INMOST_DATA_INTEGER_TYPE sw_num_nodes1[4] = {3,3,3,3};
const INMOST_DATA_INTEGER_TYPE sw_face_nodes2[12] = {0,7,4, 0,2,7, 2,4,7, 0,4,2};
const INMOST_DATA_INTEGER_TYPE sw_num_nodes2[4] = {3,3,3,3};
const INMOST_DATA_INTEGER_TYPE sw_face_nodes2[12] = {0,7,4, 0,2,7, 2,4,7, 0,4,2};
const INMOST_DATA_INTEGER_TYPE sw_num_nodes2[4] = {3,3,3,3};
const INMOST_DATA_INTEGER_TYPE sw_face_nodes3[12] = {4,6,2, 6,7,2, 4,7,6, 4,2,7};
const INMOST_DATA_INTEGER_TYPE sw_num_nodes3[4] = {3,3,3,3};
const INMOST_DATA_INTEGER_TYPE sw_face_nodes3[12] = {4,6,2, 6,7,2, 4,7,6, 4,2,7};
const INMOST_DATA_INTEGER_TYPE sw_num_nodes3[4] = {3,3,3,3};
m->CreateCell(verts,ne_face_nodes1,ne_num_nodes1,4); // Create north-east prismatic cell
m->CreateCell(verts,ne_face_nodes2,ne_num_nodes2,4); // Create north-east prismatic cell
m->CreateCell(verts,ne_face_nodes3,ne_num_nodes3,4); // Create north-east prismatic cell
m->CreateCell(verts,sw_face_nodes1,sw_num_nodes1,4); // Create south-west prismatic cell
m->CreateCell(verts,sw_face_nodes2,sw_num_nodes2,4); // Create south-west prismatic cell
m->CreateCell(verts,sw_face_nodes3,sw_num_nodes3,4); // Create south-west prismatic cell
m->CreateCell(verts,ne_face_nodes1,ne_num_nodes1,4); // Create north-east prismatic cell
m->CreateCell(verts,ne_face_nodes2,ne_num_nodes2,4); // Create north-east prismatic cell
m->CreateCell(verts,ne_face_nodes3,ne_num_nodes3,4); // Create north-east prismatic cell
m->CreateCell(verts,sw_face_nodes1,sw_num_nodes1,4); // Create south-west prismatic cell
m->CreateCell(verts,sw_face_nodes2,sw_num_nodes2,4); // Create south-west prismatic cell
m->CreateCell(verts,sw_face_nodes3,sw_num_nodes3,4); // Create south-west prismatic cell
}
Mesh * ParallelGenerator(INMOST_MPI_Comm comm, int ng, int nx, int ny, int nz)
{
int procs_per_axis[3] = {1,1,1};
int sizes[3] = {nx,ny,nz};
int rank,size;
Mesh * m = new Mesh(); // Create a mesh to be constructed
int procs_per_axis[3] = {1,1,1};
int sizes[3] = {nx,ny,nz};
int rank,size;
Mesh * m = new Mesh(); // Create a mesh to be constructed
m->SetCommunicator(comm); // Set the MPI communicator, usually MPI_COMM_WORLD
m->SetCommunicator(comm); // Set the MPI communicator, usually MPI_COMM_WORLD
#if defined(USE_MPI)
MPI_Comm_set_errhandler(comm,MPI_ERRORS_RETURN);
//MPI::COMM_WORLD.Set_errhandler(MPI::ERRORS_THROW_EXCEPTIONS);
MPI_Comm_set_errhandler(comm,MPI_ERRORS_RETURN);
//MPI::COMM_WORLD.Set_errhandler(MPI::ERRORS_THROW_EXCEPTIONS);
#endif
rank = m->GetProcessorRank(); // Get the rank of the current process
size = m->GetProcessorsNumber(); // Get the number of processors used in communicator comm
......@@ -170,11 +171,11 @@ Mesh * ParallelGenerator(INMOST_MPI_Comm comm, int ng, int nx, int ny, int nz)
int localsize[3], localstart[3], localend[3];
int avgsize[3] =
{
(int)ceil((double)sizes[0]/procs_per_axis[0]),
(int)ceil((double)sizes[1]/procs_per_axis[1]),
(int)ceil((double)sizes[2]/procs_per_axis[2])
};
{
(int)ceil((double)sizes[0]/procs_per_axis[0]),
(int)ceil((double)sizes[1]/procs_per_axis[1]),
(int)ceil((double)sizes[2]/procs_per_axis[2])
};
for(int j = 0; j < 3; j++)
{
......@@ -216,9 +217,9 @@ Mesh * ParallelGenerator(INMOST_MPI_Comm comm, int ng, int nx, int ny, int nz)
verts.push_back(newverts[V_ID(i - 0,j - 1, k - 0)]); // 5 /* |(0)*- -|- -*(2) */
verts.push_back(newverts[V_ID(i - 1,j - 0, k - 0)]); // 6 /* | / | / */
verts.push_back(newverts[V_ID(i - 0,j - 0, k - 0)]); // 7 /* | | / */
/* |/ |/ */
/* |/ |/ */
// Create cells based on parameter ng /* (1)*-------*(3) */
if (ng == 5)
if (ng == 5) // Create tetrahedral grid
{
CreateNWTetElements(m,verts);
}
......@@ -244,79 +245,86 @@ Mesh * ParallelGenerator(INMOST_MPI_Comm comm, int ng, int nx, int ny, int nz)
int main(int argc, char *argv[])
{
std::string mesh_name;
Mesh * mesh;
Mesh::Initialize(&argc,&argv);
// Specify mesh configuration
int ng = 3, nx = MESH_SIZE, ny = MESH_SIZE, nz = MESH_SIZE, args = 0;
if( argc > 4 )
{
ng = atoi(argv[1]); // 3 - Prismatic, 4 - Cubic generator
nx = atoi(argv[2]);
ny = atoi(argv[3]);
nz = atoi(argv[4]);
args = 1;
}
// Construct a mesh in parallel
//MPI_Barrier(mesh->GetCommunicator());
double tt = Timer();
mesh = ParallelGenerator(INMOST_MPI_COMM_WORLD,ng,nx,ny,nz);
tt = Timer() - tt;
tt = mesh->Integrate(tt)/mesh->GetProcessorsNumber(); //???
if( mesh->GetProcessorRank() == 0 )
{
if( args == 0 ) std::cout << "Usage: " << argv[0] << " ng nx ny nz [output.[p]vtk]" << std::endl << "ng - 3 for prismatic mesh, 4 for cubic mesh" << std::endl;
if( args == 0 ) std::cout << "Default ";
if( ng == 4 )
{
std::cout << "Cubic ";
std::stringstream str;
str << "CUBIC_" << nx << "x" << ny << "x" << nz;
mesh_name = str.str();
}
else
{
std::cout << "Prismatic ";
std::stringstream str;
str << "PRISMATIC_" << nx << "x" << ny << "x" << nz;
mesh_name = str.str();
}
std::cout << "Grid: " << nx << " x " << ny << " x " << nz << std::endl;
std::cout << "Processors: " <<mesh->GetProcessorsNumber() << std::endl;
std::cout << "Mesh generator time: " << tt << std::endl;
}
std::string filename;
if (argc > 5)
{
filename = argv[5];
}
else
{
filename = "grid";
if( mesh->GetProcessorsNumber() == 1 )
filename += ".vtk";
else
filename += ".pvtk";
}
Storage::bulk_array name = mesh->self()->BulkArray(mesh->CreateTag("GRIDNAME",DATA_BULK,MESH,NONE));
name.replace(name.begin(),name.end(),mesh_name.begin(),mesh_name.end());
std::string mesh_name;
Mesh * mesh;
Mesh::Initialize(&argc,&argv);
// Specify mesh configuration
int ng = 3, nx = MESH_SIZE, ny = MESH_SIZE, nz = MESH_SIZE, args = 0;
if( argc > 4 )
{
ng = atoi(argv[1]); // 3 - Prismatic, 4 - Cubic generator
nx = atoi(argv[2]);
ny = atoi(argv[3]);
nz = atoi(argv[4]);
args = 1;
}
// Construct a mesh in parallel
//MPI_Barrier(mesh->GetCommunicator());
double tt = Timer();
mesh = ParallelGenerator(INMOST_MPI_COMM_WORLD,ng,nx,ny,nz);
tt = Timer() - tt;
tt = mesh->Integrate(tt)/mesh->GetProcessorsNumber(); //???
if( mesh->GetProcessorRank() == 0 )
{
if( args == 0 ) std::cout << "Usage: " << argv[0] << " ng nx ny nz [output.[p]vtk]" << std::endl << "ng - 3 for prismatic mesh, 4 for cubic mesh" << std::endl;
if( args == 0 ) std::cout << "Default ";
if( ng == 5 )
{
std::cout << "Tetrahedral ";
std::stringstream str;
str << "TETRAHEDRAL_" << nx << "x" << ny << "x" << nz;
mesh_name = str.str();
}
else if( ng == 4 )
{
std::cout << "Cubic ";
std::stringstream str;
str << "CUBIC_" << nx << "x" << ny << "x" << nz;
mesh_name = str.str();
}
else
{
std::cout << "Prismatic ";
std::stringstream str;
str << "PRISMATIC_" << nx << "x" << ny << "x" << nz;
mesh_name = str.str();
}
std::cout << "Grid: " << nx << " x " << ny << " x " << nz << std::endl;
std::cout << "Processors: " <<mesh->GetProcessorsNumber() << std::endl;
std::cout << "Mesh generator time: " << tt << std::endl;
}
std::string filename;
if (argc > 5)
{
filename = argv[5];
}
else
{
filename = "grid";
if( mesh->GetProcessorsNumber() == 1 )
filename += ".vtk";
else
filename += ".pvtk";
}
Storage::bulk_array name = mesh->self()->BulkArray(mesh->CreateTag("GRIDNAME",DATA_BULK,MESH,NONE));
name.replace(name.begin(),name.end(),mesh_name.begin(),mesh_name.end());
#if defined(USE_MPI)
MPI_Barrier(mesh->GetCommunicator());
MPI_Barrier(mesh->GetCommunicator());
#endif
tt = Timer();
mesh->Save(filename); // Save constructed mesh to the file
tt = Timer();
mesh->Save(filename); // Save constructed mesh to the file
#if defined(USE_MPI)
MPI_Barrier(mesh->GetCommunicator());
MPI_Barrier(mesh->GetCommunicator());
#endif
tt = Timer() - tt;
if( mesh->GetProcessorRank() == 0 ) std::cout << "Save to file \"" << filename << "\" time: " << tt << std::endl;
tt = Timer() - tt;
if( mesh->GetProcessorRank() == 0 ) std::cout << "Save to file \"" << filename << "\" time: " << tt << std::endl;
delete mesh;
Mesh::Finalize();
return 0;
delete mesh;
Mesh::Finalize();
return 0;
}
\ No newline at end of file
Examples/OctreeCutcell/main.cpp
View file @ d2320a57
......@@ -82,6 +82,9 @@ void transformation(double xyz[3])
mat_ret_type get_material_types(double xyz[3])
{
mat_ret_type ret1;
ret1.push_back(1);
return ret1;
#if defined(MASAHIKO_MESH)
mat_ret_type ret;
//ret.push_back(1);
......@@ -274,6 +277,7 @@ Storage::real get_press_by_mat(Storage::integer mat)
void fill_K(Storage::real * center, Storage::real * Kvec, Storage::real * K)
{
return;
Storage::real Knrm[3];
if( !get_type.Normal(center,Knrm) )
{
......@@ -923,14 +927,16 @@ int main(int argc, char ** argv)
make_proj.ReadLayers(layers);
}
{
std::vector<std::string> layers;
for(int i = 1; i <= 19; i++)
for(int i = 1; i <= 1; i++)
{
std::stringstream name;
//name << "Obj/rt" << i+1 << ".obj";
name << "oil_obj2/mat/layer" << i << ".obj";
layers.push_back(name.str());
}
/*
for(int i = 1; i <= 3; i++)
{
......
Examples/OctreeCutcell/octgrid.cpp
View file @ d2320a57
......@@ -1896,7 +1896,7 @@ void cellCreateINMOST(struct grid * g, int m, bool print = false)
replace.push_back(new_edge);
}
if( replace.size() == 4 && replace[1] == replace[3] || replace.size() < 4)
if( (replace.size() == 4 && replace[1] == replace[3]) || replace.size() < 4)
{
//if( new_edge != NULL ) face_edges.push_back(new_edge);
replace.clear(); // there is only one node and one edge left - avoid this edge duplication
......
Examples/OldDrawGrid/main.cpp
View file @ d2320a57
......@@ -4948,8 +4948,9 @@ int main(int argc, char ** argv)
if(mesh->HaveTag("ADDED_ELEMENTS") )
{
Tag add = mesh->GetTag("ADDED_ELEMENTS");
for(Mesh::iteratorEdge it = mesh->BeginEdge(); it != mesh->EndEdge(); ++it) if( it->Integer(add) ) added_edges.push_back(it->self());
for(Mesh::iteratorFace it = mesh->BeginFace(); it != mesh->EndFace(); ++it) if( it->Integer(add) && !it->Boundary() ) added_faces.push_back(DrawFace(it->self()));
if( add.isDefined(EDGE) ) for(Mesh::iteratorEdge it = mesh->BeginEdge(); it != mesh->EndEdge(); ++it) if( it->Integer(add) ) added_edges.push_back(it->self());
if( add.isDefined(FACE) ) for(Mesh::iteratorFace it = mesh->BeginFace(); it != mesh->EndFace(); ++it) if( it->Integer(add) && !it->Boundary() ) added_faces.push_back(DrawFace(it->self()));
}
if(mesh->HaveTag("CONORMALS"))
{
......@@ -5034,12 +5035,12 @@ int main(int argc, char ** argv)
}
if( false )
//if( false )
//if( mesh->HaveTag("VELOCITY") && mesh->GetTag("VELOCITY").isDefined(CELL) )
{
printf("preparing octree around mesh, was sets %d\n",mesh->NumberOfSets());
Octree octsearch = Octree(mesh->CreateSet("octsearch").first);
octsearch.Construct(NODE,true);
octsearch.Construct(NODE,false); //auto-detect octree or quadtree
printf("done, sets %d\n",mesh->NumberOfSets());
printf("building streamlines\n");
Tag cell_size = mesh->CreateTag("CELL_SIZES",DATA_REAL,CELL,NONE,1);
......@@ -5055,73 +5056,76 @@ int main(int argc, char ** argv)
}
velmax = log(velmax+1.0e-25);
velmin = log(velmin+1.0e-25);
if( mesh->HaveTag("FACE_FLUX") )
{
Tag flux = mesh->GetTag("FACE_FLUX");
MarkerType visited = mesh->CreateMarker();
for(Mesh::iteratorFace f = mesh->BeginFace(); f != mesh->EndFace(); ++f) if( f->Boundary() )
{
if( f->Real(flux) > 1.0e-5 )
{
coord cntf;
f->Centroid(cntf.data());
streamlines.push_back(Streamline(octsearch,cntf,vel,cell_size,velmin,velmax,1.0,visited));
ElementArray<Edge> edges = f->getEdges();
for(ElementArray<Edge>::iterator n = edges.begin(); n != edges.end(); ++n)
{
coord cntn;
n->Centroid(cntn.data());
if( cntn[2] == cntf[2] )
{
const Storage::real coef[4] = {0.4,0.8};
for(int q = 0; q < 2; ++q)
streamlines.push_back(Streamline(octsearch,cntf*coef[q]+cntn*(1-coef[q]),vel,cell_size,velmin,velmax,1.0,visited));
}
}
}
else if( f->Real(flux) < -1.0e-5 )
{
coord cntf;
f->Centroid(cntf.data());
streamlines.push_back(Streamline(octsearch,cntf,vel,cell_size,velmin,velmax,-1.0,visited));
ElementArray<Edge> edges = f->getEdges();
for(ElementArray<Edge>::iterator n = edges.begin(); n != edges.end(); ++n)
{
coord cntn;
n->Centroid(cntn.data());
if( cntn[2] == cntf[2] )
{
const Storage::real coef[4] = {0.4,0.8};
for(int q = 0; q < 2; ++q)
streamlines.push_back(Streamline(octsearch,cntf*coef[q]+cntn*(1-coef[q]),vel,cell_size,velmin,velmax,-1.0,visited));
}
}
}
}
printf("done from boundary faces, total streamlines = %d\n",streamlines.size());
/*
for(Mesh::iteratorCell it = mesh->BeginCell(); it != mesh->EndCell(); ++it)
{
if( !it->GetMarker(visited) )
{
coord cntc;
it->Centroid(cntc.data());
if( coord(it->RealArray(vel).data()).length() > 1.0e-4 )
{
streamlines.push_back(Streamline(octsearch,cntc,vel,cell_size,velmin,velmax,1.0,0));
streamlines.push_back(Streamline(octsearch,cntc,vel,cell_size,velmin,velmax,-1.0,0));
}
}
}
printf("done from unvisited cells, total streamlines = %d\n",streamlines.size());
*/
for(Mesh::iteratorCell it = mesh->BeginCell(); it != mesh->EndCell(); ++it)
it->RemMarker(visited);
mesh->ReleaseMarker(visited);
}
/*
Tag flux = mesh->GetTag("FACE_FLUX");
MarkerType visited = mesh->CreateMarker();
if( false )
for(Mesh::iteratorFace f = mesh->BeginFace(); f != mesh->EndFace(); ++f) if( f->Boundary() )
{
if( f->Real(flux) > 1.0e-5 )
{
coord cntf;
f->Centroid(cntf.data());
streamlines.push_back(Streamline(octsearch,cntf,vel,cell_size,velmin,velmax,1.0,visited));
ElementArray<Edge> edges = f->getEdges();
for(ElementArray<Edge>::iterator n = edges.begin(); n != edges.end(); ++n)
{
coord cntn;
n->Centroid(cntn.data());
if( cntn[2] == cntf[2] )
{
const Storage::real coef[4] = {0.4,0.8};
for(int q = 0; q < 2; ++q)
streamlines.push_back(Streamline(octsearch,cntf*coef[q]+cntn*(1-coef[q]),vel,cell_size,velmin,velmax,1.0,visited));
}
}
}
else if( f->Real(flux) < -1.0e-5 )
{
coord cntf;
f->Centroid(cntf.data());
streamlines.push_back(Streamline(octsearch,cntf,vel,cell_size,velmin,velmax,-1.0,visited));
ElementArray<Edge> edges = f->getEdges();
for(ElementArray<Edge>::iterator n = edges.begin(); n != edges.end(); ++n)
{
coord cntn;
n->Centroid(cntn.data());
if( cntn[2] == cntf[2] )
{
const Storage::real coef[4] = {0.4,0.8};
for(int q = 0; q < 2; ++q)
streamlines.push_back(Streamline(octsearch,cntf*coef[q]+cntn*(1-coef[q]),vel,cell_size,velmin,velmax,-1.0,visited));
}
}
}
}
printf("done from boundary faces, total streamlines = %d\n",streamlines.size());
for(Mesh::iteratorCell it = mesh->BeginCell(); it != mesh->EndCell(); ++it)
{
if( !it->GetMarker(visited) )
{
coord cntc;
it->Centroid(cntc.data());
if( coord(it->RealArray(vel).data()).length() > 1.0e-4 )