[笔记2.1]学习AFEPack的例程

[学习例程1]：/usr/local/AFEPack/example/possion_equation/
问题：-/Delta u = f u |_{/partial /Omega} = u_b

＊＊＊＊＊ $ls
CVS D.d Makefile possion_equation.cfg possion_equation.cpp possion_equation.net README run
------------------------------------------------------------------------------------
＊＊＊＊＊ README （运行步骤）


This example will solve a Possion equation on an easymesh generated data.


The obtained solution is saved in file "u.dx".




Usage:


1. compile with make [-e debug=no];


2. generate the mesh with easymesh as:


easymesh D


3. run with ./run;


4. view the obtained numerical solution with OpenDX;


[补充］
3. 本目录的例子只需运行 $ ./run 或者 $ ./main
（后者需要先设置环境变量 $ export AFEPACK_TEMPLATE_PATH=/usr/local/AFEPack/template/triangle，否则找不到读入文件的路径）
4. OpenDX ——> Edit Visual Programs...——> 选 possion_equation.net
(possion_equation.cfg 是OpenDx用到的配置文件)
双击FileSelector，（若需要）修改输出的数值解位置，保存
Execute ——> 选Execute Once

＊＊＊＊＊ Makefile（略）

＊＊＊＊＊ possion_equation.cpp（略）
可用 $ cat -n possion_equation.cpp 加上行号显示

＊＊＊＊＊ D.d（EasyMesh输入文件，略）

＊＊＊＊＊ run（shell脚本）


#! /bin/sh




export AFEPACK_PATH="/usr/local/AFEPack"


export AFEPACK_TEMPLATE_PATH="$AFEPACK_PATH/template/triangle"




exec ./main $*

［补充］（shell编程和环境变量，略）
$*：执行shell程序时，传递给程序的所有参数组成的一个字符串。
比如：
$ ./run D 实际是运行了 ./main D
$ ./run parameters are three
此时 $* ＝ parameters are three；$1 ＝ parameters；$2 = are; $3 = three; $# = 3.
------------------------------------------------------------------------------------
＊＊＊＊＊ 关于possion_equation.cpp的补充解释（详见[笔记2.2]）
(1) Vector类（实现数值的 STL：vector）


65 Vector<double> right_hand_side;

template<typename number> class Vector;
位置：/usr/local/deal.II/lac/include/lac/vector.h 和 vector.templates.h

(2) 求解方程


73 AMGSolver solver(stiff_matrix);


74 solver.solve(solution, right_hand_side);

[运行提示信息] 其中 failed to converge...，无大碍，唯不顺眼也
AMGSolver initializing ... OK! grid levels: 3
AMGSolver begin with initial residual 0.458239 ...
failed to converge with residual 9.01090e-11 at step 21.

L2 error = 0.000629244

[原因] AMGSolver::solve是默认参数的函数，可将20放大，比如：
74行改为：solver.solve(solution, right_hand_side, 0, 0, 50);




class AMGSolver ...{


...


public:


void solve(Vector<double>& x,


const Vector<double>& r,


double tol = 0.0,


unsigned int s = 0,


unsigned int s1 = 20) const;


...


};

参数:
x the solution vector
r the right hand side vector
tol tolerence
s 0: solve mode; 1: precondition mode
s1 upper bound of iteration step

AMGSolver::solve 函数定义在：/usr/local/AFEPack/library/src/AMGSolver.cpp

[修改后结果]
73 AMGSolver solver(stiff_matrix);
74 solver.solve(solution, right_hand_side,0,0,50);

AMGSolver initializing ... OK! grid levels: 3
AMGSolver begin with initial residual 0.458239 ...
converge with residual 4.32514e-13 at step 27.

L2 error = 0.000629244

(3) 计算和精确解的L2误差（本例Dirichlet边界条件由精确解给出）

78 double error = Functional::L2Error(solution, FunctionFunction<double>(&u), 3);

(4) 计算解和近似精确解的误差（后述）
[学习例程2]：/usr/local/AFEPack/example/coefficient_possion_equation/
问题：-/nabla/cdot(A /nabla u) = f, 其中：A=(a_{ij})_{2*2}
u |_{/partial /Omega} = u_b

＊＊＊＊＊ $ls
coefficient_possion_equation.cfg coefficient_possion_equation.net D.d README
coefficient_possion_equation.cpp CVS Makefile run

如上例：.cfg .net是OpenDx文件。本例只需关注 .cpp文件
------------------------------------------------------------------------------------
＊＊＊＊＊ 比较coefficient_possion_equation.cpp和possion_equation.cpp的不同之处：
(1) 增加2*2矩阵A的定义
(2) 两问题：possion_equation 和 coefficient_possion_equation 各自的单元刚度矩阵不同
故定义刚度矩阵的一个新派生类，虚函数 getElementMatrix(...) 重载基类的函数


class Matrix : public StiffMatrix<2, double>




...{


public:


Matrix(FEMSpace<double,2>& sp) :




StiffMatrix<2,double>(sp) ...{};




virtual ~Matrix() ...{};


public:


virtual void getElementMatrix(const Element<double,2>& element0,


const Element<double,2>& element1,


const ActiveElementPairIterator<2>::State state);


};

［补充］：
1) C++允许子类的成员函数重载基类的成员函数
2) 运行时，依据类型确认调用哪个函数的能力称为多态性(polymorphism)
3) 为指明某个成员函数具有多态性，用virtual来标志其为虚函数

(3) 函数 void Matrix::getElementMatrix 的定义和 BilinearOperator.templates.h 文件中 StiffMatrix<DIM,value_type,DOW,TDIM>::getElementMatrix 定义类似，只是加进了系数A，可对比看

(4) 函数getElementMatrix是提供给stiff_matrix.build();内部调用的，所以：
possion_equation.cpp


StiffMatrix<2,double> stiff_matrix(fem_space);


stiff_matrix.algebricAccuracy() = 3;


stiff_matrix.build();

改为 coefficient_possion_equation.cpp


Matrix stiff_matrix(fem_space);


stiff_matrix.algebricAccuracy() = 3;


stiff_matrix.build();