firedrake的安装（docker环境）

firedrake安装

注意以下安装过程都是默认在docker环境下的服务器进行的，用户需要被服务器管理员添加锦docker组
1:百度搜索firedrake官网，进去，然后Ctrl + f，搜索docker，找到下面这个部分

2：点击进入上面提到的docker image，复制右边的命令：docker pull firedrakeproject/firedrake，把这段命令复制到服务器终端，然后服务器会开始下载跟firedrake相关的库，大概有两三个G.

3：下载完了以后，键入命令：
docker run -it --privileged=true firedrakeproject/firedrake /bin/bash
然后就差不多了，这个命令执行完差不多就在docker里面了，可以调用ls

说明安装好了。然后使用命令：exit可以退出docker。
4：键入命令：source firedrake/bin/activate。这样就可以使用firedrake了。
注意：使用exit退出来以后，之后每次进入docker就不是上述过程了，每次使用docker run -it – prieveleged=true firedrakeproject/firedrake /bin/bash，服务器都会开启一个新的docker
5：之后进入docker使用命令：
docker start ID
docker attach ID
这个ID是docker的ID，需要读者自己填入，可以使用命令：
docker ps -a

查看，比如博主的ID就是2d76a10c07f8
测试docker中的firedrake是否安装成功，进入docker以后，使用下面两个命令：
cd $VIRTUAL_ENV/src/firedrake
pytest tests/regression/ -k “poisson_strong or stokes_mini or dg_advection”

fireshape安装

谷歌搜索fireshape，会有一些经验，然而我们用的docker很可能用不上，下面介绍我的经验,下面假设已经进入docker
第一步：安装firedrake（参考上面）。
第二步：source firedrake/bin/activate
第三步：安装patchelf
pip install patchelf
或者sudo apt install patchelf（这一步可能会报错，可能需要先更新apt，利用sudo apt-get update）
第三步：安装ROL
pip3 install --no-cache-dir roltrilinos
pip3 install --no-cache-dir ROL
第四步：下载firedshape-master到docker，这里使用git clone可能无法克隆，所以需要先到github下载好fireshape压缩包，然后上传到docker里面，github链接
第五步：解压fireshape，解压以后一般会得到fireshape-master文件夹，进入里面然后修改里面的activate.sh第一行，把第一行修改为. .~/firedrake/bin/activate
第六步：在fireshape-master文件夹make
注意，这样make会报错，虽然在docker里面可能无法加载fireshape,不过在fireshape-master文件夹可以import fireshape，之后需要把相关的代码放到这里面运行

把文件传入传出docker

我们的代码需要在docker里面运行，为此我们需要先把文件传入docker，命令是：
cp 本地文件 ID:目的路径
比如说要把本地中/home/te 这个te文件夹复制到docker中的home里面，docker的ID 是 xx2d ,命令是
docker cp /home/te xx2d:/home
然后我们
docker start xx2d
docker attach xx2d
cd /home
就可以找到te文件夹了。
本人发现这样复制的到docker里面的文件夹，用户丧失了权限，无法使用rm -r te删除，permission denied，用sudo rm -r te可以删除
反过来，如果需要把docker里面的文件传到服务器，使用命令
docker cp xx2d:/home/te /home
这样就可以把docker里面/home目录下的te文件夹传出至本地/home目录
上面这几个cp命令都是在本地运行，在docker里面无法运行

代码运行

首先利用命令
docker start xx2d
docker attach xx2d
进入docker，然后用命令：
source firedrake/bin/activate
激活firedrake（有点像激活环境）
运行代码使用python xxx.py

neumann泊松算例

参考链接添加链接描述
然而docker里面似乎无法运行matplotlib库，无法直接画图，不过可以使用u.dat.data把网格点数值拿出来，而mesh.coordinates.dat.data表示网格点坐标，是一个$((n+1)\times (n+1),2)$的向量。

from firedrake import *
import numpy as np
import matplotlib.pyplot as plt
n = 10
mesh = UnitSquareMesh(n,n)
V = FunctionSpace(mesh,'CG',1)
u = Function(V,name = 'pred')
v = TestFunction(V)
x = SpatialCoordinate(mesh)
f = Function(V)
f.interpolate((1+8*pi*pi)*cos(x[0]*pi*2)*cos(x[1]*pi*2))

F = (inner(grad(u),grad(v)) + inner(u,v))*dx - inner(f,v)*dx
solve(F == 0,u,solver_parameters = {'ksp_type':'cg','pc_type':'none'})

u_acc = Function(V).interpolate(cos(2*pi*x[0])*cos(2*pi*x[1]))
L1_err = assemble(abs(u - u_acc)*dx)
err_norm = errornorm(u_acc,u)
print('firedrake:grid = [%d,%d],L1 error:%.3e,L2 error:%.3e'%(n,n,L1_err,err_norm))
File('helu.pvd').write(u)
#File('helacc.pvd').write(u_acc)
u_pred = u.dat.data.reshape(n + 1,n + 1).T
u_acc = u_acc.dat.data.reshape(n + 1,n + 1).T
err = u_acc - u_pred
x_train = np.linspace(0,1,n + 1,endpoint = True)
y_train = np.linspace(0,1,n + 1,endpoint = True)
X,Y = np.meshgrid(x_train,y_train)
plt.contourf(X,Y,err,40,cmap = 'rainbow')
plt.colorbar()
plt.xlabel('x')
plt.ylabel('y')
plt.title('firedrake:grid = [%d,%d],L1 error:%.3e,L2 error:%.3e'%(n,n,L1_err,err_norm))
plt.show()

Dirichlet算例

from firedrake import *

n = 10
mesh = CircleManifoldMesh(n,radius = 1)#表示圆形区域
#mesh = UnitSquareMesh(n,n)#方形区域
V = FunctionSpace(mesh,'CG',1)
u = Function(V,name = 'pred')
v = TestFunction(V)
x = SpatialCoordinate(mesh)
f = Function(V)
f.interpolate((1+8*pi*pi)*cos(x[0]*pi*2)*cos(x[1]*pi*2))
u_acc = Function(V).interpolate(cos(2*pi*x[0])*cos(2*pi*x[1]))
F = (inner(grad(u),grad(v)) + inner(u,v))*dx - inner(f,v)*dx
sub_domain = 'on_boundary'
bc = DirichletBC(V,u_acc,sub_domain)
solve(F == 0,u,solver_parameters = {'ksp_type':'cg','pc_type':'none'},bcs = [bc])
L1_err = assemble(abs(u - u_acc)*dx)
err_norm = errornorm(u_acc,u)
print('firedrake:grid = [%d,%d],L1 error:%.3e,L2 error:%.3e'%(n,n,L1_err,err_norm))

注意，上面生成的pvd文件用来可视化的，本人技术有限，暂时不知道如何把解转换为numpy进行matplotlib.pyplot画图，为此安装了一个paraview可以直接读取pvd文件，使用时需要来回传输文件在本地和服务器，docker转换
firedrake的使用一般需要和Gmesh,paraview这两个软件结合，为此读者需要事先把这两个软件也安装好

Gmesh基本入门

从官网安装好Gmesh以后，我们需要掌握Gmesh生成mesh的手法，Gmesh下载完一般会有一个tutorials文件夹，里面有一些demo，比如里面的t1.geo文件，打开Gmesh，加载t1.geo文件，选择mesh里面的2D就可以直接针对这个t1.geo设定的图形做剖分了，下面我将针对一个例子来说明Gmesh的具体用法。

第一步：点击File，选择New，自己新建一个名为house的geo文件，文件类型一直默认geo就行，然后发现在Gmesh软件上无法直接编辑house.geo，因此我们直接进入文件夹找到house.geo，直接打开（类似于记事本），这个时候会发现文件已经有了前面两行，但是似乎这两行有没有都一样，不影响

// Gmsh project created on Tue Apr 19 20:34:47 2022
SetFactory("OpenCASCADE");

第二步开始定义里面的顶点Point和剖分的边长lr

lr = 1e-1;
Point(1) = {0,0,0,lr};

上面两行表示选择原始剖分边长为lr=0.1，第一个点定义为(0,0,0)这是一个三维空间的点，如果只想定义二维空间的点，只需要把第三个坐标都设置为0就行。
第三步，定义线Line

Line(1) = {1,2};

上面这行代码表示第一条线是由第1，2两个点连接起来的。
第四步：定义面。首先需要把构成这个面的几条线都围起来

Curve Loop(1) = {1,2,3,4};
Plane Surface(1) = {1};

上面这行代码表示第一个面是由第1，2，3，4这四条线首尾连接的，注意这个有一点像向量，需要注意方向
如果是

Curve Loop(1) = {1,2,-3,4};
Plane Surface(1) = {1};

这种情况，表示第一个面遇到第三条线的时候会反方向连接。
这样就结束了，下面参考一个具体的house.geo文件。该文件写了9个点，16条边，9个面。

house.geo

// Gmsh project created on Tue Apr 19 20:34:47 2022
SetFactory("OpenCASCADE");

lr = 1e-1;
Point(1) = {0,0,0,lr};
Point(2) = {1,0,0,lr};
Point(3) = {1,1,0,lr};
Point(4) = {0,1,0,lr};

Point(5) = {0,0,1,lr};
Point(6) = {1,0,1,lr};
Point(7) = {1,1,1,lr};
Point(8) = {0,1,1,lr};

Point(9) = {0.5,0.5,1.5,lr};

//16 lines
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};

Line(5) = {1,5};
Line(6) = {2,6};
Line(7) = {3,7};
Line(8) = {4,8};

Line(9) = {5,6};
Line(10) = {6,7};
Line(11) = {7,8};
Line(12) = {8,5};

Line(13) = {5,9};
Line(14) = {6,9};
Line(15) = {7,9};
Line(16) = {8,9};

//surface
Curve Loop(1) = {1,2,3,4};
Plane Surface(1) = {1};

Curve Loop(2) = {1,6,-9,-5};
Plane Surface(2) = {2};

Curve Loop(3) = {2,7,-10,-6};
Plane Surface(3) = {3};

Curve Loop(4) = {3,8,-11,-7};
Plane Surface(4) = {4};

Curve Loop(5) = {4,5,-12,-8};
Plane Surface(5) = {5};

Curve Loop(6) = {9,14,-13};
Plane Surface(6) = {6};

Curve Loop(7) = {10,15,-14};
Plane Surface(7) = {7};

Curve Loop(8) = {11,16,-15};
Plane Surface(8) = {8};

Curve Loop(9) = {12,13,-16};
Plane Surface(9) = {9};

有了house.geo以后，我们需要利用Gmesh生成house.msh文件，表示对应的网格剖分。

选择Mesh里面的1D表示做一维剖分，2D表示二维剖分，3D表示三维剖分，然后选择save，保存好house.msh文件，然后把该文件传入docker里面就可以直接使用了。
上面这个house.geo用Gmesh得到的图像是：
使用鼠标拖拽可以看到三维效果，这个本人才发现。

mesh=Mesh(‘house.msh’)这行命令就是加载house.msh文件，不过本人发现这样得到的mesh似乎没有第三个维度，这个还有待考察，下面放一个加载t1.msh的代码

t1hel.py

from firedrake import *

mesh = Mesh('t1.msh')
V = FunctionSpace(mesh,'CG',1)
u = Function(V,name = 'pred')
v = TestFunction(V)
x = SpatialCoordinate(mesh)
f = Function(V)
f.interpolate((1+8*pi*pi)*cos(x[0]*pi*2)*cos(x[1]*pi*2))
u_acc = Function(V).interpolate(cos(2*pi*x[0])*cos(2*pi*x[1]))
F = (inner(grad(u),grad(v)) + inner(u,v))*dx - inner(f,v)*dx
sub_domain = 'on_boundary'
bc = DirichletBC(V,u_acc,sub_domain)
solve(F == 0,u,solver_parameters = {'ksp_type':'cg','pc_type':'none'},bcs = [bc])
L1_err = assemble(abs(u - u_acc)*dx)
err_norm = errornorm(u_acc,u)
print('firedrake:t1,L1 error:%.3e,L2 error:%.3e'%(L1_err,err_norm))