ffc team mailing list archive

[Bartosz Sawicki <sawickib@xxxxxxxxxxxxx>]

On 30/01/09 03:49 PM, Kristian Oelgaard wrote:
> Quoting Bartosz Sawicki <sawickib@xxxxxxxxxxxxx>:
>
> > On 30/01/09 01:04 AM, Kristian Oelgaard wrote:
> > > Quoting Bartosz Sawicki <sawickib@xxxxxxxxxxxxx>:
> > >
> > > > 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.
> > > 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.
> > 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.
>
> Embarrassing... but yes, you have to restrict all functions when computing
> interior facet integrals.
>  
> > 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:
>
> I don't know what you want to compute, so I can't help you here.
>  
> > 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. :(
>
> 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



>  
> > Have you got any other idea, how to describe separate integration for 
> > every cell?
>
> 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
> > >  
> > > > 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
> > > >
> > > >
> > > >
> > > >
> > > > _______________________________________________
> > > > FFC-dev mailing list
> > > > FFC-dev@xxxxxxxxxx
> > > > http://www.fenics.org/mailman/listinfo/ffc-dev