Issue in BoundaryCondition

[xieweidc]

魏老师，我在看Dirichlet boundary condition时，有如下几个疑问。 (1) 这个阈值threshold是用来区分混合边界的一个参数吗？ 例如在： $$ u = g_1 on \Gamma_1 u = g_2 on \Gamma_2 $$ 此时的threshold 不为 None，这样理解对吗？ https://github.com/weihuayi/fealpy/blob/1bdc6a57e3fa884e8c406b67de67255d3d50a1ff/fealpy/boundarycondition/BoundaryCondition.py#L12 (2) https://github.com/weihuayi/fealpy/blob/1bdc6a57e3fa884e8c406b67de67255d3d50a1ff/fealpy/boundarycondition/BoundaryCondition.py#L18-L47 这里返回了修改的A, b，应该外部的刚度矩阵与载荷向量都已经改变，此处的Notes是否与 https://github.com/weihuayi/fealpy/blob/1bdc6a57e3fa884e8c406b67de67255d3d50a1ff/fealpy/boundarycondition/BoundaryCondition.py#L65-L87 重叠？ 参照poisson-fem的ppt，apply这个函数是否应该理解为将边界点对应的行列全赋值为1=0，再将对角线赋为1。但这样在一开始应该对载荷向量也进行处理（减去边界点的在内部自由度所在行的值），此处是怎么处理呢？另传的参数 https://github.com/weihuayi/fealpy/blob/1bdc6a57e3fa884e8c406b67de67255d3d50a1ff/fealpy/boundarycondition/BoundaryCondition.py#L18 uh=None 对照 https://github.com/weihuayi/fealpy/blob/1bdc6a57e3fa884e8c406b67de67255d3d50a1ff/example/PoissonFEMWithQuadMesh_example.py#L56-L71 在PoissonFEMWithQuadMesh里面，uh是起了一个什么作用呢？

[weihuayi]

1. threshold 是为了区分网格的边界边的属于哪一类的边界条件，它可以是一个函数，比如下面的函数就可以做为 theshold，你把所有边界边的中点 p 传给 这个函数，它返回一个逻辑数组，是 Dirichlet 类型的边，就标记为 True，否则就标记了 False

    @cartesian
    def is_dirichlet_boundary(self, p):
        y = p[..., 1]
        return ( y == 1.0) | ( y == 0.0)

也可以直接给一个边界边的编号数组或逻辑数组。

[weihuayi]

2. 进生边界条件处理，注释中说外部的矩阵 A 没有修改，意思是它的内存没有改变，F 修改意思是它存储的部分元素发生了改变。所有边界条件的处理，都会修改 F 的内存，所以外部的 F 对应发生变化，矩阵都是重新构造返回，不影响外部矩阵的内存。

[weihuayi]

uh 传进去，边界处理程序会边界条件加到 uh 上，即修改了 uh 中的元素（默认 uh 中的元素为 0），这个只是一种充分利用已经开辟内存的小技巧。

[weihuayi]

另外， 要记住 python 函数的参数都是引用传递，意味着你在函数内部修改参数指向的内存，外部的对象也就相应修改了。

[weihuayi]

Python 的基本原则就是能不拷贝就尽量不拷贝，这样主要是为了节省内存，这对于数组化编程来说非常重要。

[xieweidc]

老师，我在运行最新fealpy版本中的fealpy/example/PoissonFEMWithQuadMesh_example.py时候，报错如下：

idc@idc-GL553VD:~/fealpy/example$ python3 PoissonFEMWithQuadMesh_example.py 
usage: PoissonFEMWithQuadMesh_example.py [-h] [-p P] [-n N] [-b] [-o O]

该示例程序采用四边形网格上的双 p(>1) 次拉格朗日有限元， 求解区域 [0, 1]^2 上的纯 Dirichlet 边界的 Poisson 方程

optional arguments:
  -h, --help  show this help message and exit
  -p P        空间次数, 默认为 1, 即双线性元
  -n N        在 x 和 y 方向上，网格的剖分段数
  -b          固定求解区域 [0, 1]^2
  -o O        把结果输出 vtu 格式的文件，可用 Paraview 打开, 默认为 None 不输出
Namespace(b=[0, 1, 0, 1], n=10, o=None, p=1)
run boxmesh2d with time: 0.00036303299930295907
Traceback (most recent call last):
  File "PoissonFEMWithQuadMesh_example.py", line 67, in <module>
    A, F = bc.apply(A, F, uh)
  File "/home/idc/fealpy/fealpy/boundarycondition/BoundaryCondition.py", line 34, in apply
    isDDof = space.set_dirichlet_bc(gD, uh, threshold=threshold)
  File "/home/idc/fealpy/fealpy/functionspace/ParametricLagrangeFiniteElementSpace.py", line 423, in set_dirichlet_bc
    uh[isDDof] = gD(ipoints[isDDof])
  File "/home/idc/fealpy/fealpy/functionspace/Function.py", line 32, in __call__
    return space.value(self, bc, index=index)
  File "/home/idc/fealpy/fealpy/decorator/coordinates.py", line 20, in add_attribute
    return func(*args, **kwargs)
  File "/home/idc/fealpy/fealpy/functionspace/ParametricLagrangeFiniteElementSpace.py", line 117, in value
    val = np.einsum(s1, phi, uh[cell2dof])
  File "<__array_function__ internals>", line 5, in einsum
  File "/usr/local/lib/python3.8/dist-packages/numpy/core/einsumfunc.py", line 1359, in einsum
    return c_einsum(*operands, **kwargs)
ValueError: einstein sum subscripts string contains too many subscripts for operand 0

此报错应该为gD参数与uh参数位置写反了，我将这一行代码进行了一下修改，文件就能正常运行。 https://github.com/weihuayi/fealpy/blob/b7645f4e3973faabd648a0b616f6fd4b09103dbd/fealpy/boundarycondition/BoundaryCondition.py#L34

isDDof = space.set_dirichlet_bc(uh, gD, threshold=threshold)

另在计算FEM真解与数值解的误差时，此处的celltype参数是指代什么呢？ https://github.com/weihuayi/fealpy/blob/b7645f4e3973faabd648a0b616f6fd4b09103dbd/fealpy/quadrature/FEMeshIntegralAlg.py#L724

[weihuayi]

celltype = True， 意思计算逐个单元的误差

[weihuayi]

边界条件里面的参数顺序是我改的，其它相应地方我还没有来得及修改

[weihuayi]

我已经修改完了，下午方便的时候会推到 github 上

[xieweidc]

我已经修改完了，下午方便的时候会推到 github 上

好的，谢谢老师

[weihuayi]

我已经更新上来了， 所有空间的设置 dirichlet 边界条件的接口都是一致的， gD 做为第一个参数：

isDDof = space.set_dirichlet_bc(gD, uh, threshold=threshold)