### 1D problem with Neumann Boundary Conditions

98
views
0
3 months ago by

I am working on applying Neumann Boundary Conditions to a 1D problem with no luck. A minimal working example is below. Any tips on what I am doing wrong would be greatly appreciated.

from fenics import *
import numpy as np

n = 200

mesh = UnitIntervalMesh(n)
V = FunctionSpace(mesh, 'CG', 1)
v = TestFunction(V)
rho = Function(V)

class LeftEdge(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0)

class RightEdge(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1)

boundaryLeft = LeftEdge()
boundaryRight = RightEdge()

boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
boundaryLeft.mark(boundaries,1)
boundaryRight.mark(boundaries,2)

ds = ds(subdomain_data = boundaries)

G = Constant(1)

F = rho * v * dx \
+ v * G * ds(1) \
- v * G * ds(2)

solve(F == 0, rho)

rho = rho.vector().get_local()

print('Resulting BC', np.diff(rho, 1./n)[[0,-1]], 'Applied BC', G.values()[0])
Community: FEniCS Project

That's right, from you weak form, it seems you want to solve $u = 0$ with Dirichlet boundary conditions. That's why you get a highly oscillatory solution – this is more obvious if you use a discontinuous function space. Also, np.diff() gives you the difference, not the derivative.