DG with Mixed FE question

It means you must write n('+') or n('-') to pick from which side the facet
normal is seen from in the dS integral.

Martin
Den 16. jan. 2013 22:31 skrev "Yaakoub El Khamra" <
<email address hidden>> følgende:

> New question #219369 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/219369
>
>
> I am learning DG through a book and thought I would try my hand at coding
> up some of the examples. The problem I am solving is heat transfer with a
> mixed formulation: (q = \nabla T and \nabla \cdot q = 0 ) and the internal
> penalty method with a coefficient of 1.0, i.e.
>
> T_hat = avg(T)
> and
> q_hat = avg(grad(T)) - jump(T)
>
> The error I am getting suggests that the facet normals must be restricted:
> raise self._exception_type(self._format_raw(*message))
> ufl.log.UFLException: FacetNormal must be restricted.
>
> Any help is greatly appreciated.
>
> Regards
> Yaakoub
>
>
> from dolfin import *
>
> # Create mesh
> mesh = UnitSquareMesh(32, 32)
>
> n = FacetNormal(mesh)
> h = CellSize(mesh)
> h_avg = (h('+') + h('-'))/2
>
> qSpace = VectorFunctionSpace(mesh, "DG",2)
> TSpace = FunctionSpace(mesh,"DG",1)
>
> W = qSpace*TSpace
>
> (q, T) = TrialFunctions(W)
> (v_q,v_T) = TestFunctions(W)
>
> zero = Constant(0.0)
>
> # for internal (dS) faces:
> # T_hat = avg(T)
> # q_hat = avg(grad(T) - C11*jump(T,n)
>
> # for external/boundary (ds) faces:
> # T_hat = T
> # q_hat = q
>
> # This does not work
> a = inner(q,v_q)*dx \
> + T*div(v_q)*dx \
> - avg(T)*inner(v_q,n)*dS \
> - T*inner(v_q,n)*ds \
> + inner(q, grad(v_T))*dx \
> - v_T*dot(avg(grad(T)),n)*dS + v_T*dot(jump(T,n),n)*dS \
> - v_T*inner(q,n)*ds
>
> L = zero*v_T*dx
>
> w = Function(W)
> solve(a==L,w)
> (u,p) = w.split()
>
> plot(u, title="Velocity u", interactive=True)
> plot(p, title="Pressure p", interactive=True)
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Thank you sir. I tried this:

a = inner(q,v_q)*dx \
  + T*div(v_q)*dx \
  - avg(T)*inner(v_q,n('+'))*dS \
  - T*inner(v_q,n('+'))*ds \
  + inner(q, grad(v_T))*dx \
  - v_T*dot(avg(grad(T)),n('+'))*dS + v_T*dot(jump(T,n),n('+'))*dS \
  - v_T*inner(q,n('+'))*ds

Now I am getting:
  File "/home/yye00/Work/FEniCS/lib/python2.7/site-packages/ufl/log.py", line 148, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Form argument must be restricted.

Can you confirm if I should/shouldn't restrict the jump(T,n) term to say jump(T,n('+')) because I got a
ufl.log.UFLException: Not expecting twice restricted expression.
error when I did.

Regards
Yaakoub

You have both spaces DG, so every term to be evaluated on internal facets (measure ds) must have its side ('+') or ('-') chosen or must be used in jump() or avg() function. On the other hand you can't use ('+'), ('-'), avg() or jump() on external facets (dS measure) because every function is defined only on one (inner) side.

> ufl.log.UFLException: Not expecting twice restricted expression.
This means that jump(T,n) does not need 2nd restriction because
jump(T,n) = inner(T('+'), n('+')) - inner(T('-'), n('+')).
So jump(T,n('+')) does not have sense.

Thank you sir. Just to confirm: ds is internal and dS is external or is it the other way around?

> Thank you sir. Just to confirm: ds is internal and dS is external or
> is it the other way around?

Other way around.

Following your instructions sir I tried again. I have this formulation now:

a = inner(q,v_q)*dx \
  + T*div(v_q)*dx \
  - avg(T)*inner(v_q,n)*dS \
  - T*inner(v_q,n)*ds \
  + inner(q, grad(v_T))*dx \
  - v_T*dot(avg(grad(T('-'))),n)*dS + v_T*dot(jump(T,n),n)*dS \
  - v_T*inner(q,n)*ds

Nothing restricted on ds or dx, and the only term that needs restriction was the avg(grad(T('-')). I am sorry but I must be missing something as I get the same error: FacetNormal must be restricted. I also tried taking out the avg(grad(T)) term. It still gave me the same error. Can you please have a look and see if anything stands out?

Please note
(q, T) = TrialFunctions(W)
(v_q,v_T) = TestFunctions(W)
I do not need to restrict the test functions as they are only locally supported right?

> a = inner(q,v_q)*dx \
> + T*div(v_q)*dx \
> - avg(T)*inner(v_q,n)*dS \
here v_q, n must be restricted!
> - T*inner(v_q,n)*ds \
> + inner(q, grad(v_T))*dx \
> - v_T*dot(avg(grad(T('-'))),n)*dS + v_T*dot(jump(T,n),n)*dS \
v_T must be restricted!

Take grad(avg(T)) rather than avg(grad(T('-'))). This is double
restriction (('-') and avg). Consider avg and grad commute.
Taking the dot product afterwards is wrong - T is scalar.

In dot(jump(T,n),n) second n must be restricted! But it is an overuse of
FacetNormal: dot(jump(T,n),n('+')) = dot(T('+') * n('+') - T('-') *
n('+'), n('+')) = T('+') - T('-')
> - v_T*inner(q,n)*ds
>
> Nothing restricted on ds or dx, and the only term that needs
> restriction was the avg(grad(T('-')). I am sorry but I must be
> missing something as I get the same error: FacetNormal must be
> restricted. I also tried taking out the avg(grad(T)) term. It still
> gave me the same error. Can you please have a look and see if
> anything stands out?
Rules are easy: every DG function (even FacetNormal) MUST be restricted
in interior facets integral. Otherwise you don't know which value take.
But apply avg() and jump() operations on nonrestriced terms because
these apply appropriate restrictions by themselves.

If this does not help, consider looking to fenics book. There are DG
examples.

> Please note
> (q, T) = TrialFunctions(W)
> (v_q,v_T) = TestFunctions(W)
> I do not need to restrict the test functions as they are only locally
> supported right?
Nope. Test function has same number of degrees of freedom as trial
function has (resulting matrix is square, not rectangular).

I guess jump would be defined as
  jump(T,n) = inner(T('+'), n('+')) + inner(T('-'), n('-'))

instead of previous
  jump(T,n) = inner(T('+'), n('+')) - inner(T('-'), n('+'))

so it could work with n=arbitrary DG vector not only for n=FacetNormal which has special property n('+') = -n('-')

But I am not sure about this. You would better check FEniCS book or UFL manual.

Thank you sir that made a lot of sense.I think I have got it (it runs without errors):

a = inner(q,v_q)*dx \
  + T*div(v_q)*dx \
  - avg(T)*inner(v_q('-'),n('+'))*dS \
  - T*inner(v_q,n)*ds \
  + inner(q, grad(v_T))*dx \
  - v_T('-')*inner(grad(avg(T)) ,n('+'))*dS \
  + eight('-')*v_T('-')*inner(jump(T,n),n('+'))*dS \
  - v_T*inner(q,n)*ds

L = Constant(0.0)*v_T*dx

While it runs without errors, I do not get the BC's enforced but I think I can manage from here. Thank you very much for your time.

Regards
Yaakoub

Very sorry, one last issue:

noslip = Constant((0.0, 0.0))

# the following does work
bc = [DirichletBC(W.sub(1), Constant(1.0), left),
      DirichletBC(W.sub(1), Constant(0.0), right),
      DirichletBC(W.sub(0), noslip, top),
      DirichletBC(W.sub(0), noslip, bottom)]

w = Function(W)
solve(a==L,w, bc)
(q,T) = w.split()

where left, right top and bottom are properly defined. I am expecting a flat T from 1 to 0 left to right and a uniform q of value (1.0,0). Not quite what I am getting. Am I applying the BC's wrong?

You're expecting q=(1,0), but it's not consistent with bc for q
> noslip = Constant((0.0, 0.0))

Also make sure that DirichletBC finds appropriate dofs. It maybe needs changing a method
  DirichletBC(W.sub(1), Constant(1.0), left, 'geometric')
for DG spaces.

Also I'm not sure if introducing dirichlet bc for T is mathematically plausible as dirichlet bc becomes natural one (set by appropriate surface terms in equation) and neumann bc becomes dirichlet bc for flux. I would reccomend reading http://fenicsproject.org/documentation/dolfin/1.1.0/python/demo/pde/mixed-poisson/python/documentation.html

Thank you sir. I changed it to remove the BC for q . It now reads:

bc = [DirichletBC(W.sub(1), Constant(1.0), left),
      DirichletBC(W.sub(1), Constant(0.0), right)]

however the q is now (0,0) everywhere and T is 0 everywhere.

You are correct sir, the boundary condition is wrong. I should have added it to the variational form. Going by the mixed-poisson example, my "u0" is not zero.

What would be the most appropriate solution to this? define an expression that is 1.0 on x[0]=1-DOLFIN_EPS and 0.0 otherwise or is there a way to break down ds into the sum of say ds(i) for i in 0,1,2,3 ?

Thank you very much for all your help sir.

On Thu, 17 Jan 2013 20:01:14 -0000
Yaakoub El Khamra <email address hidden> wrote:
> Question #219369 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/219369
>
> Yaakoub El Khamra posted a new comment:
>
> You are correct sir, the boundary condition is wrong. I should have
> added it to the variational form. Going by the mixed-poisson example,
> my "u0" is not zero.
>
> What would be the most appropriate solution to this? define an
> expression that is 1.0 on x[0]=1-DOLFIN_EPS and 0.0 otherwise or is
> there a way to break down ds into the sum of say ds(i) for i in
> 0,1,2,3 ?
Creating measures ds(i) is explained somewhere in FEniCS book.

Shortly you can do:

facet_domains = FacetFunction('size_t', mesh)
facet_domains.set_all(0) # MeshFunction needs explicit initialization!
left.mark(facet_domains, 1)
top.mark(facet_domains, 2)
right.mark(facet_domains, 3)
bottom.mark(facet_domains, 4)

where left,top,right,bottom are appropriate instances of SubDomain.
Then you end up with facet_domains MeshFunction which is 0 on interior
facets on 1,2,3,4 on exterior facets. Then you redifine define

ds = ds[facet_domains]

and you use ds(1), ds(2), ds(3), ds(4). ds means ds(0) which is empty.

Thank you very much sir. I made the changes and the boundary condition show up but I get the wrong answer. I think it is time I sat down and read the fenics book more carefully. Thank you very much for your help. I will update launchpad with the corrected code when I finish it, perhaps add it to the undocumented demos.

Thank you very much.

Regards
Yaakoub