Robin + Neumann boundary conditions with marker

7 days ago by

I'm trying to adapt what is described here:

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/ 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                                      

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


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)

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']

a = p*inner(nabla_grad(u), nabla_grad(v))*dx + sum(integrals_R_a)
L = f*v*dx + p*g*v*ds(3) + sum(integrals_R_L)

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

# Add delta
delta1 = PointSource(V, Point(in_x, in_y), M)

# 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)


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

plot(h, mesh=mesh, title = "h: Dirichlet boundary condition")

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

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

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")

plot(w, title="numerical solution u(x,0)")

plot(u.root_node(), mode='warp', title="numerical solution surface")

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

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
Please login to add an answer/comment or follow this question.

Similar posts:
Search »