计算科学 - Fenics 中的子域积分 - 吾爱随笔录

Fenics 中的子域积分

计算科学芬尼克斯

2021-12-01 11:36:42

我想计算，但在 FEniCS 中看不到真正有效的方法。我的第一个想法是为两个不同的网格组装（f*v*dx）。（try1，下面），但似乎我无法组合 1D 和 2D 网格，并且 DG 网格必须在每个方向上至少有 2 个基函数，我不确定网格拓扑如何 -> 基系数布局是可控的。 $u(y) = \int f(x,y) dx$

我的下一个想法是解决 du/dx = f 一侧的 BC 为零，在 x = 1 上累积积分。那个一直通过组装问题，然后因奇异矩阵异常而失败。

在 FEniCS 中是否有更好的方法来做到这一点？为了解决这个问题，我建议在 FEniCS 中添加通用卷积，这样可以通过适当的矩阵乘法来完成。 . $u = \int h(\vec x-\vec y) v(\vec y) dy$

from dolfin import *

m = UnitSquareMesh(12, 12)
V = FunctionSpace(m, "Lagrange", 1)

f = Function(V)
f.interpolate(Expression("x[1]-x[0]")) # int_0^1 f dx = x[1] - 0.5

def try1():
    m2 = UnitSquareMesh(1, 12)
    V2 = FunctionSpace(m2, "DG", 0)
    u = Function(V2)
    print u.vector().array().shape
    v = TestFunction(V2)
    A = assemble(f*v*dx)
    print A.array()*12.0

def try2():
    du = TrialFunction(V)
    v = TestFunction(VectorFunctionSpace(m, "Lagrange", 1))

    bc = DirichletBC(V, Constant(0.0), lambda x,b: b and x[0] < 1e-8)
    #bc2 = DirichletBC(V, Constant(0.0), lambda x,b: b and x[1] < 1e-8)

    a = inner(grad(du), v)*dx
    L = f*v[0]*dx
    #L2 = f*v[1]*dx

    u = Function(V)

    solve(a == L, u, bcs=bc)

    print u.vector().array()

2个回答

我认为没有通用的方法可以做到这一点，因为的（离散）表示必须事先定义。如果您希望它以一维有限元为基础表示，那么您可以执行以下操作 $u(y)$

定义一维 FEM 基础 $\{\phi_i(y)\}$
对于每个定义一个二维表达式 $\{\phi_i(y)\}$ $b_i(x,y) = \phi_i(y)$
根据你的测试每个以获得 $b_i$ $f$ $\int \int b_{i} (x, y) f (x, y) d x d y = \int ϕ_{i} (y) u (y) d y =: {\tilde{u}}_{i}$ $\int \int b_i(x,y)f(x,y)dxdy=\int\phi_i(y)u(y)dy =: \tilde u_i$
计算在 FEM 基础上的展开系数 - - via 其中是 FEM 质量矩阵。 $u$ $u(y) \approx \sum_i u_i \phi_i(y)$ $[\begin{matrix} ⋮ \\ u_{i} \\ ⋮ \end{matrix}] = M_{ϕ}^{- 1} [\begin{matrix} ⋮ \\ {\tilde{u}}_{i} \\ ⋮ \end{matrix}],$ $\begin{bmatrix} \vdots \\ u_i \\ \vdots \end{bmatrix} = \mathcal M_{\phi}^{-1} \begin{bmatrix} \vdots \\ \tilde u_i \\ \vdots \end{bmatrix},$ $M_\phi$

如果你有b_ifenics 和 function中的表达式f，那么你会发现

tu_i = assemble(b_i*f*dx)

作为源来求解泊松方程，然后使用高斯定律计算每个分割平面的线段上的通量。该方法仍然不会被矢量化，但至少只需要 1 个网格函数（用于标记域）。如果我做出妥协，我还不如在 python 中迭代： $f$

def slice(mesh, n, L):
    subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
    subdomains.set_all(0)
    dx = L/float(n)
    for i in range(1, n):
            b = AutoSubDomain(lambda x: x[1] > i*dx-1e-8)
            b.mark(subdomains, i)
    return subdomains

def try3():
    markers = slice(m, 12, 1.0)
    dX = dx[markers]
    v = [assemble(f*dX(i))*12.0 for i in range(12)]
    print v

其它你可能感兴趣的问题

上一篇守恒方程的细胞平均值与细胞总数的有限体积下一篇具有已知依赖关系的并行任务的最优调度