Subject: Re: [FFC-dev] [UFL-dev] [HG UFL] Added check for restriction on facet normals and form arguments

From: Harish Narayanan <[email-masked]>
Date: Mon, 08 Jun 2009 15:10:58 +0200

UFL wrote:
> One or more new changesets pushed to the primary ufl repository.
> A short summary of the last three changesets is included below.
> 
> changeset:   897:7c6c19b55ba45ffbec4bfffc68ae0dce9636d20d
> tag:         tip
> user:        "Martin Sandve AlnÃ¦s <[email-masked]>"
> date:        Mon Jun 08 10:29:42 2009 +0200
> files:       sandbox/scratch/restrictiontest.py ufl/algorithms/checks.py ufl/algorithms/propagate_restrictions.py
> description:
> Added check for restriction on facet normals and form arguments
> inside the propagate_restrictions algorithm (quick and easy fix).

With this changeset, I now see:

*** FFC: Form argument must be restricted.
*** FFC: To get more information about this error, rerun FFC with --debug.

for any form I try to compile. Must something be added to FFC for it to
work with this UFL change?

Harish

---

From: Anders Logg <[email-masked]>
Date: Mon, 8 Jun 2009 15:19:51 +0200

> work with this UFL change?

What do you get with --debug? It should then print out the traceback
for the UFL exception.

-- 
Anders

---

From: Harish Narayanan <[email-masked]>
Date: Mon, 08 Jun 2009 15:23:58 +0200

> What do you get with --debug? It should then print out the traceback

$ ffc --debug -l dolfin Poisson.ufl
This is FFC, the FEniCS Form Compiler, version 0.6.2.
For further information, visit http://www.fenics.org/ffc/.

Preprocessing form file: Poisson.ufl --> Poisson.py

Compiler stage 1: Analyzing form
--------------------------------
This form is linear in 2 arguments.
  Name:                               a
  Rank:                               2
  Cell:                               <triangle of degree 1>
  Geometric dimension:                2
  Topological dimension:              2
  Number of functions:                0
  Number of cell integrals:           1
  Number of exterior facet integrals: 0
  Number of interior facet integrals: 0
  Basis functions:                    [v_0, v_1]
  Functions:                          []
  Basis function names:               [v0, v1]
  Function names:                     []
  Unique elements:                    CG1
  Unique sub elements:                CG1

  Automatic selection of representation not implemented, defaulting to
quadrature.
  Found element in element cache: <CG1 on a <triangle of degree 1>>
  Found element in element cache: <CG1 on a <triangle of degree 1>>
  Found element in dof map cache: <CG1 on a <triangle of degree 1>>
  Found element in element cache: <CG1 on a <triangle of degree 1>>
  Found element in element cache: <CG1 on a <triangle of degree 1>>
  Found element in element cache: <CG1 on a <triangle of degree 1>>
  Found element in element cache: <CG1 on a <triangle of degree 1>>

Compiler stage 2: Computing form representation(s)
--------------------------------------------------

  QR, init, form:
  { sum_{i_1} (({ A | A_{i_0} = dv_0/dx_i_0 })[i_1]) * (({ A | A_{i_2} =
dv_1/dx_i_2 })[i_1])  } * dx0
  Computing quadrature representation
  Found element in element cache: <CG1 on a <triangle of degree 1>>
  Derivatives:
set([SpatialDerivative(BasisFunction(FiniteElement('Lagrange',
Cell('triangle', 1), 1), 0), MultiIndex((Index(0),), {Index(0): 2})),
SpatialDerivative(BasisFunction(FiniteElement('Lagrange',
Cell('triangle', 1), 1), 1), MultiIndex((Index(2),), {Index(2): 2}))])
  num_derivatives: {FiniteElement('Lagrange', Cell('triangle', 1), 1): 1}

  QR, init, psi_tables:
  {'cell': {1: {FiniteElement('Lagrange', Cell('triangle', 1), 1):
{None: [{(0, 0): array([[ 0.33333333],
         [ 0.33333333],
         [ 0.33333333]])}, {(0, 1): array([[ -1.00000000e+00],
         [  7.85046229e-17],
         [  1.00000000e+00]]), (1, 0): array([[-1.],
         [ 1.],
         [ 0.]])}]}}}, 'exterior_facet': {}, 'interior_facet': {}}

  QR, init, quadrature_weights:
  {'cell': {1: array([ 0.5])}, 'exterior_facet': {}, 'interior_facet': {}}

Compiler stage 3: Optimizing form representation
------------------------------------------------
  Optimization of tensor contraction representation currently broken (to
be fixed).

Compiler stage 4: Generating code
---------------------------------
  Generating code for finite elements...
  Removing unused variable: Jinv_11
  Removing unused variable: Jinv_10
  Removing unused variable: Jinv_01
  Removing unused variable: Jinv_00
  Removing unused variable: Jinv_11
  Removing unused variable: Jinv_10
  Removing unused variable: Jinv_01
  Removing unused variable: Jinv_00
  Removing unused variable: Jinv_11
  Removing unused variable: Jinv_10
  Removing unused variable: Jinv_01
  Removing unused variable: Jinv_00
  Removing unused variable: Jinv_11
  Removing unused variable: Jinv_10
  Removing unused variable: Jinv_01
  Removing unused variable: Jinv_00
  done
  Generating code for finite dof maps...
  done
  Generating code for form...
  done
  Generating code using quadrature representation

  QG-utils, psi_tables:
  {1: {FiniteElement('Lagrange', Cell('triangle', 1), 1): {None: [{(0,
0): array([[ 0.33333333],
         [ 0.33333333],
         [ 0.33333333]])}, {(0, 1): array([[ -1.00000000e+00],
         [  7.85046229e-17],
         [  1.00000000e+00]]), (1, 0): array([[-1.],
         [ 1.],
         [ 0.]])}]}}}

  QG-utils, flatten_tables, points:
  1

  QG-utils, flatten_tables, elem_dict:
  {FiniteElement('Lagrange', Cell('triangle', 1), 1): {None: [{(0, 0):
array([[ 0.33333333],
         [ 0.33333333],
         [ 0.33333333]])}, {(0, 1): array([[ -1.00000000e+00],
         [  7.85046229e-17],
         [  1.00000000e+00]]), (1, 0): array([[-1.],
         [ 1.],
         [ 0.]])}]}}

  QG-utils, flatten_tables, elem:
  <CG1 on a <triangle of degree 1>>

  QG-utils, flatten_tables, facet_tables:
  {None: [{(0, 0): array([[ 0.33333333],
         [ 0.33333333],
         [ 0.33333333]])}, {(0, 1): array([[ -1.00000000e+00],
         [  7.85046229e-17],
         [  1.00000000e+00]]), (1, 0): array([[-1.],
         [ 1.],
         [ 0.]])}]}

  QG-utils, flatten_tables, derivs:
  (0, 0)

  QG-utils, flatten_tables, psi_table:
  [[ 0.33333333]
   [ 0.33333333]
   [ 0.33333333]]
  Table name: FE0

  QG-utils, flatten_tables, derivs:
  (0, 1)

  QG-utils, flatten_tables, psi_table:
  [[ -1.00000000e+00]
   [  7.85046229e-17]
   [  1.00000000e+00]]
  Table name: FE0_D01

  QG-utils, flatten_tables, derivs:
  (1, 0)

  QG-utils, flatten_tables, psi_table:
  [[-1.]
   [ 1.]
   [ 0.]]
  Table name: FE0_D10

  QG-utils, psi_tables, flat_tables:
  {'FE0_D10': array([[-1.,  1.,  0.]]), 'FE0_D01': array([[
-1.00000000e+00,   7.85046229e-17,   1.00000000e+00]]), 'FE0': array([[
0.33333333,  0.33333333,  0.33333333]])}

  tables: {'FE0_D10': array([[-1.,  1.,  0.]]), 'FE0_D01': array([[
-1.00000000e+00,   7.85046229e-17,   1.00000000e+00]]), 'FE0': array([[
0.33333333,  0.33333333,  0.33333333]])}

  name_map: {}

  inv_name_map: {'FE0_D10': 'FE0_D10', 'FE0_D01': 'FE0_D01', 'FE0': 'FE0'}

  QG-utils, psi_tables, unique_tables:
  {'FE0_D10': array([[-1.,  1.,  0.]]), 'FE0_D01': array([[-1.,  0.,
1.]]), 'FE0': array([[ 0.33333333,  0.33333333,  0.33333333]])}

  QG-utils, psi_tables, name_map:
  {'FE0_D10': ('FE0_D10', (), False, False), 'FE0_D01': ('FE0_D01', (),
False, False), 'FE0': ('FE0', (), False, False)}
  Generating code for cell integrals using quadrature representation.

  QG, cell_integral, integrals:
  {1:
Integral(IndexSum(Product(Indexed(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange',
Cell('triangle', 1), 1), 0), MultiIndex((Index(0),), {Index(0): 2})),
MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),),
{Index(1): 2})),
Indexed(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange',
Cell('triangle', 1), 1), 1), MultiIndex((Index(2),), {Index(2): 2})),
MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),),
{Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2})),
Measure('cell', 0, {}))}
  Looping points: 1
  integral: { sum_{i_1} (({ A | A_{i_0} = dv_0/dx_i_0 })[i_1]) * (({ A |
A_{i_2} = dv_1/dx_i_2 })[i_1])  } * dx0<{}>

  Integral tree_format: Integral:
      domain type: cell
      domain id: 0
      integrand:
          IndexSum
          (
              Product
              (
                  Indexed
                  (
                      ComponentTensor
                      (
                          SpatialDerivative
                          (
                              BasisFunction(FiniteElement('Lagrange',
Cell('triangle', 1), 1), 0)
                              MultiIndex((Index(0),), {Index(0): 2})
                          )
                          MultiIndex((Index(0),), {Index(0): 2})
                      )
                      MultiIndex((Index(1),), {Index(1): 2})
                  )
                  Indexed
                  (
                      ComponentTensor
                      (
                          SpatialDerivative
                          (
                              BasisFunction(FiniteElement('Lagrange',
Cell('triangle', 1), 1), 1)
                              MultiIndex((Index(2),), {Index(2): 2})
                          )
                          MultiIndex((Index(2),), {Index(2): 2})
                      )
                      MultiIndex((Index(1),), {Index(1): 2})
                  )
              )
              MultiIndex((Index(1),), {Index(1): 2})
          )

  QG, Using Transformer
Form argument must be restricted.
Traceback (most recent call last):
  File "/Users/harish/Work/FEniCS/unstable/build/bin/ffc", line 186, in
<module>
    sys.exit(main(sys.argv[1:]))
  File "/Users/harish/Work/FEniCS/unstable/build/bin/ffc", line 127, in main
    execfile(script, {})
  File "Poisson.py", line 34, in <module>
    compile([a, L, M, element], "Poisson", {'log_level': 10, 'format':
'dolfin', 'form_postfix': True, 'quadrature_order': 'auto', 'precision':
'15', 'cpp optimize': False, 'cache_dir': None, 'split': False,
'representation': 'auto', 'optimize': False, 'output_dir': '.'}, globals())
  File
"/Users/harish/Work/FEniCS/unstable/build/lib/python2.5/site-packages/ffc/compiler/compiler.py",
line 99, in compile
    form_code = generate_form_code(form_data, representations, prefix,
format.format, options)
  File
"/Users/harish/Work/FEniCS/unstable/build/lib/python2.5/site-packages/ffc/compiler/compiler.py",
line 191, in generate_form_code
    codes.append(code_generator.generate_integrals(representations[i],
format))
  File
"/unstable/lib/python2.5/site-packages/ffc/compiler/codegeneration/quadrature/quadraturegenerator.py",
line 74, in generate_integrals
    code.update(self.generate_cell_integrals(form_representation, format))
  File
"/unstable/lib/python2.5/site-packages/ffc/compiler/codegeneration/quadrature/quadraturegenerator.py",
line 102, in generate_cell_integrals
    self.generate_cell_integral(form_representation, transformer,
integrals, format)
  File
"/unstable/lib/python2.5/site-packages/ffc/compiler/codegeneration/quadrature/quadraturegenerator.py",
line 165, in generate_cell_integral
    integrals, Indent, format)
  File
"/unstable/lib/python2.5/site-packages/ffc/compiler/codegeneration/quadrature/quadraturegenerator.py",
line 352, in __generate_element_tensor
    generate_code(integral.integrand(), transformer, Indent, format)
  File
"/unstable/lib/python2.5/site-packages/ffc/compiler/codegeneration/quadrature/quadraturetransformer.py",
line 1132, in generate_code
    new_integrand = propagate_restrictions(new_integrand)
  File
"/unstable/lib/python2.5/site-packages/ufl/algorithms/propagate_restrictions.py",
line 38, in propagate_restrictions
    return RestrictionPropagator().visit(expression)
  File
"/unstable/lib/python2.5/site-packages/ufl/algorithms/transformations.py",
line 128, in visit
    r = h(o, *map(self.visit, o.operands()))
  File
"/unstable/lib/python2.5/site-packages/ufl/algorithms/transformations.py",
line 128, in visit
    r = h(o, *map(self.visit, o.operands()))
  File
"/unstable/lib/python2.5/site-packages/ufl/algorithms/transformations.py",
line 128, in visit
    r = h(o, *map(self.visit, o.operands()))
  File
"/unstable/lib/python2.5/site-packages/ufl/algorithms/transformations.py",
line 132, in visit
    r = h(o)
  File
"/unstable/lib/python2.5/site-packages/ufl/algorithms/propagate_restrictions.py",
line 31, in form_argument
    ufl_assert(self.current_restriction is not None, "Form argument must
be restricted.")
  File "/unstable/lib/python2.5/site-packages/ufl/assertions.py", line
20, in ufl_assert
    if not condition: error(*message)
  File "/unstable/lib/python2.5/site-packages/ufl/log.py", line 97, in error
    raise UFLException(self._format_raw(*message))
ufl.log.UFLException: Form argument must be restricted.

Harish

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 08 Jun 2009 16:20:17 +0200

>     new_integrand = propagate_restrictions(new_integrand)

I guess the error occurs here. I apply a series of UFL algorithms:

new_integrand = expand_indices(integrand)
new_integrand = purge_list_tensors(new_integrand)
new_integrand = propagate_restrictions(new_integrand)

before starting the transformation.

Should this not be legal/work?

Kristian

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

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 8 Jun 2009 17:29:55 +0200

> Should this not be legal/work?

I figured propagate_restrictions should only be applied
to the integrands of interior facet integrals.

Martin

---

From: Anders Logg <[email-masked]>
Date: Mon, 8 Jun 2009 17:45:39 +0200

> >> Compiler stage 1: Analyzing form
5F> >>   Cell:                               <triangle of degree 1>
> to the integrands of interior facet integrals.

Is propagate_restrictions performing the test?

In that case, it would be natural to have it check that terms on
interior facets *are* restricted and terms not on interior facet are
*not* restricted.

Also, would it be a good idea to wrap expand_indices,
purge_list_tensors and propagate_restrictions into a common function
to be applied for massaging a form before using it?

-- 
Anders

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 8 Jun 2009 17:54:09 +0200

> to be applied for massaging a form before using it?

No. But propagate_restrictions can probably be merged with
renumber_indices and some other stuff. However, it's best to
wait until things are more streamlined before deciding that.

(purge_list_tensors is currently implemented by calling expand_indices,
which does much more than purge_list_tensors should do)

Martin