Subject: [FFC-dev] defining sudbomain integrals

From: Jake Ostien <[email-masked]>
Date: Wed, 30 Jan 2008 21:23:42 -0500

Hi,

Something doesn't seem right when I try and define multiple integrals 
over subdomains in a form file.

See the attached test.form and test.txt.  Basically I copied the element 
tensor section of the generated header file from two different 
compilations of test.form, one where

a = v*u*dx  (for a reference)

and one with

dx0 = Integral("cell", 0)
dx1 = Integral("cell", 1)

a = v*u*dx0 + 10*v*u*dx1

Now I would expect that the cell_integral_0 from the second compilation 
to be identical to the reference compilation, and that cell_integral_1 
would be just scaled by 10, but that is not what is happening.  It 
appears that in the second compilation the two subdomain integrals are 
just added together (see test.txt).

I didn't get too much out of running ffc with the -d1 flag, but I really 
don't know what I am looking for.  Any help on where to look?

Thanks,
Jake

---

From: Jake Ostien <[email-masked]>
Date: Wed, 30 Jan 2008 21:26:18 -0500

And for anyone paying attention, here are the missing files.

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

element = FiniteElement("CG", "tetrahedron", 1)

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

#a = v*u*dx
a = v*u*dx0 + 10.0*v*u*dx1
reference
a = v*u*dx

    // Compute element tensor
    A[0] = 0.0166666666666666*G0_;
    A[1] = 0.00833333333333331*G0_;
    A[2] = 0.00833333333333331*G0_;
    A[3] = 0.00833333333333331*G0_;
    A[4] = 0.00833333333333331*G0_;
    A[5] = 0.0166666666666666*G0_;
    A[6] = 0.00833333333333331*G0_;
    A[7] = 0.00833333333333331*G0_;
    A[8] = 0.00833333333333331*G0_;
    A[9] = 0.00833333333333331*G0_;
    A[10] = 0.0166666666666666*G0_;
    A[11] = 0.00833333333333331*G0_;
    A[12] = 0.00833333333333331*G0_;
    A[13] = 0.00833333333333331*G0_;
    A[14] = 0.00833333333333331*G0_;
    A[15] = 0.0166666666666666*G0_;

Using subdomains
a = v*u*dx0 + 10.0*v*u*dx1
cell_integral_0

    // Compute element tensor
    A[0] = 0.0166666666666666*G0_ + 0.166666666666666*G1_;
    A[1] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[2] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[3] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[4] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[5] = 0.0166666666666666*G0_ + 0.166666666666666*G1_;
    A[6] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[7] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[8] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[9] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[10] = 0.0166666666666666*G0_ + 0.166666666666666*G1_;
    A[11] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[12] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[13] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[14] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[15] = 0.0166666666666666*G0_ + 0.166666666666666*G1_;

cell_integral_1

    // Compute element tensor
    A[0] = 0.0166666666666666*G0_ + 0.166666666666666*G1_;
    A[1] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[2] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[3] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[4] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[5] = 0.0166666666666666*G0_ + 0.166666666666666*G1_;
    A[6] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[7] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[8] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[9] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[10] = 0.0166666666666666*G0_ + 0.166666666666666*G1_;
    A[11] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[12] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[13] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[14] = 0.00833333333333331*G0_ + 0.0833333333333331*G1_;
    A[15] = 0.0166666666666666*G0_ + 0.166666666666666*G1_;

---

From: Kristian Oelgaard <[email-masked]>
Date: Thu, 31 Jan 2008 09:40:10 +0100

> > don't know what I am looking for.  Any help on where to look?

I'm not quite sure about what is going on, or what should be going on. It does,
however, look a bit strange that cell_integral_0 also contains the dx1 part of
the form.

The values look OK to me so it looks like a bug in determining which monomial
terms go to which integral? I don't really know where this takes place.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Thu, 31 Jan 2008 11:27:36 +0100

> don't know what I am looking for.  Any help on where to look?

Looks like an obvious bug that should be easy to fix. I will take a
look.

-- 
Anders

---

From: Anders Logg <[email-masked]>
Date: Thu, 31 Jan 2008 12:09:24 +0100

> Looks like an obvious bug that should be easy to fix. I will take a

ok, try again and see if it works now. It turned out that the
integrals for the different sub domains were computed correctly but
then completely ignored when generating code.

-- 
Anders

---

From: Jake Ostien <[email-masked]>
Date: Thu, 31 Jan 2008 08:56:34 -0500

> then completely ignored when generating code.
Looks good.

Jake