AFEPack
|
00001 00011 #ifndef __MovingMesh2D_h__ 00012 #define __MovingMesh2D_h__ 00013 00014 #include <lac/sparsity_pattern.h> 00015 #include <lac/sparse_matrix.h> 00016 #include "AMGSolver.h" 00017 #include "Geometry.h" 00018 #include "EasyMesh.h" 00019 00020 class MovingMesh2D : public EasyMesh 00021 { 00022 public: 00023 typedef GeometryBM::bmark_t bound_t; 00024 struct Vertex : public afepack::Point<2> { 00025 int index; 00026 int boundary_mark; 00027 }; 00028 struct Edge { 00029 int index; 00030 int vertex[2]; 00031 int boundary_mark; 00032 double normal[2]; 00033 double logical_normal[2]; 00034 }; 00035 struct Domain { 00036 int n_vertex; 00037 int n_edge; 00038 std::vector<Vertex> physical_domain_vertex; 00039 std::vector<Vertex> logical_domain_vertex; 00040 std::vector<Edge> edge; 00041 }; 00042 class Solver { 00043 private: 00044 bool is_initialized; 00045 int n_project; 00046 std::vector<SparseMatrix<double> *> project_matrix; 00047 std::vector<SparseMatrix<double> *> project_matrix_r; 00048 std::vector<SparseMatrix<double> *> projected_matrix; 00049 std::vector<const std::vector<int> *> boundary_mark; 00050 u_int smooth_step; 00051 u_int min_order; 00052 const Domain * domain; 00053 public: 00054 Solver(); 00055 ~Solver(); 00056 public: 00057 u_int smoothStep() const {return smooth_step;} 00058 u_int& smoothStep() {return smooth_step;} 00059 u_int minimalOrder() const {return min_order;} 00060 u_int& minimalOrder() {return min_order;} 00061 00062 void reinit(SparseMatrix<double>&, 00063 const std::vector<int>&, 00064 const Domain&); 00065 void clear(); 00066 void solve(std::vector<Vector<double> >& x, 00067 const std::vector<Vector<double> >& r, 00068 u_int steps = 5) const; 00069 private: 00070 void Project(const SparseMatrix<double>& M, 00071 const std::vector<int>& bm, 00072 SparseMatrix<double> *& P, 00073 SparseMatrix<double> *& PMPT, 00074 SparseMatrix<double> *& Pt, 00075 std::vector<bound_t> *& pbm); 00076 void Init(SparseMatrix<double>&, 00077 const std::vector<int>&, 00078 const Domain&); 00079 void GaussSidel(const SparseMatrix<double>& M, 00080 std::vector<Vector<double> >& x, 00081 const std::vector<Vector<double> >& r, 00082 const std::vector<int>& bm, 00083 const u_int& s) const; 00084 SparseMatrix<double> * getPMPT(const SparseMatrix<double> & P, 00085 const SparseMatrix<double> & M, 00086 const SparseMatrix<double> & Pt) const; 00087 void lazyPMPT(const SparseMatrix<double>& P, 00088 const SparseMatrix<double>& M, 00089 const SparseMatrix<double>& Pt, 00090 SparseMatrix<double>& A) const; 00091 }; 00092 00093 private: 00094 Domain domain; 00095 std::vector<afepack::Point<2> > logical_node; 00096 double move_step_length; 00097 int n_move_step; 00098 std::vector<afepack::Point<2> > move_direction; 00099 std::vector<afepack::Point<2> > logical_move_direction; 00100 std::vector<double> mon; 00101 00102 std::vector<std::vector<int> > mb_node; 00103 std::vector<int> boundary_mark; 00104 00105 SparsityPattern spM; 00106 SparseMatrix<double> M; 00107 Solver solver; 00108 00109 static int primes[]; 00110 u_int solve_step; 00111 00112 u_int max_step; 00114 double tol; 00116 double energy_reduction; 00117 00118 public: 00119 MovingMesh2D(); 00120 virtual ~MovingMesh2D(); 00121 00122 public: 00123 const std::vector<double>& monitor() const {return mon;}; 00124 std::vector<double>& monitor() {return mon;}; 00125 const double& monitor(const int& i) const {return mon[i];}; 00126 double& monitor(const int& i) {return mon[i];} 00127 00128 const std::vector<afepack::Point<2> >& moveDirection() const {return move_direction;}; 00129 std::vector<afepack::Point<2> >& moveDirection() {return move_direction;}; 00130 const afepack::Point<2>& moveDirection(const int& i) const {return move_direction[i];}; 00131 afepack::Point<2>& moveDirection(const int& i) {return move_direction[i];}; 00132 00133 std::vector<double> moveDirection(const afepack::Point<2>& point, 00134 const int& element) const; 00135 std::vector<std::vector<double> > moveDirection(const std::vector<afepack::Point<2> >& point, 00136 const int& element) const; 00137 double moveDirectionDivergence(const int& element) const; 00138 00139 double moveStepLength() const {return move_step_length;}; 00140 double& moveStepLength() {return move_step_length;}; 00141 00142 int moveStep() const {return n_move_step;}; 00143 int& moveStep() {return n_move_step;}; 00144 00145 u_int solveStep() const {return solve_step;} 00146 u_int& solveStep() {return solve_step;} 00147 00148 double tolerence() const {return tol;} 00149 double& tolerence() {return tol;} 00150 00151 u_int maxStep() const {return max_step;} 00152 u_int& maxStep() {return max_step;} 00153 00154 void moveMesh(); 00155 void outputPhysicalMesh(const std::string& file); 00156 void outputLogicalMesh(const std::string& file); 00157 virtual void getMonitor(); 00158 virtual void smoothMonitor(int step = 1); 00159 virtual void updateMesh(); 00160 virtual void updateSolution() = 0; 00161 virtual void outputSolution() = 0; 00162 virtual void getMoveStepLength(); 00163 void readDomain(const std::string& file); 00164 00165 private: 00166 void getLogicalMesh(); 00167 void getMoveDirection(); 00168 void readDummy(std::ifstream& is); 00169 void parseBoundary(); 00170 00171 void outputMoveDirection(const std::string&); 00172 }; 00173 00174 #endif