Problems with more than one value for g

variational formulation
1

Hello everybody,

i have a problem with my specific values for g (Neumann condition).

from dolfin import *
from mshr import *
import matplotlib.pylab as plt
import math

# Gebiet und Mesh erstellen------------------------------------------------
domain_vertices = [Point(2.0, 2.0),
                  Point(1.0, 3.0),
                  Point(-1.0, 1.0),
                  Point(4.0, -4.0),
                  Point(5.0, -3.0),
                  Point(2.0, 0.0),
                  Point(0.001, 0.0)]


domain = Polygon(domain_vertices)
mesh = generate_mesh(domain,10)
coordinates = mesh.coordinates()


#Functionspace-----------------------------------------------------------
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)  
v = TestFunction(V)
uh = Function(V)

x = SpatialCoordinate(mesh)
f = sin(2.0*pi*x[0])*cos(2.0*pi*x[1])
n = FacetNormal(mesh)

# Tangential gradient operator:
def grad_t(u):
    return grad(u) - dot(grad(u),n)*n

# Randbedingungen-----------------------------------------------------------
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (between(x[1], (-1.0, 1.0)) and between(x[0], (4.0, -4.0)))

class Right(SubDomain):
    def inside(self, x, on_boundary):
        #return near(x[1], 1.0)
        return on_boundary and ((between(x[1], (2.0, 2.0)) and between(x[0], (1.0, 3.0))) and (between(x[1], (0.1, 0.0)) and between(x[0], (2.0, 2.0))) and (between(x[1], (2.0, 0.0)) and between(x[0], (0.1, 0.0))) and (between(x[1], (5.0, -3.0)) and between(x[0], (2.0, 0.0))))

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (between(x[1], (5.0, -3.0)) and between(x[0], (4.0, -4.0)))

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (between(x[1], (1.0, 3.0)) and between(x[0], (-1.0, 1.0)))

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundaries.set_all(4)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)


# Define Dirichlet boundary condition at top
bc2 = DirichletBC(V, 0.0, boundaries, 2)

#Take a look at the solution:
ds=Measure('ds', domain=mesh, subdomain_data=boundaries)
area_gamma1 = assemble(1*ds(4))
g = Expression("1/area_gamma1 * x4")
leitfaehigkeit=8.6e6
# Residual of Poisson problem with the specified BC:----------------------------------------------
res2 = inner(grad(u),grad(v))*dx - leitfaehigkeit*inner(g,grad_t(v))*ds(4) - inner(f,v)*dx
solve(lhs(res2)==rhs(res2),uh, bc2)

from matplotlib import pyplot as plt
plot(uh)
plt.show()
plt.savefig("solution2.png")

With x4 are my values for the current at several time steps.

I get the following error

Blockquote
g = Expression(“1/area_gamma1 * x4”)
File “/usr/local/lib/python3.6/dist-packages/dolfin/function/expression.py”, line 400, in init
self._cpp_object = jit.compile_expression(cpp_code, params)
File “/usr/local/lib/python3.6/dist-packages/dolfin/function/jit.py”, line 158, in compile_expression
expression = compile_class(cpp_data, mpi_comm=mpi_comm)
File “/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py”, line 170, in compile_class
raise RuntimeError(“Unable to compile C++ code with dijitso”)
RuntimeError: Unable to compile C++ code with dijitso
Blockquote

Could somebody please help me with this? Currently i am a bit confused.

Greetings

2

Please make this a minimal example, and describe what your issue with having multiple Neumann boundaries gives rise to, to make it easier for people to help you.

3

Sure. Thank you for the hint. I hope it is fine now.

4

There are several issues with how you specify g.
First of all, g has to be a vector expression, as you are doing the inner product with grad_t(v).
Secondly, you need to define the variables used in the string in the Expression, and supply a degree to specify which function space it should be interpolated into.
Here is a line that works:

g = Expression(("1/area_gamma1 * x4","2"), x4=1, area_gamma1=area_gamma1, degree=2)

However, in the code above, it is not clear what the second component should be (so I set it to 2), and x4 is not defined, so I set it to one.

5

Thank you very much. I was so stupid. Now I understand using g and I was able to correct the form.