Assign mixed element function values on subdomain
fenics.FunctionSpacemade 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?
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 = \ fenics.inner, fenics.dot, fenics.grad, fenics.div, fenics.sym F = (dot(v, dot(grad(u), u)) - div(v)*p \ + 2.*mu*inner(sym(grad(v)), sym(grad(u))) - 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.) class FixedWalls(fenics.SubDomain): def inside(self, x, on_boundary): return on_boundary and (near(x,0.) | near(x,1.) | near(x,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 > 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.
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:
To avoid hacking; perhaps the better answer lies within the
Edit: Removed comment about this not working with adapted meshes. One just has to call
leaf_node() in all of the right places.