Open Boundary Problem


126
views
1
9 weeks ago by
Dear Community,

to solve electrostatic problems, it would be perfect to have a feature for open boundaries. Could anybody guide me to an example to implement open boundaries in fenics?

Regards, Markus
Community: FEniCS Project

2 Answers


1
9 weeks ago by
There is a large literature on this, which I haven't read, so the following is likely neither novel nor ideal, but the approach of mapping a finite domain to an infinite one is made relatively convenient by UFL and appears to work surprisingly well for the Poisson problem:

from dolfin import *
mesh = UnitSquareMesh(50,50)

# (The mapping tends to pump up the estimated quadrature degree to an
# impractical level.)
dx = dx(metadata={"quadrature_degree":4})

# Convergence appears to hold for messy unstructured meshes too
#from mshr import *
#mesh = generate_mesh(Rectangle(Point(0,0),Point(1,1)),128)

# Use to define an exact solution decaying to zero at infinity
def gaussian(x,sigma):
    r = x[0]**2 + x[1]**2
    return exp(-0.5*((r/sigma)**2))/(sigma*sqrt(2.0*pi))

# Map the unit interval to the entire real line
def I_to_R(x):
    return tan(pi*(x-0.5))

# Map the unit square to R^2
def I2_to_R2(xhat):
    return as_vector([I_to_R(xhat[0]),I_to_R(xhat[1])])

# Map mesh coordinates to infinite domain
xhat = SpatialCoordinate(mesh)
x = I2_to_R2(xhat)
Dx = grad(x)

# Define grad and div w.r.t. x
def gradx(f):
    return dot(grad(f),inv(Dx))
def divx(f):
    n = rank(f)
    ii = indices(n)
    return as_tensor(gradx(f)[ii+(ii[n-1],)],ii[0:n-1])

# Integral change-of-variables
J = det(Dx)

# Manufacture source term for the Gaussian exact solution
u_exact = gaussian(x,0.2)
f = -divx(gradx(u_exact))

# Define weak Poisson problem
V = FunctionSpace(mesh,"Lagrange",2)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(gradx(u),gradx(v))*J*dx
L = inner(f,v)*J*dx

# Homogeneous Dirichlet BCs on parametric domain (so, at infinity in
# physical domain); this only really makes sense for a 2D electric potential
# if the net charge is zero, since the free-space Green's function does
# not decay to zero until at least three dimensions.  
bc = DirichletBC(V,0.0,"on_boundary")

# Solve
u = Function(V)
solve(a==L,u,bcs=[bc,])

# Check L^2 error on bi-unit square in physical domain; some quick
# mesh-doubling tests indicate that this converges optimally.
import math
chi = conditional(gt(Max(abs(x[0]),abs(x[1])),1.0),Constant(0.0),Constant(1.0))

# Error in full domain also appears to converge despite J blowing up at edges
#chi = Constant(1.0)

print("L^2 error: "+str(math.sqrt(assemble(chi*((u-u_exact)**2)*J*dx))))

print("Net charge: "+str(assemble(f*J*dx)))

####### Visualization #######

# Output for visualizing the solution
u.rename("u","u")
File("u.pvd") << u

# Using the standard project() function is oscillatory for singular functions
def lumpedProject(f):
    V = FunctionSpace(mesh,"Lagrange",1)
    v = TestFunction(V)
    lhs = assemble(inner(Constant(1.0),v)*dx)
    rhs = assemble(inner(f,v)*dx)
    u = Function(V)
    as_backend_type(u.vector())\
        .vec().pointwiseDivide(as_backend_type(rhs).vec(),
                               as_backend_type(lhs).vec())
    return u

# Regularize mapping used for visualization slightly, to avoid too much
# craziness near the edges:
cx = Constant((0.5,0.5))
x_reg = I2_to_R2(0.9999*(xhat-cx)+cx)
dispx = lumpedProject(x_reg[0]-xhat[0])
dispy = lumpedProject(x_reg[1]-xhat[1])
dispx.rename("dispx","dispx")
dispy.rename("dispy","dispy")
File("dispx.pvd") << dispx
File("dispy.pvd") << dispy

# To plot, load all three files, apply an Append Attributes filter to them, use
# a Calculator filter to define the vector field
#
#   dispx*iHat + dispy*jHat
#
# and use this field in a Warp by Vector filter.  Set the warping scale to 1,
# color by u, and zoom in to the domain of interest.  Using the
# "Surface with Edges" mode to visualize shows the effect of the
# geometrical mapping clearly.

​

You might also come up with a more complicated piecewise mapping that is identity within the region of interest, in which case the solution in that region would not require mesh warping to visualize correctly.  (It would also give you more direct control over the resolution there.)
0
8 weeks ago by
Thanks for the input. A good literature for electro statics would be this paper:

Imhoff, J. F.; Meunier G.; Sabonnadière, J. C: Finite element method of open boundary problems
IEEE Transactions on Magnetics, Vol. MAG-26 No. 2, 1990,
S. 588 - 591
Please login to add an answer/comment or follow this question.

Similar posts:
Search »