fealpy
fealpy copied to clipboard
weihuayi
the FEMeshIntegralAlg bug in ScaledMonomialSpace2d
魏老师, 可以先把 FEMeshIntegralAlg 从 ScaledMonomialSpace2d 中去掉吗?
对于 tri mesh, 做积分时 `FEMeshIntegralAlg' 仍然存在 bug, 魏老师可以尝试运行一下下面的测试程序
from fealpy.functionspace.ScaledMonomialSpace2d import ScaledMonomialSpace2d
import numpy as np
from fealpy.mesh import TriangleMesh
# --- mesh --- #
node = np.array([
(0, 0),
(1, 0),
(1, 1),
(0, 1)], dtype=np.float)
cell = np.array([
(1, 2, 0),
(3, 0, 2)], dtype=np.int)
mesh = TriangleMesh(node, cell)
p = 1
smspace = ScaledMonomialSpace2d(mesh, p)
# --- integral --- #
def f(x, index=None):
gphi = smspace.grad_basis(x, index=index)
gpphi = smspace.grad_basis(x, index=index, p=p+1)
return np.einsum('...mn, ...kn->...km', gphi, gpphi)
S = smspace.integralalg.integral(f, celltype=True)
问题还是出现在 function f 的 index 的取值上, 因为在 FEMeshIntegralAlg 中并不需要 index, 所以在 ScaledMonomialSpace2d --- line 274 and line 275
index = index if index is not None else np.s_[:]
phi[..., 1:3] = (point - self.cellbarycenter[index])/h[index].reshape(-1, 1)
这里会有 bug.
这里不是程序的 Bug, 是少输入了一个参数, FEMeshIntegralAlg 中默认处理的是 三角形, 四边形, 四面体, 六面体这样规则的单元, 它的输入点, 都是重心坐标的形式, 所以 barycenter 默认设置为 True, 你要处理输入是笛卡尔坐标点时, 就要设置为 False.
S = smspace.integralalg.integral(f, celltype=True, barycenter=False)
我正在尝试一种新坐标系 Decorator, 用户要积分的函数, 加上一个坐标系的 Decorator, 就会自动在你的函数对象里加一个属性, 这样数值积分程序, 可以自动判断你输入的函数需要什么样的坐标点了. 但还在尝试中, 但还有一些问题没有完全解决掉.
def cartesian(func):
@wraps(func)
def add_attribute(*args, **kwargs):
return func(*args, **kwargs)
add_attribute.__dict__['coordtype'] = 'cartesian'
return add_attribute
def barycentric(func):
@wraps(func)
def add_attribute(*args, **kwargs):
return func(*args, **kwargs)
add_attribute.__dict__['coordtype'] = 'barycentric'
return add_attribute
@cartesian
def fp(ps):
print('This function need a cartesian type input!')
@barycentric
def fb(bs):
print('This function need a barycentric type input!')
你可以查看在 StackOverflow 上的提问和回答 https://stackoverflow.com/questions/62550259/how-to-distinguish-two-kinds-of-similar-functions-or-methods-in-python
ScaledMonomialSpace2d 中引入 FEMeshIntegralAlg 是为混合元, 棱元等标准的向量有限元空间服务的, 是不能去掉的.
这里不是程序的 Bug, 是少输入了一个参数, FEMeshIntegralAlg 中默认处理的是 三角形, 四边形, 四面体, 六面体这样规则的单元, 它的输入点, 都是重心坐标的形式, 所以 barycenter 默认设置为 True, 你要处理输入是笛卡尔坐标点时, 就要设置为 False.
这里确实不是 FEMeshIntegralAlg 本身的 bug, 而是放在 ScaledMonomialSpace2d 中会对统一的程序有影响. 我猜想 ScaledMonomialSpace2d 的设计之初是为任意多边形分片间断 (DG) 空间考虑的 (VEM 空间也是要用到的), 所以 FEMeshIntegralAlg 在函数的参数的赋值上是和 PolygonMeshIntegralAlg 有部分冲突的, 这使得适用于多边形网格的程序不兼容三角形.
ScaledMonomialSpace2d 中引入 FEMeshIntegralAlg 是为混合元, 棱元等标准的向量有限元空间服务的, 是不能去掉的.
目前只是简单在 ScaledMonomialSpace2d 中通过 meshtype 来判断调用 PolygonMeshIntegralAlg 还是 FEMeshIntegralAlg, 这并不合理. 如果仅仅是为了标准的向量空间, 重新 import FEMeshIntegralAlg 不行吗. 或者应该提供另外的一个接口, 起码可让用户自己在程序中选择使用哪种积分.
我正在尝试一种新坐标系 Decorator, 用户要积分的函数, 加上一个坐标系的 Decorator, 就会自动在你的函数对象里加一个属性, 这样数值积分程序, 可以自动判断你输入的函数需要什么样的坐标点了. 但还在尝试中, 但还有一些问题没有完全解决掉.
期待.
从本质上讲, 一类数值离散方法有默认的网格类型要求. 而你写的是 HG 的程序, 输入网格类型应是多边形网格. 你的 HG 程序不应去适应所有的网格. 我前面实现的 VEM 的程序, 都是默认网格都是 PolygonMesh, 其它类型的网格都要做转化 PolygonMesh 才能使用. 这是符合数值离散方法的预期的.
而 ScaledMonomialSpace2d 空间开始确实是专门为多边形网格设计的, 但后来我意识到这个空间其实和拉格朗日空间一样基础(甚至更基础和灵活), 它们应得到同样的重视, 它们只是多项式空间的两种不同的表达形式.
嗯嗯 好的~.
不过我还是觉得在 ScaledMonomialSpace2d 中用 meshtype 来区分不同的空间不太合理 😂, 并且 fealpy 下 tri mesh 是可以直接用 PolygonMeshIntegralAlg 做积分的. 也就是程序本来不用考虑 meshtype, 现在需要区别一下, 影响也不大.