How to make custom tensor expression on subdomains?


145
views
0
4 months ago by
Dario  
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  

1 Answer


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!

Please login to add an answer/comment or follow this question.

Similar posts:
Search »