From a1a756efe61d8318efd7830ef47be0b3c8a78a08 Mon Sep 17 00:00:00 2001 From: wpx <1143615697@qq.com> Date: Wed, 25 Dec 2024 16:44:36 +0800 Subject: [PATCH] update --- app/tssim/3d_NS/main_ipcs.py | 49 +++++-------- app/tssim/NS-CH-GNBC/main.py | 15 +--- app/tssim/NS-CH-GNBC/pde.py | 15 ++-- app/tssim/NS-CH-GNBC/solver.py | 7 +- fealpy/cfd/ns_fem_solver.py | 43 ++++++++++- fealpy/fem/__init__.py | 3 +- .../fem/fluid_boundary_friction_integrator.py | 59 +-------------- fealpy/fem/scalar_diffusion_integrator.py | 33 +++++++-- fealpy/fem/viscous_work_integrator.py | 66 +++++++++++++++++ fealpy/functionspace/lagrange_fe_space.py | 73 +++++++++++++++++++ fealpy/functionspace/tensor_space.py | 13 ++++ fealpy/old/fvm/scalar_diffusion_integrator.py | 2 - 12 files changed, 250 insertions(+), 128 deletions(-) create mode 100644 fealpy/fem/viscous_work_integrator.py diff --git a/app/tssim/3d_NS/main_ipcs.py b/app/tssim/3d_NS/main_ipcs.py index 18ca5b1a4..e46848cc6 100644 --- a/app/tssim/3d_NS/main_ipcs.py +++ b/app/tssim/3d_NS/main_ipcs.py @@ -17,7 +17,7 @@ from fealpy.solver import spsolve output = './' -T = 10 +T = 1 nt = 500 n = 16 @@ -31,58 +31,43 @@ space = LagrangeFESpace(mesh, p=2) uspace = TensorFunctionSpace(space, (2,-1)) -solver = NSFEMSolver(pde, mesh, pspace, space, dt, q=5) +solver = NSFEMSolver(pde, mesh, pspace, uspace, dt, q=5) ugdof = uspace.number_of_global_dofs() pgdof = pspace.number_of_global_dofs() +u0 = uspace.function() +us = uspace.function() u1 = uspace.function() +p0 = pspace.function() p1 = pspace.function() -''' -ipoint = space.interpolation_points() -import matplotlib.pylab as plt -fig = plt.figure() -axes = fig.gca() -mesh.add_plot(axes) -#mesh.find_edge(axes,fontsize=20,showindex=True) -mesh.find_node(axes,node=ipoint,fontsize=20,showindex=True) -plt.show() -''' + fname = output + 'test_'+ str(0).zfill(10) + '.vtu' mesh.nodedata['u'] = u1.reshape(2,-1).T mesh.nodedata['p'] = p1 mesh.to_vtk(fname=fname) -''' + BC = DirichletBC(space=uspace, gd=pde.velocity, threshold=pde.is_u_boundary, method='interp') -''' -BForm = solver.IPCS_BForm_0(None) -A = BForm.assembly() -print(A.to_dense()) -#A = BC.apply_matrix(A) -print(bm.sum(bm.abs(A.to_dense()))) - - - -exit() -LForm = solver.Ossen_LForm() +BForm0 = solver.IPCS_BForm_0(None) +LForm0 = solver.IPCS_LForm_0() +A0 = BForm0.assembly() -for i in range(10): +for i in range(1): t = timeline.next_time_level() print(f"第{i+1}步") print("time=", t) - solver.NS_update(u1) - A = BForm.assembly() - b = LForm.assembly() - A,b = BC.apply(A,b) + solver.update_ipcs_0(u0, p0) + print(bm.sum(bm.abs(A0.to_dense()))) + b0 = LForm0.assembly() + A0,b0 = BC.apply(A0,b0) + print(bm.sum(bm.abs(b0))) - x = spsolve(A, b, 'mumps') - u1[:] = x[:ugdof] - p1[:] = x[ugdof:] + us[:] = spsolve(A0, b0, 'mumps') fname = output + 'test_'+ str(i+1).zfill(10) + '.vtu' mesh.nodedata['u'] = u1.reshape(2,-1).T diff --git a/app/tssim/NS-CH-GNBC/main.py b/app/tssim/NS-CH-GNBC/main.py index 008cf6ae2..c9bd0e6ad 100644 --- a/app/tssim/NS-CH-GNBC/main.py +++ b/app/tssim/NS-CH-GNBC/main.py @@ -28,7 +28,8 @@ #bm.set_default_device('cuda') output = './' -h = 1/256 +#h = 1/256 +h = 1/10 T = 2 nt = int(T/(0.1*h)) @@ -47,20 +48,8 @@ space = LagrangeFESpace(mesh, p=2) uspace = TensorFunctionSpace(space, (2,-1)) -''' -ipoint = space.interpolation_points() -import matplotlib.pylab as plt -fig = plt.figure() -axes = fig.gca() -mesh.add_plot(axes) -#mesh.find_edge(axes,fontsize=20,showindex=True) -mesh.find_node(axes,node=ipoint,fontsize=20,showindex=True) -plt.show() -''' - solver = Solver(pde, mesh, pspace, phispace, uspace, dt, q=5) - u0 = uspace.function() u1 = uspace.function() u2 = uspace.function() diff --git a/app/tssim/NS-CH-GNBC/pde.py b/app/tssim/NS-CH-GNBC/pde.py index ff6fe1460..5bede971a 100644 --- a/app/tssim/NS-CH-GNBC/pde.py +++ b/app/tssim/NS-CH-GNBC/pde.py @@ -12,25 +12,21 @@ from fealpy.mesh import TriangleMesh class CouetteFlow: - ''' - @brief: CouetteFlow - ''' - def __init__(self, eps=1e-10, h=1/256): self.eps = eps ## init the parameter self.R = 5.0 ##dimensionless self.l_s = 0.0025 ##dimensionless slip length - self.L_s = self.l_s / 100 + self.L_s = self.l_s self.epsilon = 0.004 ## the thickness of interface self.L_d = 0.0005 ##phenomenological mobility cofficient self.lam = 12.0 ##dimensionless self.V_s = 200.0 ##dimensionless self.s = 2.5 ##stablilizing parameter - #self.theta_s = bm.array(bm.pi/2) - self.theta_s = bm.array(77.6/180 * bm.pi) + self.theta_s = bm.array(bm.pi/2) + #self.theta_s = bm.array(77.6/180 * bm.pi) self.h = h def mesh(self): @@ -56,6 +52,11 @@ def is_uy_Dirichlet(self,p): return (bm.abs(p[..., 1] - 0.125) < self.eps) | \ (bm.abs(p[..., 1] + 0.125) < self.eps) + @cartesian + def is_slip_boundary(self,p): + return self.is_wall_boundary(p) + #return bm.logical_not(self.is_wall_boundary(p))) + @cartesian def init_phi(self,p): x = p[..., 0] diff --git a/app/tssim/NS-CH-GNBC/solver.py b/app/tssim/NS-CH-GNBC/solver.py index 3f564a526..6276368f9 100644 --- a/app/tssim/NS-CH-GNBC/solver.py +++ b/app/tssim/NS-CH-GNBC/solver.py @@ -32,7 +32,7 @@ def __init__(self, pde, mesh, pspace, phispace, uspace, dt, q=5): self.pde = pde self.dt = dt self.q = q - + def CH_BForm(self): phispace = self.phispace dt = self.dt @@ -54,7 +54,8 @@ def CH_BForm(self): A10 = BilinearForm(phispace) A10.add_integrator(ScalarDiffusionIntegrator(coef=-epsilon, q=q)) A10.add_integrator(ScalarMassIntegrator(coef=-s/epsilon, q=q)) - A10.add_integrator(BoundaryFaceMassIntegrator(coef=-3/(2*dt*V_s), q=q, threshold=self.pde.is_wall_boundary)) + A10.add_integrator(BoundaryFaceMassIntegrator(coef=-3/(2*dt*V_s), + q=q, threshold=self.pde.is_slip_boundary)) A11 = BilinearForm(phispace) A11.add_integrator(ScalarMassIntegrator(coef=1, q=q)) @@ -191,6 +192,7 @@ def NS_update(self, u_0, u_1, mu_2, phi_2, phi_1): lam = self.pde.lam epsilon = self.pde.epsilon normal = self.mesh.edge_unit_normal() + normal[..., 1] = 1 tangent = self.mesh.edge_unit_tangent() tangent[..., 0] = 1 @@ -214,6 +216,7 @@ def u_SI_coef(bcs, index): def u_BF_SI_coef(bcs, index): L_phi = epsilon*bm.einsum('eld, ed -> el', phi_2.grad_value(bcs, index), normal[index,:]) + print(phi_2.grad_value(bcs, index)) L_phi -= 2*(bm.sqrt(bm.array(2))/6)*bm.pi*bm.cos(theta_s)*bm.cos((bm.pi/2)*phi_2(bcs, index)) L_phi += (bm.sqrt(bm.array(2))/6)*bm.pi*bm.cos(theta_s)*bm.cos((bm.pi/2)*phi_1(bcs, index)) diff --git a/fealpy/cfd/ns_fem_solver.py b/fealpy/cfd/ns_fem_solver.py index 8b853d41f..5df1dd26e 100644 --- a/fealpy/cfd/ns_fem_solver.py +++ b/fealpy/cfd/ns_fem_solver.py @@ -16,7 +16,7 @@ SourceIntegrator, PressWorkIntegrator, FluidBoundaryFrictionIntegrator, - FluidBoundaryFrictionIntegratorP) + ViscousWorkIntegrator) from fealpy.fem import (BoundaryFaceMassIntegrator, BoundaryFaceSourceIntegrator) @@ -94,8 +94,43 @@ def IPCS_BForm_0(self, threshold): q = self.q Bform = BilinearForm(uspace) - M = ScalarMassIntegrator(coef=R , q=q) - S = ScalarDiffusionIntegrator(coef=R , q=q) - F = FluidBoundaryFrictionIntegrator(coef=1, q=5, threshold=threshold) + M = ScalarMassIntegrator(coef=R/dt, q=q) + F = FluidBoundaryFrictionIntegrator(coef=-1, q=q, threshold=threshold) + VW = ViscousWorkIntegrator(coef=2, q=q) + Bform.add_integrator(VW) Bform.add_integrator(F) + Bform.add_integrator(M) return Bform + + def IPCS_LForm_0(self, pthreshold=None): + pspace = self.pspace + uspace = self.uspace + dt = self.dt + R = self.pde.R + q = self.q + + Lform = LinearForm(uspace) + self.ipcs_lform_SI = SourceIntegrator(q=q) + self.ipcs_lform_SI.keep_data() + self.ipcs_lform_BSI = BoundaryFaceSourceIntegrator(q=q, threshold=pthreshold) + self.ipcs_lform_BSI.keep_data() + return Lform + + def update_ipcs_0(self, u0, p0): + dt = self.dt + R = self.pde.R + + def coef(bcs, index): + result = 1/dt*u0(bcs, index) + result += np.einsum('...j, ....,ij -> ...i', u0(bcs, index), u0.grad_value(bcs, index)) + result += np.repeat(p0(bcs,index)[...,np.newaxis], 2, axis=-1) + return + + def B_coef(bcs, index): + result = np.einsum('..ij, ....j->...ij', p(bcs, index), self.mesh.edge_unit_normal(bcs, index)) + return + + self.ipcs_lform_SI.source = coef + self.ipcs_lform_BSI.source = B_coef + + diff --git a/fealpy/fem/__init__.py b/fealpy/fem/__init__.py index 989a15d9e..dcafcb367 100644 --- a/fealpy/fem/__init__.py +++ b/fealpy/fem/__init__.py @@ -20,7 +20,7 @@ from .curlcurl_integrator import CurlCurlIntegrator from .nonlinear_elastic_integrator import NonlinearElasticIntegrator from .div_integrator import DivIntegrator - +from .viscous_work_integrator import ViscousWorkIntegrator ### Cell Source from .cell_source_integrator import CellSourceIntegrator @@ -33,7 +33,6 @@ from .scalar_robin_bc_integrator import ScalarRobinBCIntegrator from .face_mass_integrator import BoundaryFaceMassIntegrator, InterFaceMassIntegrator from .fluid_boundary_friction_integrator import FluidBoundaryFrictionIntegrator -from .fluid_boundary_friction_integrator import FluidBoundaryFrictionIntegratorP ### Face Source from .scalar_neumann_bc_integrator import ScalarNeumannBCIntegrator, ScalarRobinSourceIntegrator diff --git a/fealpy/fem/fluid_boundary_friction_integrator.py b/fealpy/fem/fluid_boundary_friction_integrator.py index e5544281f..30b9a1e68 100644 --- a/fealpy/fem/fluid_boundary_friction_integrator.py +++ b/fealpy/fem/fluid_boundary_friction_integrator.py @@ -81,64 +81,7 @@ def assembly(self, space: _FS): mesh = getattr(space, 'mesh', None) bcs, ws, phi, gphi, fm, index, n = self.fetch(space) val = process_coef_func(coef, bcs=bcs, mesh=mesh, etype='face', index=index) - gphin = bm.einsum('ej, eqlj->eql', n, gphi) + gphin = bm.einsum('e...i, eql...ij->eql...j', n, gphi) result = bilinear_integral(phi, gphin, ws, fm, val, batched=self.batched) return result -class FluidBoundaryFrictionIntegratorP(LinearInt, OpInt, FaceInt): - def __init__(self, coef: Optional[CoefLike]=None, q: Optional[int]=None, *, - threshold: Optional[Threshold]=None, - batched: bool=False, mesh): - super().__init__() - self.coef = coef - self.q = q - self.threshold = threshold - self.batched = batched - self.mesh = mesh - - def make_index(self, space: _FS) -> TensorLike: - threshold = self.threshold - mesh = self.mesh - if isinstance(threshold, TensorLike): - index = threshold - else: - index = mesh.boundary_face_index() - if callable(threshold): - bc = mesh.entity_barycenter('face', index=index) - index = index[threshold(bc)] - return index - - @enable_cache - def to_global_dof(self, space: _FS) -> TensorLike: - index = self.make_index(space) - return space.face_to_dof(index=index) - - @enable_cache - def fetch(self, space: _FS): - space0 = space[0] - space1 = space[1] - index = self.make_index(space) - mesh = self.mesh - - if not isinstance(mesh, HomogeneousMesh): - raise RuntimeError("The ScalarRobinBCIntegrator only support spaces on" - f"homogeneous meshes, but {type(mesh).__name__} is" - "not a subclass of HomoMesh.") - - facemeasure = mesh.entity_measure('face', index=index) - q = space.p+3 if self.q is None else self.q - qf = mesh.quadrature_formula(q, 'face') - bcs, ws = qf.get_quadrature_points_and_weights() - phi1 = space[1].face_basis(bcs, index) - phi0 = space[0].face_basis(bcs, index) - n = mesh.edge_normal(index) - return bcs, ws, phi0, phi1, facemeasure, index, n - - def assembly(self, space: _FS): - coef = self.coef - mesh = getattr(space, 'mesh', None) - bcs, ws, phi0, phi1, fm, index, n = self.fetch(space) - val = process_coef_func(coef, bcs=bcs, mesh=mesh, etype='face', index=index) - phi0n = bm.einsum('ej, e...j->e...j', n, phi0) - result = bilinear_integral(phi1, phi0n, ws, fm, val, batched=self.batched) - return result diff --git a/fealpy/fem/scalar_diffusion_integrator.py b/fealpy/fem/scalar_diffusion_integrator.py index 28104370b..e1bf99000 100644 --- a/fealpy/fem/scalar_diffusion_integrator.py +++ b/fealpy/fem/scalar_diffusion_integrator.py @@ -1,4 +1,5 @@ from typing import Optional, Literal +from ..mesh.mesh_base import SimplexMesh from ..backend import backend_manager as bm from ..typing import TensorLike, Index, _S @@ -62,7 +63,6 @@ def assembly(self, space: _FS, /, indices=None) -> TensorLike: index = self.entity_selection(indices) coef = process_coef_func(coef, bcs=bcs, mesh=mesh, etype='cell', index=index) gphi = self.fetch_gphix(space, indices) - return bilinear_integral(gphi, gphi, ws, cm, coef, batched=self.batched) @assemblymethod('fast') @@ -72,13 +72,30 @@ def fast_assembly(self, space: _FS, /, indices=None) -> TensorLike: TODO: 加入 assert """ mesh = space.mesh - _, ws = self.fetch_qf(space) - gphi = self.fetch_gphiu(space, indices) - M = bm.einsum('q, qik, qjl -> ijkl', ws, gphi, gphi) - cm = self.fetch_measure(space, indices) - glambda = mesh.grad_lambda(index=self.entity_selection(indices)) - A = bm.einsum('ijkl, ckm, clm, c -> cij', M, glambda, glambda, cm) - return A + if isinstance(mesh, SimplexMesh): + + _, ws = self.fetch_qf(space) + gphi = self.fetch_gphiu(space, indices) + M = bm.einsum('q, qik, qjl -> ijkl', ws, gphi, gphi) + cm = self.fetch_measure(space, indices) + glambda = mesh.grad_lambda(index=self.entity_selection(indices)) + result = bm.einsum('ijkl, ckm, clm, c -> cij', M, glambda, glambda, cm) + else: + coef = self.coef + mesh = space.mesh + index = self.entity_selection(indices) + cm = self.fetch_measure(space, indices) + bcs,ws = self.fetch_qf(space) + coef = process_coef_func(coef, bcs=bcs, mesh=mesh, etype='cell', index=index) + + gphiu = self.fetch_gphiu(space, indices) + M = bm.einsum('qim, qjn, q -> qijmn', gphiu, gphiu, ws) + J = mesh.jacobi_matrix(bcs, index) + G = mesh.first_fundamental_form(J) + G = bm.linalg.inv(G) + JG = bm.einsum("cqkm, cqmn -> cqkn", J, G) + result = bm.einsum('cqkn, qijmn, cqkm, c -> cij', JG, M, JG, cm) # (NC, NQ, ldof, GD) + return result @assemblymethod('nonlinear') def nonlinear_assembly(self, space: _FS, /, indices=None) -> TensorLike: diff --git a/fealpy/fem/viscous_work_integrator.py b/fealpy/fem/viscous_work_integrator.py new file mode 100644 index 000000000..274b1c826 --- /dev/null +++ b/fealpy/fem/viscous_work_integrator.py @@ -0,0 +1,66 @@ +#!/usr/bin/python3 +'''! + @Author: wpx + @File Name: vector_viscous_work_integrator.py + @Mail: wpx15673207315@gmail.com + @Created Time: Tue 17 Dec 2024 02:50:59 PM CST + @bref + @ref +''' +from typing import Optional, Literal + +from ..backend import backend_manager as bm +from ..typing import TensorLike, Index, _S + +from ..mesh import HomogeneousMesh +from ..functionspace.space import FunctionSpace as _FS +from ..utils import process_coef_func +from ..functional import bilinear_integral, linear_integral, get_semilinear_coef +from .integrator import ( + LinearInt, OpInt, CellInt, + enable_cache, + assemblymethod, + CoefLike +) + + +class ViscousWorkIntegrator(LinearInt, OpInt, CellInt): + """ + construct the mu * (epslion(u), epslion(v)) fem matrix + epsion(u) = 1/2*(\\nabla u+ (\\nabla u).T) + """ + def __init__(self, coef: Optional[CoefLike] = None, q: Optional[int] = None, *, + region: Optional[TensorLike] = None, + batched: bool = False, + method: Optional[str]=None) -> None: + super().__init__(method=method if method else 'assembly') + self.coef = coef + self.q = q + self.set_region(region) + self.batched = batched + + @enable_cache + def to_global_dof(self, space: _FS, /, indices=None) -> TensorLike: + return space.cell_to_dof(index=self.entity_selection(indices)) + + @enable_cache + def fetch_assemble(self, space: _FS, /, indices=None): + mesh = space.mesh + index=self.entity_selection(indices) + q = space.p+3 if self.q is None else self.q + qf = mesh.quadrature_formula(q, 'cell') + bcs, ws = qf.get_quadrature_points_and_weights() + gphi = space.grad_basis(bcs, index) + cm = mesh.entity_measure('cell', index=index) + return bcs, ws, gphi, cm + + def assembly(self, space: _FS, /, indices=None) -> TensorLike: + coef = self.coef + mesh = space.mesh + bcs, ws, gphi, cm = self.fetch_assemble(space) + index=self.entity_selection(indices) + coef = process_coef_func(coef, bcs=bcs, mesh=mesh, etype='cell', index=index) + sym_gphi = 0.5*(bm.swapaxes(gphi, -2, -1) + gphi) + result = bilinear_integral(sym_gphi, sym_gphi, ws, cm, coef, batched=self.batched) + return result + diff --git a/fealpy/functionspace/lagrange_fe_space.py b/fealpy/functionspace/lagrange_fe_space.py index 76194a7ca..e92d95912 100644 --- a/fealpy/functionspace/lagrange_fe_space.py +++ b/fealpy/functionspace/lagrange_fe_space.py @@ -198,3 +198,76 @@ def grad_value(self, uh: TensorLike, bc: TensorLike, index: Index=_S) -> TensorL e2dof = self.dof.entity_to_dof(TD, index=index) val = bm.einsum('cilm, cl -> cim', gphi, uh[e2dof]) return val + + def grad_recovery(self, uh: TensorLike, method: str='simple'): + """ + + Notes + ----- + + uh 是线性有限元函数，该程序把 uh 的梯度(分片常数）恢复到分片线性连续空间 + 中。 + + """ + GD = self.GD + cell2dof = self.cell_to_dof() + gdof = self.number_of_global_dofs() + ldof = self.number_of_local_dofs() + p = self.p + bc = self.dof.multiIndex/p + guh = uh.grad_value(bc) + guh = guh.swapaxes(0, 1) + rguh = self.function(dim=GD) + + if method == 'simple': + deg = bm.bincount(cell2dof.flat, minlength = gdof) + if GD > 1: + bm.add.at(rguh, (cell2dof, bm.s_[:]), guh) + else: + bm.add.at(rguh, cell2dof, guh) + + elif method == 'area': + measure = self.mesh.entity_measure('cell') + ws = bm.einsum('i, j->ij', measure,bm.ones(ldof)) + deg = bm.bincount(cell2dof.flat,weights = ws.flat, minlength = gdof) + guh = bm.einsum('ij..., i->ij...', guh, measure) + if GD > 1: + bm.add.at(rguh, (cell2dof, bm.s_[:]), guh) + else: + bm.add.at(rguh, cell2dof, guh) + + elif method == 'distance': + ipoints = self.interpolation_points() + bp = self.mesh.entity_barycenter('cell') + v = bp[:, bm.newaxis, :] - ipoints[cell2dof, :] + d = bm.sqrt(bm.sum(v**2, axis=-1)) + deg = bm.bincount(cell2dof.flat,weights = d.flat, minlength = gdof) + guh = bm.einsum('ij..., ij->ij...', guh, d) + if GD > 1: + bm.add.at(rguh, (cell2dof, bm.s_[:]), guh) + else: + bm.add.at(rguh, cell2dof, guh) + + elif method == 'area_harmonic': + measure = 1/self.mesh.entity_measure('cell') + ws = bm.einsum('i, j->ij', measure,bm.ones(ldof)) + deg = bm.bincount(cell2dof.flat,weights = ws.flat, minlength = gdof) + guh = bm.einsum('ij..., i->ij...', guh, measure) + if GD > 1: + bm.add.at(rguh, (cell2dof, bm.s_[:]), guh) + else: + bm.add.at(rguh, cell2dof, guh) + + elif method == 'distance_harmonic': + ipoints = self.interpolation_points() + bp = self.mesh.entity_barycenter('cell') + v = bp[:, bm.newaxis, :] - ipoints[cell2dof, :] + d = 1/bm.sqrt(bm.sum(v**2, axis=-1)) + deg = bm.bincount(cell2dof.flat,weights = d.flat, minlength = gdof) + guh = bm.einsum('ij..., ij->ij...',guh,d) + if GD > 1: + bm.add.at(rguh, (cell2dof, bm.s_[:]), guh) + else: + bm.add.at(rguh, cell2dof, guh) + rguh /= deg.reshape(-1, 1) + return rguh diff --git a/fealpy/functionspace/tensor_space.py b/fealpy/functionspace/tensor_space.py index 888618c7e..9483ee788 100644 --- a/fealpy/functionspace/tensor_space.py +++ b/fealpy/functionspace/tensor_space.py @@ -74,6 +74,19 @@ def face_basis(self, p: TensorLike, index: Index=_S, **kwargs) -> TensorLike: def grad_basis(self, p: TensorLike, index: Index=_S, **kwargs) -> TensorLike: gphi = self.scalar_space.grad_basis(p, index, **kwargs) return generate_tensor_grad_basis(gphi, self.dof_shape, self.dof_priority) + + + @barycentric + def cell_basis_on_edge(self, bc: TensorLike, eindex: TensorLike) -> TensorLike: + result = self.scalar_space.cell_basis_on_edge(bc, eindex) + return generate_tensor_basis(result, self.dof_shape, self.dof_priority) + + @barycentric + def cell_grad_basis_on_edge(self, bc: TensorLike, eindex: TensorLike, + isleft = True) -> TensorLike: + result = self.scalar_space.cell_grad_basis_on_edge(bc, eindex, isleft) + return generate_tensor_grad_basis(result, self.dof_shape, self.dof_priority) + def cell_to_dof(self, index: Index=_S) -> TensorLike: """Get the cell to dof mapping. diff --git a/fealpy/old/fvm/scalar_diffusion_integrator.py b/fealpy/old/fvm/scalar_diffusion_integrator.py index 6544e78b7..93fc62c01 100644 --- a/fealpy/old/fvm/scalar_diffusion_integrator.py +++ b/fealpy/old/fvm/scalar_diffusion_integrator.py @@ -39,8 +39,6 @@ def cell_center_matrix(self, bf, is_bf_boundary): val = -h/h*np.ones(I.shape) A += csr_matrix((val,(I, J)), shape=(NC, NC)) - - b = np.zeros(NC) index = np.where(flag)[0] index1 = cell2edge[flag]