the `integral` function in `class FEMeshIntegralAlg()`

[yoczhang]

从 line-621 --- line-640

# old api 
def integral(self, u, celltype=False, barycenter=True):
    """
        """
    qf = self.integrator
    bcs = qf.quadpts  # 积分点 (NQ, 3)
    ws = qf.weights  # 积分点对应的权重 (NQ, )
    if barycenter or u.coordtype == 'barycentric':
        val = u(bcs)
    else:
        ps = self.mesh.bc_to_point(bcs)  # (NQ, NC, 2)
        val = u(ps)
    dim = len(ws.shape)
    s0 = 'abcde'
    s1 = '{}, {}j..., j->j...'.format(s0[0:dim], s0[0:dim])
    e = np.einsum(s1, ws, val, self.cellmeasure)
    if celltype is True:
        return e
    else:
        return e.sum()

• 1. 首先 barycenter=True 是否改为 barycenter=False 更合理点? 因为 if 的判断 if barycenter or u.coordtype == 'barycentric': 已经默认可以执行 barycentric 的情况. 并且改为 barycenter=False 时, 将有可能统一多边形和三角形上的体积分.
• 2. 在执行 val = u(ps) 前, 能否定义一个 index = np.arange(NC), 然后将 val = u(ps) 改为 val = u(ps, index), 目的也是统一多边形和三角形上的体积分. 但是不知道这样更改对其他程序有没有影响.

---

魏老师能否做下下面的测试, 会出现上面的问题

from fealpy.functionspace.ScaledMonomialSpace2d import ScaledMonomialSpace2d
import numpy as np
from fealpy.mesh import TriangleMesh
from fealpy.decorator import cartesian

# --- mesh
n = 0
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)
mesh.uniform_refine(n)

p = 1
smspace = ScaledMonomialSpace2d(mesh, p)

uh = smspace.function()

@cartesian
def f1(x, index=None):
    return uh.value(x, index)

S = smspace.integralalg.integral(f1, celltype=True, barycenter=False)

[weihuayi]

integral 这个 api 是没有引入 coordtype 以前的 API， 目前为兼容以前的例子，就留在这里。单元上的积分， 建议使用 cell_integral.

或者，integral 保留， barycenter 这个参数也去掉

    def integral(self, u, etype='cell', q=None):
        """

        Notes
        -----
            计算函数 u 在指定网格实体上的整体积分。
        """
        if etype == 'cell':
            e = np.sum(self.cell_integral(u, q=q))
        elif etype == 'face':
            e = np.sum(self.face_integral(u, q=q))
        elif etype == 'edge':
            e = np.sum(self.edge_integral(u, q=q))

        return e

你看怎么样？

[weihuayi]

或者命名为 mesh_integral(...)?

[weihuayi]

指明是在整个网格上的积分

[weihuayi]

原始的 integral, 仍然保留， 只是不再判断 coordtype, 这样还和一些老的程序兼容。 而 PolygonMesh 上的积分算法也会考虑和 FEM 网格积分算法的接口调整为一致的。

[yoczhang]

单元上的积分， 建议使用 cell_integral.

FEMeshIntegralAlg() 下的 cell_integral 在 ps = mesh.bc_to_point(bcs, 'cell') 调用 f = f(ps), 仍然需要 index=np.arange(NC) 参数, 即 f = f(ps, index), 原因在于 ScaledMonomialSpace2d() 的 value() 函数

    @cartesian
    def value(self, uh, point, index=np.s_[:]):
        phi = self.basis(point, index=index)
        cell2dof = self.dof.cell2dof
        dim = len(uh.shape) - 1
        s0 = 'abcdefg'
        s1 = '...ij, ij{}->...i{}'.format(s0[:dim], s0[:dim])
        return np.einsum(s1, phi, uh[cell2dof[index]])

当 index 没指定时, uh[cell2dof[index]] 会多出一个维度, 或者用 np.squeeze(uh[cell2dof[index]]) 消去冗余的维度??

这里在某个地方的逻辑上肯定还是有点问题的.

魏老师可以简单的调试一下我上面提到的那个测试.

---

原始的 integral, 仍然保留，

保留还是挺好的

而 PolygonMesh 上的积分算法也会考虑和 FEM 网格积分算法的接口调整为一致的

是的, 目前, 比如 PolygonMeshIntegralAlg() 和 FEMeshIntegralAlg() 中的 cell_integral 接口还是不一样.

[weihuayi]

好的， 我再测试一下。     ------------------ Original ------------------ From:  "YcZhang"<[email protected]>; Date:  Wed, Aug 5, 2020 09:38 AM To:  "weihuayi/fealpy"<[email protected]>; Cc:  "魏华祎"<[email protected]>; "Comment"<[email protected]>; Subject:  Re: [weihuayi/fealpy] the integral function in class FEMeshIntegralAlg() (#26)

 

单元上的积分， 建议使用 cell_integral.

FEMeshIntegralAlg() 下的 cell_integral 在 ps = mesh.bc_to_point(bcs, 'cell') 调用 f = f(ps), 仍然需要 index=np.arange(NC) 参数, 即 f = f(ps, index), 原因在于 ScaledMonomialSpace2d() 的 value() 函数 @cartesian def value(self, uh, point, index=np.s_[:]): phi = self.basis(point, index=index) cell2dof = self.dof.cell2dof dim = len(uh.shape) - 1 s0 = 'abcdefg' s1 = '...ij, ij{}->...i{}'.format(s0[:dim], s0[:dim]) return np.einsum(s1, phi, uh[cell2dof[index]])

当 index 没指定时, uh[cell2dof[index]] 会多出一个维度, 或者用 np.squeeze(uh[cell2dof[index]]) 消去冗余的维度??

这里在某个地方的逻辑上肯定还是有点问题的.

魏老师可以简单的调试一下我上面提到的那个测试.

原始的 integral, 仍然保留，

保留还是挺好的

而 PolygonMesh 上的积分算法也会考虑和 FEM 网格积分算法的接口调整为一致的

是的, 目前, 比如 PolygonMeshIntegralAlg() 和 FEMeshIntegralAlg() 中的 cell_integral 接口还是不一样.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[weihuayi]

你定义函数的时候， index 不能默认为 None， 应该是 np.s_[:]

@cartesian
def f1(x, index=np.s_[:]):
    return uh.value(x, index)

[yoczhang]

好的, 魏老师, 还有一个问题, ScaledMonomialSpace2d() 中, 是否应将网格 'quad' 类型也加入 if mtype in {'polygon', 'halfedge', 'halfedge2d'}: 中 ? 即 if mtype in {'quad', 'polygon', 'halfedge', 'halfedge2d'}: ?

[weihuayi]

quad 的积分和三角形的模式是一样的。     ------------------ Original ------------------ From:  "YcZhang"<[email protected]>; Date:  Thu, Aug 6, 2020 10:09 AM To:  "weihuayi/fealpy"<[email protected]>; Cc:  "魏华祎"<[email protected]>; "Comment"<[email protected]>; Subject:  Re: [weihuayi/fealpy] the integral function in class FEMeshIntegralAlg() (#26)

 

好的, 魏老师, 还有一个问题, ScaledMonomialSpace2d() 中, 是否应将网格 'quad' 类型也加入 if mtype in {'polygon', 'halfedge', 'halfedge2d'}: 中 ? 即 if mtype in {'quad', 'polygon', 'halfedge', 'halfedge2d'}: ?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[yoczhang]

嗯, 但是 ScaledMonomialSpace2d() 中是不是并没有对应 'quad' mesh 的 integralalg 呢?

[weihuayi]

是的， 我需要把逻辑修正一下，就可以了。