### Robin + Neumann boundary conditions with marker

38
views
0
7 days ago by
Hello,

I'm trying to adapt what is described here: https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1005.html#setting-multiple-dirichlet-neumann-and-robin-conditions

I successfully did it in a rectangular domain, Newmann condition    $\frac{\left\{\partial u\right\}}{\partial n}=h${u}n =h   at the top, and Dirichlet conditions   $u=g$u=g  at the rest of  the boundary.

Now I'm trying over the same rectangular domain, with the same Newmann condition at the top, but Robin's condition at the rest : $\alpha u+\beta\frac{\partial u}{\partial n}=h$αu+βun =h
To start, I'm trying with  $g=h=0$g=h=0 and f a dirac's delta. But got this error:

/usr/local/lib/python3.6/dist-packages/matplotlib/colors.py:504: RuntimeWarning: invalid value encountered in less
xa[xa < 0] = -1​

Here is my code, with boundary integrals commented:

from fenics import *
from dolfin import *
import matplotlib.pyplot as plt
# plot = lambda *args, **kwargs: None
import numpy as np

#parameters
x0 = 0.0#Constant(0)
xN = 8.0#Constant(8)
y0 = -4.0#Constant(-4)
yM = 0.0#Constant(0)
nx = 50#Constant(10)
ny = 15#Constant(10)
hx = (xN-x0)/(nx-1)
hy = (yM-y0)/(ny-1)

# fuente y salida
in_x = 4.0
in_y = 0.0
out_x = 7.0
out_y = 0.0

d=1

#mesh
mesh = RectangleMesh(Point(x0,y0), Point(xN,yM), nx, ny, "right")
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
V = FunctionSpace(mesh, 'Lagrange', d)

#boundaries
class BoundaryX0(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[0], x0, self.tol)

class BoundaryXN(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[0], xN, self.tol)

class BoundaryY0(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[1], y0, self.tol)

class BoundaryYM(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[1], yM, self.tol)

bx0 = BoundaryX0()
bx0.mark(boundary_markers, 4)
bx1 = BoundaryXN()
bx1.mark(boundary_markers, 1)
by0 = BoundaryY0()
by0.mark(boundary_markers, 2)
byM = BoundaryYM()
byM.mark(boundary_markers, 3)

# dirichlet
h = Constant(0) #Expression("100*cos(x[0]*x[1])", degree=d) #Constant(0)
# newmann
g = Constant(0) #Expression("cos(x[0])", degree=1) #Constant(0)

boundary_conditions = {4: {'Robin': h},
1: {'Robin': h},
2: {'Robin': h},
3: {'Neumann': g}}

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

p = Expression("0.015+0.005*cos(3.5+x[1])", degree=d) #Constant(0.01) # #
I = 1.0
DV = hx*hy
M = I/DV
alpha = Expression("x[0]/(x[0]*x[0]+x[1]*x[1])", degree=d)
beta = Constant(1.0)

# source
f =  Constant(0) #Expression("cos(x[0]+x[1])", degree=d) #Constant(0) #source term

integrals_R_L = []
integrals_R_a = []
for i in boundary_conditions:
if 'Robin' in boundary_conditions[i]:
if boundary_conditions[i]['Robin'] != 0:
hi = boundary_conditions[i]['Robin']
integrals_R_L.append((p*hi*v/beta)*ds(i))
integrals_R_a.append((alpha*p*u*v/beta)*ds(i))

L = f*v*dx + p*g*v*ds(3) + sum(integrals_R_L)

# Compute solution
A, b = assemble_system(a, L)

delta1 = PointSource(V, Point(in_x, in_y), M)
delta1.apply(b)

# delta2 = PointSource(V, Point(out_x, out_y), -M)
# delta2.apply(b)

u = Function(V)
solve(A, u.vector(), b)

#plot y=k
y = 0.0
meshx = IntervalMesh(nx, x0, xN)
W = FunctionSpace(meshx,"Lagrange",d)
w = Function(W)
coor = meshx.coordinates()
w_array = w.vector()#.array()

if meshx.num_vertices() == len(w_array):
for i in range(meshx.num_vertices()):
w_array[i] = u(coor[i][0],y)

w.vector().set_local(w_array)

# plots
plt.subplot(2,3,1)
plot(f, mesh=mesh, title = "f: source term")
# plt.savefig('vc_poisson/f.png')

plt.subplot(2,3,2)
plot(h, mesh=mesh, title = "h: Dirichlet boundary condition")

plt.subplot(2,3,3)
plot(g, mesh=mesh, title = "g: Newmman boundary condition")
# plt.savefig('vc_poisson/g.png')

plt.subplot(2,3,4)
plot(p, mesh=mesh, title = "\sigma: resistivity distribution")
# plt.savefig('vc_poisson/g.png')

plt.subplot(2,3,5)
plot(u, title = "numerical solution")
# plt.savefig('vc_poisson/sol.png')

# plt.subplot(2,3,5)
# plot(mesh)
# plot(u, interactive=True, title = "numerical solution + mesh")

plt.subplot(2,3,6)
plot(w, title="numerical solution u(x,0)")
plt.savefig('vc_poisson/solution.png')

plt.figure()
plot(u.root_node(), mode='warp', title="numerical solution surface")
plt.savefig('vc_poisson/solution_surface.png')

# Save solution to file in VTK format
vtkfile = File('vc_poisson/solution.pvd')
vtkfile << u

plt.show()  ​

These are my operators:

$a(u,v)=\int_\Omega\sigma\nabla u\cdot\nabla v+\int_{\Gamma_R}\frac{\alpha\sigma}{\beta}uv$
$L(v)=\int_\Omega fv+\int_{\Gamma_N}\sigma\,g\,v+\int_{\Gamma_R}\frac{\sigma}{\beta}hv$

If I comment the boundary integral on $a$, no error appears but the solution isn't what i'm looking for. I must add I'm 2 weeks old in fenics

Community: FEniCS Project