Commit ffe8de78 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

update N-S turek test

parent f3cc075f
......@@ -44,7 +44,9 @@ int main(int argc,char ** argv)
Mesh * m = new Mesh(); // Create an empty mesh
{ // Load the mesh
ttt = Timer();
m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
//~ m->SetParallelFileStrategy(0);
//~ m->SetParallelStrategy(1);
//~ m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
if( m->GetProcessorRank() == 0 ) std::cout << "Processors: " << m->GetProcessorsNumber() << std::endl;
if( m->isParallelFileFormat(argv[1]) ) //The format is
{
......@@ -341,8 +343,10 @@ int main(int argc,char ** argv)
if( R.Norm() < 1.0e-4 ) break;
//Solver S(Solver::INNER_ILU2);
Solver S(Solver::INNER_MPTILUC);
//~ Solver S(Solver::INNER_MPTILUC);
Solver S(Solver::K3BIILU2);
//Solver S("superlu");
S.SetParameter("verbosity","1");
S.SetParameter("relative_tolerance", "1.0e-14");
S.SetParameter("absolute_tolerance", "1.0e-12");
S.SetParameter("drop_tolerance", "1.0e-2");
......
......@@ -32,7 +32,7 @@ int main(int argc, char ** argv)
if( argc < 2 )
{
std::cout << "Usage: " << argv[0] << " mesh [mesh_out=grid_out.pmf] [Umax=2.25]" << std::endl;
std::cout << "Usage: " << argv[0] << " mesh [mesh_out=grid_out.pmf] [Umax=2.25] [fix_cylinder=false]" << std::endl;
return 0;
}
......@@ -50,8 +50,10 @@ int main(int argc, char ** argv)
std::string fout = "grid_out.pmf";
double Umax = 2.25;
bool fix_cylinder = false;
if( argc > 2 ) fout = std::string(argv[2]);
if( argc > 3 ) Umax = atof(argv[3]);
if( argc > 4 ) fix_cylinder = atoi(argv[4]);
......@@ -84,14 +86,111 @@ int main(int argc, char ** argv)
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->Boundary() )
it->FixNormalOrientation();
{ // prepare geometrical data on the mesh
Mesh::GeomParam table;
table[CENTROID] = CELL | FACE; //Compute averaged center of mass
table[NORMAL] = FACE; //Compute normals
table[MEASURE] = CELL | FACE; //Compute volumes and areas
table[BARYCENTER] = CELL | FACE; //Compute volumetric center of mass
m->RemoveGeometricData(table); //Ask to precompute the data
}
if( fix_cylinder )
{
MarkerType cylinder = m->CreateMarker();
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->Boundary() )
{
double n[3], c[3];
it->UnitNormal(n);
it->Centroid(c);
if( fabs(n[0]-1) < 1.0e-3 && c[0] > 2.5-eps) // outflow
{
}
else if( fabs(n[0]+1) < 1.0e-3 && c[0] < 0.0+eps) //inflow
{
}
else //no-slip walls
{
if( fabs(n[2]) < 0.7 && c[0] > 0.45-eps && c[0] < 0.55+eps && c[1] > 0.15-eps && c[1] < 0.25+eps )
it->SetMarker(cylinder);
}
}
for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it)
if( it->nbAdjElements(FACE,cylinder) ) it->SetMarker(cylinder);
TagRealArray vec_t = m->CreateTag("vec_t",DATA_REAL,FACE,FACE,2);
std::cout << "project boundary nodes onto cylinder " << std::endl;
for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) if( it->GetMarker(cylinder) )
{
double x = it->Coords()[0], y = it->Coords()[1], dx, dy;
double r = sqrt((x-0.5)*(x-0.5) + (y-0.2)*(y-0.2));
if( r )
{
dx = (0.05/r-1.0)*(x-0.5);
dy = (0.05/r-1.0)*(y-0.2);
//std::cout << "at " << x << "," << y << " r " << r << " dx " << dx << " dy " << dy << std::endl;
it->Coords()[0] += dx;
it->Coords()[1] += dy;
}
else std::cout << "node at center of cylinder: " << x << "," << y << std::endl;
}
std::cout << "project centers of boundary faces onto cylinder " << std::endl;
int iter = 0;
while(iter < 1000)
{
double A = 0, err = 0, fA;
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->GetMarker(cylinder) )
{
double cnt[3], dx, dy;
it->Centroid(cnt);
double r = sqrt((cnt[0]-0.5)*(cnt[0]-0.5) + (cnt[1]-0.2)*(cnt[1]-0.2));
if( r )
{
dx = (0.05/r-1.0)*(cnt[0]-0.5);
dy = (0.05/r-1.0)*(cnt[1]-0.2);
//std::cout << "at " << cnt[0] << "," << cnt[1] << " r " << r << " dx " << dx << " dy " << dy << std::endl;
vec_t[*it][0] = dx;
vec_t[*it][1] = dy;
fA = it->Area();
err += sqrt(dx*dx+dy*dy)*fA;
A += fA;
}
else std::cout << " face center is at center of cylinder: " << cnt[0] << "," << cnt[1] << std::endl;
}
err = err/A;
if( iter % 50 == 0 )
std::cout << "iter " << iter << " error: " << err << " area " << A << std::endl;
if( err < 1.0e-8 ) break;
for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) if( it->GetMarker(cylinder) )
{
ElementArray<Face> faces = it->getFaces(cylinder);
double dxy[2] = {0,0}, dxyA = 0;
for(ElementArray<Face>::iterator jt = faces.begin(); jt != faces.end(); ++jt)
{
fA = jt->Area();
dxy[0] += vec_t[*jt][0]*fA;
dxy[1] += vec_t[*jt][1]*fA;
dxyA += fA;
}
dxy[0] /= dxyA;
dxy[1] /= dxyA;
it->Coords()[0] += dxy[0];
it->Coords()[1] += dxy[1];
}
iter++;
}
m->DeleteTag(vec_t);
m->ReleaseMarker(cylinder,FACE);
}
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->Boundary() )
{
double n[3], c[3];
it->UnitNormal(n);
it->Barycenter(c);
it->Centroid(c);
bcphi[*it][0] = 1;
bcphi[*it][1] = 0;
bcphi[*it][2] = 0;
......
#ifdef _MSC_VER //kill some warnings
#define _CRT_SECURE_NO_WARNINGS
#endif
#include "inmost.h"
#include <cfloat>
#if defined(USE_MESH)
//vtk states
#define R_VERSION 0
#define R_USERDATA 1
#define R_DATATYPE 2
#define R_WAITDATA 3
#define R_SGRID 4
#define R_RGRID 5
#define R_UGRID 6
#define R_POLYDATA 7
#define R_FIELD 8
#define R_SPOINTS 9
#define R_ATTRIBUTES 10
#define R_ATTRDATA 11
#define R_QUIT 100
//static int __isnan__(double x) { return x != x; }
//static int isinf(double x) { return !isnan(x) && isnan(x - x); }
//static int __isinf__(double x) { return fabs(x) > DBL_MAX; }
//static int __isbad(double x) { return __isnan__(x) || __isinf__(x); }
template<typename T>
void ReadCoords(FILE * f,INMOST_DATA_REAL_TYPE c[3])
{
T temp[3];
if( fread(temp,sizeof(T),3,f) == 3 )
{
c[0] = temp[0];
c[1] = temp[1];
c[2] = temp[2];
}
else
{
std::cout << __FILE__ << ":" << __LINE__ << " Cannot read 3 coordinates from file " << std::endl;
throw INMOST::BadFile;
}
}
namespace INMOST
{
std::string Space2Underscore(const std::string & inp)
{
std::string ret = inp;
for (size_t k = 0; k < ret.size(); ++k)
{
if (ret[k] == ' ') ret[k] = '_';
}
return ret;
}
int VtkElementType(ElementType t)
{
switch(t)
{
case Element::Line: return 3;
case Element::Tri: return 5;
case Element::Quad: return 9;
case Element::MultiLine: return 4;
case Element::Polygon: return 7;
case Element::Tet: return 10;
case Element::Hex: return 12;
case Element::Prism: return 13;
case Element::Pyramid: return 14;
case Element::Polyhedron: return 42;
case Element::MultiPolygon: return 42;
}
assert(false);
return -1;
}
INMOST_DATA_ENUM_TYPE VtkElementNodes(ElementType t)
{
switch(t)
{
case Element::Line: return 2;
case Element::Tri: return 3;
case Element::Quad: return 4;
case Element::MultiLine: return ENUMUNDEF;
case Element::Polygon: return ENUMUNDEF;
case Element::Tet: return 4;
case Element::Hex: return 8;
case Element::Prism: return 6;
case Element::Pyramid: return 5;
case Element::Polyhedron: return ENUMUNDEF;
case Element::MultiPolygon: return ENUMUNDEF;
}
assert(false);
return ENUMUNDEF;
}
void Mesh::SaveVTK(std::string File)
{
integer dim = GetDimensions();
if( dim > 3 )
{
printf("VTK file supports 3 dimensions max\n");
return;
}
FILE * f = fopen(File.c_str(),"w");
if( !f ) throw BadFileName;
fprintf(f,"# vtk DataFile Version 3.0\n");
if (this->GetFileOption("VTK_OUTPUT_FACES") == "1") fprintf(f, "VTK_OUTPUT_FACES file is written by INMOST\n");
else fprintf(f,"file is written by INMOST\n");
fprintf(f,"ASCII\n");
fprintf(f,"DATASET UNSTRUCTURED_GRID\n");
bool output_faces = false;
for (INMOST_DATA_ENUM_TYPE k = 0; k < file_options.size(); ++k)
{
if (file_options[k].first == "VTK_OUTPUT_FACES")
{
output_faces = true;
}
}
//ReorderEmpty(CELL | NODE);
Tag set_id = CreateTag("TEMPORARY_ELEMENT_ID",DATA_INTEGER,CELL |FACE| NODE,NONE,1);
Storage::integer cur_num = 0;
for(Mesh::iteratorNode it = BeginNode(); it != EndNode(); ++it) it->IntegerDF(set_id) = cur_num++;
cur_num = 0;
for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); ++it) it->IntegerDF(set_id) = cur_num++;
cur_num = 0;
if (output_faces) for (Mesh::iteratorFace it = BeginFace(); it != EndFace(); ++it) it->IntegerDF(set_id) = cur_num++;
fprintf(f, "POINTS %u double\n", NumberOfNodes());
for(Mesh::iteratorNode it = BeginNode(); it != EndNode(); it++)
{
Storage::real_array coords = it->RealArray(CoordsTag());
for(integer i = 0; i < dim; i++)
{
double temp = coords[i];
fprintf(f,"%.10f ",temp);
}
for(integer i = dim; i < 3; i++)
fprintf(f,"0 ");
fprintf(f,"\n");
}
{
dynarray<int,64> values;
for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++)
{
switch(it->GetGeometricType())
{
case Element::Tri:
case Element::Quad:
//case Element::MultiLine:
//case Element::Polygon:
/*
case Element::Tet:
case Element::Hex:
case Element::Pyramid:
case Element::Prism:
*/
{
ElementArray<Node> nodes = it->getNodes();
if( nodes.size() != VtkElementNodes(it->GetGeometricType()) ) goto safe_output;
values.push_back(static_cast<integer>(nodes.size()));
for(ElementArray<Node>::iterator jt = nodes.begin(); jt != nodes.end(); jt++)
values.push_back(jt->IntegerDF(set_id));
break;
}
case Element::MultiLine:
case Element::Polygon:
{
ElementArray<Node> nodes = it->getNodes();
values.push_back(static_cast<integer>(nodes.size()));
for(ElementArray<Node>::iterator jt = nodes.begin(); jt != nodes.end(); jt++)
values.push_back(jt->IntegerDF(set_id));
break;
}
/*
case Element::Prism:
{
ElementArray<Node> nodes = it->getNodes();
if( nodes.size() != 6 ) goto safe_output;
values.push_back(static_cast<integer>(nodes.size()));
values.push_back(nodes[0].IntegerDF(set_id));
values.push_back(nodes[2].IntegerDF(set_id));
values.push_back(nodes[1].IntegerDF(set_id));
values.push_back(nodes[3].IntegerDF(set_id));
values.push_back(nodes[5].IntegerDF(set_id));
values.push_back(nodes[4].IntegerDF(set_id));
break;
}
*/
/*
case Element::Hex:
{
ElementArray<Node> nodes = it->getNodes();
if( nodes.size() != 8 ) goto safe_output;
values.push_back(static_cast<integer>(nodes.size()));
values.push_back(nodes[0].IntegerDF(set_id));
values.push_back(nodes[3].IntegerDF(set_id));
values.push_back(nodes[2].IntegerDF(set_id));
values.push_back(nodes[1].IntegerDF(set_id));
values.push_back(nodes[4].IntegerDF(set_id));
values.push_back(nodes[7].IntegerDF(set_id));
values.push_back(nodes[6].IntegerDF(set_id));
values.push_back(nodes[5].IntegerDF(set_id));
break;
}
case Element::Pyramid:
{
ElementArray<Node> nodes = it->getNodes();
if( nodes.size() != 5 ) goto safe_output;
values.push_back(static_cast<integer>(nodes.size()));
values.push_back(nodes[0].IntegerDF(set_id));
values.push_back(nodes[3].IntegerDF(set_id));
values.push_back(nodes[2].IntegerDF(set_id));
values.push_back(nodes[1].IntegerDF(set_id));
values.push_back(nodes[4].IntegerDF(set_id));
break;
}
*/
case Element::Tet:
case Element::Hex:
case Element::Pyramid:
case Element::Prism:
case Element::Polyhedron:
case Element::MultiPolygon:
{
safe_output:
//printf("polyhedron!!!\n");
ElementArray<Face> faces = it->getFaces();
integer totalNum = 1 + static_cast<integer>(faces.size());
for(ElementArray<Face>::iterator jt = faces.begin(); jt != faces.end(); jt++)
totalNum += jt->nbAdjElements(NODE);
values.push_back(totalNum);
values.push_back(static_cast<integer>(faces.size()));
for(ElementArray<Face>::iterator jt = faces.begin(); jt != faces.end(); jt++)
{
ElementArray<Node> nodes = jt->getNodes();
values.push_back(static_cast<integer>(nodes.size()));
if( jt->FaceOrientedOutside(it->self()) )
for(ElementArray<Node>::iterator kt = nodes.begin(); kt != nodes.end(); kt++)
values.push_back(kt->IntegerDF(set_id));
else
for(ElementArray<Node>::reverse_iterator kt = nodes.rbegin(); kt != nodes.rend(); kt++)
values.push_back(kt->IntegerDF(set_id));
}
break;
}
default: printf("This should not happen %s\n",Element::GeometricTypeName(it->GetGeometricType()));
}
}
if (output_faces)
{
for (Mesh::iteratorFace it = BeginFace(); it != EndFace(); it++)
{
switch (it->GetGeometricType())
{
case Element::Line:
{
ElementArray<Node> nodes = it->getNodes();
values.push_back(static_cast<integer>(nodes.size()));
for (ElementArray<Node>::iterator jt = nodes.begin(); jt != nodes.end(); jt++)
values.push_back(jt->IntegerDF(set_id));
break;
}
case Element::Tri:
case Element::Quad:
case Element::MultiLine:
case Element::Polygon:
{
ElementArray<Node> nodes = it->getNodes();
values.push_back(static_cast<integer>(nodes.size()));
for (ElementArray<Node>::iterator jt = nodes.begin(); jt != nodes.end(); jt++)
values.push_back(jt->IntegerDF(set_id));
break;
}
default: printf("This should not happen %s\n", Element::GeometricTypeName(it->GetGeometricType()));
}
}
fprintf(f, "CELLS %u %d\n", NumberOfCells() + NumberOfFaces(), (int)values.size());
}
else
fprintf(f, "CELLS %u %d\n", NumberOfCells(), (int)values.size());
for(dynarray<Storage::integer,64>::size_type i = 0; i < values.size(); i++)
{
fprintf(f,"%d ",values[i]);
if( (i+1) % 20 == 0) fprintf(f,"\n");
}
fprintf(f,"\n");
}
if (output_faces){
fprintf(f, "CELL_TYPES %u\n", NumberOfCells() + NumberOfFaces());
}
else
fprintf(f, "CELL_TYPES %u\n", NumberOfCells());
for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++)
{
INMOST_DATA_ENUM_TYPE nnodes = VtkElementNodes(it->GetGeometricType());
//if( nnodes == ENUMUNDEF || nnodes == it->nbAdjElements(NODE) ) //nodes match - output correct type
// fprintf(f,"%d\n",VtkElementType(it->GetGeometricType()));
//else //number of nodes mismatch with expected - some topology checks must be off
fprintf(f,"%d\n",VtkElementType(Element::MultiPolygon));
}
if (output_faces){
for (Mesh::iteratorFace it = BeginFace(); it != EndFace(); it++)
{
INMOST_DATA_ENUM_TYPE nnodes = VtkElementNodes(it->GetGeometricType());
if (nnodes == ENUMUNDEF || nnodes == it->nbAdjElements(NODE)) //nodes match - output correct type
fprintf(f, "%d\n", VtkElementType(it->GetGeometricType()));
else //number of nodes mismatch with expected - some topology checks must be off
fprintf(f, "%d\n", VtkElementType(Element::MultiPolygon));
}
}
DeleteTag(set_id);
{
std::vector<std::string> tag_names;
std::vector<Tag> tags;
ListTagNames(tag_names);
for(unsigned int i = 0; i < tag_names.size(); i++)
{
Tag t = GetTag(tag_names[i]);
//printf("%s %d %d %d\n",tag_names[i].c_str(),t.isDefined(CELL),!t.isSparse(CELL),t.GetDataType() != DATA_BULK);
if (((t.isDefined(CELL) /*&& !t.isSparse(CELL)*/)
|| (t.isDefined(FACE) && output_faces)) &&
t.GetPrint() && //Temporary solution: @see Mesh::file_option
t.GetDataType() != DATA_BULK &&
t.GetDataType() != DATA_REFERENCE &&
t.GetDataType() != DATA_REMOTE_REFERENCE &&
t != CoordsTag() &&
t != SharedTag() &&
t != SendtoTag() &&
t != ProcessorsTag())
{
//printf("added!\n");
tags.push_back(t);
}
}
if (!tags.empty() && output_faces)
{
fprintf(f, "CELL_DATA %u\n", NumberOfCells() + NumberOfFaces());
}
else if (!tags.empty())fprintf(f, "CELL_DATA %u\n", NumberOfCells());
for(unsigned int i = 0; i < tags.size(); i++)
{
unsigned int comps = tags[i].GetSize();
if( comps == ENUMUNDEF )
{
//printf("Warning: vtk don't support arrays of variable size (tag name: %s)\n",tags[i].GetTagName().c_str());
continue;
}
else
{
{
std::string type_str = "int";
if( tags[i].GetDataType() == DATA_REAL
#if defined(USE_AUTODIFF)
|| tags[i].GetDataType() == DATA_VARIABLE
#endif
) type_str = "double";
fprintf(f,"SCALARS %s %s %d\n",Space2Underscore(tags[i].GetTagName()).c_str(),type_str.c_str(),comps);
fprintf(f,"LOOKUP_TABLE default\n");
for (Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++)
{
switch (tags[i].GetDataType())
{
case DATA_REAL:
{
if (tags[i].isDefined(CELL) && it->HaveData(tags[i]))
{
Storage::real_array arr = it->RealArray(tags[i]);
for (unsigned int m = 0; m < comps; m++)
{
double val = static_cast<double>(arr[m]);
fprintf(f, "%14e ", __isbad(val) ? -0.9999E30 : val);
}
}
else for (unsigned int m = 0; m < comps; m++) fprintf(f, "%14e ", -0.9999E30);
fprintf(f, "\n");
}
break;
case DATA_INTEGER:
{
if (tags[i].isDefined(CELL) && it->HaveData(tags[i]))
{
Storage::integer_array arr = it->IntegerArray(tags[i]);
for (unsigned int m = 0; m < comps; m++) fprintf(f, "%d ", arr[m]);
}
else for (unsigned int m = 0; m < comps; m++) fprintf(f, "%d ",INT_MIN);
fprintf(f, "\n");
}
break;
#if defined(USE_AUTODIFF)
case DATA_VARIABLE:
{
if (tags[i].isDefined(CELL) && it->HaveData(tags[i]))
{
Storage::var_array arr = it->VariableArray(tags[i]);
for (unsigned int m = 0; m < comps; m++)
{
double val = static_cast<double>(arr[m].GetValue());
fprintf(f, "%14e ", __isbad(val) ? -0.9999E30 : val);
}
}
else for (unsigned int m = 0; m < comps; m++) fprintf(f, "%14e ", -0.9999E30);
fprintf(f, "\n");
}
break;
#endif
default: continue;
}
}
if (output_faces)
{
for (Mesh::iteratorFace it = BeginFace(); it != EndFace(); it++)
{
switch (tags[i].GetDataType())
{
case DATA_REAL: