### How to make custom tensor expression on subdomains?

145
views
0
4 months ago by
I'm trying to make custom tensor expression defined on subdomains.

from fenics import *
from mshr import *

# Domain definition and meshing
domain = Rectangle(Point(-2, -2), Point(2, 2))
patch = Rectangle(Point(0, 0), Point(1, 1))

domain.set_subdomain(1, patch)
mesh = generate_mesh(domain, n_mesh)
V = FunctionSpace(mesh, 'P', 1)

# Subdomain parameters
markers = MeshFunction('size_t', mesh, 2, mesh.domains())
class koef(Expression):
def __init__(self, mesh, **kwargs):
self.markers = markers
def eval_cell(self, values, x, cell):
if markers[cell.index] == 0:
values[0] = 1.0
values[1] = 0.0
values[2] = 0.0
values[3] = 1.0
else:
values[0] = 0.8
values[1] = 0.3
values[2] = 0.1
values[3] = 0.6
def value_shape(self):
return (2,2)

A = koef(mesh, degree=0)​

# Boundary condition
def boundary(x, on_boundary):
return on_boundary

bc = DirichletBC(V, u_0, boundary)

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

# Source
f = Constant(1)
a = (-inner(A*grad(u), grad(v)) + dot(u, v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

Error i get is:

Traceback (most recent call last):
File "helmholtz2.py", line 65, in <module>
solve(a == L, u, bc)
File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 300, in solve
_solve_varproblem(*args, **kwargs)
File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 330, in _solve_varproblem
solver.solve()
File "helmholtz2.py", line 31, in eval_cell
values[1] = 0.0
IndexError: index 1 is out of bounds for axis 0 with size 1
Community: FEniCS Project
1
I haven't tested it, but I guess it just could be the indentation of the lines
def value_shape(self):
return (2,2)​

which in your code is inside your eval_cell method.

written 4 months ago by klunkean

0
4 months ago by
I'm not sure why it isn't reading the shape, but you can also give an element for the construction of the element which automatically sets its shape:

A = koef(mesh, element = TensorElement('CG', triangle, 3))​

This runs in my computer, so please confirm as well in yours. I didn't check the solution, but perhaps it makes sense. Let me know. Best regards!