How to evaluate an expression only on a subdomain?


293
views
0
6 months ago by
Marcin  
Hello all,
What is the proper syntax for evaluating an expression on a subdomain of a mesh? To be more precise, I have an output of a simulation (E) that I would like to process further by evaluating f(E), where the form of f(E) is varies between subdomains.

There is no integration involved so multiplying by dx(i) will not work.


from dolfin import *

# mesh with two domains dx(0) and dx(1)

class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (0.5, 0.7)) and between(x[0], (0.2, 1.0)))

obstacle = Obstacle()
# Define mesh
mesh = UnitSquareMesh(64, 64)
# Initialize mesh function for interior domains
domains = CellFunction("size_t", mesh)
domains.set_all(0)
obstacle.mark(domains, 1)

# in my case this will be an output of some simulation
V = VectorElement("Lagrange", mesh.ufl_cell(), 2, dim = 2)
em_space  = FunctionSpace(mesh, V)
# define some field in the em space
E_analytical = Expression(('cos(pi*x[0])', 'sin(pi*x[1])'),
                          element = el_space.ufl_element())

E = interpolate(E_analytical, em_space)

# quesion? how to define the following
# some f(E) that depends on domains
a = 1.0
b = 2.0
f = grad(c1*E) # on dx(0)
f = grad(c2*E) # on dx(1)
# I do not want to define a variable coefficent c on the whole mesh,
# because c*E will be discontinous across boundaries and thus grad(c*E) will
# have infinities

# and plot
plot(f, interactive = True)
​
Community: FEniCS Project

1 Answer


3
6 months ago by
Hi, consider this:

from dolfin import *


mesh = UnitSquareMesh(64, 64)

class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (0.5, 0.7)) and between(x[0], (0.2, 1.0)))

obstacle = Obstacle()

# Initialize mesh function for interior domains
domains = CellFunction("size_t", mesh)
domains.set_all(0)
obstacle.mark(domains, 1)

V = VectorElement("Lagrange", mesh.ufl_cell(), 2, dim = 2)

E_analytical = Expression(('cos(pi*x[0])', 'sin(pi*x[1])'), degree=2)

E1 = SubMesh(mesh, domains, 0)
E2 = SubMesh(mesh, domains, 1)

E1space = FunctionSpace(E1, V)
E2space = FunctionSpace(E2, V)

E_1 = project(E_analytical, E1space)
E_2 = project(E_analytical, E2space)
plot(E_1, interactive=True)
plot(E_2, interactive=True)

dx = Measure("dx")(subdomain_data = domains)



Thank you, this will help.
written 6 months ago by Marcin  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »