Subject: [FFC-dev] Internal integration in discontinuous vector element

From: Bartosz Sawicki <[email-masked]>
Date: Thu, 29 Jan 2009 18:20:04 -0700

I have a problem how to define form for integration over internal facets 
inside VectorElement("Discontinuous Lagrange", "tetrahedron", 0).
Having vector function "j" which is constant over element, I want to 
integrate over all internal facets term dot(j,n), where n is facet 
normal vector.

Consider form:
=====
element = FiniteElement("Discontinuous Lagrange", "tetrahedron", 0)
vectorelement = VectorElement("Discontinuous Lagrange", "tetrahedron", 0)

v = TestFunction(element)
e = TrialFunction(element)
j = Function(vectorelement)
n = FacetNormal("tetrahedron")

a = v*e*dx
L = v*dot(j, n)*dS
=====

Unfortunately it doesn't work. FFC communicates that "Integrand must be 
restricted ('+') or ('-') in an interior facet integral". That is 
confusing for me.

Maybe interior integration operator works differently?
So, would it be possible to use form language to describe such simple 
algorithm?

iterate for every cell c:
   e_i = 0
   integrate for every facet in cell:
     e_i += dot(j,n)
   e(c) = e_i

regrds.
BArtek

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 30 Jan 2009 09:04:00 +0100

> restricted ('+') or ('-') in an interior facet integral". That is 

Functions on interior facets can potentially be double valued (in your case they
really are). So you need to tell FFC which value you want.

Maybe you want something like:

L = v*jump(j, n)*dS = v*(dot(j('+'), n('+')) + dot(j('-'), n('-')))*dS

The discontinuous operators are explained in more detail in the FFC manual.

Kristian
 
> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev
>

---

From: Bartosz Sawicki <[email-masked]>
Date: Fri, 30 Jan 2009 13:03:45 -0700

> The discontinuous operators are explained in more detail in the FFC manual.

Thanks for your answer.

I've read FCC manual carefully, but there is "..., one _may_ restrict
expressions to the positive or negative side of the facet". What sugest
that I don't have to use restrictions.
It seems that, when using dS operator we _have to_ restrict every
function. Even your example doesn't compile, because v in not restricted.

Something what compiles correctly is:

a = v*e*dx
L = v('+')*dot(j('+'),n('+'))*dS + v('-')*dot(j('-'),n('-'))*dS

but I'm not sure if it is what I need. Moreover when I've tried to a 
solve the problem:

Function e;
VariationalProblem problem(a,L);
problem.solve(e);

I get runtime error "*** Error: Coefficient 1 has not been defined."
So something is definitely wrong with this form. :(

Have you got any other idea, how to describe separate integration for 
every cell?

regrs.
BArtek

---

From: Bartosz Sawicki <[email-masked]>
Date: Fri, 30 Jan 2009 16:52:15 -0700

On 30/01/09 03:49 PM, Kristian Oelgaard wrote:
> Embarrassing... but yes, you have to restrict all functions when computing
> interior facet integrals.
>  
> I don't know what you want to compute, so I can't help you here.
>  
>> Function e;
> Did you initialise the normal vector? and the functin j in you linear form?
> Have a look at the DG demos demo/pde/dg/

Thanks again. This time your advice helped perfectly ! :)
I did not initialize FacetNormal. I fix it, and now everything compiles,
and results look reasonable.

Have a nice weekend :)
BArtek

> The linear form that you have compiling should do just that.
> You can shorten it a bit by doing:
> 
> L = mult(v, dot(j, n))('+')*dS + mult(v, dot(j, n))('-')*dS
> 
> which is equivalent, I tested it :)
> 
> Kristian
> 
>> regrs.
>> BArtek
>>
>>
>>
>>
>>
>>> Kristian
>>>