Commit 843e751c authored by Dmitry Bagaev's avatar Dmitry Bagaev
Browse files

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

parents f16a9d90 bcca5599
...@@ -4896,17 +4896,21 @@ void svg_draw(std::ostream & file) ...@@ -4896,17 +4896,21 @@ void svg_draw(std::ostream & file)
if( oclipper ) if( oclipper )
{ {
//std::vector<face2gl> temp_boundary; /*
//oclipper->gen_clip(temp_boundary,n); std::vector<face2gl> temp_boundary;
//face2gl::radix_sort_dist(temp_boundary); oclipper->gen_clip(temp_boundary,n);
//svg_draw_faces(file,temp_boundary,modelview,projection,viewport); for(INMOST_DATA_ENUM_TYPE q = 0; q < temp_boundary.size() ; q++)
temp_boundary[q].compute_dist(campos);
face2gl::radix_sort_dist(temp_boundary);
svg_draw_faces(file,temp_boundary,modelview,projection,viewport);
*/
std::vector<face2gl> sorted_clip_boundary(clip_boundary); std::vector<face2gl> sorted_clip_boundary(clip_boundary);
for(INMOST_DATA_ENUM_TYPE q = 0; q < sorted_clip_boundary.size() ; q++) for(INMOST_DATA_ENUM_TYPE q = 0; q < sorted_clip_boundary.size() ; q++)
sorted_clip_boundary[q].compute_dist(campos); sorted_clip_boundary[q].compute_dist(campos);
face2gl::radix_sort_dist(sorted_clip_boundary); face2gl::radix_sort_dist(sorted_clip_boundary);
svg_draw_faces(file,sorted_clip_boundary,modelview,projection,viewport); svg_draw_faces(file,sorted_clip_boundary,modelview,projection,viewport);
//file << "<g stroke=\"black\">" << std::endl; //file << "<g stroke=\"black\">" << std::endl;
//svg_draw_edges(file,sorted_clip_boundary,modelview,projection,viewport); //svg_draw_edges(file,sorted_clip_boundary,modelview,projection,viewport);
//file << "</g>" << std::endl; //file << "</g>" << std::endl;
......
...@@ -298,10 +298,11 @@ namespace INMOST ...@@ -298,10 +298,11 @@ namespace INMOST
flag = 0; flag = 0;
break; break;
} }
if( l == 0 ) break;
if (fabs(get_value(Sigma(nm,nm))) + anorm == anorm) if (fabs(get_value(Sigma(nm,nm))) + anorm == anorm)
break; break;
} }
if (flag) if (flag && l != 0)
{ {
c = 0.0; c = 0.0;
s = 1.0; s = 1.0;
...@@ -377,7 +378,7 @@ namespace INMOST ...@@ -377,7 +378,7 @@ namespace INMOST
} }
z = pythag(f, h); z = pythag(f, h);
Sigma(j,j) = z; Sigma(j,j) = z;
if (z) if (get_value(z))
{ {
z = 1.0 / z; z = 1.0 / z;
c = f * z; c = f * z;
...@@ -1115,6 +1116,19 @@ namespace INMOST ...@@ -1115,6 +1116,19 @@ namespace INMOST
n = _n; n = _n;
m = _m; m = _m;
} }
Matrix PseudoInvert(INMOST_DATA_REAL_TYPE tol = 0)
{
Matrix U,S,V;
SVD(U,S,V);
for(int k = 0; k < S.Cols(); ++k)
{
if( S(k,k) > tol )
S(k,k) = 1.0/S(k,k);
else
S(k,k) = 0.0;
}
return V*S*U.Transpose();
}
}; };
typedef Matrix<INMOST_DATA_REAL_TYPE> rMatrix; //shortcut for real matrix typedef Matrix<INMOST_DATA_REAL_TYPE> rMatrix; //shortcut for real matrix
......
...@@ -1501,35 +1501,44 @@ namespace INMOST ...@@ -1501,35 +1501,44 @@ namespace INMOST
template<class A> template<class A>
class pow_const_expression : public shell_expression<pow_const_expression<A> > class pow_const_expression : public shell_expression<pow_const_expression<A> >
{ {
const A & left; const A & left;
INMOST_DATA_REAL_TYPE value, ldmult; INMOST_DATA_REAL_TYPE value, ldmult, ldmult2;
public: public:
pow_const_expression(const shell_expression<A> & pleft, INMOST_DATA_REAL_TYPE pright) : left(pleft) pow_const_expression(const shell_expression<A> & pleft, INMOST_DATA_REAL_TYPE pright) : left(pleft)
{ {
INMOST_DATA_REAL_TYPE lval = left.GetValue(); INMOST_DATA_REAL_TYPE lval = left.GetValue();
value = ::pow(lval,pright); value = ::pow(lval,pright);
if( lval != 0 ) if( lval != 0 )
ldmult = value * pright / lval; {
else ldmult = value * pright / lval;
ldmult = 0; ldmult2 = ldmult * (pright-1) / lval;
} }
pow_const_expression(const pow_const_expression & other) else
:left(other.left), value(other.value), ldmult(other.ldmult) {} ldmult = ldmult2 = 0;
__INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; } }
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const pow_const_expression(const pow_const_expression & other)
{ :left(other.left), value(other.value), ldmult(other.ldmult) {}
left.GetJacobian(mult*ldmult,r); __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
} __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const {
{ left.GetJacobian(mult*ldmult,r);
left.GetJacobian(mult*ldmult,r); }
} __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
__INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const {
{ left.GetJacobian(mult*ldmult,r);
throw NotImplemented; }
} __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
}; {
//(F(G))'' = (F'(G)G')' = (F''(G)G'+F'(G)G'')
//(g(x)^n)'' = n g(x)^(n - 2) (g(x) g''(x) + (n - 1) g'(x)^2)
Sparse::Row JL;
Sparse::HessianRow HL;
left.GetHessian(1,JL,1,HL); //retrive jacobian row and hessian matrix of the left expression
Sparse::HessianRow::MergeJacobianHessian(multH*ldmult2,JL,JL,multH*ldmult,HL,H);
for(Sparse::Row::iterator it = JL.Begin(); it != JL.End(); ++it) it->second *= ldmult*multJ;
}
};
template<class A> template<class A>
class const_pow_expression : public shell_expression<const_pow_expression<A> > class const_pow_expression : public shell_expression<const_pow_expression<A> >
......
...@@ -22,6 +22,26 @@ namespace INMOST ...@@ -22,6 +22,26 @@ namespace INMOST
/// @see Sparse::Solve /// @see Sparse::Solve
class Solver class Solver
{ {
public:
//Backward-compatibility
typedef std::string Type;
static const Type INNER_ILU2; ///< inner Solver based on BiCGStab(L) solver with second order ILU factorization as preconditioner.
static const Type INNER_DDPQILUC; ///< inner Solver based on BiCGStab(L) solver with second order Crout-ILU with inversed-based condition estimation and unsymmetric reordering for diagonal dominance as preconditioner.
static const Type INNER_MPTILUC; ///< inner Solver based on BiCGStab(L) solver with second order Crout-ILU with inversed-based condition estimation and maximum product transversal reordering as preconditioner.
static const Type INNER_MPTILU2; ///< inner Solver based on BiCGStab(L) solver with second order ILU and maximum product transversal reordering as preconditioner.
static const Type Trilinos_Aztec; ///< external Solver AztecOO from Trilinos package.
static const Type Trilinos_Belos; ///< external Solver Belos from Trilinos package, currently without preconditioner.
static const Type Trilinos_ML; ///< external Solver AztecOO with ML preconditioner.
static const Type Trilinos_Ifpack;///< external Solver AztecOO with Ifpack preconditioner.
static const Type PETSc; ///< external Solver PETSc, @see http://www.mcs.anl.gov/petsc/
static const Type ANI; ///< external Solver from ANI3D based on ILU2 (sequential Fortran version), @see http://ani3d.sourceforge.net/
static const Type FCBIILU2; ///< external FCBIILU2 Solver (BIILU2 parallel F2C version).
static const Type K3BIILU2; ///< inner K3BIILU2 Solver (BIILU2 parallel version).
static const Type SUPERLU; ///< external Solver SuperLU @see https://github.com/starseeker/SuperLU
void SetParameterEnum(std::string name, INMOST_DATA_ENUM_TYPE value);
void SetParameterReal(std::string name, INMOST_DATA_REAL_TYPE value);
std::string GetReason() {return ReturnReason();}
//Backward-compatibility
private: private:
static std::vector<SolverParameters> parameters; static std::vector<SolverParameters> parameters;
static int *argc; static int *argc;
...@@ -196,7 +216,7 @@ namespace INMOST ...@@ -196,7 +216,7 @@ namespace INMOST
/// Trilinos_ML: trilinos_ml_options.xml /// Trilinos_ML: trilinos_ml_options.xml
/// Trilinos_Aztec: trilinos_aztec_options.xml /// Trilinos_Aztec: trilinos_aztec_options.xml
/// Trilinos_Belos: trilinos_belos_options.xml /// Trilinos_Belos: trilinos_belos_options.xml
static void Initialize(int *argc, char ***argv, const char *database); static void Initialize(int *argc, char ***argv, const char *database = NULL);
/// Finalize the stage of parallel solution. /// Finalize the stage of parallel solution.
/// If MPI was initialized in Solver::Initialize, then it will be finalized. /// If MPI was initialized in Solver::Initialize, then it will be finalized.
......
...@@ -104,7 +104,8 @@ namespace INMOST ...@@ -104,7 +104,8 @@ namespace INMOST
FILE * f = fopen(File.c_str(),"w"); FILE * f = fopen(File.c_str(),"w");
if( !f ) throw BadFileName; if( !f ) throw BadFileName;
fprintf(f,"# vtk DataFile Version 3.0\n"); fprintf(f,"# vtk DataFile Version 3.0\n");
fprintf(f,"file is written by INMOST\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,"ASCII\n");
fprintf(f,"DATASET UNSTRUCTURED_GRID\n"); fprintf(f,"DATASET UNSTRUCTURED_GRID\n");
...@@ -633,6 +634,7 @@ safe_output: ...@@ -633,6 +634,7 @@ safe_output:
case R_USERDATA: case R_USERDATA:
{ {
//printf("Attached info: %s\n",readline); //printf("Attached info: %s\n",readline);
if (!strncmp(readline, "VTK_OUTPUT_FACES", 16)) this->SetFileOption("VTK_OUTPUT_FACES", "1");
state = R_DATATYPE; state = R_DATATYPE;
break; break;
} }
......
...@@ -99,7 +99,7 @@ namespace INMOST ...@@ -99,7 +99,7 @@ namespace INMOST
dynarray< char , 256 > hide_column; dynarray< char , 256 > hide_column;
dynarray< char , 256 > hide_row; dynarray< char , 256 > hide_row;
dynarray< char , 256 > stub_row; dynarray< char , 256 > stub_row;
//std::vector< std::pair<std::vector<int>,double> > remember; std::vector< std::pair<std::vector<int>,double> > remember;
double min_loop_measure; double min_loop_measure;
bool print; bool print;
AbstractCoords * coords; AbstractCoords * coords;
...@@ -293,37 +293,61 @@ namespace INMOST ...@@ -293,37 +293,61 @@ namespace INMOST
} }
} while(true); } while(true);
data.RemPrivateMarker(mrk); data.RemPrivateMarker(mrk);
mesh->ReleasePrivateMarker(mrk); Storage::real fnrm[3], fcnt[3], ccnt[3] = {0,0,0}, nnodes = 0;
Storage::real nrm[3], cnt[3]; Storage::real_array x0,a,b;
Storage::real_array v0,v1,v2; //find cell centroid
for(typename ElementArray<T>::size_type j = 0; j < data.size(); j++)
{
ElementArray<Node> nodes = data[j]->getAsFace()->getNodes(mrk,true);
nodes.SetPrivateMarker(mrk);
for(ElementArray<Node>::iterator qt = nodes.begin(); qt != nodes.end(); ++qt)
{
x0 = coords->Get(qt->self());
ccnt[0] += x0[0];
ccnt[1] += x0[1];
ccnt[2] += x0[2];
nnodes += 1;
}
}
ccnt[0] /= nnodes;
ccnt[1] /= nnodes;
ccnt[2] /= nnodes;
for(typename ElementArray<T>::size_type j = 0; j < data.size(); j++) for(typename ElementArray<T>::size_type j = 0; j < data.size(); j++)
{ {
//compute normal to face //compute normal to face
ElementArray<Node> nodes = data[j]->getAsFace()->getNodes(); ElementArray<Node> nodes = data[j]->getAsFace()->getNodes();
nrm[0] = nrm[1] = nrm[2] = 0; nodes.RemPrivateMarker(mrk);
v0 = v1 = coords->Get(nodes[0]); fnrm[0] = fnrm[1] = fnrm[2] = 0;
cnt[0] = v1[0]; fcnt[0] = fcnt[1] = fcnt[2] = 0;
cnt[1] = v1[1]; nnodes = 0;
cnt[2] = v1[2]; x0 = coords->Get(nodes[0]);
a = x0;
for(unsigned i = 0; i < nodes.size(); i++) for(unsigned i = 0; i < nodes.size(); i++)
{ {
v2 = coords->Get(nodes[(i+1)%nodes.size()]); b = coords->Get(nodes[(i+1)%nodes.size()]);
cnt[0] += v2[0]; fnrm[0] += (a[1]-x0[1])*(b[2]-x0[2]) - (a[2]-x0[2])*(b[1]-x0[1]);
cnt[1] += v2[1]; fnrm[1] += (a[2]-x0[2])*(b[0]-x0[0]) - (a[0]-x0[0])*(b[2]-x0[2]);
cnt[2] += v2[2]; fnrm[2] += (a[0]-x0[0])*(b[1]-x0[1]) - (a[1]-x0[1])*(b[0]-x0[0]);
nrm[0] += (v1[1]-v0[1])*(v2[2]-v0[2]) - (v1[2]-v0[2])*(v2[1]-v0[1]);
nrm[1] += (v1[2]-v0[2])*(v2[0]-v0[0]) - (v1[0]-v0[0])*(v2[2]-v0[2]); fcnt[0] += a[0];
nrm[2] += (v1[0]-v0[0])*(v2[1]-v0[1]) - (v1[1]-v0[1])*(v2[0]-v0[0]); fcnt[1] += a[1];
v1.swap(v2); fcnt[2] += a[2];
a.swap(b);
nnodes += 1;
} }
nrm[0] *= 0.5; fnrm[0] *= 0.5;
nrm[1] *= 0.5; fnrm[1] *= 0.5;
nrm[2] *= 0.5; fnrm[2] *= 0.5;
cnt[0] /= (Storage::real)nodes.size(); fcnt[0] /= nnodes;
cnt[1] /= (Storage::real)nodes.size(); fcnt[1] /= nnodes;
cnt[2] /= (Storage::real)nodes.size(); fcnt[2] /= nnodes;
measure += (data[j]->GetPrivateMarker(rev) ? -1.0 : 1.0)*dot_prod(cnt,nrm); for(int r = 0; r < 3; ++r)
fcnt[r] = fcnt[r]-ccnt[r];
measure += (data[j]->GetPrivateMarker(rev) ? -1.0 : 1.0)*dot_prod(fcnt,fnrm);
} }
mesh->ReleasePrivateMarker(mrk);
data.RemPrivateMarker(rev); data.RemPrivateMarker(rev);
mesh->ReleasePrivateMarker(rev); mesh->ReleasePrivateMarker(rev);
measure /= 3.0; measure /= 3.0;
...@@ -447,7 +471,7 @@ namespace INMOST ...@@ -447,7 +471,7 @@ namespace INMOST
std::cout << " " << cnt[0] << " " << cnt[1] << " " << cnt[2]; std::cout << " " << cnt[0] << " " << cnt[1] << " " << cnt[2];
std::cout << std::endl; std::cout << std::endl;
} }
/*
std::cout << "loops [" << remember.size() << "]:" << std::endl; std::cout << "loops [" << remember.size() << "]:" << std::endl;
for(size_t k = 0; k < remember.size(); ++k) for(size_t k = 0; k < remember.size(); ++k)
{ {
...@@ -456,7 +480,7 @@ namespace INMOST ...@@ -456,7 +480,7 @@ namespace INMOST
std::cout << remember[k].first[l] << ", "; std::cout << remember[k].first[l] << ", ";
std::cout << remember[k].first.back() << " measure " << remember[k].second << std::endl; std::cout << remember[k].first.back() << " measure " << remember[k].second << std::endl;
} }
*/
if( GetHandleElementType(head_column[0]) == EDGE ) if( GetHandleElementType(head_column[0]) == EDGE )
{ {
std::cout << "edges [" << head_column.size() << "]" << std::endl; std::cout << "edges [" << head_column.size() << "]" << std::endl;
...@@ -582,7 +606,7 @@ namespace INMOST ...@@ -582,7 +606,7 @@ namespace INMOST
if( min_loop.empty() ) if( min_loop.empty() )
{ {
if( print ) std::cout << "abandon " << first << std::endl; if( print ) std::cout << "abandon " << first << std::endl;
//remember.push_back(std::make_pair(std::vector<int>(1,first),-1.0)); remember.push_back(std::make_pair(std::vector<int>(1,first),-1.0));
visits[first]--; //don't start again from this element visits[first]--; //don't start again from this element
} }
} }
...@@ -595,8 +619,8 @@ namespace INMOST ...@@ -595,8 +619,8 @@ namespace INMOST
if( !ret.empty() ) if( !ret.empty() )
{ {
//std::vector<int> add_remember; std::vector<int> add_remember;
//add_remember.reserve(min_loop.size()); add_remember.reserve(min_loop.size());
MarkerType hide_marker = mesh->CreatePrivateMarker(); MarkerType hide_marker = mesh->CreatePrivateMarker();
for(typename ElementArray<T>::size_type k = 0; k < ret.size(); k++) mesh->SetPrivateMarker(ret.at(k),hide_marker); for(typename ElementArray<T>::size_type k = 0; k < ret.size(); k++) mesh->SetPrivateMarker(ret.at(k),hide_marker);
if( print ) std::cout << "return loop [" << ret.size() << "]:"; if( print ) std::cout << "return loop [" << ret.size() << "]:";
...@@ -604,13 +628,13 @@ namespace INMOST ...@@ -604,13 +628,13 @@ namespace INMOST
if( mesh->GetPrivateMarker(head_column[k],hide_marker) ) if( mesh->GetPrivateMarker(head_column[k],hide_marker) )
{ {
visits[k]--; visits[k]--;
//add_remember.push_back((int)k); add_remember.push_back((int)k);
if( print ) std::cout << k << " "; if( print ) std::cout << k << " ";
} }
if( print ) std::cout << std::endl; if( print ) std::cout << std::endl;
for(typename ElementArray<T>::size_type k = 0; k < ret.size(); k++) mesh->RemPrivateMarker(ret.at(k),hide_marker); for(typename ElementArray<T>::size_type k = 0; k < ret.size(); k++) mesh->RemPrivateMarker(ret.at(k),hide_marker);
mesh->ReleasePrivateMarker(hide_marker); mesh->ReleasePrivateMarker(hide_marker);
//remember.push_back(std::make_pair(add_remember,min_loop_measure)); remember.push_back(std::make_pair(add_remember,min_loop_measure));
return true; return true;
} }
return false; return false;
......
...@@ -1251,6 +1251,20 @@ namespace INMOST ...@@ -1251,6 +1251,20 @@ namespace INMOST
if (!loop.empty()) ret.push_back(m->CreateCell(loop).first); if (!loop.empty()) ret.push_back(m->CreateCell(loop).first);
else break; else break;
} while( true ); } while( true );
/* //debug
if(!mat.all_visited())
{
mat.print_matrix();
incident_matrix<Face> mat(m,temp.begin(),temp.end(),ninner,GridCoords(),true);
do
{
mat.find_shortest_loop(loop);
if (!loop.empty()) ret.push_back(m->CreateCell(loop).first);
else break;
} while( true );
}
*/
return ret; return ret;
} }
......
...@@ -46,6 +46,21 @@ ...@@ -46,6 +46,21 @@
#endif #endif
namespace INMOST { namespace INMOST {
//initialization of backward-compatibility constants
const Solver::Type Solver::INNER_ILU2 = "inner_ilu2";
const Solver::Type Solver::INNER_DDPQILUC = "inner_ddpqiluc2";
const Solver::Type Solver::INNER_MPTILUC = "inner_mptiluc";
const Solver::Type Solver::INNER_MPTILU2 = "inner_mptilu2";
const Solver::Type Solver::Trilinos_Aztec = "trilinos_aztec";
const Solver::Type Solver::Trilinos_Belos = "trilinos_belos";
const Solver::Type Solver::Trilinos_ML = "trilinos_ml";
const Solver::Type Solver::Trilinos_Ifpack = "trilinos_ifpack";
const Solver::Type Solver::PETSc = "petsc";
const Solver::Type Solver::ANI = "inner_ilu2";
const Solver::Type Solver::FCBIILU2 = "fcbiilu2";
const Solver::Type Solver::K3BIILU2 = "k3biilu2";
const Solver::Type Solver::SUPERLU = "superlu";
SolverInterface *SolverMaster::getSolver(std::string name) { SolverInterface *SolverMaster::getSolver(std::string name) {
if (name == "inner_ilu2") return new SolverILU2(); if (name == "inner_ilu2") return new SolverILU2();
......
...@@ -10,6 +10,9 @@ namespace INMOST { ...@@ -10,6 +10,9 @@ namespace INMOST {
bool Solver::is_initialized = false; bool Solver::is_initialized = false;
bool Solver::is_finalized = false; bool Solver::is_finalized = false;
std::vector<SolverParameters> Solver::parameters = std::vector<SolverParameters>(); std::vector<SolverParameters> Solver::parameters = std::vector<SolverParameters>();
void Solver::SetParameterEnum(std::string name, INMOST_DATA_ENUM_TYPE value) {SetParameter(name,to_string(value));}
void Solver::SetParameterReal(std::string name, INMOST_DATA_REAL_TYPE value) {SetParameter(name,to_string(value));}
Solver::Solver(std::string solverName, std::string prefix, INMOST_MPI_Comm _comm) { Solver::Solver(std::string solverName, std::string prefix, INMOST_MPI_Comm _comm) {
std::string lowerName = string_to_lower(solverName); std::string lowerName = string_to_lower(solverName);
......
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