fealpy
fealpy copied to clipboard
weihuayi
Bug in QuadrangleMesh
我在6月22日课上讲的stiffmatrix.py将三角形网格换成矩形网格。 代码如下：
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import csr_matrix
from fealpy.mesh import MeshFactory as MF
from fealpy.functionspace import LagrangeFiniteElementSpace
p = 1 # 一次元
q = 3 # Gauss Legendre积分阶数
box = [0, 1, 0, 1]
mesh = MF.boxmesh2d(box, nx=1, ny=1, meshtype='quad')
# mesh = MF.boxmesh2d(box, nx=1, ny=1, meshtype='tri')
NN = mesh.number_of_nodes()
space = LagrangeFiniteElementSpace(mesh, p=p)
gdof = space.number_of_global_dofs()
qf = mesh.integrator(q, 'cell')
bcs, ws = qf.get_quadrature_points_and_weights()
cellmeasure = mesh.entity_measure('cell') #
gphi = space.grad_basis(bcs) # (NQ, NC, ldof, 2)
# (NC, ldof, ldof)
A = np.einsum('q, qcid, qcjd, c->cij', ws, gphi, gphi, cellmeasure)
cell2dof = space.cell_to_dof() # (NC, ldof)
print('cell2dof:\n', cell2dof)
# (NC, ldof) --> (NC, ldof, 1) --> (NC, ldof, ldof)
I = np.broadcast_to(cell2dof[:, :, None], shape=A.shape)
print('I:\n', I)
# (NC, ldof) --> (NC, 1, ldof) --> (NC, ldof, ldof)
J = np.broadcast_to(cell2dof[:, None, :], shape=A.shape)
print('J:\n', J)
A = csr_matrix((A.flat, (I.flat, J.flat)), shape=(gdof, gdof))
phi = space.basis(bcs) # (NQ, NC, ldof)
# (NC, ldof, ldof)
M = np.einsum('q, qci, qcj, c->cij', ws, phi, phi, cellmeasure)
M = csr_matrix((M.flat, (I.flat, J.flat)), shape=(gdof, gdof))
def f(p):
pi = np.pi
x = p[..., 0]
y = p[..., 1]
return np.sin(pi*p[..., 0])*np.sin(pi*p[..., 1])
ps = mesh.bc_to_point(bcs) # (NQ, NC, GD)
val = f(ps) # (NQ, NC)
# (NC, ldof)
bb = np.einsum('q, qc, qci, c->ci', ws, val, phi, cellmeasure)
print('bb:\n', bb)
b = np.zeros(gdof, dtype=np.float64)
np.add.at(b, cell2dof, bb)
print('b:\n', b)
ipoints = space.interpolation_points()
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)
mesh.find_node(axes, node=ipoints, showindex=True, color='r', fontsize=24)
plt.show()
print('bcs:\n', bcs)
print('ws:\n', ws)
报错如下：
run boxmesh2d with time: 0.0027847709999946346
Traceback (most recent call last):
File "/home/idc/桌面/python/fealpy_why/tutorial/stiff_matrix.py", line 16, in <module>
space = LagrangeFiniteElementSpace(mesh, p=p)
File "/home/idc/anaconda3/lib/python3.8/site-packages/fealpy/functionspace/LagrangeFiniteElementSpace.py", line 68, in __init__
self.integralalg = FEMeshIntegralAlg(
File "/home/idc/anaconda3/lib/python3.8/site-packages/fealpy/quadrature/FEMeshIntegralAlg.py", line 20, in __init__
self.edgeintegrator = mesh.integrator(q, 'edge')
File "/home/idc/anaconda3/lib/python3.8/site-packages/fealpy/mesh/QuadrangleMesh.py", line 58, in integrator
return GaussLegendreQuadrature(k)
NameError: name 'GaussLegendreQuadrature' is not defined
我将mesh/QuadrangleMesh添加了
from ..GaussLegendreQuadrature import GaussLegendreQuadrature
GaussLegendreQuadrature(k)没有定义的错误被修正，但又出现了：
run boxmesh2d with time: 0.0027926709999519517
Traceback (most recent call last):
File "/home/idc/桌面/python/fealpy_why/tutorial/stiff_matrix.py", line 18, in <module>
gdof = space.number_of_global_dofs()
File "/home/idc/anaconda3/lib/python3.8/site-packages/fealpy/functionspace/LagrangeFiniteElementSpace.py", line 80, in number_of_global_dofs
return self.dof.number_of_global_dofs()
AttributeError: 'LagrangeFiniteElementSpace' object has no attribute 'dof'
报错表明拉格朗日有限元空间没有‘dof’这个属性，但 https://github.com/weihuayi/fealpy/blob/42425503eff3e19f7965ef28af42a83169e9c84b/fealpy/functionspace/LagrangeFiniteElementSpace.py#L22-L56
但functionspace/LagrangeFiniteElementSpace.py的这一段不是已经把'dof'属性给LagrangeFiniteElementSpace这个类了吗？ 另从我的运行信息中可以观察到fealpy是安装在site-packages下面的，希望老师能解答一下fealpy在Windows和Ubuntu系统下软硬安装的办法。
首先， LagrangeFiniteElementSpace 只适用于区间，三角形，四面体网格，不适用于四边形网格. 四边形网格的例子，见 https://github.com/weihuayi/fealpy/blob/master/example/PoissonFEMWithQuadMesh_example.py
老师我在学习四边形网格的例子时，它采用了ParametricLagrangeFiniteElementSpace这个函数空间， https://github.com/weihuayi/fealpy/blob/42425503eff3e19f7965ef28af42a83169e9c84b/fealpy/functionspace/ParametricLagrangeFiniteElementSpace.py#L12-L13
https://github.com/weihuayi/fealpy/blob/42425503eff3e19f7965ef28af42a83169e9c84b/fealpy/functionspace/ParametricLagrangeFiniteElementSpace.py#L130
这两行代码中的参数c各自代表的意思是什么呢？
老师我在学习四边形网格的例子时，它采用了ParametricLagrangeFiniteElementSpace这个函数空间， https://github.com/weihuayi/fealpy/blob/42425503eff3e19f7965ef28af42a83169e9c84b/fealpy/functionspace/ParametricLagrangeFiniteElementSpace.py#L12-L13
https://github.com/weihuayi/fealpy/blob/42425503eff3e19f7965ef28af42a83169e9c84b/fealpy/functionspace/ParametricLagrangeFiniteElementSpace.py#L130
这两行代码中的参数c各自代表的意思是什么呢？
第一个spacetype='C'的C代表continuous。第二个c是装矩阵时候的系数，None就是1，也可以传一个函数或者矩阵进去。
https://github.com/weihuayi/fealpy/blob/029882c7e366d4623f1c79a4323edad191499623/fealpy/functionspace/ParametricLagrangeFiniteElementSpace.py#L174-L181 老师，在调用ParametricLagrangeFiniteElementSpace时，这三个变量可以理解为计算每单元的面积吗？这里的第一基本形式是什么意思呢？
老师，在调用
ParametricLagrangeFiniteElementSpace时，这三个变量可以理解为计算每单元的面积吗？这里的第一基本形式是什么意思呢？
不是单元的，因为这个地方要用quadrature在参数曲面上F(u,v)算积分，这个是算在积分点处的参数曲面的Jacobi矩阵相乘F_u叉乘F_v。
最简单的例子就是https://en.wikipedia.org/wiki/Surface_integral 里面第一个公式里面那个叉乘。
这个空间开始设置了mesh之后(simplicial或者tensorial)，就会直接调用mesh的first_fundamental_form函数，比如这个：
https://github.com/weihuayi/fealpy/blob/029882c7e366d4623f1c79a4323edad191499623/fealpy/mesh/LagrangeTriangleMesh.py#L336-L365
老师，在调用
ParametricLagrangeFiniteElementSpace时，这三个变量可以理解为计算每单元的面积吗？这里的第一基本形式是什么意思呢？不是单元的，因为这个地方要用quadrature在参数曲面上
F(u,v)算积分，这个是算在积分点处的参数曲面的Jacobi矩阵相乘F_u叉乘F_v。最简单的例子就是https://en.wikipedia.org/wiki/Surface_integral 里面第一个公式里面那个叉乘。
这个空间开始设置了
mesh之后(simplicial或者tensorial)，就会直接调用mesh的first_fundamental_form函数，比如这个：https://github.com/weihuayi/fealpy/blob/029882c7e366d4623f1c79a4323edad191499623/fealpy/mesh/LagrangeTriangleMesh.py#L336-L365
好的，谢谢