# different quadature degree

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?

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Martin Sandve Alnæs

- Solved:
- 2013-04-16

- Last query:
- 2013-04-16

- Last reply:
- 2013-04-11

Martin Sandve Alnæs (martinal) said : | #1 |

Yes, there is a ufl demo file in ffc showing how. Grep for quad.

Den 10. apr. 2013 18:45 skrev "Jan Blechta" <

<email address hidden>> følgende:

> New question #226398 on DOLFIN:

> https:/

>

> 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.

>

Jan Blechta (blechta) said : | #2 |

On Wed, 10 Apr 2013 17:31:07 -0000

Martin Sandve Alnæs <email address hidden> wrote:

> Your question #226398 on DOLFIN changed:

> https:/

>

> 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(

u = TrialFunction(

v = TestFunction(

f = Coefficient(

dx0_1 = dx(0, {"quadrature_

dx0_2 = dx(0, {"quadrature_

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_

b = u**4*dx(0, {'quadrature_

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_

b = u**4*dx(0, {'quadrature_

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:/

> >

> > 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.

> >

>

Martin Sandve Alnæs (martinal) said : | #3 |

Is it using tensor representation? If so try adding "representation

Jan Blechta (blechta) said : | #4 |

On Wed, 10 Apr 2013 19:16:15 -0000

Martin Sandve Alnæs <email address hidden> wrote:

> Your question #226398 on DOLFIN changed:

> https:/

>

> Status: Open => Answered

>

> Martin Sandve Alnæs proposed the following answer:

> Is it using tensor representation? If so try adding

> "representation

> 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

Garth Wells (garth-wells) said : | #5 |

On 11 April 2013 07:26, Jan Blechta

<email address hidden> wrote:

> Question #226398 on DOLFIN changed:

> https:/

>

> 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:/

>>

>> Status: Open => Answered

>>

>> Martin Sandve Alnæs proposed the following answer:

>> Is it using tensor representation? If so try adding

>> "representation

>> 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_

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:/

>

> 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:/

> >

> > 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(

>

> u = TrialFunction(

> v = TestFunction(

> f = Coefficient(

>

> dx0_1 = dx(0, {"quadrature_

> dx0_2 = dx(0, {"quadrature_

>

> 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_

> 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_

> b = u**4*dx(0, {'quadrature_

> 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(

f = Coefficient(

m0 = f**3*dx(0, {'quadrature_

m1 = f**4*dx(0, {'quadrature_

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*

}// 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*

}// 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_

> b = u**4*dx(0, {'quadrature_

> 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:/

> > >

> > > 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.

>

Jan Blechta (blechta) said : | #7 |

On Thu, 11 Apr 2013 06:41:39 -0000

Kristian B. Ølgaard <email address hidden> wrote:

> Your question #226398 on DOLFIN changed:

> https:/

>

> 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:/

> >

> > 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:/

> > >

> > > 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(

> >

> > u = TrialFunction(

> > v = TestFunction(

> > f = Coefficient(

> >

> > dx0_1 = dx(0, {"quadrature_

> > dx0_2 = dx(0, {"quadrature_

> >

> > 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_

>

>

> > 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_

> > b = u**4*dx(0, {'quadrature_

> > 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(

>

> f = Coefficient(

>

> m0 = f**3*dx(0, {'quadrature_

> m1 = f**4*dx(0, {'quadrature_

> 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*

> }// 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*

> }// 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_

> > b = u**4*dx(0, {'quadrature_

> > 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:/

> > > >

> > > > 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.

> >

>

Martin Sandve Alnæs (martinal) said : | #8 |

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:/

>

> 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:/

>>

>> 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:/

>> >

>> > 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:/

>> > >

>> > > 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(

>> >

>> > u = TrialFunction(

>> > v = TestFunction(

>> > f = Coefficient(

>> >

>> > dx0_1 = dx(0, {"quadrature_

>> > dx0_2 = dx(0, {"quadrature_

>> >

>> > 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_

>>

>>

>> > 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_

>> > b = u**4*dx(0, {'quadrature_

>> > 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(

>>

>> f = Coefficient(

>>

>> m0 = f**3*dx(0, {'quadrature_

>> m1 = f**4*dx(0, {'quadrature_

>> 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*

>> }// 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*

>> }// 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_

>> > b = u**4*dx(0, {'quadrature_

>> > 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:/

>> > > >

>> > > > 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.

Jan Blechta (blechta) said : | #9 |

On Thu, 11 Apr 2013 14:56:15 -0000

Martin Sandve Alnæs <email address hidden> wrote:

> Your question #226398 on DOLFIN changed:

> https:/

>

> 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:/

> >

> > 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:/

> >>

> >> 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:/

> >> >

> >> > 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:/

> >> > >

> >> > > 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(

> >> >

> >> > u = TrialFunction(

> >> > v = TestFunction(

> >> > f = Coefficient(

> >> >

> >> > dx0_1 = dx(0, {"quadrature_

> >> > dx0_2 = dx(0, {"quadrature_

> >> >

> >> > 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_

> >>

> >>

> >> > 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_

> >> > b = u**4*dx(0, {'quadrature_

> >> > 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(

> >>

> >> f = Coefficient(

> >>

> >> m0 = f**3*dx(0, {'quadrature_

> >> m1 = f**4*dx(0, {'quadrature_

> >> 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*

> >> }// 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*

> >> }// 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_

> >> > b = u**4*dx(0, {'quadrature_

> >> > 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:/

> >> > > >

> >> > > > 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.

>

Johan Hake (johan-hake) said : | #10 |

On 04/11/2013 04:11 PM, Jan Blechta wrote:

> Question #226398 on DOLFIN changed:

> https:/

>

> 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:/

>>

>> 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:/

>>>

>>> 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:/

>>>>

>>>> 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(

>>>

>>> u = TrialFunction(

>>> v = TestFunction(

>>> f = Coefficient(

>>>

>>> dx0_1 = dx(0, {"quadrature_

>>> dx0_2 = dx(0, {"quadrature_

>>>

>>> 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_

>>

>>

>>> 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_

>>> b = u**4*dx(0, {'quadrature_

>>> 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(

>>

>> f = Coefficient(

>>

>> m0 = f**3*dx(0, {'quadrature_

>> m1 = f**4*dx(0, {'quadrature_

>> 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*

>> }// 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*

>> }// 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.

or

parameters.

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_

>>> b = u**4*dx(0, {'quadrature_

>>> 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:/

>>>>>

>>>>> 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.

>>>

>>

>

Jan Blechta (blechta) said : | #11 |

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