Projection of a 2D function on parts of a 3D domain (extrusion in one direction)

Question

Dear All,

I have the following problem: I have a 3D domain Click here,left, in which I want to solve an advection-diffusion problem. The bottom bit of the domain is supposed to have a fluid velocity of zero (0,0,0), whilst the top part (the part extending the bottom in z-direction) has a flow field (0,0,w(x,y)). Therefore, the flow field is constant in time and z-direction.

I know w(x,y) from a previous finite element simulation in 2D, which is shown here here, right.. I know in particular that at the narrow end of the domain x=0 w(x=0,y)=0 due to the boundary condition I specified then.

The question: How can I project it on the new domain?, so that the flow field in the top third of the domain becomes (0,0,w(x,y)) for all z?

What I have tried:
I defined a function u(x,y,z) = (0,0,w(x,y)) if I am in the top region of the domain, other wise it evaluates to zero. The projection of this function on the new mesh has two problems:
(1). It is not constant in z-direction, and, (2) and the projected flow field is not (0,0,0) for x=0 as I would expect since w(x=1,y)=0.

I also tried to adapt this idea, which only led to much worse results.

I would be grateful for any ideas, and please be gentle since I am at the beginning of my FEM journey.

Please let me know if you need further information.

Best wishes,
Merlin :)

Answer 1

Hi, based on how I understand your problem I suggest you consider the following

from dolfin import *

mesh2d = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh2d, 'CG', 1)
f = interpolate(Expression('sin(pi*(x[0]+x[1]))', degree=4), V)

class Foo(Expression):
    def __init__(self, f): self.f = f 

    def value_shape(self): return (3, )

    def eval(self, values, x):
        values[0] = 0.0
        values[1] = 0.0
        values[2] = 0.0 if x[2] < 0.5 else self.f(x[:2])

mesh3d = UnitCubeMesh(10, 10, 10)
V = VectorFunctionSpace(mesh3d, 'CG', 1)
v = interpolate(Foo(f), V)
plot(v)
interactive()