fealpy icon indicating copy to clipboard operation
fealpy copied to clipboard

Issue in ParametricLagrangeFiniteElementSpace

Open xieweidc opened this issue 3 years ago • 3 comments

在调用此空间中进行变系数问题的质量矩阵组装时,会报错 ps 未定义,报错原因是因为重心坐标下的数组没有定义,错误地址如下:

https://github.com/weihuayi/fealpy/blob/bc39a9248b62dde7959f3ef8d47519bd17ca0dc7/fealpy/functionspace/ParametricLagrangeFiniteElementSpace.py#L296-L300

ps = self.mesh.bc_to_point(bcs)

如果是添加完此行以后,又会发现c数组的einsum维数错误,一个是四个数组,却要对五个数组进行运算

https://github.com/weihuayi/fealpy/blob/bc39a9248b62dde7959f3ef8d47519bd17ca0dc7/fealpy/functionspace/ParametricLagrangeFiniteElementSpace.py#L304-L311

这一个问题是由于第308行, 暂时还没有想清楚比较好的解决方法,是不是应该仿照刚度矩阵的组装,对c进行一个分类处理? https://github.com/weihuayi/fealpy/blob/bc39a9248b62dde7959f3ef8d47519bd17ca0dc7/fealpy/functionspace/ParametricLagrangeFiniteElementSpace.py#L206-L237

xieweidc avatar Nov 22 '21 10:11 xieweidc

"如果是添加完此行以后,又会发现c数组的einsum维数错误,一个是四个数组,却要对五个数组进行运算" 这句话, 我不太理解. 提问最好的方式, 你应该提供一个最小可用的测试例子, 可以复现你的问题. 这样就可以更快定位问题.

weihuayi avatar Nov 25 '21 00:11 weihuayi

测试代码如下:

import numpy as np

from fealpy.mesh import MeshFactory as MF
from fealpy.functionspace import ParametricLagrangeFiniteElementSpace
from fealpy.decorator import cartesian

pi = np.pi
sin, cos = np.sin, np.cos

@cartesian
def c(p):
    x = p[..., 0]
    y = p[..., 1]
    
    val = (2+1.5*sin(2*pi*x))*(2+1.5*sin(2*pi))
    val = 1/val
    
    return val

box = np.array([0, 1, 0, 1])
mesh = MF.boxmesh2d(box, nx=3, ny=3, meshtype='quad', p=1)
space = ParametricLagrangeFiniteElementSpace(mesh, p=1)

M = space.mass_matrix(c=c, q=20)

报错如下:

Traceback (most recent call last):

  File "C:\Users\IDC\.spyder-py3\temp.py", line 24, in <module>
    M = space.mass_matrix(c=c, q=20)

  File "c:\users\idc\desktop\fealpy\fealpy\functionspace\ParametricLagrangeFiniteElementSpace.py", line 300, in mass_matrix
    c = c(ps)

NameError: name 'ps' is not defined

错误原因是 ps 数组未定义,地址如下: https://github.com/weihuayi/fealpy/blob/bc39a9248b62dde7959f3ef8d47519bd17ca0dc7/fealpy/functionspace/ParametricLagrangeFiniteElementSpace.py#L299-L300 添加

elif c.coordtype='cartesian':
     ps = self.mesh.bc_to_point(bcs)
     c = c(ps)

之后报错为

Traceback (most recent call last):

  File "C:\Users\IDC\.spyder-py3\temp.py", line 24, in <module>
    M = space.mass_matrix(c=c, q=20)

  File "c:\users\idc\desktop\fealpy\fealpy\functionspace\ParametricLagrangeFiniteElementSpace.py", line 309, in mass_matrix
    M = np.einsum('q, qci, qcj, qc, c->cij', ws*rm, phi, phi, d) # (NC, ldof, ldof)

  File "<__array_function__ internals>", line 5, in einsum

  File "C:\Users\IDC\anaconda3\lib\site-packages\numpy\core\einsumfunc.py", line 1361, in einsum
    return c_einsum(*operands, **kwargs)

ValueError: more operands provided to einstein sum function than specified in the subscripts string

错误原因是 https://github.com/weihuayi/fealpy/blob/bc39a9248b62dde7959f3ef8d47519bd17ca0dc7/fealpy/functionspace/ParametricLagrangeFiniteElementSpace.py#L308

xieweidc avatar Nov 25 '21 01:11 xieweidc

已经把 bug 改过来了, 我暂时推不到 github 上. gitlab 的上的仓库已经更新。

weihuayi avatar Dec 04 '21 01:12 weihuayi