Mixed boundary conditions

Also, I was suspecting that the mixed formulation may be needed (solve q, and u) instead of removing q for k grad u. Would this help the sovler in that q's bc would be a function of u, instead of u's bc being a function of u? It would seem this would help the solver?

Such fluxes, Robin condition, can be applied just as you think they would ;)

For a linear diffusion equation written on nonlinear form (because it is
easier ;) )

  ...
  v, u = TestFunction(V), Function(V)
  u0 = Constant(0)
  F = inner(grad(u), grad(v))*dx + (u-u0)*ds
  J = derivative(F, u)
  bcs = ...

  VariationalProblem(F, J, bcs)

Note that the syntax for VariationProblem has changed in the development
branch.

Johan

On Friday June 24 2011 06:26:05 Charles Cook wrote:
> New question #162603 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/162603
>
>
> I'm new to the fenics project and am going through the tutorials and test
> cases (in C++). I'm trying to run a simple test on a unit-square with
> three boundary conditions (Essential/Direchlet, Natural, and mixed) for
> Newton's law of cooling.
>
> I'm having trouble with the mixed boundary condition, where q dot n = H (u
> - u_inf), specifically the normal flux is a function of u. How may I
> apply a normal flux that is a function of the solution?
>
> I was also looking through the NS examples hoping to find a far-field
> boundary condition that might address this.
>
> On a side note, is the fenics-book the go-to source for the project?
>
> Thanks!

Wow, maybe its that I'm new... But dx is applied within the domain and ds has the meaning at the boundary? If I may ask another question, how would I constrain (u-u0)*ds to just one boundary? My guess is that a new member for the operator (F) will be added to which I can assign a domain?

I've been working my way up in complexity. I started with the Poisson example, made it 1-d, and removed the source term to get a very simple 1d diffusion problem. When adding a linear convection term I was amazed at just having to add

... +ac*inner(grad(u), v)*dx

By the way, I am using the PPA packages in Ubuntu (I'm new to Ubuntu). Would it be better to move to the development branch? Do I obtain them from launchpad with some version control system? (cvs/svn/mercurial/git?)

Thank you very much!

On Friday June 24 2011 08:51:07 Charles Cook wrote:
> Question #162603 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/162603
>
> Charles Cook posted a new comment:
> Wow, maybe its that I'm new... But dx is applied within the domain and
> ds has the meaning at the boundary?

Yes!

> If I may ask another question, how
> would I constrain (u-u0)*ds to just one boundary?

New syntax:

  markers = SomeFacetFunction()
  dss = ds[markers]
  F = ... + (u-u0)*dss(3)
  VariationalProblem(F, J, bcs)

Stable syntax:

  F = ... + (u-u0)*dss(3)
  VariationalProblem(F, J, bcs, exterior_facet_domains=markers)

> My guess is that a
> new member for the operator (F) will be added to which I can assign a
> domain?

For the new C++ interface yes. But in the Python interface you can just follow
the above syntax. If you wonder how to create a particular FacetFunction with
markers at different boundaries. Have a look in:

  demo/undocumented/subdomains

> I've been working my way up in complexity. I started with the Poisson
> example, made it 1-d, and removed the source term to get a very simple 1d
> diffusion problem. When adding a linear convection term I was amazed at
> just having to add
>
> ... +ac*inner(grad(u), v)*dx

:)

> By the way, I am using the PPA packages in Ubuntu (I'm new to Ubuntu).
> Would it be better to move to the development branch? Do I obtain them
> from launchpad with some version control system? (cvs/svn/mercurial/git?)

Instal bazaar and:

  bzr branch lp:foo

where foo is your fenics project you want to grab code for.

You might want to have a look at dorsal:

   http://bit.ly/mddNER

Johan

> Thank you very much!

I tried installing with dorsal but had no luck with ubuntu 11.04 (64). It appears 11.04 is not supported as of yet?

Thank for the feedback, I am still going through a lot of it.

Do I need a nonlinear version of ffc? I am having some problems adding the nonlinear term

To add a little more detail, I'm getting errors about terms not having the same rank. This makes me assume I have my problem setup wrong (tensor rank disagreement).

I can say nothing if there are no minimal code attached to the problem

Johan

On Friday June 24 2011 11:31:12 Charles Cook wrote:
> Question #162603 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/162603
>
> Charles Cook posted a new comment:
> To add a little more detail, I'm getting errors about terms not having
> the same rank. This makes me assume I have my problem setup wrong
> (tensor rank disagreement).

Of course. I am working on my UFL file now:

a = inner(grad(u), grad(v))*dx + (u-u0)*v*ds

from the partial integration i expect the test function, v, to be there. I have substituted the derivative for my BC of u-u0

Full listing :)

element = FiniteElement("Lagrange", triangle, 1)

u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)
g = Coefficient(element)
u0 = Constant(element)

a = inner(grad(u), grad(v))*dx + inner((u-u0),v)*ds
L = f*v*dx + g*v*ds

My apologies for so many posts. The above includes a change I was trying (inner product on second to last line)

You need to shift the linear (in v) term in a to the RHS (into L),

    a = inner(grad(u), grad(v))*dx +u*v*ds
    L = f*v*dx + g*v*ds + u0*v*ds

Thank you! As it always seems to happen I found the part of the fenics-book [2011-06-03] that discusses this as you answered the question.

For others, section 2.6.3 discusses this in detail.

Thanks Garth Wells, that solved my question.

On Friday June 24 2011 11:50:56 Charles Cook wrote:
> Question #162603 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/162603
>
> Charles Cook posted a new comment:
> Full listing :)
>
> element = FiniteElement("Lagrange", triangle, 1)
>
> u = TrialFunction(element)
> v = TestFunction(element)
> f = Coefficient(element)
> g = Coefficient(element)
> u0 = Constant(element)
>
> a = inner(grad(u), grad(v))*dx + inner((u-u0),v)*ds
> L = f*v*dx + g*v*ds

So you are correct wrt rank being mixed here. As you have a rank 1 expression
for in the bilinear form.

In the way you set up the problem you need to move split the mixed boundary
part:

   a = inner(grad(u), grad(v))*dx + u*v*ds
   L = f*v*dx + g*v*ds + u0*v*ds

You can also use lhs and rhs functions:

  form = inner(grad(u), grad(v))*dx + inner((u-u0),v)*ds - f*v*dx - g*v*ds
  a = lhs(form)
  L = rhs(form)

That is why I described the original problem as a nonlinear one. It is easier
and clearer. It uses a coefficient Function for u instead of a TrialFunction.

Johan

Thank you Johan for the additional clarification, much appreciated as always!

If I may tag on one last question about marking domains... I've been going through the text/demos and the python translation for the robin chapter is not clear to me. I feel like I am almost there, just missing some concept

I believe I have assembled my boundary parts, and DirichletBC, but what type of boundary condition do I use for the natural and robin BC?

    MeshFunction<unsigned int> boundary_parts(mesh, mesh.topology().dim() - 1);

    // Define dirichlet boundary condition
    Constant ug(0.0);
    Gamma_g gamma_g;
    gamma_g.mark(boundary_parts,2);
    DirichletBC bc_g(V, ug, gamma_g);

    Gamma_h gamma_h;
    gamma_h.mark(boundary_parts, 1);

    Gamma_H gamma_H;
    gamma_H.mark(boundary_parts, 0);

 std::vector<const BoundaryCondition*> bcs;
 bcs.push_back(&bc_g);
 bcs.push_back(&gamma_h); // this of course does not work
 bcs.push_back(&gamma_H);

    // Setup the problem
    VariationalProblem problem(a, L, bcs);

Found through an undocumented demo.

    a.exterior_facet_domains = boundary_parts;

Thanks guys!

Thanks again guys, I matched my analytical solution with all three boundary types :)