00001 #ifndef _MY_NumericMethods_
00002 #define _MY_NumericMethods_
00003 #include "globals.h"
00004 #include "sfunctions.h"
00005 #include <lac/sparse_matrix.h>
00006 #include <lac/sparse_direct.h>
00007
00008
00009 namespace NumericMethods {
00010
00011 void fillVector(VecDouble &v,double number);
00012 double interpolateLinear1D(double X,double Px1,double Py1,double Px2,double Py2);
00013 double interpolateLinear1D(double X,double Px1,double Py1,double Px2,double Py2,double Px3,double Py3);
00014
00015 double interpolateRectLin2D(Point2D &X,Point2D &BL,Point2D &UR,double Z1,double Z2, double Z3, double Z4);
00016 inline bool isInRectangle(const Point2D &X,const Point2D &P1,const Point2D &P2){return (X(0) <= P2(0) && X(0)>=P1(0) && X(1) <= P2(1) && X(1)>= P1(1));}
00017 inline bool isInRectangle(const Point1D &X,const Point1D &P1,const Point1D &P2){return (X(0) <= P2(0) && X(0)>=P1(0));}
00018 bool isInCube(const Point3D &X,const Point3D &P1,const Point3D &P2);
00019 bool isInCube(const VecDouble &X,const VecDouble &P1,const VecDouble &P2);
00020
00021 void rotate90(Point2D &X);
00022 void rotate270(Point2D &X);
00023 double intThreePointsLinear(double P1,double Sw1,double P2,double Sw2, double P3,double Sw3, double Xstart,double Xend);
00024
00025
00026
00027 extern unsigned CellDirectionToCellNeighborsIndexTbl[4][2];
00028 void solve3x3(double M[3][4]);
00029
00030 void multiplyEachElement(VecDouble &values,const VecDouble &v1,const VecDouble &v2);
00031 void divideEachElement(VecDouble &values,const VecDouble &v1,const VecDouble &v2);
00032
00033
00034
00035 double _div(double a,double b);
00036
00037 double minMod(double a,double b,double c);
00038 double minMod(double a,double b);
00039 double maxMod(double a,double b,double c);
00040 double maxMod(double a,double b);
00041 double min(double a,double b,double c);
00042
00043 inline double insideInterval(double x,double a,double b){return ( (x >= a && x <= b) || (x <= a && x >= b));}
00044
00045
00046 bool isNearZero(double dd);
00047
00048
00049
00050
00051 void IterateOrder1EDOEuler(double &u,unsigned numIteracoes,double dTau,Function1D &f,double c);
00052 void IterateOrder1EDOSystem(double &u,double &x,unsigned numIteracoes,double dTau,Function1D &fV,double por,double v,Function1D &fR,double r);
00053 void normalizePoint2D (Point2D &p);
00054 void StlVectorDoubleToVecDouble(const std::vector<double> &vV,VecDouble &vValues);
00055 void readVecDouble(std::ifstream &fIn,VecDouble &vValues);
00056 void writeVecDouble(std::ofstream &fOut,VecDouble &vValues);
00057 void getSparseMatrixWithoutZeroLines(SparseMatrix<double> &M,SparseMatrix<double> &C,SparsityPattern &spStripped);
00058 void getTranspose(SparseMatrix<double> &M,SparseMatrix<double> &MT,SparsityPattern &patternMT);
00059 void getInvMatrix(unsigned dim,SparseDirectUMFPACK &umf_solver,SparseMatrix<double> &A,SparsityPattern &patternA);
00060 void makeSumPattern(SparseMatrix<double> &A,SparseMatrix<double> &B,SparsityPattern &patternS);
00061 bool equal(const VecDouble &V1,const VecDouble &V2,double tol=TOLERANCE);
00062
00063 double vectorMinValue(const VecDouble &v);
00064 double vectorMaxValue(const VecDouble &v);
00065 void vectorMinMaxValues(const VecDouble &v,double &min,double &max);
00066 double vectorMaxAbsValue(const VecDouble &v);
00067 double vectorMinAbsValue(const VecDouble &v);
00068 double vectorSum(const VecDouble &v);
00069 void vectorBounds(const VecDouble &values,double &min,double &max);
00070 void vectorModBounds(const VecDouble &values,double &min,double &media,double &max);
00071 void vectorFill(VecDouble &vec,double start,double inc);
00072 void vectorDivideEntries(VecDouble &V1,const VecDouble &vQuot);
00073 void vectorDivideEntriesAvoidingNAN(VecDouble &V1,const VecDouble &vQuot);
00074 double vectorMaxDiff(const VecDouble &v1,const VecDouble &v2);
00075 double vectorRelError(const VecDouble &v1,const VecDouble &v2);
00076 double vectorMaxRelDiff(const VecDouble &v1,const VecDouble &v2);
00077
00078 double matrixEfficience(SparseMatrix<double> &A,double tol=1.0e-11);
00079 void matrixMultiply(SparseMatrix<double> &A,SparseMatrix<double> &B,SparseMatrix<double> &C,SparsityPattern &patternC);
00080 void matrixSum(SparseMatrix<double> &A,SparseMatrix<double> &B,SparseMatrix<double> &C);
00081 void matrixFill(SparseMatrix<double> &A,double start,double inc) ;
00082 void matrixFill(Matrix &A,double value);
00083 void matrixGetColumn(VecDouble &v, const Matrix &M, int col);
00084
00085
00086 void multiplySubMatrix(const SparseMatrix<double> &C,std::vector<unsigned> &equations,std::vector<unsigned> °s,VecDouble &dst,const VecDouble &src);
00087 void multiplyBlockSubMatrix(const SparseMatrix<double> &C,unsigned sRow,unsigned eRow,unsigned sCol,unsigned eCols,VecDouble &dst,const VecDouble &src);
00088
00089
00090 bool IsSymetric(const SparseMatrix<double> &C,double tol=1E-08);
00091 void checkLinearSolution(const SparseMatrix<double> &C,const VecDouble &vSol,const VecDouble &vB,double &min,double &media,double &max);
00092 void getMinMaxValueByColumn(Matrix &M,VecDouble &vMinValues,VecDouble &vMaxValues);
00093 void getMinMaxValueByColumn(Matrix &M,VecDouble &vMinValues,VecDouble &vMaxValues,VecDouble &vIndexMin,VecDouble &vIndexMax);
00094
00095
00096 void printPoint(Point3D p);
00097 void printVectorWithIndices(const VecDouble &vV,double tol=0.0,bool printZeros=false);
00098 void printVectorWithIndices(const BlockVecDouble &vV,double tol=0.0);
00099 void printVectorWithIndices(const VecTag &vV,double tol=0.0,bool printZeros=false);
00100 void printVectorWithIndices(const VecIndex &vI,std::ostream &out=std::cout);
00101
00102 void printVector(const VecDouble &vV,bool bASCII);
00103 void printVectorWithIndices(const VecTag &vV);
00104
00105
00106
00107 void printVector(const VecDouble &vV,double tol=0.0);
00108 void printVector(const BlockVecDouble &vV,double tol);
00109 void printMatrix(const Matrix &M);
00110 void printMatrix(std::string file,const SparseMatrix<double> &M);
00111
00112 void printSparseMatrixOctave(std::ostream &out,const SparseMatrix<double> &M);
00113 void printNonZeroEntriesPerLine(const SparseMatrix<double> &M);
00114 double triangleMeasure(const Point3D &u,const Point3D &v);
00115 double squareMeasure(const Point3D &p1,const Point3D &p2,const Point3D &p3,const Point3D &p4);
00116 bool a_equal(double a, double b,double tol);
00117 void adjustBounds(double &c,const double a,const double b);
00118 void printPrescBoundMapping(MapIntDouble &prescVelMap);
00119 void getSubVector(const VecDouble &values,const VecIndex &vI,VecDouble &subV);
00120
00121 void VecDoubleToSTL(const VecDouble &v1,std::vector<double> &v2);
00122 void STLToVecDouble(VecDouble &v1,const std::vector<double> &v2);
00123
00124 void getCompactRowFormat(SparseMatrix<double> &A,std::vector<int> &cnt,std::vector<int> &col, std::vector<double> &ele);
00125
00126 void addMapValues(MapIntDouble &map,VecDouble &values);
00127
00128
00129
00130 double det3x3(const Matrix &M);
00131 double abs_vector_product(const VecDouble &v1,const VecDouble &v2);
00132 };
00133
00134 #endif
