AFEPack 使用 Tutorial（二）：解带系数二维泊松方程

最新推荐文章于 2024-04-24 09:41:51 发布

努力的老周

最新推荐文章于 2024-04-24 09:41:51 发布

阅读量315

点赞数

分类专栏： AFEPack 文章标签： PoissonEquation AFEPack PDE

本文链接：https://blog.csdn.net/justidle/article/details/116299257

版权

AFEPack 专栏收录该内容

3 篇文章 1 订阅

订阅专栏

AFEPack Step by Step Tutorial 2: solve Coefficient Poisson Equation

Coefficient Poisson Equation 的源码在 example/coefficient_possion_equation 目录下。

运行样例

编译 EasyMesh 网格文件

~/AFEPack/example/coefficient_possion_equation$ easymesh D
Working.  Please wait !
Improving the grid quality.  Please wait !
Renumerating nodes, elements and sides. Please wait !
Processing material marks. Please wait !

编译程序

使用 make 进行编译，具体编译信息可以自己查看 Makefile。

~/AFEPack/example/coefficient_possion_equation$ make
mpicxx -c -o coefficient_possion_equation.o coefficient_possion_equation.cpp -I. -O2 -fno-delete-null-pointer-checks -I/usr/include/trilinos -fPIC -D__SERIALIZATION__ -DMULTITHREAD -pthread -ftemplate-depth-256 -I/home/yizhou/AFEPack/library/include -I/home/yizhou/local/include -I/home/yizhou/local/include/deal.II -I/usr/include/trilinos
mpicxx -o main  coefficient_possion_equation.o -L/home/yizhou/local/lib -L../../library/lib -lAFEPack -ldeal_II  -ltbb -ldl -lm  -pthread

运行结果

~/AFEPack/example/coefficient_possion_equation$ ./run D
Reading easymesh data file ...
        reading node data ... OK!
        reading side data ... OK!
        reading element data ... OK!
AFEPack library file found: /home/yizhou/AFEPack/template/triangle/triangle.tmp_geo
AFEPack library file opened: /home/yizhou/AFEPack/template/triangle/triangle.tmp_geo
Opening shared library /home/yizhou/AFEPack/template/triangle/./triangle.geometry.so ...
        tried /home/yizhou/AFEPack/template/triangle/./triangle.geometry.so: success
AFEPack library file found: /home/yizhou/AFEPack/template/triangle/triangle.crd_trs
AFEPack library file opened: /home/yizhou/AFEPack/template/triangle/triangle.crd_trs
Opening shared library /home/yizhou/AFEPack/template/triangle/./triangle.geometry.so ...
        tried /home/yizhou/AFEPack/template/triangle/./triangle.geometry.so: success
AFEPack library file found: /home/yizhou/AFEPack/template/triangle/triangle.1.tmp_dof
AFEPack library file opened: /home/yizhou/AFEPack/template/triangle/triangle.1.tmp_dof
AFEPack library file found: /home/yizhou/AFEPack/template/triangle/triangle.1.bas_fun
AFEPack library file opened: /home/yizhou/AFEPack/template/triangle/triangle.1.bas_fun
Opening shared library /home/yizhou/AFEPack/template/triangle/./triangle.1.bas_fun.so ...
        tried /home/yizhou/AFEPack/template/triangle/./triangle.1.bas_fun.so: success
Opening shared library /home/yizhou/AFEPack/template/triangle/./triangle.1.bas_fun.so ...
        tried /home/yizhou/AFEPack/template/triangle/./triangle.1.bas_fun.so: success
Opening shared library /home/yizhou/AFEPack/template/triangle/./triangle.1.bas_fun.so ...
        tried /home/yizhou/AFEPack/template/triangle/./triangle.1.bas_fun.so: success
Building degree of freedom for the FEM space ...
        total 3009 degree of freedom found.
AMGSolver initializing ... OK! grid levels: 5
AMGSolver begin with initial residual 880.211 ...
        failed to converge with residual 0.633937 at step 21.

从上面可以看出运行结果出错了，failed to converge，第 21 步没法收敛。问题我看看是不是出在网格。

Coefficient Poisson Equation

这里，我们的二维带系数泊松方程为
$\left\{\begin{matrix} -K \cdot \Delta u = f\\ u |_{\partial \Omega} = u_b\\ u_b(x,y)=sin(x)*sin(y), \ if \ (x,y) \in \partial \Omega\\ f(x,y) = 1, \ for \ (x,y) \in \Omega \end{matrix}\right.$
其中 $u_b$ 是边界条件， $f$ 是右端项， $K$ 是一个 $\times 2$ 矩阵，定义如下：
$K=\begin{bmatrix} k_{11}(x,y) & k_{12}(x,y)\\ k_{21}(x,y) & k_{22}(x,y) \end{bmatrix}= \begin{bmatrix} \left\{\begin{matrix} 10.0+x \ \ \forall x<y\\ 2.0+y \ \ \forall x \geq y \end{matrix}\right. & 1.0 \\ 2.0 & \left\{\begin{matrix} 10.0+x \ \ \forall x+y>1\\ 1.0+y \ \ \forall x+y \leq 1 \end{matrix}\right. \end{bmatrix}$ 矩阵 $K$ 是 positive definite。

核心代码分析

Coefficient Poisson Equation 代码大部分和 2D Poisson Equation 代码是一样的。差别只是在构造基函数方面。2D Poisson Equation 在构造刚度矩阵（Stiff Matrix）的时候，直接使用 AFEPack 提供的模板类 StiffMatrix 即可。Coefficient Poisson Equation 需要自己构造基函数，所以我们需要从 StiffMatrix 类中派出一个新类。

BilinearOperator 类

BilinearOperator is the discretized version of the bilinear form on two linear space
$\qquad u \in U, v \in V]$ , thus for $U$ and $V$ to be finite dimensional space with basis function as $\phi^i, 1\leq i\leq M$ and $\psi^j, 1\leq j \leq N$ , respectively, the bilinear form can be given as a matrix with its entry as $f[a^{ij} = a(\phi^i, \psi^j)]$ .
The linear finite dimensional spaces here are two finite element space and the resulted matrix is a sparse matrix. The method getElementMatrix() is a pure virtual function that you must specify this method in derived class to contruct an object.

从上面的代码注释中可以看到，我们需要实现该类的纯虚函数 getElementMatrix()。该模板类在 AFEPack/library/include/BilinearOperator.h 文件中。定义如下：

template <int DIM, class value_type0, class value_type1=value_type0,
  int DOW=DIM, int TDIM0=DIM, int TDIM1 = DIM>
  class BilinearOperator : public SparseMatrix<double>

StiffMatrix 类

头文件位置

该模板类在 AFEPack/library/include/BilinearOperator.h 文件中。定义如下：

template <int DIM, class value_type, int DOW=DIM, int TDIM=DIM>
  class StiffMatrix : public BilinearOperator<DIM, value_type, value_type, DOW, TDIM, TDIM>
  {
  public:
  StiffMatrix() {};
  StiffMatrix(FEMSpace<typename StiffMatrix::value_type,DIM,DOW,TDIM>& sp) :
  BilinearOperator<DIM, typename StiffMatrix::value_type, 
  typename StiffMatrix::value_type,DOW,TDIM,TDIM>(sp, sp) {};
  virtual ~StiffMatrix() {};
  public:
  void reinit(FEMSpace<typename StiffMatrix::value_type,DIM,DOW,TDIM>& sp);
  virtual void getElementMatrix(const Element<typename StiffMatrix::value_type,DIM,DOW,TDIM>&, 
                                const Element<typename StiffMatrix::value_type,DIM,DOW,TDIM>&,
                                const typename ActiveElementPairIterator<DIM>::State state = ActiveElementPairIterator<DIM>::EQUAL);
  };

关键函数

getElementMatrix()

如上小节所示，我们需要实现 StiffMatrix 类中的虚函数 getElementMatrix()。这个虚函数将在父类 BilinearOperator 的 buildSparseMatrix() 函数中调用。函数 buildSparseMatrix() 将在父类 BilinearOperator 的 build() 中调用。
具体参看文件 AFEPack/library/include/BilinearOperator.templates.h。

类实现

我们从 StiffMatrix 类中派生，并实现纯虚函数 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);
};

getElementMatrix() 实现

根据带系数二维泊松方程定义，我们函数实现如下：

void Matrix::getElementMatrix(
		const Element<double,2>& element0,
		const Element<double,2>& element1,
		const ActiveElementPairIterator<2>::State state)
{
	int n_element_dof0 = elementDof0().size();
	int n_element_dof1 = elementDof1().size();
	double volume = element0.templateElement().volume();
	const QuadratureInfo<2>& quad_info = element0.findQuadratureInfo(algebricAccuracy());
	std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
	int n_quadrature_point = quad_info.n_quadraturePoint();
	std::vector<AFEPack::Point<2> > q_point = element0.local_to_global(quad_info.quadraturePoint());
	std::vector<std::vector<std::vector<double> > > basis_gradient = element0.basis_function_gradient(q_point);
	for (int l = 0;l < n_quadrature_point;l ++) {
		double Jxw = quad_info.weight(l)*jacobian[l]*volume;
		AFEPack::Point<2> q_point = element0.local_to_global(quad_info.quadraturePoint(l));
		double a11 = a11_(q_point);
		double a12 = a12_(q_point);
		double a21 = a21_(q_point);
		double a22 = a22_(q_point);
		for (int j = 0; j < n_element_dof0; j++) {
			for (int k = 0; k < n_element_dof1; k++) {
				elementMatrix(j,k) += Jxw*(
					a11 * basis_gradient[j][l][0] * basis_gradient[k][l][0]
					+ a12 * basis_gradient[j][l][0] * basis_gradient[k][l][1]
					+ a21 * basis_gradient[j][l][1] * basis_gradient[k][l][0]
					+ a22 * basis_gradient[j][l][1] * basis_gradient[k][l][1]
				);
			}
		}
	}
};

这样，我们就完成了核心代码实现。其他和二维泊松方程就一样了。

完整实现

Step 1: 网格文件

	EasyMesh mesh;
	mesh.readData("D");

Step 2: 构建模板单元

        // construct template geometry
        TemplateGeometry<2> triangle_template_geometry;
        triangle_template_geometry.readData("triangle.tmp_geo");

        // construct coordinate transformation
        CoordTransform<2,2> triangle_coord_transform;
        triangle_coord_transform.readData("triangle.crd_trs");

        // construct template degree of freedom
        TemplateDOF<2> triangle_template_dof(triangle_template_geometry);
        triangle_template_dof.readData("triangle.1.tmp_dof");

        // construct basis functions
        BasisFunctionAdmin<double,2,2> triangle_basis_function(triangle_template_dof);
        triangle_basis_function.readData("triangle.1.bas_fun");

        // construct template element, 
        // because to construct a finite element space, we general need a list of
        // template elements that we store the template elements in a vector and
        // provide it to a finite element space
        std::vector<TemplateElement<double,2,2> > template_element(1);
        template_element.reinit(triangle_template_geometry,
                triangle_template_dof,
                triangle_coord_transform,
                triangle_basis_function);

Step 3: 构建有限元空间

        // construct the finite element space on the mesh and the template element
        FEMSpace<double,2> fem_space(mesh, template_element);

        // this is the number of the geometry elements
        int n_element = mesh.n_geometry(2);
        // and set the number of elements in the finite element space the same as geometry elements
        fem_space.element().resize(n_element);
        // appoint that for every element, its template element is the 0-th one in the vector template_element
        for (int i = 0;i < n_element;i ++)
                fem_space.element(i).reinit(fem_space,i,0);

        // build the element information in the finite element space
        fem_space.buildElement();
        fem_space.buildDof();
        fem_space.buildDofBoundaryMark();

Step 4: 构建刚度矩阵

前面三步和二维泊松方程是一样的，这里和二维泊松方程不同。

	Matrix stiff_matrix(fem_space);
	stiff_matrix.algebricAccuracy() = 3;
	stiff_matrix.build();

Step 5: 解拉普拉斯方程

	FEMFunction<double,2> solution(fem_space);
	Vector<double> right_hand_side;
	Operator::L2Discretize(&f, fem_space, right_hand_side, 3);

	BoundaryFunction<double,2> boundary(BoundaryConditionInfo::DIRICHLET, 1, &u);
	BoundaryConditionAdmin<double,2> boundary_admin(fem_space);
	boundary_admin.add(boundary);
	boundary_admin.apply(stiff_matrix, solution, right_hand_side);

	AMGSolver solver(stiff_matrix);
	solver.solve(solution, right_hand_side);	

	solution.writeOpenDXData("u.dx");

定义 $u_b$ 函数

double u(const double * p)
{
	return sin(p[0]) * sin(p[1]);
};

定义 $f$ 函数

double f(const double * p)
{
	return 1.0;
};

这样，我们就完成了带系数二维泊松方程。

努力的老周

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
打赏
0
评论
AFEPack 使用 Tutorial（二）：解带系数二维泊松方程

AFEPack Step by Step Tutorial 2: solve Coefficient Poisson EquationCoefficient Poisson Equation 的源码在 example/coefficient_possion_equation 目录下。运行样例编译 EasyMesh 网格文件~/AFEPack/example/coefficient_possion_equation$ easymesh DWorking. Please wait !Improvi
复制链接

扫一扫