n次拉格朗日元

[guihunkun]

想请问一下老师，老师开发的支持n次拉格朗日元，一般像mfem也是，具体在程序包里，都是内置写了自由度，以及基函数，是用switch控制选择，还是，自由度节点，以及基函数也是程序自己去求的？

[weihuayi]

我们的实现是与维度和次数无关的实现方式， 没有switch 或者 if-else。 自由度节点和基函数都是程序根据网格和次数，自动计算出的。

[weihuayi]

@barycentric
    def basis(self, bc, index=np.s_[:], p=None):
        """
        compute the basis function values at barycentric point bc

        Parameters
        ----------
        bc : numpy.ndarray
            the shape of `bc` can be `(TD+1,)` or `(NQ, TD+1)`
        Returns
        -------
        phi : numpy.ndarray
            the shape of 'phi' can be `(1, ldof)` or `(NQ, 1, ldof)`

        See Also
        --------

        Notes
        -----

        """
        if p is None:
            p = self.p

        if p == 0 and self.spacetype == 'D':
            if len(bc.shape) == 1:
                return np.ones(1, dtype=self.ftype)
            else:
                return np.ones((bc.shape[0], 1), dtype=self.ftype)

        TD = bc.shape[-1] - 1 
        multiIndex = self.multi_index_matrix[TD](p)

        c = np.arange(1, p+1, dtype=np.int)
        P = 1.0/np.multiply.accumulate(c)
        t = np.arange(0, p)
        shape = bc.shape[:-1]+(p+1, TD+1)
        A = np.ones(shape, dtype=self.ftype)
        A[..., 1:, :] = p*bc[..., np.newaxis, :] - t.reshape(-1, 1)
        np.cumprod(A, axis=-2, out=A)
        A[..., 1:, :] *= P.reshape(-1, 1)
        idx = np.arange(TD+1)
        phi = np.prod(A[..., multiIndex, idx], axis=-1)
        return phi[..., np.newaxis, :] # (..., 1, ldof)

[weihuayi]

上面的代码对一维，二维和三维的网格，还有任意的 p 都是可以 work 的

[guihunkun]

谢谢老师，我学习一下。

[guihunkun]

老师方便说一下自由度节点根据网格和次数，自动计算的代码在那个文件里吗？

[weihuayi]

https://github.com/weihuayi/fealpy/blob/master/fealpy/functionspace/LagrangeFiniteElementSpace.py

[scaomath]

魏老师能否讲讲这里面的einsum是啥意思，每次看到einsum就头大： https://github.com/weihuayi/fealpy/blob/c9ca79a2976904fa5046f988ab9b79cf3381f913/fealpy/functionspace/LagrangeFiniteElementSpace.py#L219-L227

[weihuayi]

第一个 einsum 想把 measure 变成一个 （gdof, ldof) 的数组， 就是按列重复 ldof （这是插值点的个数）次，其实不需要这样做的。 guh 存储的是每个单元上每个插值点处梯度值，第二个 einsum 意思是，每个单元上每个插值点上的梯度值都对应乘以对应单元上的面积。

            measure = self.mesh.entity_measure('cell')
            NC = len(measure)
            ws = np.broadcast_to(measure[:, None], shape=(gdof, ldof))
            deg = np.zeros(gdof, dtype=self.ftype)
            np.add.at(deg, cell2dof, ws)
            guh = np.einsum('ij..., i->ij...', guh, measure) # 省略号表示对 1， 2,  3 维的都可以用
            if GD > 1:
                np.add.at(rguh, (cell2dof, np.s_[:]), guh)
            else:
                np.add.at(rguh, cell2dof, guh)

这段代码对高次元也可以用，但只有线性元的梯度恢复可以起到改善精度的效果。

[weihuayi]

第二个 einsum 其实也可以不用, 下面这个做法也可以

guh *= measure[:, None, None]

[weihuayi]

这个梯度恢复的函数从效率上讲有很多可以改进的地方，我抽时间改一个代码。einsum 函数功能很强大，灵活运用可以代替很多复杂的运算，但 numpy 的 einsum 并不支持多线程的运算，效率上有些问题。