Skip to content

I
INMOST
Commit a03370ea authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files
Options

remove Examples/TestSuite (moved to separate reposetory); fix some more... 

remove Examples/TestSuite (moved to separate reposetory); fix some more warnings; processing of data for webgl in DrawGrid
parent 0d56213d
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 268 additions and 7548 deletions
Examples/ADFVDiscr/main.cpp
View file @ a03370ea
......@@ -225,7 +225,7 @@ int main(int argc,char ** argv)
for(Storage::integer iface = 0; iface < m->FaceLastLocalID(); ++iface ) if( m->isValidFace(iface) )
{
Face face = Face(m,ComposeFaceHandle(iface));
Element::Status s1,s2;
Element::Status s1 = Element::Any, s2 = Element::Any;
Cell r1 = face->BackCell();
Cell r2 = face->FrontCell();
if( ((!r1->isValid() || (s1 = r1->GetStatus()) == Element::Ghost)?0:1) +
......
Examples/ADMFD/elastic.cpp
View file @ a03370ea
......@@ -247,7 +247,11 @@ void KTensor(const real_array & Cv, rMatrix & K)
void GetBC(const real_array & bc, const rMatrix & n, rMatrix & Ra, rMatrix & Rb, rMatrix & r)
{
const rMatrix I = rMatrix::Unit(3);
Storage::real alpha_perp, alpha_parallel, beta_perp, beta_parallel;
Storage::real
alpha_perp = 0,
alpha_parallel = 0,
beta_perp = 1,
beta_parallel = 1;
if( bc.size() == 6 )
{
......
Examples/ADVDIFF/conv_diff.cpp
View file @ a03370ea
......@@ -117,18 +117,18 @@ ConvectionDiffusion::ConvectionDiffusion(Mesh * _m, Tag _tag_U, Tag _tag_K, Tag
real A, //area of the face
U, //normal component of the velocity
C, //coefficient for upstream cell
T, //two-point transmissibility
R, //right hand side
dK, //distance from center to interface at cell K
dL, //distance from center to interface at cell L
lambdaK, //projection of co-normal onto normal at cell K
lambdaL, //projection of co-normal onto normal at cell L
T = 0, //two-point transmissibility
R = 0, //right hand side
dK = 0, //distance from center to interface at cell K
dL = 0, //distance from center to interface at cell L
lambdaK = 0, //projection of co-normal onto normal at cell K
lambdaL = 0, //projection of co-normal onto normal at cell L
div //divider
;
const real eps = degenerate_diffusion_regularization;
Cell cK, cL, cU;
Face fKL;
bulk flag_DIFF, flag_CONV;
bulk flag_DIFF = 0, flag_CONV = 0;
KK.Zero();
KL.Zero();
KD.Zero();
......
Examples/AdaptiveMesh/amesh.cpp
View file @ a03370ea
......@@ -386,6 +386,7 @@ namespace INMOST
void AdaptiveMesh::CheckParentSet(std::string file, int line)//, TagInteger indicator)
{
(void)file,(void)line;
ENTER_FUNC();
#if !defined(NDEBUG)
Storage::integer err = 0;
......@@ -1665,6 +1666,7 @@ namespace INMOST
{
//one (or both) of the adjacent cells were coarsened and has lower level
bool visited = false;
(void)visited;
ElementArray<Cell> cells = f.getCells();
for(ElementArray<Cell>::size_type kt = 0; kt < cells.size(); ++kt)
{
......@@ -1737,6 +1739,7 @@ namespace INMOST
{
//at least one face must have lower level
bool visited = false;
(void)visited;
ElementArray<Face> faces = e.getFaces();
for(ElementArray<Face>::size_type kt = 0; kt < faces.size(); ++kt)
{
......
Examples/AdaptiveMesh/main.cpp
View file @ a03370ea
......@@ -123,7 +123,7 @@ int main(int argc, char ** argv)
cnt[0] = cnt0[0] + 0.25*r0*sin(k/20.0*M_PI);
cnt[1] = cnt0[1] + 0.25*r0*cos(k/20.0*M_PI);
m.ClearFile();
//m.ClearFile();
std::fill(nc.begin(),nc.end(),0); nc[m.GetProcessorRank()] = m.NumberOfCells(); m.Integrate(&nc[0],nc.size()); if( !m.GetProcessorRank() ) {std::cout << "start "; for(unsigned q = 0; q < nc.size(); ++q) std::cout << nc[q] << " "; std::cout << std::endl;}
......
Examples/CMakeLists.txt
View file @ a03370ea
......@@ -23,7 +23,6 @@ if (USE_AUTODIFF AND USE_SOLVER AND USE_MESH)
add_subdirectory(ADFVDiscr)
add_subdirectory(ADMFD)
add_subdirectory(ADVDIFF)
add_subdirectory(TestSuite)
add_subdirectory(MFD-ES)
endif (USE_AUTODIFF AND USE_SOLVER AND USE_MESH)
#add_subdirectory(OctreeCutcell)
......
Examples/DrawGrid/face2gl.cpp
View file @ a03370ea
......@@ -371,6 +371,89 @@ namespace INMOST
glVertex3dv(&verts[(k + 3) % verts.size()]);
}
}
__INLINE void vfill4(const double vi[3], double tc, float v[4])
{
v[0] = vi[0], v[1] = vi[1], v[2] = vi[2], v[3] = tc;
}
__INLINE void vfill3(const double vi[3], float v[3])
{
v[0] = vi[0], v[1] = vi[1], v[2] = vi[2];
}
void face2gl::wgl_draw(std::ostream & file) const
{
size_t k;
for (k = 0; k < verts.size()-3; k += 3)
{
file << cnt[0] << " " << cnt[1] << " " << cnt[2] << " " << (texcoords.empty() ? 1.0 : cnttexcoord) << std::endl;
file << verts[k+0] << " " << verts[k+1] << " " << verts[k+2] << " " << (texcoords.empty() ? 1.0 : texcoords[k/3+0]) << std::endl;
file << verts[k+3] << " " << verts[k+4] << " " << verts[k+5] << " " << (texcoords.empty() ? 1.0 : texcoords[k/3+1]) << std::endl;
}
file << cnt[0] << " " << cnt[1] << " " << cnt[2] << " " << (texcoords.empty() ? 1.0 : cnttexcoord) << std::endl;
k = verts.size()-3;
file << verts[k+0] << " " << verts[k+1] << " " << verts[k+2] << " " << (texcoords.empty() ? 1.0 : texcoords[k/3+0]) << std::endl;
k = 0;
file << verts[k+0] << " " << verts[k+1] << " " << verts[k+2] << " " << (texcoords.empty() ? 1.0 : texcoords[k/3+0]) << std::endl;
}
void face2gl::wgl_draw_bin(std::ostream & file) const
{
size_t k;
float v[4];
for (k = 0; k < verts.size()-3; k += 3)
{
vfill4(cnt ,texcoords.empty() ? 1.0 : cnttexcoord,v);
file.write(reinterpret_cast<char *>(v),sizeof(float)*4);
vfill4(&verts[k+0],texcoords.empty() ? 1.0 : texcoords[k/3+0],v);
file.write(reinterpret_cast<char *>(v),sizeof(float)*4);
vfill4(&verts[k+3],texcoords.empty() ? 1.0 : texcoords[k/3+1],v);
file.write(reinterpret_cast<char *>(v),sizeof(float)*4);
}
vfill4(cnt ,texcoords.empty() ? 1.0 : cnttexcoord,v);
file.write(reinterpret_cast<char *>(v),sizeof(float)*4);
k = verts.size()-3;
vfill4(&verts[k+0],texcoords.empty() ? 1.0 : texcoords[k/3+0],v);
file.write(reinterpret_cast<char *>(v),sizeof(float)*4);
k = 0;
vfill4(&verts[k+0],texcoords.empty() ? 1.0 : texcoords[k/3+0],v);
file.write(reinterpret_cast<char *>(v),sizeof(float)*4);
}
void face2gl::wgl_draw_edges(std::ostream & file) const
{
size_t k;
for (k = 0; k < verts.size()-3; k += 3)
{
file << verts[k+0] << " " << verts[k+1] << " " << verts[k+2] << std::endl;
file << verts[k+3] << " " << verts[k+4] << " " << verts[k+5] << std::endl;
}
k = verts.size()-3;
file << verts[k+0] << " " << verts[k+1] << " " << verts[k+2] << std::endl;
k = 0;
file << verts[k+0] << " " << verts[k+1] << " " << verts[k+2] << std::endl;
}
void face2gl::wgl_draw_edges_bin(std::ostream & file) const
{
size_t k;
float v[3];
for (k = 0; k < verts.size()-3; k += 3)
{
vfill3(&verts[k+0],v);
file.write(reinterpret_cast<char *>(v),sizeof(float)*3);
vfill3(&verts[k+3],v);
file.write(reinterpret_cast<char *>(v),sizeof(float)*3);
}
k = verts.size()-3;
vfill3(&verts[k+0],v);
file.write(reinterpret_cast<char *>(v),sizeof(float)*3);
k = 0;
vfill3(&verts[k+0],v);
file.write(reinterpret_cast<char *>(v),sizeof(float)*3);
}
void face2gl::svg_draw(std::ostream & file, bool drawedges, double modelview[16], double projection[16], int viewport[4]) const
{
......@@ -528,6 +611,31 @@ namespace INMOST
glPrintError();
}
void wgl_draw_faces(std::ostream & file, const std::vector<face2gl> & set)
{
for (size_t q = 0; q < set.size(); q++)
set[q].wgl_draw(file);
}
void wgl_draw_edges(std::ostream & file, const std::vector<face2gl> & set)
{
for (size_t q = 0; q < set.size(); q++)
set[q].wgl_draw_edges(file);
}
void wgl_draw_faces_bin(std::ostream & file, const std::vector<face2gl> & set)
{
for (size_t q = 0; q < set.size(); q++)
set[q].wgl_draw_bin(file);
}
void wgl_draw_edges_bin(std::ostream & file, const std::vector<face2gl> & set)
{
for (size_t q = 0; q < set.size(); q++)
set[q].wgl_draw_edges_bin(file);
}
void svg_draw_faces(std::ostream & file, std::vector<face2gl> & set, bool drawedges, double modelview[16], double projection[16], int viewport[4], int highlight)
{
......
Examples/DrawGrid/face2gl.h
View file @ a03370ea
......@@ -33,6 +33,10 @@ namespace INMOST
void draw_colour_alpha(double alpha) const;
void draw() const;
void svg_draw(std::ostream & file, bool drawedges, double modelview[16], double projection[16], int viewport[4]) const;
void wgl_draw_bin(std::ostream & file) const;
void wgl_draw_edges_bin(std::ostream & file) const;
void wgl_draw(std::ostream & file) const;
void wgl_draw_edges(std::ostream & file) const;
void drawedges() const;
void svg_drawedges(std::ostream & file, double modelview[16], double projection[16], int viewport[4]) const;
bool operator <(const face2gl & other) const { return dist < other.dist; }
......@@ -65,6 +69,10 @@ namespace INMOST
void svg_draw_faces_nc(std::ostream & file, std::vector<face2gl> & set, bool drawedges, double modelview[16], double projection[16], int viewport[4], int highlight = -1);
void svg_draw_faces(std::ostream & file, std::vector<face2gl> & set, bool drawedges, double modelview[16], double projection[16], int viewport[4], int highlight = -1);
void svg_draw_edges(std::ostream & file, std::vector<face2gl> & set, double modelview[16], double projection[16], int viewport[4], int highlight = -1);
void wgl_draw_faces_bin(std::ostream & file, const std::vector<face2gl> & set);
void wgl_draw_edges_bin(std::ostream & file, const std::vector<face2gl> & set);
void wgl_draw_faces(std::ostream & file, const std::vector<face2gl> & set);
void wgl_draw_edges(std::ostream & file, const std::vector<face2gl> & set);
void draw_faces_nc(std::vector<face2gl> & set, int highlight = -1);
void draw_faces(std::vector<face2gl> & set, int highlight = -1);
......
Examples/DrawGrid/main.cpp
View file @ a03370ea
......@@ -852,6 +852,25 @@ void keyboard(unsigned char key, int x, int y)
{
write_state("state.txt");
}
else if( key == 'g' )
{
{
std::ofstream file("wgl.bin",std::ios::binary);
wgl_draw_faces_bin(file,clip_boundary);
}
{
std::ofstream file("wgle.bin",std::ios::binary);
wgl_draw_edges_bin(file,clip_boundary);
}
{
std::ofstream file("wgl.txt",std::ios::binary);
wgl_draw_faces(file,clip_boundary);
}
{
std::ofstream file("wgle.txt",std::ios::binary);
wgl_draw_edges(file,clip_boundary);
}
}
}
void keyboard2(unsigned char key, int x, int y)
......@@ -1126,7 +1145,7 @@ double display_elem_info(Element e, double top, double left, double interval)
if( e->HaveData(*t) )
{
std::stringstream sstr;
int dsize;
int dsize = 0;
switch(t->GetDataType())
{
case DATA_INTEGER:
......@@ -1217,7 +1236,7 @@ double display_elem_info_cout(Element e)
if( e->HaveData(*t) )
{
std::stringstream sstr;
int dsize;
int dsize = 0;
switch(t->GetDataType())
{
case DATA_INTEGER:
......@@ -1322,7 +1341,7 @@ void ProcessCommonInput(char inpstr[8192], int inptype)
}
else if( inptype == 3 )
{
int k1 = 0, k2, slen = (int)strlen(inpstr);
int k1 = 0, slen = (int)strlen(inpstr), k2 = slen;
for(int k = 0; k < slen; ++k)
{
if( inpstr[k] == ',' )
......
Examples/FVDiscr/main.cpp
View file @ a03370ea
......@@ -211,7 +211,7 @@ int main(int argc,char ** argv)
{
Face face(m,ComposeHandle(FACE,fi));
//~ std::cout << face->LocalID() << " / " << m->NumberOfFaces() << std::endl;
Element::Status s1,s2;
Element::Status s1 = Element::Any,s2 = Element::Any;
Cell r1 = face->BackCell();
Cell r2 = face->FrontCell();
if( ((!r1->isValid() || (s1 = r1->GetStatus()) == Element::Ghost)?0:1) +
......
Examples/GridTools/kershaw_grid.cpp
View file @ a03370ea
......@@ -66,7 +66,7 @@ int main(int argc, char *argv[])
for(int i = 0; i < n; ++i)
{
double x = i * 1.0 / (n - 1);
double ymid;
double ymid = 0;
if( i <= nsub[0] )
ymid = alpha;
else if( i <= nsub[1] )
......
Examples/GridTools/slice_func.cpp
View file @ a03370ea
......@@ -6,7 +6,7 @@ using namespace INMOST;
Storage::real Slice::Search(Storage::real r0, Storage::real r1, Storage::real c0[3], Storage::real c1[3], Storage::real p[3], bool binary) const
{
Storage::real rp = 1.0e20, rp_min = 1.0e20, p_min[3];
Storage::real rp = 1.0e20, rp_min = 1.0e20, p_min[3] = {0,0,0};
int iters = 0;
do
{
......
Examples/GridTools/split_nonplanar.cpp
View file @ a03370ea
......@@ -276,7 +276,7 @@ int main(int argc, char ** argv)
Tag proj_coords = A.CreateTag("PROJECTED_COORDS",DATA_REAL,NODE,NONE,2);
MarkerType myedges = A.CreateMarker();
real nrm[3], cnt[3], ncnt[3], scnt[3], pcnt[3], orthx[3], orthy[3], d, nd;
real nrm[3], cnt[3], ncnt[3], scnt[3], pcnt[3], orthx[3] = {0,0,0}, orthy[3] = {0,0,0}, d, nd;
(void)nd;
std::cout << "Start splitting faces" << std::endl;
int nsplit = 0, had_faces = A.NumberOfFaces();
......
Examples/GridTools/tetra.cpp
View file @ a03370ea
......@@ -345,7 +345,7 @@ int main(int argc, char ** argv)
std::vector<segment> segments, joints;
std::vector<node> loop;
std::map<HandleType,real> hits;
real nrm[3], cnt[3], ncnt[3], scnt[3], pcnt[3], fcnt[3], ray[3], orthx[3], orthy[3], d, nd;
real nrm[3], cnt[3], ncnt[3], scnt[3], pcnt[3], fcnt[3], ray[3], orthx[3] = {0,0,0}, orthy[3] = {0,0,0}, d, nd;
//for(Mesh::iteratorNode it = m.BeginNode(); it != m.EndNode(); ++it)
for(ElementArray<Node>::iterator it = nodes_rcm.begin(); it != nodes_rcm.end(); ++it)
{
......
Examples/TestSuite/CMakeLists.txt deleted 100644 → 0
View file @ 0d56213d
This diff is collapsed.
Examples/TestSuite/adv00.cpp deleted 100644 → 0
View file @ 0d56213d
#include "inmost.h"
using namespace INMOST;
const Storage::real pi = 3.1415926535897932384626433832795;
const Storage::real eps = 1.0e-6;
std::string problem_name = "knobloch_transverse_diffusion";
int main(int argc, char ** argv)
{
if( argc < 3 )
{
std::cout << "Usage: " << argv[0] << " mesh mesh_out" << std::endl;
return 0;
}
Mesh * m = new Mesh;
m->SetFileOption("VERBOSITY","2");
try{m->Load(argv[1]);} catch(...) { std::cout << "Cannot load the mesh " << argv[1] << std::endl; return -1;}
Mesh::GeomParam t;
t[MEASURE] = CELL | FACE;
t[ORIENTATION] = FACE;
t[NORMAL] = FACE;
//t[BARYCENTER] = CELL;
t[CENTROID] = CELL | FACE;
m->AssignGlobalID(CELL|FACE);
m->PrepareGeometricData(t);
double max[3] = {-1.0e20, -1.0e20, -1.0e20}, min[3] = {1.0e20,1.0e20,1.0e20};
Storage::real c[3] = {0,0,0}, nrm[3] = {0,0,0};
for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it)
{
it->Centroid(c);
if( c[0] > max[0] ) max[0] = c[0];
if( c[1] > max[1] ) max[1] = c[1];
if( c[2] > max[2] ) max[2] = c[2];
if( c[0] < min[0] ) min[0] = c[0];
if( c[1] < min[1] ) min[1] = c[1];
if( c[2] < min[2] ) min[2] = c[2];
}
if( max[0] <= min[0] ) {std::cout << "strange X " << min[0] << ":" << max[0] << std::endl; return -1;}
if( max[1] <= min[1] ) {std::cout << "strange Y " << min[1] << ":" << max[1] << std::endl; return -1;}
if( max[2] <= min[2] )
{
//2d mesh
if( m->GetDimensions() == 3 )
{
//offset from z-plane
min[2] -= 0.0001;
max[2] += 0.0001;
}
else
{
min[2] = -0.0001;
max[2] = +0.0001;
}
}
std::cout << "Mesh bounds: " << min[0] << ":" << max[0] << " " << min[1] << ":" << max[1] << " " << min[2] << ":" << max[2] << std::endl;
Tag material;
//if( m->HaveTag("PERM") ) m->DeleteTag(m->GetTag("PERM"));
//Tag force = m->CreateTag("FORCE",DATA_REAL,CELL,NONE,1);
//Tag tensor = m->CreateTag("PERM",DATA_REAL,CELL,NONE,6);
//Tag solution_val = m->CreateTag("REFERENCE_SOLUTION",DATA_REAL,CELL|FACE|EDGE|NODE,NONE,1);
//Tag solution_flux = m->CreateTag("REFERENCE_FLUX",DATA_REAL,FACE,FACE,1);
Tag bndcond = m->CreateTag("BOUNDARY_CONDITION",DATA_REAL,FACE,FACE,3);
Tag velocity = m->CreateTag("VELOCITY",DATA_REAL,FACE,NONE,1);
{
Storage::bulk_array name = m->self()->BulkArray(m->CreateTag("PROBLEMNAME",DATA_BULK,MESH,NONE));
name.replace(name.begin(),name.end(),problem_name.begin(),problem_name.end());
}
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it)
{
it->Centroid(c);
it->UnitNormal(nrm);
//fill in normal component of the velocity
it->Real(velocity) = nrm[0]*cos(pi/3.0) - nrm[1]*sin(pi/3.0);
Storage::real x = c[0];//(c[0]-min[0])/(max[0]-min[0]);
Storage::real y = c[1];//(c[1]-min[1])/(max[1]-min[1]);
Storage::real z = c[2];//(c[2]-min[2])/(max[2]-min[2]);
if( fabs(x) < 1.0e-7 || fabs(y-1) < 1.0e-7 ) //top of the domain or left domain
{
Storage::real_array bc = it->RealArray(bndcond); //create bc
//zero dirichlet bc
bc[0] = 1;
bc[1] = 0;
bc[2] = 1;
if( y <= 0.7 ) bc[2] = 0;
if( it->Real(velocity) > 0.0 ) std::cout << "Dirichlet condition next to outlet cell " << x << "," << y << "," << z << std::endl;
}
else if( it->Boundary() && it->Real(velocity) < 0.0 ) std::cout << "Neumann condition next to inlet cell " << x << "," << y << "," << z << std::endl;
}
std::cout << "Saving output to " << argv[2] << std::endl;
m->Save(argv[2]);
delete m;
}
Examples/TestSuite/adv01.cpp deleted 100644 → 0
View file @ 0d56213d
#include "inmost.h"
using namespace INMOST;
typedef Storage::real real;
typedef Storage::real_array real_array;
typedef Storage::bulk_array bulk_array;
const real pi = 3.1415926535897932384626433832795;
const real eps = 1.0e-6;
std::string problem_name = "zalesak_disk";
real max[3] = {-1.0e20, -1.0e20, -1.0e20}, min[3] = {1.0e20,1.0e20,1.0e20};
real disc_cx = 0.5;
real disc_cy = 0.75;
real disc_rad = 15.0/100.0;
real slot_top = disc_cy + 0.1;
real slot_w = 5.0/100.0;
real slot_l = disc_cx - slot_w*0.5;
real slot_r = disc_cx + slot_w*0.5;
real zalesak_shape(real * coord, real time)
{
(void)time;
real x = (coord[0]-min[0])/(max[0]-min[0]); //normalize x into [0,1]
real y = (coord[1]-min[1])/(max[1]-min[1]); //normalize y into [0,1]
if( sqrt((x-disc_cx)*(x-disc_cx) + (y-disc_cy)*(y-disc_cy)) < disc_rad )
{
if( y < slot_top && x > slot_l && x < slot_r )
return 1.0;
return 3.0;
}
return 1.0;
}
int main(int argc, char ** argv)
{
if( argc < 3 )
{
std::cout << "Usage: " << argv[0] << " mesh mesh_out cfl" << std::endl;
return 0;
}
Mesh * m = new Mesh;
m->SetFileOption("VERBOSITY","2");
try{m->Load(argv[1]);} catch(...) { std::cout << "Cannot load the mesh " << argv[1] << std::endl; return -1;}
Mesh::GeomParam t;
t[MEASURE] = CELL | FACE;
t[ORIENTATION] = FACE;
t[NORMAL] = FACE;
//t[BARYCENTER] = CELL;
t[CENTROID] = CELL | FACE;
m->AssignGlobalID(CELL|FACE);
m->PrepareGeometricData(t);
real c[3] = {0,0,0}, nrm[3] = {0,0,0};
for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it)
{
it->Centroid(c);
if( c[0] > max[0] ) max[0] = c[0];
if( c[1] > max[1] ) max[1] = c[1];
if( c[2] > max[2] ) max[2] = c[2];
if( c[0] < min[0] ) min[0] = c[0];
if( c[1] < min[1] ) min[1] = c[1];
if( c[2] < min[2] ) min[2] = c[2];
}
if( max[0] <= min[0] ) {std::cout << "strange X " << min[0] << ":" << max[0] << std::endl; return -1;}
if( max[1] <= min[1] ) {std::cout << "strange Y " << min[1] << ":" << max[1] << std::endl; return -1;}
if( max[2] <= min[2] )
{
//2d mesh
if( m->GetDimensions() == 3 )
{
//offset from z-plane
min[2] -= 0.0001;
max[2] += 0.0001;
}
else
{
min[2] = -0.0001;
max[2] = +0.0001;
}
}
real mesh_radius = 0;
for(Mesh::iteratorCell e = m->BeginCell(); e != m->EndCell(); ++e)
{
real maxmin[6];
maxmin[0] = -1e20;
maxmin[1] = 1e20;
maxmin[2] = -1e20;
maxmin[3] = 1e20;
maxmin[4] = -1e20;
maxmin[5] = 1e20;
ElementArray<Node> nodes = e->getNodes();
for (ElementArray<Node>::iterator it = nodes.begin(); it != nodes.end(); it++)
{
real_array c = it->Coords();
for (int i = 0; i < (int)c.size(); i++)
{
if (maxmin[2 * i + 0] < c[i]) maxmin[2 * i + 0] = c[i]; //max
if (maxmin[2 * i + 1] > c[i]) maxmin[2 * i + 1] = c[i]; //min
}
if( c.size() < 3 )
{
for(int i = c.size(); i < 3; i++)
{
maxmin[2*i+0] = maxmin[2*i+1] = 0;
}
}
}
mesh_radius = std::max(mesh_radius,maxmin[0]-maxmin[1]);
mesh_radius = std::max(mesh_radius,maxmin[2]-maxmin[3]);
//mesh_radius = std::max(mesh_radius,maxmin[4]-maxmin[5]);
}
mesh_radius *= 1;
std::cout << "Mesh bounds: " << min[0] << ":" << max[0] << " " << min[1] << ":" << max[1] << " " << min[2] << ":" << max[2] << " radius " << mesh_radius << std::endl;
Tag material;
//if( m->HaveTag("PERM") ) m->DeleteTag(m->GetTag("PERM"));
Tag time = m->CreateTag("TIME_INFO",DATA_REAL,MESH,NONE,2); //first entry - time step, second entry - total time
Tag bndcond = m->CreateTag("BOUNDARY_CONDITION",DATA_REAL,FACE,FACE,3); //boundary conditions
Tag velocity = m->CreateTag("VELOCITY",DATA_REAL,FACE,NONE,1); //normal component of the velocity
Tag conc = m->CreateTag("CONCENTRATION",DATA_REAL,CELL,NONE,1); //initial concentration
{
bulk_array name = m->self()->BulkArray(m->CreateTag("PROBLEMNAME",DATA_BULK,MESH,NONE));
name.replace(name.begin(),name.end(),problem_name.begin(),problem_name.end());
}
real_array time_info = m->self()->RealArray(time);
time_info[0] = 1.0/628.0;// * 14.277739148561303614695220241265 / 2.5; // cfl 10 on cartesian mesh 1/100
time_info[1] = 1;
real dt = time_info[0], cfl = 0;
real origin_x = (max[0]+min[0])*0.5;
real origin_y = (max[1]+min[1])*0.5;
real omega = 2*pi; //one circle in one second
std::cout << "Origin axis: (" << origin_x << "," << origin_y << ") angular velocity " << omega << std::endl;
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it)
{
it->Centroid(c);
it->UnitNormal(nrm);
real v = omega*(c[0] - origin_x);
real u = omega*(origin_y - c[1]);
//fill in normal component of the velocity
it->Real(velocity) = nrm[0]*u + nrm[1]*v;
if( it->Boundary() )
{
real_array bc = it->RealArray(bndcond); //create bc
//dirichlet bc with concentration 1
bc[0] = 1;
bc[1] = 0;
bc[2] = 1;
}
}
Tag velrec = m->CreateTag("VELREC",DATA_REAL,CELL,NONE,3);
real Vmax = 0.0;
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
{
real_array V = it->RealArray(velrec);
real xK[3], xF[3], sgn, A, U, Vol = it->Volume();
ElementArray<Face> faces = it->getFaces();
it->Centroid(xK);
for(int k = 0; k < (int)faces.size(); ++k)
{