velocity field over subdomains in 1 plot

On Thu, Oct 11, 2012 at 11:26:04AM -0000, Sander Rhebergen wrote:
> New question #210928 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/210928
>
> Dear,
>
> I am trying to solve an advection diffusion equation in which the given velocity field is defined over 3 subdomains. I have my own mesh files:
>
> mesh = Mesh("domain.xml")
> subdomains = MeshFunction("uint", mesh, "domain_physical_region.xml")
> boundaries = MeshFunction("uint", mesh, "domain_facet_region.xml")
>
> and on each subdomain I have a different velocity, for example:
>
> class uone(Expression):
> def eval(self, values, x):
> xx=x[0]
> yy=x[1]
> cc=numpy.sin(xx)
> cd=numpy.cos(yy)
> ux=costheta*cc + sintheta*cd
> uy=sintheta*cc - costheta*cd
> values[0]=ux
> values[1]=uy
> def value_shape(self):
> return (2,)
>
> u0=Expression(("1","-1"))
> u1=uone()
> u2=Expression(("0","0"))
>
> With these different velocities I can solve an advection diffusion equation for T:
>
> dx=Measure("dx")[subdomains]
> ds=Measure("ds")[boundaries]
> T = TrialFunction(TT)
> w = TestFunction(TT)
> f = Constant(0.0)
> at = (u0[i]*T.dx(i)*w + w.dx(i)*T.dx(i))*dx(0) + \
> (u1[i]*T.dx(i)*w + w.dx(i)*T.dx(i))*dx(1) + \
> (u2[i]*T.dx(i)*w + w.dx(i)*T.dx(i))*dx(2)
> Lt = f*w*dx(0) + f*w*dx(1) + f*w*dx(2)
>
> A, bb = assemble_system(at, Lt, bcs)
> T = Function(TT)
> solve(A, T.vector(), bb)
>
> I have two questions about this code. First of all, as also mentioned in the FEniCS book, defining class uone as I did above is slow. Is there a nice way to speed this up? The second question is, is there a way to plot the velocity field on the whole domain? I have been able to plot the solution on subdomains, but then I get 3 different plots:
>
> sub_mesh0 = SubMesh(mesh,subdomains,0)
> sub_mesh1 = SubMesh(mesh,subdomains,1)
> sub_mesh2 = SubMesh(mesh,subdomains,2)
> sV0 = VectorFunctionSpace(sub_mesh0,"CG",2)
> sV1 = VectorFunctionSpace(sub_mesh1,"CG",2)
> sV2 = VectorFunctionSpace(sub_mesh2,"CG",2)
> Iv0=interpolate(u0,sV0)
> Iv1=interpolate(u1,sV1)
> Iv2=interpolate(u2,sV2)
> plot(Iv0)
> plot(Iv1)
> plot(Iv2)
>
> I was thinking, maybe there is a way to interpolate u0, u1, u2 to one vector U on the whole domain? The velocity, however, may be discontinuous across the subdomain boundaries, so maybe a DG space? Any help is appreciated.

1. If you want to speed things up, try help(Expression) and see how to
implement your expression as a C++ code snippet from Python.

2. Use the eval_cell version of the callback for your Expression (see
again the docstring for Expression) and use that to select between the
different domains based on the cell number and indexing into your
subdomains MeshFunction.

--
Anders

Thanks! Part 1 worked. Part 2 I think I have working as well, but it is very slow. I'll have a go at writing that as well as a C++ snippet.

Thanks

On Fri, Oct 12, 2012 at 2:16 PM, Sander Rhebergen
<email address hidden> wrote:
> Question #210928 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/210928
>
> Sander Rhebergen posted a new comment:
> Thanks! Part 1 worked. Part 2 I think I have working as well, but it is
> very slow. I'll have a go at writing that as well as a C++ snippet.
>

Take a look at

    demo/undocumented/tensor-weighted-poisson/python/demo_tensor-weighted-poisson.py

for an example of (2).

Garth

> --
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.

Thanks!