From 1692c13bf215e5c41b22fc69a629145d76893c5a Mon Sep 17 00:00:00 2001 From: "2655493031@qq.com" <2655493031@qq.com> Date: Thu, 19 Dec 2024 12:55:35 +0800 Subject: [PATCH] update --- example/fem/surface_poisson_lfem_example.py | 1 + fealpy/mesh/lagrange_quadrangle_mesh.py | 90 +++++++++++++++------ fealpy/mesh/lagrange_triangle_mesh.py | 20 ++--- test/mesh/test_lagrange_quadrangle_mesh.py | 20 +++-- test/mesh/test_lagrange_triangle_mesh.py | 6 +- 5 files changed, 96 insertions(+), 41 deletions(-) diff --git a/example/fem/surface_poisson_lfem_example.py b/example/fem/surface_poisson_lfem_example.py index c9694e516..e1dfabfd1 100644 --- a/example/fem/surface_poisson_lfem_example.py +++ b/example/fem/surface_poisson_lfem_example.py @@ -97,6 +97,7 @@ def coo(A): data = A._values indices = A._indices return coo_array((data, indices), shape=A.shape) + print(coo(A).shape, C.reshape(-1, 1).shape, C.shape) A = bmat([[coo(A), C.reshape(-1,1)], [C, None]], format='coo') A = COOTensor(bm.stack([A.row, A.col], axis=0), A.data, spshape=A.shape) diff --git a/fealpy/mesh/lagrange_quadrangle_mesh.py b/fealpy/mesh/lagrange_quadrangle_mesh.py index 7a35a5ffe..0344a02e1 100644 --- a/fealpy/mesh/lagrange_quadrangle_mesh.py +++ b/fealpy/mesh/lagrange_quadrangle_mesh.py @@ -26,13 +26,16 @@ def __init__(self, node: TensorLike, cell: TensorLike, p=1, surface=None, self.construct() self.meshtype = 'lquad' - self.linearmesh = None # 网格的顶点必须在球面上 + self.linearmesh = None self.nodedata = {} self.edgedata = {} self.celldata = {} self.meshdata = {} + #self.cell_shape_function = self.shape_function + #self.cell_grad_shape_function = self.grad_shape_function + def reference_cell_measure(self): return 1 @@ -43,11 +46,11 @@ def interpolation_points(self, p: int, index: Index=_S): @classmethod def from_quadrangle_mesh(cls, mesh, p: int, surface=None): - bnode = mesh.entity('node') + init_node = mesh.entity('node') node = mesh.interpolation_points(p) cell = mesh.cell_to_ipoint(p) if surface is not None: - bnode[:],_ = surface.project(bnode) + init_node[:],_ = surface.project(init_node) node,_ = surface.project(node) lmesh = cls(node, cell, p=p, construct=False) @@ -79,43 +82,82 @@ def bc_to_point(self, bc: TensorLike, index: Index=_S, etype='cell'): return p # shape function - def shape_function(self, bc: TensorLike, p=None): + def shape_function(self, bc: TensorLike, p: int=1): """ @berif bc 是一个长度为 TD 的 tuple 数组 bc[i] 是一个一维积分公式的重心坐标数组 假设 bc[0]==bc[1]== ... ==bc[TD-1] """ - return self.linearmesh.cell_shape_function(bc, p) + return self.shape_function(bc, p=p) - def grad_shape_function(self, bc: TensorLike, p=None, + def grad_shape_function(self, bc: TensorLike, p: int=1, index: Index=_S, variables='x'): - return self.linearmesh.cell_grad_shape_function(bc, p, index, variables) + return self.grad_shape_function(bc, p=p, variables=variables) # ipoints def number_of_local_ipoints(self, p:int, iptype:Union[int, str]='cell'): """ - @berif 每个ltri单元上插值点的个数 + @berif 每个lquad单元上插值点的个数 """ - #TODO - pass + if isinstance(iptype, str): + iptype = estr2dim(self, iptype) + return tensor_ldof(p, iptype) - def number_of_global_ipoints(self, p:int): + def number_of_global_ipoints(self, p:int) -> int: """ - @berif ltri网格上插值点总数 + @berif lquad网格上插值点总数 """ - # TODO - pass + num = [self.entity(i).shape[0] for i in range(self.TD+1)] + return tensor_gdof(p, num) def cell_to_ipoint(self, p:int, index:Index=_S): """ - @berif 获取单元与插值点的对应关系 + @brief 获取单元上的双 p 次插值点 """ - #TODO 复制粘贴 QuadrangleMesh 的代码并测试 - pass + sp = self.p + cell = self.cell[:, [0, -sp-1, 1]] # 取角点 + + if p == 0: + return bm.arange(len(cell)).reshape((-1, 1))[index] + if p == 1: + return cell[index, [0, 3, 1, 2]] + + cell2ipoint = bm.zeros((NC, (p + 1) * (p + 1)), dtype=self.itype, device=bm.get_device(cell)) + c2p = cell2ipoint.reshape((NC, p + 1, p + 1)) + e2p = self.edge_to_ipoint(p) + + flag = edge2cell[:, 2] == 0 + c2p = bm.set_at(c2p, (edge2cell[flag, 0], slice(None), 0), e2p[flag]) + + flag = edge2cell[:, 2] == 1 + c2p = bm.set_at(c2p, (edge2cell[flag, 0], -1, slice(None)), e2p[flag]) + + flag = edge2cell[:, 2] == 2 + c2p = bm.set_at(c2p, (edge2cell[flag, 0], slice(None), -1), bm.flip(e2p[flag], axis=-1)) + + flag = edge2cell[:, 2] == 3 + c2p = bm.set_at(c2p, (edge2cell[flag, 0], 0, slice(None)), bm.flip(e2p[flag], axis=-1)) + + iflag = edge2cell[:, 0] != edge2cell[:, 1] + flag = iflag & (edge2cell[:, 3] == 0) + c2p = bm.set_at(c2p, (edge2cell[flag, 1], slice(None), 0), bm.flip(e2p[flag], axis=-1)) + + flag = iflag & (edge2cell[:, 3] == 1) + c2p = bm.set_at(c2p, (edge2cell[flag, 1], -1, slice(None)), bm.flip(e2p[flag], axis=-1)) + + flag = iflag & (edge2cell[:, 3] == 2) + c2p = bm.set_at(c2p, (edge2cell[flag, 1], slice(None), -1), e2p[flag]) + + flag = iflag & (edge2cell[:, 3] == 3) + c2p = bm.set_at(c2p, (edge2cell[flag, 1], 0, slice(None)), e2p[flag]) + + c2p = bm.set_at(c2p, (slice(None), slice(1, -1), slice(1, -1)), NN + NE * (p - 1) + + bm.arange(NC * (p - 1) * (p - 1)).reshape(NC, p - 1, p - 1)) + + return cell2ipoint[index] def face_to_ipoint(self, p: int, index: Index=_S): - #TODO test return self.edge_to_ipoint(p, index) def entity_measure(self, etype=2, index:Index=_S): @@ -167,7 +209,7 @@ def cell_unit_normal(self, bc: TensorLike, index: Index=_S): n /= l return n - def jacobi_matrix(self, bc: TensorLike, p=None, index: Index=_S, return_grad=False): + def jacobi_matrix(self, bc: TensorLike, index: Index=_S, return_grad=False): """ @berif 计算参考单元 （xi, eta) 到实际 Lagrange 四边形(x) 之间映射的 Jacobi 矩阵。 @@ -175,9 +217,10 @@ def jacobi_matrix(self, bc: TensorLike, p=None, index: Index=_S, return_grad=Fal """ TD = len(bc) entity = self.entity(TD, index) - gphi = self.grad_shape_function(bc, p=p, variables='u') - J = bm.einsum( - 'cin, cqim -> cqnm', + gphi = self.grad_shape_function(bc) + print('bc',bc.shape) + print('gphi', gphi.shape) + J = bm.einsum('cim, cqin -> cqmn', self.node[entity[index], :], gphi) #(NC,ldof,GD),(NC,NQ,ldof,TD) if return_grad is False: return J #(NC,NQ,GD,TD) @@ -191,8 +234,7 @@ def first_fundamental_form(self, bc: Union[TensorLike, Tuple[TensorLike]], Compute the first fundamental form of a mesh surface at integration points. """ TD = len(bc) - - J = self.jacobi_matrix(bc, index=index, + J = self.jacobi_matrix(bc, index=index, return_grad=return_grad) if return_grad: diff --git a/fealpy/mesh/lagrange_triangle_mesh.py b/fealpy/mesh/lagrange_triangle_mesh.py index c4639b627..a69acef7c 100644 --- a/fealpy/mesh/lagrange_triangle_mesh.py +++ b/fealpy/mesh/lagrange_triangle_mesh.py @@ -7,6 +7,7 @@ from .utils import simplex_gdof, simplex_ldof from .mesh_base import HomogeneousMesh, estr2dim from .triangle_mesh import TriangleMesh +#from fealpy.functionspace.dofs import LinearMeshCFEDof class LagrangeTriangleMesh(HomogeneousMesh): def __init__(self, node: TensorLike, cell: TensorLike, p=1, surface=None, @@ -32,7 +33,7 @@ def __init__(self, node: TensorLike, cell: TensorLike, p=1, surface=None, self.construct() self.meshtype = 'ltri' - self.linearmesh = None # 网格的顶点必须在球面上 + self.linearmesh = None self.nodedata = {} self.edgedata = {} @@ -52,7 +53,6 @@ def generate_local_lagrange_edges(self, p: int) -> TensorLike: localEdge = bm.zeros((3, p+1), dtype=bm.int32) localEdge[2, :], = bm.where(multiIndex[:, 2] == 0) localEdge[1,:] = bm.flip(bm.where(multiIndex[:, 1] == 0)[0]) - #localEdge[1, -1::-1], = bm.where(multiIndex[:, 1] == 0) localEdge[0, :], = bm.where(multiIndex[:, 0] == 0) return localEdge @@ -96,11 +96,13 @@ def interpolation_points(self, p: int, index: Index=_S): @classmethod def from_triangle_mesh(cls, mesh, p: int, surface=None): - bnode = mesh.entity('node') + init_node = mesh.entity('node') + #cls.dof = LinearMeshCFEDof(mesh, p=p) + node = mesh.interpolation_points(p) cell = mesh.cell_to_ipoint(p) if surface is not None: - bnode[:], _ = surface.project(bnode) + init_node[:], _ = surface.project(init_node) node, _ = surface.project(node) lmesh = cls(node, cell, p=p, construct=True) @@ -143,7 +145,8 @@ def shape_function(self, bc: TensorLike, p=None): phi = bm.simplex_shape_function(bc, p=p) return phi[None, :, :] - def grad_shape_function(self, bc: TensorLike, p=None, index: Index=_S, variables='x'): + def grad_shape_function(self, bc: TensorLike, p=None, + index: Index=_S, variables='x'): """ @berif 计算单元形函数关于参考单元变量 u=(xi, eta) 或者实际变量 x 梯度。 @@ -152,7 +155,7 @@ def grad_shape_function(self, bc: TensorLike, p=None, index: Index=_S, variables lambda_2 = eta """ - p = self.p if p is None else p + p = self.p if p is None else p TD = bc.shape[-1] - 1 if TD == 2: Dlambda = bm.array([[-1, -1], [1, 0], [0, 1]], dtype=bm.float64) @@ -298,16 +301,15 @@ def cell_unit_normal(self, bc: TensorLike, index: Index=_S): n /= l return n - def jacobi_matrix(self, bc: TensorLike, p=None, index: Index=_S, return_grad=False): + def jacobi_matrix(self, bc: TensorLike, index: Index=_S, return_grad=False): """ @berif 计算参考单元 （xi, eta) 到实际 Lagrange 三角形(x) 之间映射的 Jacobi 矩阵。 x(xi, eta) = phi_0 x_0 + phi_1 x_1 + ... + phi_{ldof-1} x_{ldof-1} """ - TD = bc.shape[-1] - 1 entity = self.entity(TD, index) - gphi = self.grad_shape_function(bc, p=p, variables='u') + gphi = self.grad_shape_function(bc, variables='u') J = bm.einsum( 'cin, cqim -> cqnm', self.node[entity[index], :], gphi) #(NC,ldof,GD),(NC,NQ,ldof,TD) diff --git a/test/mesh/test_lagrange_quadrangle_mesh.py b/test/mesh/test_lagrange_quadrangle_mesh.py index 14596bc1b..fc9482f56 100644 --- a/test/mesh/test_lagrange_quadrangle_mesh.py +++ b/test/mesh/test_lagrange_quadrangle_mesh.py @@ -33,12 +33,21 @@ def test_cell_area(self, data, backend): surface = SphereSurface() #以原点为球心，1 为半径的球 mesh = QuadrangleMesh.from_unit_sphere_surface() - # 计算收敛阶 - maxit = 4 + # 计算收敛阶 + p = 1 + lmesh = LagrangeQuadrangleMesh.from_quadrangle_mesh(mesh, p=p, surface=surface) + + cm[i] = np.sum(lmesh.cell_area()) + + x = bm.to_numpy(cm[i]) + y = data["sphere_cm"] + em[i] = np.abs(x - y) # absolute error + """ + maxit = 1 cm = np.zeros(maxit, dtype=np.float64) em = np.zeros(maxit, dtype=np.float64) for i in range(maxit): - lmesh = LagrangeQuadrangleMesh.from_quadrangle_mesh(mesh, p=3, surface=surface) + lmesh = LagrangeQuadrangleMesh.from_quadrangle_mesh(mesh, p=p, surface=surface) cm[i] = np.sum(lmesh.cell_area()) @@ -51,9 +60,10 @@ def test_cell_area(self, data, backend): em_ratio = em[0:-1] / em[1:] print("unit_sphere:", em_ratio) + """ if __name__ == "__main__": a = TestLagrangeQuadrangleMeshInterfaces() - a.test_surface_mesh('numpy') - #a.test_cell_area(cell_area_data[0], 'numpy') + #a.test_surface_mesh('numpy') + a.test_cell_area(cell_area_data[0], 'numpy') #pytest.main(["./test_lagrange_triangle_mesh.py"]) diff --git a/test/mesh/test_lagrange_triangle_mesh.py b/test/mesh/test_lagrange_triangle_mesh.py index 3e8576671..93e23aa85 100644 --- a/test/mesh/test_lagrange_triangle_mesh.py +++ b/test/mesh/test_lagrange_triangle_mesh.py @@ -126,11 +126,11 @@ def test_error(self, backend): if __name__ == "__main__": - a = TestLagrangeTriangleMeshInterfaces() + #a = TestLagrangeTriangleMeshInterfaces() #a.test_init_mesh(init_data[0], 'numpy') #a.test_from_triangle_mesh(from_triangle_mesh_data[0], 'numpy') - a.test_surface_mesh('numpy') + #a.test_surface_mesh('numpy') #a.test_cell_area(cell_area_data[0], 'numpy') #a.test_(cell_[0], 'numpy') #a.test_error('numpy') - #pytest.main(["./test_lagrange_triangle_mesh.py"]) + pytest.main(["./test_lagrange_triangle_mesh.py"])