Dear all,
I am trying to solve a linear elasticity problem on domains with different materials in 2D. The whole domain is a half disk divided into a smaller half disk and its complement. The model has both Dirichlet and Neumann boundary conditions. The Dirichlet bc corresponds to "Bottom" and Newmann "Top". Here is my code:
from __future__ import print_function
from dolfin import *
from mshr import *
import sys, math, numpy
# Create classes for defining parts of the boundaries and the interior
# of the domain
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0.0)
class Top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (between(x[1], (0.5, 1.1)))
class Obstacle(SubDomain):
def inside(self, x, on_boundary):
return (x[0]*x[0] + x[1]*x[1] <= 0.25)
# Define parameters
r = 1.0
MeshDensity = 20
# Initialize sub-domain instances
top = Top()
bottom = Bottom()
obstacle = Obstacle()
# Define the domain and generate the mesh
domain = Circle(Point(0, 0), r) - Rectangle(Point(-1.0, -2.0), Point(1.0, 0.0))
mesh = generate_mesh(domain,MeshDensity,"cgal")
#plot(mesh, 'Mesh')
# Initialize mesh function for interior domains
domains = CellFunction("size_t", mesh)
domains.set_all(0)
obstacle.mark(domains, 1)
#print(domains.array())
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
bottom.mark(boundaries, 0)
top.mark(boundaries, 1)
#print(boundaries.array())
# Define input data
a1 = Constant(1.0)
a0 = Constant(0.01)
g = Constant((0.0, -1.0))
f = Constant((0.0, 1.0))
# Define function space and basis functions
V = VectorFunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
# Define Dirichlet boundary ( x=0, y=0 )
c = Expression(("0.0", "0.0"))
# Define Dirichlet boundary conditions at top and bottom boundaries
bcs = DirichletBC(V, c, boundaries, 0)
# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
# Define variational form
F = (dot(a0*grad(u), grad(v))*dx(0) + (dot(a1*grad(u), grad(v)))*dx(1)
- inner(g, v)*ds(1)
- inner(f, v)*dx(0) - inner(f, v)*dx(1))
# Separate left and right hand sides of equation
a, L = lhs(F), rhs(F)
# Solve problem
u_h = Function(V)
solve(a == L, u_h.vector(), bcs)
# Plot solution
plot(u, title="u")
interactive()
The error message goes as follows:
Can only integrate scalar expressions. The integrand is a tensor expression with value rank 2 and free indices ().
Traceback (most recent call last):
File "/home/jason/Programs/Python/Questions/08-Feb-2016.py", line 72, in
- inner(f, v)dx(0) - inner(f, v)dx(1))
File "/usr/lib/python2.7/dist-packages/ufl/measure.py", line 374, in rmul
error(msg % (integrand.rank(), integrand.ufl_free_indices))
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 151, in error
raise self._exception_type(self._format_raw(*message))
UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value rank 2 and free indices ().
The way I see it, after taking the dot product (inner product), all of my integrands are scalars. Is it because the way I define the expression or the function space is wrong, maybe both? It would be really appreciated if you could tell me how to modify it to make it work.
