different quadature degree

On Wed, 10 Apr 2013 17:31:07 -0000
Martin Sandve Alnæs <email address hidden> wrote:
> Your question #226398 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/226398
>
> Status: Open => Answered
>
> Martin Sandve Alnæs proposed the following answer:
> Yes, there is a ufl demo file in ffc showing how. Grep for quad.

Ok, there is Poisson.ufl:
_________________________________________________
element = FiniteElement("Lagrange", triangle, 1)

u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)

dx0_1 = dx(0, {"quadrature_order": 1})
dx0_2 = dx(0, {"quadrature_order": 2})

a = dot(grad(v), grad(u))*dx0_1
L = v*f*dx0_2
_________________________________________________

Is there a reason why form a has degree 1? I would expect that degree 0
suffices for exact quadrature as gradients of both u and v have
polynomial degree 0.

It seems that this code
________________________________________
from dolfin import *
mesh = UnitSquareMesh(2,2)
V = FunctionSpace(mesh, 'CG', 2)
u = Function(V)
a = u**3*dx(0, {'quadrature_order': 2})
b = u**4*dx(0, {'quadrature_order': 8})
M = assemble(a+b)
________________________________________
compiles form a+b together with only one quadrature rule. Maybe I don't
read its ufc code correctly.

It seems for me that only solution is this:
________________________________________
from dolfin import *
mesh = UnitSquareMesh(2,2)
V = FunctionSpace(mesh, 'CG', 2)
u = Function(V)
a = u**3*dx(0, {'quadrature_order': 2})
b = u**4*dx(0, {'quadrature_order': 8})
M = assemble(a)
M += assemble(b) # for forms of rank = 0
#M = assemble(b, tensor=M, add_values=True) # for forms of rank > 0
________________________________________
Am I right?

Jan

> Den 10. apr. 2013 18:45 skrev "Jan Blechta" <
> <email address hidden>> følgende:
>
> > New question #226398 on DOLFIN:
> > https://answers.launchpad.net/dolfin/+question/226398
> >
> > Is it possible to JIT compile a form by FFC whose parts would have
> > different qudrature degrees? Consider for example
> >
> > u = Function(V)
> > v = TestFunction(V)
> > phi = Function(Q)
> > F = u*v*dx + phi*v*dx
> >
> > Is there the way to tell FFC to compile first term with some (say 2
> > or automatic) degree and second with another (say 4) quadrature
> > degree? As far as I understandt it I would say that compiled F has
> > common degree corresponding to term with the highest degree
> > estimated by UFL, doesn't it?
> >
> > Or is only solution define forms part by part, assemble each part
> > separately and add tensors together?
> >
> > --
> > You received this question notification because you are a member of
> > DOLFIN Team, which is an answer contact for DOLFIN.
> >
>

Is it using tensor representation? If so try adding "representation":"quadrature" to the dicts. Or maybe this is not completely implemented in ffc; I know there is infrastructure for it but I have never used it myself.

On Wed, 10 Apr 2013 19:16:15 -0000
Martin Sandve Alnæs <email address hidden> wrote:
> Your question #226398 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/226398
>
> Status: Open => Answered
>
> Martin Sandve Alnæs proposed the following answer:
> Is it using tensor representation? If so try adding
> "representation":"quadrature" to the dicts. Or maybe this is not
> completely implemented in ffc; I know there is infrastructure for it
> but I have never used it myself.
>

Thanks, Martin. I checked again and it seems that compiled code for form
a+b contains geometry tensors / quadrature rules of both form a and
form b. So it seems working on side of UFL, FFC, UFC for both
representations. But I would like to know that DOLFIN assembler takes
this into account correctly (i.e it assembles one part of form with
one degree and other part with other degree) and I've no clue how to
investigate this. Do you? Or can you suspect a developer who is
aware of the thing?

Jan

On 11 April 2013 07:26, Jan Blechta
<email address hidden> wrote:
> Question #226398 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/226398
>
> Status: Answered => Open
>
> Jan Blechta is still having a problem:
> On Wed, 10 Apr 2013 19:16:15 -0000
> Martin Sandve Alnæs <email address hidden> wrote:
>> Your question #226398 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/226398
>>
>> Status: Open => Answered
>>
>> Martin Sandve Alnæs proposed the following answer:
>> Is it using tensor representation? If so try adding
>> "representation":"quadrature" to the dicts. Or maybe this is not
>> completely implemented in ffc; I know there is infrastructure for it
>> but I have never used it myself.
>>
>
> Thanks, Martin. I checked again and it seems that compiled code for form
> a+b contains geometry tensors / quadrature rules of both form a and
> form b. So it seems working on side of UFL, FFC, UFC for both
> representations. But I would like to know that DOLFIN assembler takes
> this into account correctly (i.e it assembles one part of form with
> one degree and other part with other degree) and I've no clue how to
> investigate this.

DOLFIN just makes a call to ufc::tabulate_tensor(.....), so it's
oblivious to the integration scheme. To see what's happening you'll
need to look at the generated tabulate_tenor function.

Kristian Oelgaard knows most about what goes on with the quadrature
representation.

Garth

> Do you? Or can you suspect a developer who is
> aware of the thing?
>
> Jan
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

On 10 April 2013 20:41, Jan Blechta <email address hidden>wrote:

> Question #226398 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/226398
>
> Status: Answered => Open
>
> Jan Blechta is still having a problem:
> On Wed, 10 Apr 2013 17:31:07 -0000
> Martin Sandve Alnæs <email address hidden> wrote:
> > Your question #226398 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/226398
> >
> > Status: Open => Answered
> >
> > Martin Sandve Alnæs proposed the following answer:
> > Yes, there is a ufl demo file in ffc showing how. Grep for quad.
>
> Ok, there is Poisson.ufl:
> _________________________________________________
> element = FiniteElement("Lagrange", triangle, 1)
>
> u = TrialFunction(element)
> v = TestFunction(element)
> f = Coefficient(element)
>
> dx0_1 = dx(0, {"quadrature_order": 1})
> dx0_2 = dx(0, {"quadrature_order": 2})
>
> a = dot(grad(v), grad(u))*dx0_1
> L = v*f*dx0_2
> _________________________________________________
>
> Is there a reason why form a has degree 1? I would expect that degree 0
> suffices for exact quadrature as gradients of both u and v have
> polynomial degree 0.
>
>
Using degree 0 or 1 will result in a one point integration scheme, so it
doesn't matter.
Note that you should use 'quadrature_order' --> 'quadrature_degree'.

> It seems that this code
> ________________________________________
> from dolfin import *
> mesh = UnitSquareMesh(2,2)
> V = FunctionSpace(mesh, 'CG', 2)
> u = Function(V)
> a = u**3*dx(0, {'quadrature_order': 2})
> b = u**4*dx(0, {'quadrature_order': 8})
> M = assemble(a+b)
> ________________________________________
> compiles form a+b together with only one quadrature rule. Maybe I don't
> read its ufc code correctly.
>
>
Compiling the following form with FFC 1.0 and quadrature representation:

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

f = Coefficient(element)

m0 = f**3*dx(0, {'quadrature_degree': 2})
m1 = f**4*dx(0, {'quadrature_degree': 8})
M = (m0+m1)

results in the code:

[snip]

    for (unsigned int ip = 0; ip < 25; ip++)
    {

      // Coefficient declarations.
      double F0 = 0.0;

      // Total number of operations to compute function values = 12
      for (unsigned int r = 0; r < 6; r++)
      {
        F0 += FE1[ip][r]*w[0][r];
      }// end loop over 'r'

      // Number of operations for primary indices: 6
      // Number of operations to compute entry: 6
      A[0] += F0*F0*F0*F0*W25[ip]*det;
    }// end loop over 'ip'

    // Loop quadrature points for integral.
    // Number of operations to compute element tensor for following IP loop
= 51
    for (unsigned int ip = 0; ip < 3; ip++)
    {

      // Coefficient declarations.
      double F0 = 0.0;

      // Total number of operations to compute function values = 12
      for (unsigned int r = 0; r < 6; r++)
      {
        F0 += FE0[ip][r]*w[0][r];
      }// end loop over 'r'

      // Number of operations for primary indices: 5
      // Number of operations to compute entry: 5
      A[0] += F0*F0*F0*W3[ip]*det;
    }// end loop over 'ip'
  }

As you see two loops over integration points, as one would expect.
I do, however, receive the warning 'Quadrature degree must be equal within
each sub domain, using degree 8.' from FFC,
but I guess that's related to the tensor representation.
And as Garth says, Dolfin knows nothing about the integration schemes which
are used inside tabulate_tensor(), it just calls the function.

> It seems for me that only solution is this:
> ________________________________________
> from dolfin import *
> mesh = UnitSquareMesh(2,2)
> V = FunctionSpace(mesh, 'CG', 2)
> u = Function(V)
> a = u**3*dx(0, {'quadrature_order': 2})
> b = u**4*dx(0, {'quadrature_order': 8})
> M = assemble(a)
> M += assemble(b) # for forms of rank = 0
> #M = assemble(b, tensor=M, add_values=True) # for forms of rank > 0
> ________________________________________
> Am I right?
>
> Jan
>
>
> > Den 10. apr. 2013 18:45 skrev "Jan Blechta" <
> > <email address hidden>> følgende:
> >
> > > New question #226398 on DOLFIN:
> > > https://answers.launchpad.net/dolfin/+question/226398
> > >
> > > Is it possible to JIT compile a form by FFC whose parts would have
> > > different qudrature degrees? Consider for example
> > >
> > > u = Function(V)
> > > v = TestFunction(V)
> > > phi = Function(Q)
> > > F = u*v*dx + phi*v*dx
> > >
> > > Is there the way to tell FFC to compile first term with some (say 2
> > > or automatic) degree and second with another (say 4) quadrature
> > > degree? As far as I understandt it I would say that compiled F has
> > > common degree corresponding to term with the highest degree
> > > estimated by UFL, doesn't it?
> > >
> > > Or is only solution define forms part by part, assemble each part
> > > separately and add tensors together?
> > >
> > > --
> > > You received this question notification because you are a member of
> > > DOLFIN Team, which is an answer contact for DOLFIN.
> > >
> >
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

On Thu, 11 Apr 2013 06:41:39 -0000
Kristian B. Ølgaard <email address hidden> wrote:
> Your question #226398 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/226398
>
> Kristian B. Ølgaard proposed the following answer:
> On 10 April 2013 20:41, Jan Blechta
> <email address hidden>wrote:
>
> > Question #226398 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/226398
> >
> > Status: Answered => Open
> >
> > Jan Blechta is still having a problem:
> > On Wed, 10 Apr 2013 17:31:07 -0000
> > Martin Sandve Alnæs <email address hidden> wrote:
> > > Your question #226398 on DOLFIN changed:
> > > https://answers.launchpad.net/dolfin/+question/226398
> > >
> > > Status: Open => Answered
> > >
> > > Martin Sandve Alnæs proposed the following answer:
> > > Yes, there is a ufl demo file in ffc showing how. Grep for quad.
> >
> > Ok, there is Poisson.ufl:
> > _________________________________________________
> > element = FiniteElement("Lagrange", triangle, 1)
> >
> > u = TrialFunction(element)
> > v = TestFunction(element)
> > f = Coefficient(element)
> >
> > dx0_1 = dx(0, {"quadrature_order": 1})
> > dx0_2 = dx(0, {"quadrature_order": 2})
> >
> > a = dot(grad(v), grad(u))*dx0_1
> > L = v*f*dx0_2
> > _________________________________________________
> >
> > Is there a reason why form a has degree 1? I would expect that
> > degree 0 suffices for exact quadrature as gradients of both u and v
> > have polynomial degree 0.
> >
> >
> Using degree 0 or 1 will result in a one point integration scheme, so
> it doesn't matter.
> Note that you should use 'quadrature_order' --> 'quadrature_degree'.
>
>
> > It seems that this code
> > ________________________________________
> > from dolfin import *
> > mesh = UnitSquareMesh(2,2)
> > V = FunctionSpace(mesh, 'CG', 2)
> > u = Function(V)
> > a = u**3*dx(0, {'quadrature_order': 2})
> > b = u**4*dx(0, {'quadrature_order': 8})
> > M = assemble(a+b)
> > ________________________________________
> > compiles form a+b together with only one quadrature rule. Maybe I
> > don't read its ufc code correctly.
> >
> >
> Compiling the following form with FFC 1.0 and quadrature
> representation:
>
> element = FiniteElement("Lagrange", triangle, 2)
>
> f = Coefficient(element)
>
> m0 = f**3*dx(0, {'quadrature_degree': 2})
> m1 = f**4*dx(0, {'quadrature_degree': 8})
> M = (m0+m1)
>
> results in the code:
>
>
> [snip]
>
> for (unsigned int ip = 0; ip < 25; ip++)
> {
>
> // Coefficient declarations.
> double F0 = 0.0;
>
> // Total number of operations to compute function values = 12
> for (unsigned int r = 0; r < 6; r++)
> {
> F0 += FE1[ip][r]*w[0][r];
> }// end loop over 'r'
>
> // Number of operations for primary indices: 6
> // Number of operations to compute entry: 6
> A[0] += F0*F0*F0*F0*W25[ip]*det;
> }// end loop over 'ip'
>
> // Loop quadrature points for integral.
> // Number of operations to compute element tensor for following
> IP loop = 51
> for (unsigned int ip = 0; ip < 3; ip++)
> {
>
> // Coefficient declarations.
> double F0 = 0.0;
>
> // Total number of operations to compute function values = 12
> for (unsigned int r = 0; r < 6; r++)
> {
> F0 += FE0[ip][r]*w[0][r];
> }// end loop over 'r'
>
> // Number of operations for primary indices: 5
> // Number of operations to compute entry: 5
> A[0] += F0*F0*F0*W3[ip]*det;
> }// end loop over 'ip'
> }
>
> As you see two loops over integration points, as one would expect.
> I do, however, receive the warning 'Quadrature degree must be equal
> within each sub domain, using degree 8.' from FFC,
> but I guess that's related to the tensor representation.
> And as Garth says, Dolfin knows nothing about the integration schemes
> which are used inside tabulate_tensor(), it just calls the function.

Ok. It seems that
"Quadrature degree must be equal within each sub domain, using degree n"
warning does not apply for quadrature representation but is displayed.
Is it a bug?

Is stdout of FFC accessible somehow when JIT compiling from DOLFIN?
Perhapse Johan could know?

Jan

PS: Here is only simple demo of what can happen. I wrote it before I
realized that everything works quite satisfactory. But it could be
useful to somebody. Just skip it if not insterested. I checked codes
generated by FFC trunk for this code
_________________________________
U = FiniteElement('CG', 'triangle', 2)
V = FiniteElement('CG', 'triangle', 3)
u = Coefficient(U)
v = Coefficient(V)
aT = u**2*dx(0, {'representation': 'tensor'})
bT = v**4*dx(0, {'representation': 'tensor'})
aQ = u**2*dx(0, {'representation': 'quadrature'})
bQ = v**4*dx(0, {'representation': 'quadrature'})
aD = u**2*dx(0)
bD = v**4*dx(0)
________________________________
and found that

form | representaion
----------------------
aT + bT | G
aQ + bQ | W6 + W49
aD + bD | W49
aT + bQ | W6 + W49
aQ + bT | W6 + W49
aQ + bD | W6 + W49
aD + bQ | W49
aT + bD | W6 + W49
aD + bT | G
----------------------

where Wnp is quadrature with np = number of points and G is tensor. The
moral is that one has to check the generated code if it does what's
intended. It is also somehow readable from FFC stdout "Analyzing forms"
step.

>
>
> > It seems for me that only solution is this:
> > ________________________________________
> > from dolfin import *
> > mesh = UnitSquareMesh(2,2)
> > V = FunctionSpace(mesh, 'CG', 2)
> > u = Function(V)
> > a = u**3*dx(0, {'quadrature_order': 2})
> > b = u**4*dx(0, {'quadrature_order': 8})
> > M = assemble(a)
> > M += assemble(b) # for forms of rank = 0
> > #M = assemble(b, tensor=M, add_values=True) # for forms of rank > 0
> > ________________________________________
> > Am I right?
> >
> > Jan
> >
> >
> > > Den 10. apr. 2013 18:45 skrev "Jan Blechta" <
> > > <email address hidden>> følgende:
> > >
> > > > New question #226398 on DOLFIN:
> > > > https://answers.launchpad.net/dolfin/+question/226398
> > > >
> > > > Is it possible to JIT compile a form by FFC whose parts would
> > > > have different qudrature degrees? Consider for example
> > > >
> > > > u = Function(V)
> > > > v = TestFunction(V)
> > > > phi = Function(Q)
> > > > F = u*v*dx + phi*v*dx
> > > >
> > > > Is there the way to tell FFC to compile first term with some
> > > > (say 2 or automatic) degree and second with another (say 4)
> > > > quadrature degree? As far as I understandt it I would say that
> > > > compiled F has common degree corresponding to term with the
> > > > highest degree estimated by UFL, doesn't it?
> > > >
> > > > Or is only solution define forms part by part, assemble each
> > > > part separately and add tensors together?
> > > >
> > > > --
> > > > You received this question notification because you are a
> > > > member of DOLFIN Team, which is an answer contact for DOLFIN.
> > > >
> > >
> >
> > --
> > You received this question notification because you are a member of
> > DOLFIN Team, which is an answer contact for DOLFIN.
> >
>

I think that... when representation is "auto", ffc tries to build the
tensor representation. If that fails at some point it throws an
exception and the quadrature representation is used instead. Probably
not the most efficient or elegant approach but it works. So the
warning probably occurs while tensor representation is tried.

On 11 April 2013 16:11, Jan Blechta
<email address hidden> wrote:
> Question #226398 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/226398
>
> Status: Answered => Open
>
> Jan Blechta is still having a problem:
> On Thu, 11 Apr 2013 06:41:39 -0000
> Kristian B. Ølgaard <email address hidden> wrote:
>> Your question #226398 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/226398
>>
>> Kristian B. Ølgaard proposed the following answer:
>> On 10 April 2013 20:41, Jan Blechta
>> <email address hidden>wrote:
>>
>> > Question #226398 on DOLFIN changed:
>> > https://answers.launchpad.net/dolfin/+question/226398
>> >
>> > Status: Answered => Open
>> >
>> > Jan Blechta is still having a problem:
>> > On Wed, 10 Apr 2013 17:31:07 -0000
>> > Martin Sandve Alnæs <email address hidden> wrote:
>> > > Your question #226398 on DOLFIN changed:
>> > > https://answers.launchpad.net/dolfin/+question/226398
>> > >
>> > > Status: Open => Answered
>> > >
>> > > Martin Sandve Alnæs proposed the following answer:
>> > > Yes, there is a ufl demo file in ffc showing how. Grep for quad.
>> >
>> > Ok, there is Poisson.ufl:
>> > _________________________________________________
>> > element = FiniteElement("Lagrange", triangle, 1)
>> >
>> > u = TrialFunction(element)
>> > v = TestFunction(element)
>> > f = Coefficient(element)
>> >
>> > dx0_1 = dx(0, {"quadrature_order": 1})
>> > dx0_2 = dx(0, {"quadrature_order": 2})
>> >
>> > a = dot(grad(v), grad(u))*dx0_1
>> > L = v*f*dx0_2
>> > _________________________________________________
>> >
>> > Is there a reason why form a has degree 1? I would expect that
>> > degree 0 suffices for exact quadrature as gradients of both u and v
>> > have polynomial degree 0.
>> >
>> >
>> Using degree 0 or 1 will result in a one point integration scheme, so
>> it doesn't matter.
>> Note that you should use 'quadrature_order' --> 'quadrature_degree'.
>>
>>
>> > It seems that this code
>> > ________________________________________
>> > from dolfin import *
>> > mesh = UnitSquareMesh(2,2)
>> > V = FunctionSpace(mesh, 'CG', 2)
>> > u = Function(V)
>> > a = u**3*dx(0, {'quadrature_order': 2})
>> > b = u**4*dx(0, {'quadrature_order': 8})
>> > M = assemble(a+b)
>> > ________________________________________
>> > compiles form a+b together with only one quadrature rule. Maybe I
>> > don't read its ufc code correctly.
>> >
>> >
>> Compiling the following form with FFC 1.0 and quadrature
>> representation:
>>
>> element = FiniteElement("Lagrange", triangle, 2)
>>
>> f = Coefficient(element)
>>
>> m0 = f**3*dx(0, {'quadrature_degree': 2})
>> m1 = f**4*dx(0, {'quadrature_degree': 8})
>> M = (m0+m1)
>>
>> results in the code:
>>
>>
>> [snip]
>>
>> for (unsigned int ip = 0; ip < 25; ip++)
>> {
>>
>> // Coefficient declarations.
>> double F0 = 0.0;
>>
>> // Total number of operations to compute function values = 12
>> for (unsigned int r = 0; r < 6; r++)
>> {
>> F0 += FE1[ip][r]*w[0][r];
>> }// end loop over 'r'
>>
>> // Number of operations for primary indices: 6
>> // Number of operations to compute entry: 6
>> A[0] += F0*F0*F0*F0*W25[ip]*det;
>> }// end loop over 'ip'
>>
>> // Loop quadrature points for integral.
>> // Number of operations to compute element tensor for following
>> IP loop = 51
>> for (unsigned int ip = 0; ip < 3; ip++)
>> {
>>
>> // Coefficient declarations.
>> double F0 = 0.0;
>>
>> // Total number of operations to compute function values = 12
>> for (unsigned int r = 0; r < 6; r++)
>> {
>> F0 += FE0[ip][r]*w[0][r];
>> }// end loop over 'r'
>>
>> // Number of operations for primary indices: 5
>> // Number of operations to compute entry: 5
>> A[0] += F0*F0*F0*W3[ip]*det;
>> }// end loop over 'ip'
>> }
>>
>> As you see two loops over integration points, as one would expect.
>> I do, however, receive the warning 'Quadrature degree must be equal
>> within each sub domain, using degree 8.' from FFC,
>> but I guess that's related to the tensor representation.
>> And as Garth says, Dolfin knows nothing about the integration schemes
>> which are used inside tabulate_tensor(), it just calls the function.
>
> Ok. It seems that
> "Quadrature degree must be equal within each sub domain, using degree n"
> warning does not apply for quadrature representation but is displayed.
> Is it a bug?
>
> Is stdout of FFC accessible somehow when JIT compiling from DOLFIN?
> Perhapse Johan could know?
>
> Jan
>
>
> PS: Here is only simple demo of what can happen. I wrote it before I
> realized that everything works quite satisfactory. But it could be
> useful to somebody. Just skip it if not insterested. I checked codes
> generated by FFC trunk for this code
> _________________________________
> U = FiniteElement('CG', 'triangle', 2)
> V = FiniteElement('CG', 'triangle', 3)
> u = Coefficient(U)
> v = Coefficient(V)
> aT = u**2*dx(0, {'representation': 'tensor'})
> bT = v**4*dx(0, {'representation': 'tensor'})
> aQ = u**2*dx(0, {'representation': 'quadrature'})
> bQ = v**4*dx(0, {'representation': 'quadrature'})
> aD = u**2*dx(0)
> bD = v**4*dx(0)
> ________________________________
> and found that
>
> form | representaion
> ----------------------
> aT + bT | G
> aQ + bQ | W6 + W49
> aD + bD | W49
> aT + bQ | W6 + W49
> aQ + bT | W6 + W49
> aQ + bD | W6 + W49
> aD + bQ | W49
> aT + bD | W6 + W49
> aD + bT | G
> ----------------------
>
> where Wnp is quadrature with np = number of points and G is tensor. The
> moral is that one has to check the generated code if it does what's
> intended. It is also somehow readable from FFC stdout "Analyzing forms"
> step.
>
>>
>>
>> > It seems for me that only solution is this:
>> > ________________________________________
>> > from dolfin import *
>> > mesh = UnitSquareMesh(2,2)
>> > V = FunctionSpace(mesh, 'CG', 2)
>> > u = Function(V)
>> > a = u**3*dx(0, {'quadrature_order': 2})
>> > b = u**4*dx(0, {'quadrature_order': 8})
>> > M = assemble(a)
>> > M += assemble(b) # for forms of rank = 0
>> > #M = assemble(b, tensor=M, add_values=True) # for forms of rank > 0
>> > ________________________________________
>> > Am I right?
>> >
>> > Jan
>> >
>> >
>> > > Den 10. apr. 2013 18:45 skrev "Jan Blechta" <
>> > > <email address hidden>> følgende:
>> > >
>> > > > New question #226398 on DOLFIN:
>> > > > https://answers.launchpad.net/dolfin/+question/226398
>> > > >
>> > > > Is it possible to JIT compile a form by FFC whose parts would
>> > > > have different qudrature degrees? Consider for example
>> > > >
>> > > > u = Function(V)
>> > > > v = TestFunction(V)
>> > > > phi = Function(Q)
>> > > > F = u*v*dx + phi*v*dx
>> > > >
>> > > > Is there the way to tell FFC to compile first term with some
>> > > > (say 2 or automatic) degree and second with another (say 4)
>> > > > quadrature degree? As far as I understandt it I would say that
>> > > > compiled F has common degree corresponding to term with the
>> > > > highest degree estimated by UFL, doesn't it?
>> > > >
>> > > > Or is only solution define forms part by part, assemble each
>> > > > part separately and add tensors together?
>> > > >
>> > > > --
>> > > > You received this question notification because you are a
>> > > > member of DOLFIN Team, which is an answer contact for DOLFIN.
>> > > >
>> > >
>> >
>> > --
>> > You received this question notification because you are a member of
>> > DOLFIN Team, which is an answer contact for DOLFIN.
>> >
>>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

On Thu, 11 Apr 2013 14:56:15 -0000
Martin Sandve Alnæs <email address hidden> wrote:
> Your question #226398 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/226398
>
> Status: Open => Answered
>
> Martin Sandve Alnæs proposed the following answer:
> I think that... when representation is "auto", ffc tries to build the
> tensor representation. If that fails at some point it throws an
> exception and the quadrature representation is used instead. Probably
> not the most efficient or elegant approach but it works. So the
> warning probably occurs while tensor representation is tried.

Nice guess, but no:) It is also printed if you select quadrature
representation manually.

>
> On 11 April 2013 16:11, Jan Blechta
> <email address hidden> wrote:
> > Question #226398 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/226398
> >
> > Status: Answered => Open
> >
> > Jan Blechta is still having a problem:
> > On Thu, 11 Apr 2013 06:41:39 -0000
> > Kristian B. Ølgaard <email address hidden> wrote:
> >> Your question #226398 on DOLFIN changed:
> >> https://answers.launchpad.net/dolfin/+question/226398
> >>
> >> Kristian B. Ølgaard proposed the following answer:
> >> On 10 April 2013 20:41, Jan Blechta
> >> <email address hidden>wrote:
> >>
> >> > Question #226398 on DOLFIN changed:
> >> > https://answers.launchpad.net/dolfin/+question/226398
> >> >
> >> > Status: Answered => Open
> >> >
> >> > Jan Blechta is still having a problem:
> >> > On Wed, 10 Apr 2013 17:31:07 -0000
> >> > Martin Sandve Alnæs <email address hidden> wrote:
> >> > > Your question #226398 on DOLFIN changed:
> >> > > https://answers.launchpad.net/dolfin/+question/226398
> >> > >
> >> > > Status: Open => Answered
> >> > >
> >> > > Martin Sandve Alnæs proposed the following answer:
> >> > > Yes, there is a ufl demo file in ffc showing how. Grep for
> >> > > quad.
> >> >
> >> > Ok, there is Poisson.ufl:
> >> > _________________________________________________
> >> > element = FiniteElement("Lagrange", triangle, 1)
> >> >
> >> > u = TrialFunction(element)
> >> > v = TestFunction(element)
> >> > f = Coefficient(element)
> >> >
> >> > dx0_1 = dx(0, {"quadrature_order": 1})
> >> > dx0_2 = dx(0, {"quadrature_order": 2})
> >> >
> >> > a = dot(grad(v), grad(u))*dx0_1
> >> > L = v*f*dx0_2
> >> > _________________________________________________
> >> >
> >> > Is there a reason why form a has degree 1? I would expect that
> >> > degree 0 suffices for exact quadrature as gradients of both u
> >> > and v have polynomial degree 0.
> >> >
> >> >
> >> Using degree 0 or 1 will result in a one point integration scheme,
> >> so it doesn't matter.
> >> Note that you should use 'quadrature_order' -->
> >> 'quadrature_degree'.
> >>
> >>
> >> > It seems that this code
> >> > ________________________________________
> >> > from dolfin import *
> >> > mesh = UnitSquareMesh(2,2)
> >> > V = FunctionSpace(mesh, 'CG', 2)
> >> > u = Function(V)
> >> > a = u**3*dx(0, {'quadrature_order': 2})
> >> > b = u**4*dx(0, {'quadrature_order': 8})
> >> > M = assemble(a+b)
> >> > ________________________________________
> >> > compiles form a+b together with only one quadrature rule. Maybe I
> >> > don't read its ufc code correctly.
> >> >
> >> >
> >> Compiling the following form with FFC 1.0 and quadrature
> >> representation:
> >>
> >> element = FiniteElement("Lagrange", triangle, 2)
> >>
> >> f = Coefficient(element)
> >>
> >> m0 = f**3*dx(0, {'quadrature_degree': 2})
> >> m1 = f**4*dx(0, {'quadrature_degree': 8})
> >> M = (m0+m1)
> >>
> >> results in the code:
> >>
> >>
> >> [snip]
> >>
> >> for (unsigned int ip = 0; ip < 25; ip++)
> >> {
> >>
> >> // Coefficient declarations.
> >> double F0 = 0.0;
> >>
> >> // Total number of operations to compute function values = 12
> >> for (unsigned int r = 0; r < 6; r++)
> >> {
> >> F0 += FE1[ip][r]*w[0][r];
> >> }// end loop over 'r'
> >>
> >> // Number of operations for primary indices: 6
> >> // Number of operations to compute entry: 6
> >> A[0] += F0*F0*F0*F0*W25[ip]*det;
> >> }// end loop over 'ip'
> >>
> >> // Loop quadrature points for integral.
> >> // Number of operations to compute element tensor for following
> >> IP loop = 51
> >> for (unsigned int ip = 0; ip < 3; ip++)
> >> {
> >>
> >> // Coefficient declarations.
> >> double F0 = 0.0;
> >>
> >> // Total number of operations to compute function values = 12
> >> for (unsigned int r = 0; r < 6; r++)
> >> {
> >> F0 += FE0[ip][r]*w[0][r];
> >> }// end loop over 'r'
> >>
> >> // Number of operations for primary indices: 5
> >> // Number of operations to compute entry: 5
> >> A[0] += F0*F0*F0*W3[ip]*det;
> >> }// end loop over 'ip'
> >> }
> >>
> >> As you see two loops over integration points, as one would expect.
> >> I do, however, receive the warning 'Quadrature degree must be equal
> >> within each sub domain, using degree 8.' from FFC,
> >> but I guess that's related to the tensor representation.
> >> And as Garth says, Dolfin knows nothing about the integration
> >> schemes which are used inside tabulate_tensor(), it just calls the
> >> function.
> >
> > Ok. It seems that
> > "Quadrature degree must be equal within each sub domain, using
> > degree n" warning does not apply for quadrature representation but
> > is displayed. Is it a bug?
> >
> > Is stdout of FFC accessible somehow when JIT compiling from DOLFIN?
> > Perhapse Johan could know?
> >
> > Jan
> >
> >
> > PS: Here is only simple demo of what can happen. I wrote it before I
> > realized that everything works quite satisfactory. But it could be
> > useful to somebody. Just skip it if not insterested. I checked codes
> > generated by FFC trunk for this code
> > _________________________________
> > U = FiniteElement('CG', 'triangle', 2)
> > V = FiniteElement('CG', 'triangle', 3)
> > u = Coefficient(U)
> > v = Coefficient(V)
> > aT = u**2*dx(0, {'representation': 'tensor'})
> > bT = v**4*dx(0, {'representation': 'tensor'})
> > aQ = u**2*dx(0, {'representation': 'quadrature'})
> > bQ = v**4*dx(0, {'representation': 'quadrature'})
> > aD = u**2*dx(0)
> > bD = v**4*dx(0)
> > ________________________________
> > and found that
> >
> > form | representaion
> > ----------------------
> > aT + bT | G
> > aQ + bQ | W6 + W49
> > aD + bD | W49
> > aT + bQ | W6 + W49
> > aQ + bT | W6 + W49
> > aQ + bD | W6 + W49
> > aD + bQ | W49
> > aT + bD | W6 + W49
> > aD + bT | G
> > ----------------------
> >
> > where Wnp is quadrature with np = number of points and G is tensor.
> > The moral is that one has to check the generated code if it does
> > what's intended. It is also somehow readable from FFC stdout
> > "Analyzing forms" step.
> >
> >>
> >>
> >> > It seems for me that only solution is this:
> >> > ________________________________________
> >> > from dolfin import *
> >> > mesh = UnitSquareMesh(2,2)
> >> > V = FunctionSpace(mesh, 'CG', 2)
> >> > u = Function(V)
> >> > a = u**3*dx(0, {'quadrature_order': 2})
> >> > b = u**4*dx(0, {'quadrature_order': 8})
> >> > M = assemble(a)
> >> > M += assemble(b) # for forms of rank
> >> > = 0 #M = assemble(b, tensor=M, add_values=True) # for forms of
> >> > rank > 0 ________________________________________
> >> > Am I right?
> >> >
> >> > Jan
> >> >
> >> >
> >> > > Den 10. apr. 2013 18:45 skrev "Jan Blechta" <
> >> > > <email address hidden>> følgende:
> >> > >
> >> > > > New question #226398 on DOLFIN:
> >> > > > https://answers.launchpad.net/dolfin/+question/226398
> >> > > >
> >> > > > Is it possible to JIT compile a form by FFC whose parts would
> >> > > > have different qudrature degrees? Consider for example
> >> > > >
> >> > > > u = Function(V)
> >> > > > v = TestFunction(V)
> >> > > > phi = Function(Q)
> >> > > > F = u*v*dx + phi*v*dx
> >> > > >
> >> > > > Is there the way to tell FFC to compile first term with some
> >> > > > (say 2 or automatic) degree and second with another (say 4)
> >> > > > quadrature degree? As far as I understandt it I would say
> >> > > > that compiled F has common degree corresponding to term with
> >> > > > the highest degree estimated by UFL, doesn't it?
> >> > > >
> >> > > > Or is only solution define forms part by part, assemble each
> >> > > > part separately and add tensors together?
> >> > > >
> >> > > > --
> >> > > > You received this question notification because you are a
> >> > > > member of DOLFIN Team, which is an answer contact for DOLFIN.
> >> > > >
> >> > >
> >> >
> >> > --
> >> > You received this question notification because you are a member
> >> > of DOLFIN Team, which is an answer contact for DOLFIN.
> >> >
> >>
> >
> > --
> > You received this question notification because you are a member of
> > DOLFIN Team, which is an answer contact for DOLFIN.
>

On 04/11/2013 04:11 PM, Jan Blechta wrote:
> Question #226398 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/226398
>
> Status: Answered => Open
>
> Jan Blechta is still having a problem:
> On Thu, 11 Apr 2013 06:41:39 -0000
> Kristian B. Ølgaard <email address hidden> wrote:
>> Your question #226398 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/226398
>>
>> Kristian B. Ølgaard proposed the following answer:
>> On 10 April 2013 20:41, Jan Blechta
>> <email address hidden>wrote:
>>
>>> Question #226398 on DOLFIN changed:
>>> https://answers.launchpad.net/dolfin/+question/226398
>>>
>>> Status: Answered => Open
>>>
>>> Jan Blechta is still having a problem:
>>> On Wed, 10 Apr 2013 17:31:07 -0000
>>> Martin Sandve Alnæs <email address hidden> wrote:
>>>> Your question #226398 on DOLFIN changed:
>>>> https://answers.launchpad.net/dolfin/+question/226398
>>>>
>>>> Status: Open => Answered
>>>>
>>>> Martin Sandve Alnæs proposed the following answer:
>>>> Yes, there is a ufl demo file in ffc showing how. Grep for quad.
>>>
>>> Ok, there is Poisson.ufl:
>>> _________________________________________________
>>> element = FiniteElement("Lagrange", triangle, 1)
>>>
>>> u = TrialFunction(element)
>>> v = TestFunction(element)
>>> f = Coefficient(element)
>>>
>>> dx0_1 = dx(0, {"quadrature_order": 1})
>>> dx0_2 = dx(0, {"quadrature_order": 2})
>>>
>>> a = dot(grad(v), grad(u))*dx0_1
>>> L = v*f*dx0_2
>>> _________________________________________________
>>>
>>> Is there a reason why form a has degree 1? I would expect that
>>> degree 0 suffices for exact quadrature as gradients of both u and v
>>> have polynomial degree 0.
>>>
>>>
>> Using degree 0 or 1 will result in a one point integration scheme, so
>> it doesn't matter.
>> Note that you should use 'quadrature_order' --> 'quadrature_degree'.
>>
>>
>>> It seems that this code
>>> ________________________________________
>>> from dolfin import *
>>> mesh = UnitSquareMesh(2,2)
>>> V = FunctionSpace(mesh, 'CG', 2)
>>> u = Function(V)
>>> a = u**3*dx(0, {'quadrature_order': 2})
>>> b = u**4*dx(0, {'quadrature_order': 8})
>>> M = assemble(a+b)
>>> ________________________________________
>>> compiles form a+b together with only one quadrature rule. Maybe I
>>> don't read its ufc code correctly.
>>>
>>>
>> Compiling the following form with FFC 1.0 and quadrature
>> representation:
>>
>> element = FiniteElement("Lagrange", triangle, 2)
>>
>> f = Coefficient(element)
>>
>> m0 = f**3*dx(0, {'quadrature_degree': 2})
>> m1 = f**4*dx(0, {'quadrature_degree': 8})
>> M = (m0+m1)
>>
>> results in the code:
>>
>>
>> [snip]
>>
>> for (unsigned int ip = 0; ip < 25; ip++)
>> {
>>
>> // Coefficient declarations.
>> double F0 = 0.0;
>>
>> // Total number of operations to compute function values = 12
>> for (unsigned int r = 0; r < 6; r++)
>> {
>> F0 += FE1[ip][r]*w[0][r];
>> }// end loop over 'r'
>>
>> // Number of operations for primary indices: 6
>> // Number of operations to compute entry: 6
>> A[0] += F0*F0*F0*F0*W25[ip]*det;
>> }// end loop over 'ip'
>>
>> // Loop quadrature points for integral.
>> // Number of operations to compute element tensor for following
>> IP loop = 51
>> for (unsigned int ip = 0; ip < 3; ip++)
>> {
>>
>> // Coefficient declarations.
>> double F0 = 0.0;
>>
>> // Total number of operations to compute function values = 12
>> for (unsigned int r = 0; r < 6; r++)
>> {
>> F0 += FE0[ip][r]*w[0][r];
>> }// end loop over 'r'
>>
>> // Number of operations for primary indices: 5
>> // Number of operations to compute entry: 5
>> A[0] += F0*F0*F0*W3[ip]*det;
>> }// end loop over 'ip'
>> }
>>
>> As you see two loops over integration points, as one would expect.
>> I do, however, receive the warning 'Quadrature degree must be equal
>> within each sub domain, using degree 8.' from FFC,
>> but I guess that's related to the tensor representation.
>> And as Garth says, Dolfin knows nothing about the integration schemes
>> which are used inside tabulate_tensor(), it just calls the function.
>
> Ok. It seems that
> "Quadrature degree must be equal within each sub domain, using degree n"
> warning does not apply for quadrature representation but is displayed.
> Is it a bug?
>
> Is stdout of FFC accessible somehow when JIT compiling from DOLFIN?
> Perhapse Johan could know?

Try:

parameters.form_compiler.log_level = INFO

or

parameters.form_compiler.log_level = DEBUG

Johan

>
> Jan
>
>
> PS: Here is only simple demo of what can happen. I wrote it before I
> realized that everything works quite satisfactory. But it could be
> useful to somebody. Just skip it if not insterested. I checked codes
> generated by FFC trunk for this code
> _________________________________
> U = FiniteElement('CG', 'triangle', 2)
> V = FiniteElement('CG', 'triangle', 3)
> u = Coefficient(U)
> v = Coefficient(V)
> aT = u**2*dx(0, {'representation': 'tensor'})
> bT = v**4*dx(0, {'representation': 'tensor'})
> aQ = u**2*dx(0, {'representation': 'quadrature'})
> bQ = v**4*dx(0, {'representation': 'quadrature'})
> aD = u**2*dx(0)
> bD = v**4*dx(0)
> ________________________________
> and found that
>
> form | representaion
> ----------------------
> aT + bT | G
> aQ + bQ | W6 + W49
> aD + bD | W49
> aT + bQ | W6 + W49
> aQ + bT | W6 + W49
> aQ + bD | W6 + W49
> aD + bQ | W49
> aT + bD | W6 + W49
> aD + bT | G
> ----------------------
>
> where Wnp is quadrature with np = number of points and G is tensor. The
> moral is that one has to check the generated code if it does what's
> intended. It is also somehow readable from FFC stdout "Analyzing forms"
> step.
>
>>
>>
>>> It seems for me that only solution is this:
>>> ________________________________________
>>> from dolfin import *
>>> mesh = UnitSquareMesh(2,2)
>>> V = FunctionSpace(mesh, 'CG', 2)
>>> u = Function(V)
>>> a = u**3*dx(0, {'quadrature_order': 2})
>>> b = u**4*dx(0, {'quadrature_order': 8})
>>> M = assemble(a)
>>> M += assemble(b) # for forms of rank = 0
>>> #M = assemble(b, tensor=M, add_values=True) # for forms of rank > 0
>>> ________________________________________
>>> Am I right?
>>>
>>> Jan
>>>
>>>
>>>> Den 10. apr. 2013 18:45 skrev "Jan Blechta" <
>>>> <email address hidden>> følgende:
>>>>
>>>>> New question #226398 on DOLFIN:
>>>>> https://answers.launchpad.net/dolfin/+question/226398
>>>>>
>>>>> Is it possible to JIT compile a form by FFC whose parts would
>>>>> have different qudrature degrees? Consider for example
>>>>>
>>>>> u = Function(V)
>>>>> v = TestFunction(V)
>>>>> phi = Function(Q)
>>>>> F = u*v*dx + phi*v*dx
>>>>>
>>>>> Is there the way to tell FFC to compile first term with some
>>>>> (say 2 or automatic) degree and second with another (say 4)
>>>>> quadrature degree? As far as I understandt it I would say that
>>>>> compiled F has common degree corresponding to term with the
>>>>> highest degree estimated by UFL, doesn't it?
>>>>>
>>>>> Or is only solution define forms part by part, assemble each
>>>>> part separately and add tensors together?
>>>>>
>>>>> --
>>>>> You received this question notification because you are a
>>>>> member of DOLFIN Team, which is an answer contact for DOLFIN.
>>>>>
>>>>
>>>
>>> --
>>> You received this question notification because you are a member of
>>> DOLFIN Team, which is an answer contact for DOLFIN.
>>>
>>
>

Thanks Martin Sandve Alnæs, that solved my question.