Solving a 1d linear problem


288
views
0
7 months ago by

Good evening every one

I try to volve the very simple problem v"=0, v(0)=0, v(1) = 1 on [0, 1]

I write the following fenics code

from fenics import *
#
# 1d variationnal linear problem
#
# Create mesh and define fuction space
#
mesh = IntervalMesh(20, 0.0, 1.0)

V = FunctionSpace(mesh, "CG", 1)
#
# Define boundary condition
#
tol = 1E-14

def x0_boundary(x, on_boundary):
return on_boundary and x[0] < tol

def x1_boundary(x, on_boundary):
return on_boundary and abs(x[0]-1.0) < tol

bc0 = DirichletBC(V, Constant(0), x0_boundary)
bc1 = DirichletBC(V, Constant(1), x1_boundary)

bc = [bc0, bc1]
#
# Define variational problem
#
u = TrialFunction(V)
v = TestFunction(V)
a = -inner(grad(u)*grad(v))*dx
T = Constant(0)
L = T*v*dx
#
# Compute solution
#
u = Function(V)
solve(a == L, u, bc)

and I obtain the following error message
Invalid ranks 1 and 1 in product.

Traceback (most recent call last):

  File "drying.py", line 30, in <module>

    a = -inner(grad(u)*grad(v))*dx

  File "/usr/local/lib/python2.7/dist-packages/ufl/exproperators.py", line 193, in _mul

    return _mult(self, o)

  File "/usr/local/lib/python2.7/dist-packages/ufl/exproperators.py", line 167, in _mult

    error("Invalid ranks {0} and {1} in product.".format(r1, r2))

  File "/usr/local/lib/python2.7/dist-packages/ufl/log.py", line 172, in error

    raise self._exception_type(self._format_raw(*message))

ufl.log.UFLException: Invalid ranks 1 and 1 in product.

Aborted

I really do not understand where I am wrong.
Second querstion: How can I replace the Dirichlet boundary condition v(1)=0 by a Neuman condition (v'(1)=1 for example).

Thank you very much for help.

Regards

Xavier

Community: FEniCS Project

2 Answers


3
7 months ago by
The error occurs because you are trying to multiply two vectors together by writing grad(u)*grad(v). You should change that corresponding line to a = -inner(grad(u), grad(v))*dx.

Thank you very much

I do appreciate your taking time to answer

written 7 months ago by xavier chateau  
2
7 months ago by
And, to answer your second question, to put the Neumann condition v'(1)=2 on the right end point (I put a value of 2, to avoid confusion), you would leave out the DirichletBC at that point, and you would subtract the "integral" of  the Neumann boundary data times the test function v over the Neumann boundary to the functional L.  (Since the Neumann boundary is just the single point x=1, the thing you need to subtract from L is Constant(2.)*v*ds.  The ds means "integration" over both boundary points, but it doesn't matter since the test functions v will all vanish at the left end point.

The entire code (to solve u'' = 0, u(0) = 0, u'(1) = 2), is thus:
from fenics import *
mesh = IntervalMesh(20, 0.0, 1.0)
V = FunctionSpace(mesh, "CG", 1)
def x0_boundary(x, on_boundary):
    return on_boundary and near(x[0], 0.)
bc = DirichletBC(V, Constant(0.), x0_boundary)
u = TrialFunction(V)
v = TestFunction(V)
a = -inner(grad(u), grad(v))*dx
T = Constant(0.) 
L = T*v*dx - Constant(2.)*v*ds
u = Function(V)
solve(a == L, u, bc)
​
Thank you very much

I do appreciate your taking time to answer
written 7 months ago by xavier chateau  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »