How to make custom tensor expression on subdomains?
4 months ago by
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 = 1.0 values = 0.0 values = 0.0 values = 1.0 else: values = 0.8 values = 0.3 values = 0.1 values = 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
File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 330, in _solve_varproblem
File "helmholtz2.py", line 31, in eval_cell
values = 0.0
IndexError: index 1 is out of bounds for axis 0 with size 1
Community: FEniCS Project
4 months ago by
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.