Poisson's Equation

[adtzlr]

Poisson's Equation in FElupe

this should be added to the docs as an example using https://felupe.readthedocs.io/en/latest/howto/forms.html.

Poisson's equation with fixed boundaries

$$ \text{div}(\boldsymbol{\nabla} v) + f = 0 \quad \text{in} \quad \Omega$$

$$ v = 0 \quad \text{on} \quad \Gamma_v$$

is transformed into integral form representation by the divergence (Gauss's) theorem.

$$ \int_\Omega \boldsymbol{\nabla} v \cdot \boldsymbol{\nabla} u \ d\Omega = \int_\Omega f \cdot v \ d\Omega$$

import felupe as fem

mesh = fem.Rectangle(n=2**6).triangulate()
region = fem.RegionTriangle(mesh)
field = fem.FieldContainer([fem.Field(region, dim=1)])

@fem.Form(v=field, u=field, grad_v=[True], grad_u=[True])
def bilinearform():
    return [lambda gradv, gradu: fem.math.ddot(gradv, gradu)]

# rhs of poisson's equation with unit load
@fem.Form(v=field, grad_v=[False])
def linearform():
    return [lambda v: 1.0 * v]

b = linearform.assemble(v=field)
A = bilinearform.assemble(v=field, u=field)

bounds = {
    "bottom-left":fem.Boundary(field[0], fx=0, fy=0, mode="or"),
    "top-right":fem.Boundary(field[0], fx=1, fy=1, mode="or"),
}

dof0, dof1 = fem.dof.partition(field, bounds)
system = fem.solve.partition(field, A,  dof1, dof0, r=-b)
field += fem.solve.solve(*system)

fem.ViewMesh(
    mesh=mesh, 
    point_data={"Field": field[0].values}
).plot("Field", show_undeformed=False, show_edges=False).show()

[adtzlr]

For the Newton-method to converge (to enable the "force"-balance check) it is necessary to include a linear form of the Poisson equation. This is also true for linear elasticity.

Added in dfd18c7.

[adtzlr]

This one uses an alternative approach.

import felupe as fem
import numpy as np

mesh = fem.Rectangle(n=21)
region = fem.RegionQuad(mesh, uniform=True)
field = fem.FieldContainer([fem.Field(region, dim=1)])

boundaries = dict(
    bottom=fem.Boundary(field[0], fy=0),
    top=fem.Boundary(field[0], fy=1),
    left=fem.Boundary(field[0], fx=0),
    right=fem.Boundary(field[0], fx=1),
)

# K = ∫ ∇_j(δu_i) ∇_j(u_k) dΩ
# K = ∫ ∇_j(δu_i) : (δ_ik δ_jl) : ∇_l(u_k) dΩ
lhs = fem.IntegralForm(
    fun=[fem.math.cdya_ik(np.eye(1), np.eye(2)).reshape(1, 2, 1, 2, 1, 1)],
    v=field,
    u=field,
    dV=region.dV,
    grad_v=[True],
    grad_u=[True],
)

# f = ∫ δu_i dΩ
rhs = fem.IntegralForm(fun=[np.ones((1, 1, 1))], v=field, dV=region.dV, grad_v=[False])

dof0, dof1 = fem.dof.partition(field, boundaries)
system = fem.solve.partition(field, lhs.assemble(), dof1, dof0, r=-rhs.assemble())
u = fem.solve.solve(*system)

view = mesh.view(point_data={"T": u})
view.plot("T").show()