### Assign mixed element function values on subdomain

157
views
0
10 weeks ago by
Short question: Given a fenics.Function from a fenics.FunctionSpace made from a fenics.MixedElement, how does one assign a subset of the function values restricted to a subdomain, i.e. without changing values outside of the subdomain?

Details:

The following solves the lid-driven cavity with mixed finite elements (tested with FEniCS 2017.2.0):

import fenics

mesh = fenics.UnitSquareMesh(20, 20)

P2P1 = fenics.MixedElement([
fenics.VectorElement('P', mesh.ufl_cell(), 2),
fenics.FiniteElement('P', mesh.ufl_cell(), 1)])

W = fenics.FunctionSpace(mesh, P2P1)

solution = fenics.Function(W)

u, p = fenics.split(solution)

v, q = fenics.TestFunctions(W)

mu = fenics.Constant(0.01)

inner, dot, grad, div, sym = \

F = (dot(v, dot(grad(u), u)) - div(v)*p \
- q*div(u))*fenics.dx

near = fenics.near

class Lid(fenics.SubDomain):

def inside(self, x, on_boundary):

return on_boundary and near(x[1], 1.)

class FixedWalls(fenics.SubDomain):

def inside(self, x, on_boundary):

return on_boundary and (near(x[0],0.) | near(x[0],1.) | near(x[1],0.))

bcs = [
fenics.DirichletBC(W.sub(0), (1., 0.), Lid()),
fenics.DirichletBC(W.sub(0), (0., 0.), FixedWalls())]

fenics.solve(F == 0., solution, bcs)​

Here's the resulting velocity vector field for reference:
import matplotlib
%matplotlib inline

u, p = fenics.split(solution)

fenics.plot(u)​

Now I want to use most of this solution as the initial guess for a different boundary value problem. For example, I could change the lid velocity.

Assume that I need to assign the velocity local to some subdomain by the lid, while preserving the other guess values from our previous solution. I could for example define the subdomain with

class TopLayer(fenics.SubDomain):

def inside(self, x, on_boundary):

return x[1] > 0.9​

How do I change the values of the solution only in this subdomain? For now I only need to assign a constant value, and it's okay to create a new function (from the same function space) rather than assigning the values in place.

I remember seeing some related posts that play with degree of freedom maps for simpler cases; but I haven't seen a mixed finite element example, and asking the user to figure out the mapping for this case seems a bit complicated. I sure hope I'm missing something more straightforward.
Community: FEniCS Project

3
10 weeks ago by
@klunkean came up with a hack that solves this problem in the discussion here: https://www.allanswered.com/post/mqpbr/#bamzr

So continuing from the script in my question, we can set the velocity and pressure in the subdomain with

hack = fenics.DirichletBC(W.sub(0), fenics.Constant((0., 0.)), TopLayer())

hack.apply(solution.vector())

hack = fenics.DirichletBC(W.sub(1), fenics.Constant(0.), TopLayer())

hack.apply(solution.vector())​

which appears to give the right result:
fenics.plot(u)​

fenics.plot(p)

fenics.plot(mesh)

To avoid hacking; perhaps the better answer lies within the fenics.DirichletBC.apply code.

Edit: Removed comment about this not working with adapted meshes. One just has to call leaf_node() in all of the right places.

At least when actually instantiating a fenics.NonlinearVariationalProblem, rather than calling fenics.solve as done above, this seems to somehow corrupt the boundary conditions of the problem. This can be avoided by creating new temporary copies of the function space and solution on which to call DirichletBC.apply, e.g.

function_space_copy = fenics.FunctionSpace(mesh, P2P1)

new_solution = fenics.Function(function_space_copy)

new_solution.vector()[:] = sim.solution.vector()

hack = fenics.DirichletBC(
function_space_copy.sub(0),
(0., 0.),
TopLayer())

hack.apply(new_solution.vector())

solution.vector()[:] = new_solution.vector()​
written 10 weeks ago by Alexander G. Zimmerman
The BCs are homogenised if you solve a nonlinear problem.
written 10 weeks ago by Nate
Perhaps I shouldn't have added my comment about fenics.NonlinearVariationalProblem without showing exactly what I meant by "corrupt". The same issue should arise with linear problems. It's not the values of the BC's that get "corrupted", but rather the data structure gets ruined somehow and an error is thrown about the BC's being in the wrong space.
written 10 weeks ago by Alexander G. Zimmerman