including facet-based variable coefficients in bilinear form

Question

Hello,

I would like to implement in FEniCS a nonlinear finite element method, for example one of its parts is of the form \sum_{facets F} a_F(u_h) (\nabla u_h, \nabla v_h) dS

where a_F(u_h) is a function defined on facets F and its values are real numbers which depend on some conditions for \nabla u_h.

Is that possible to implement this in FEniCS in python level?

Thank you in advance, Fotini

Answer 1

Yes. $\phantom{ } $

Comments

Hi Jan,

can you please give me some hints? I have been looking on the web and I haven't found
a way to implement that ...

thank you,
Fotini

---

Check nonlinear-poisson demo for a solution of nonlinear problems. For your a_F(u_h) term you would need to invent this by yourself probably - I guess that this can be naturally defined by appropriate choice of finite element.

---

I am not sure that you can defined elements on facet on FEniCS.
I do not see how the file of nonlinear poisson can help me ...

what I need is the following:
To define a function
a_F(u_h) = tau_F if \nabla u_h^+, \nabla u_h^- on F satisfy certain conditions
and the values tau_F depend each time on the current values of \nabla u_h^+ and \nabla u_h^-.

Maybe I did not understand your suggestion ...

thanks again
Fotini

---

Your formula
$$ \sum_{facets F} a_F(u_h) (\nabla u_h, \nabla v_h) dS $$
is basically interior facet integral and $a_F(u_h)$ is merely a function of $u_h^\pm$. Where is the problem?

---

Hi Jan,

thanks for you answer. This bilinear form is included in a fixed point iteration.
a_F(u_h) can be zero if u_h has a minimum to a vertex connected to F, 1 in some other areas, can be a real number which is based on the local l1-norm of \nabla u_h+ and \nabla u_h^- etc

So, I have to have a loop over the facets, to check if some certain conditions are
satisfied and then to decide the value of a_F(u_h).
I have a code that the values of a_F(u_h) are 0 or 1 and then I use a FacetFunction.

Maybe it is more clear what I need to do.

best,
fotini

---

If I understand it correctly, then it is not probably possible to use Newton solver. You can do that if you can express a dependence of your integrand on unknowns using UFL (i.e. simple operations). Note that multi-branch function is possible using conditional but it should be continuous in order to Jacobian computed by automatic differentiation was correct.

Fixed-point iteration is probably possible. You already have the procedure for computation of $a_F(u_h)$ term and you have to figure out the FE space which preserves its values and project/interpolate onto it. Isn't it the RT space? Or some Quadrature space?

---

Hi again,

I haven't managed yet to implement the above method in FEniCS. I have computed
the coefficients based on facets and I have created a CR function a_uh whose
coefficient corresponding to F is exactly the value a_F(uh).
I need to include now in the bilinear form quantities of the form

V = FunctionSpace(mesh, "CG", 1)
VCR = FunctionSpace(mesh, "CR", 1)
uh = TrialFunction(V)
vh = TestFunction(V)
a_uh = Function(VCR)  #a_uh(facet F dof) = a_F(u_h)

1/h('+')*a_uh('+')*dot(grad(uh)('+'), grad(vh)('+'))*dS 
   1/h('-')* + a_uh('-')*dot(grad(uh)('-'), grad(vh)('-'))*dS

but the above code will not give the correct result since 1/h('+') \int_F a_uh('+') dS \ne a_F(uh).

If I also tried the restriction on facet
VCR = FunctionSpace(mesh, "CR", 1, facet)
but I am not sure what does it mean ...

One solution would be to have the midpoint rule to evaluate the integral on F. Does it exist?
Which are the DOFs for the space FunctionSpace(mesh, "Quadrature", 1, facet)?

Thank you in advance for your help.
Fotini