Subject: Re: [Ffc] [Branch ~ffc-core/ffc/dev] Rev 1429: Finished support for 3D Lagrange in evaluate_basis.

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 8 Jan 2010 17:35:17 +0100 (CET)

2010/1/8  <[email-masked]>:
> ------------------------------------------------------------
> revno: 1429
> committer: Kristian B. Ølgaard <[email-masked]>
> branch nick: dev
> timestamp: Fri 2010-01-08 16:54:41 +0100
> message:
>  Finished support for 3D Lagrange in evaluate_basis.
> modified:
>  ffc/evaluatebasis.py

2D and 3D works now for Lagrange elements, 1D is not implemented because of a bug in FIAT.
Before fixing BDM, Nedelec and MixedElements we should fix the following bugs:

In FIAT.FiniteElement: set the correct mapping (self._mapping) for all elements

in UFL (for MixedElements):
Traceback (most recent call last):
  File "/home/oelgaard/fenics/branches/dev/scripts/ffc", line 167, in <module>
    sys.exit(main(sys.argv[1:]))
  File "/home/oelgaard/fenics/branches/dev/scripts/ffc", line 155, in main
    compile_form(ufd.forms, prefix, options)
  File "/home/oelgaard/fenics/branches/dev/ffc/compiler.py", line 133, in compile_form
    ir = compute_ir(preprocessed_form, form_data, options)
  File "/home/oelgaard/fenics/branches/dev/ffc/representation.py", line 52, in compute_ir
    ir_integrals = compute_integrals_ir(form, form_data, options)
  File "/home/oelgaard/fenics/branches/dev/ffc/representation.py", line 132, in compute_integrals_ir
    ir = TensorRepresentation(form, form_data)
  File "/home/oelgaard/fenics/branches/dev/ffc/tensor/tensorrepresentation.py", line 80, in __init__
    transform_monomial_form(monomial_form)
  File "/home/oelgaard/fenics/branches/dev/ffc/tensor/monomialtransformation.py", line 36, in transform_monomial_form
    integrand.monomials[i] = TransformedMonomial(monomial)
  File "/home/oelgaard/fenics/branches/dev/ffc/tensor/monomialtransformation.py", line 278, in __init__
    same = ufl_element.component_element(component.index_range[0])[0].mapping()
AttributeError: 'VectorElement' object has no attribute 'component_element'

I will start on evaluate_basis_derivatives for Lagrange elements.

Kristian

>
> --
> lp:~ffc-core/ffc/dev
> https://code.launchpad.net/~ffc-core/ffc/dev
>
> You are subscribed to branch lp:~ffc-core/ffc/dev.
> To unsubscribe from this branch go to https://code.launchpad.net/~ffc-core/ffc/dev/+edit-subscription.
>
> === modified file 'ffc/evaluatebasis.py'
> --- ffc/evaluatebasis.py        2010-01-08 11:48:02 +0000
> +++ ffc/evaluatebasis.py        2010-01-08 15:54:41 +0000
> @@ -607,16 +607,30 @@
>
>  #    return code + [""]
>
> -# FIAT_NEW code (compute index function)
> +# FIAT_NEW code (compute index function) TriangleExpansionSet
>  # def idx(p,q):
>  #    return (p+q)*(p+q+1)/2 + q
> -def _idx(p, q):
> +def _idx2D(p, q):
>     f1 = create_float(1)
>     f2 = create_float(2)
>     idx = create_fraction(create_product([(p+q).expand(), (p+q+f1).expand()]), f2) + q
>     return idx
>
> -# FIAT_NEW code (helper variables)
> +# FIAT_NEW code (compute index function) TetrahedronExpansionSet
> +# def idx(p,q,r):
> +#     return (p+q+r)*(p+q+r+1)*(p+q+r+2)/6 + (q+r)*(q+r+1)/2 + r
> +def _idx3D(p, q, r):
> +    f1 = create_float(1)
> +    f2 = create_float(2)
> +    f6 = create_float(6)
> +    fac1 = create_fraction( (p + q + r + f2).expand(), f6)
> +    fac2 = create_fraction( (q + r + f1).expand(), f2)
> +    fac3 = create_product([(p + q + r).expand(), (p + q + r + f1).expand(), fac1])
> +    fac4 = create_product([(q + r).expand(), (q + r + f1).expand(), fac2])
> +    idx = fac3 + fac4 + r
> +    return idx
> +
> +# FIAT_NEW code (helper variables) TriangleExpansionSet and TetrahedronExpansionSet
>  # def jrc( a , b , n ):
>  #    an = float( ( 2*n+1+a+b)*(2*n+2+a+b)) \
>  #        / float( 2*(n+1)*(n+1+a+b))
> @@ -644,7 +658,6 @@
>
>     return (an, bn, cn)
>
> -
>  def _compute_basisvalues(data, Indent, format):
>     """From FIAT_NEW.expansions."""
>
> @@ -667,17 +680,19 @@
>     format_free_indices = format["free secondary indices"]
>     format_r            = format_free_indices[0]
>     format_s            = format_free_indices[1]
> +    format_t            = format_free_indices[2]
>
>     idx0, idx1, idx2    = [format["evaluate_basis aux index"](i) for i in range(1,4)]
> -    f1, f2, f3          = [create_symbol(format["evaluate_basis aux factor"](i), CONST) for i in range(1,4)]
> +    f1, f2, f3, f4, f5  = [create_symbol(format["evaluate_basis aux factor"](i), CONST) for i in range(1,6)]
>     an, bn, cn          = [create_symbol(format["evaluate_basis aux value"](i), CONST) for i in range(3)]
>
>     # Get embedded degree
>     embedded_degree = data["embedded_degree"]
>
>     # Create helper symbols
> -    symbol_r = create_symbol(format_r, CONST)
> -    symbol_s = create_symbol(format_s, CONST)
> +    symbol_p = create_symbol(format_r, CONST)
> +    symbol_q = create_symbol(format_s, CONST)
> +    symbol_r = create_symbol(format_t, CONST)
>     symbol_x = create_symbol(format_x, CONST)
>     symbol_y = create_symbol(format_y, CONST)
>     symbol_z = create_symbol(format_z, CONST)
> @@ -690,7 +705,10 @@
>     float_3 = create_float(3)
>     float_n = create_float(embedded_degree)
>     float_n1 = create_float(embedded_degree + 1)
> +    float_nm1 = create_float(embedded_degree - 1)
> +    float_1_5 = create_float(1.5)
>     float_0_5 = create_float(0.5)
> +    float_0_25 = create_float(0.25)
>
>     # Init return code
>     code = []
> @@ -704,7 +722,7 @@
>     code += [(name, value)]
>
>     # Declare helper variables, will be removed if not used.
> -    code += [Indent.indent(format["comment"]("Declare helper variables"))]
> +    code += ["", Indent.indent(format["comment"]("Declare helper variables"))]
>     code += [(format_uint + idx0, 0)]
>     code += [(format_uint + idx1, 0)]
>     code += [(format_uint + idx2, 0)]
> @@ -719,6 +737,7 @@
>     # 1D
>     if (element_cell_domain == "interval"):
>         count = 0
> +        error("1D is not implemented")
>  #        for i in range(0, data.degree() + 1):
>
>  #            factor = math.sqrt(1.0*i + 0.5)
> @@ -738,28 +757,28 @@
>     # 2D
>     elif (element_cell_domain == "triangle"):
>         # FIAT_NEW.expansions.TriangleExpansionSet
> -        # FIXME: KBO: Move common stuff to general functions
> +
> +        # Compute helper factors
> +        # FIAT_NEW code
> +        # f1 = (1.0+2*x+y)/2.0
> +        # f2 = (1.0 - y) / 2.0
> +        # f3 = f2**2
> +        fac1 = create_fraction(float_1 + float_2*symbol_x + symbol_y, float_2)
> +        fac2 = create_fraction(float_1 - symbol_y, float_2)
> +        code += [(format_float_decl + str(f1), fac1)]
> +        code += [(format_float_decl + str(f2), fac2)]
> +        code += [(format_float_decl + str(f3), f2*f2)]
> +
> +        code += ["", Indent.indent(format["comment"]("Compute basisvalues"))]
> +        # The initial value basisvalue 0 is always 1.0
> +        # FIAT_NEW code
> +        # for ii in range( results.shape[1] ):
> +        #    results[0,ii] = 1.0 + apts[ii,0]-apts[ii,0]+apts[ii,1]-apts[ii,1]
> +        code += [(format_basisvalue(0), format_float(1.0))]
>
>         # Only continue if the embedded degree is larger than zero
>         if embedded_degree > 0:
>
> -            # FIAT_NEW code
> -            # f1 = (1.0+2*x+y)/2.0
> -            # f2 = (1.0 - y) / 2.0
> -            # f3 = f2**2
> -            fac1 = create_fraction(float_1 + float_2*symbol_x + symbol_y, float_2)
> -            fac2 = create_fraction(float_1 - symbol_y, float_2)
> -            code += [(format_float_decl + str(f1), fac1)]
> -            code += [(format_float_decl + str(f2), fac2)]
> -            code += [(format_float_decl + str(f3), f2*f2)]
> -
> -            code += ["", Indent.indent(format["comment"]("Compute basisvalues"))]
> -            # The initial value basisvalue 0 is always 1.0
> -            # FIAT_NEW code
> -            # for ii in range( results.shape[1] ):
> -            #    results[0,ii] = 1.0 + apts[ii,0]-apts[ii,0]+apts[ii,1]-apts[ii,1]
> -            code += [(format_basisvalue(0), format_float(1.0))]
> -
>             # The initial value of basisfunction 1 is equal to f1
>             # FIAT_NEW code
>             # results[idx(1,0),:] = f1
> @@ -775,16 +794,16 @@
>                 #        - p/(1.0+p) * f3 *results[idx(p-1,0),:]
>                 # FIXME: KBO: Is there an error in FIAT? why is b not used?
>                 lines = []
> -                loop_vars = [(format_r, 1, embedded_degree)]
> +                loop_vars = [(str(symbol_p), 1, embedded_degree)]
>                 # Compute indices
> -                lines.append((idx0, _idx(symbol_r + float_1, float_0)))
> -                lines.append((idx1, _idx(symbol_r, float_0)))
> -                lines.append((idx2, _idx(symbol_r - float_1, float_0)))
> +                lines.append((idx0, _idx2D(symbol_p + float_1, float_0)))
> +                lines.append((idx1, _idx2D(symbol_p, float_0)))
> +                lines.append((idx2, _idx2D(symbol_p - float_1, float_0)))
>                 # Compute single helper variable an
> -                lines.append((str(an), create_fraction(float_2*symbol_r + float_1, symbol_r + float_1)))
> +                lines.append((str(an), create_fraction(float_2*symbol_p + float_1, symbol_p + float_1)))
>                 # Compute value
>                 fac0 = create_product([an, f1, basis_idx1])
> -                fac1 = create_product([create_fraction(symbol_r, float_1 + symbol_r), f3, basis_idx2])
> +                fac1 = create_product([create_fraction(symbol_p, float_1 + symbol_p), f3, basis_idx2])
>                 lines.append((str(basis_idx0), fac0 - fac1))
>                 # Create loop (block of lines)
>                 code += generate_loop(lines, loop_vars, Indent, format)
> @@ -797,21 +816,21 @@
>                 #            = ( a1 * y + a2 ) * results[idx(p,q)] \
>                 #            - a3 * results[idx(p,q-1)]
>                 lines = []
> -                loop_vars = [(format_r, 0, embedded_degree - 1),\
> -                             (format_s, 1, float_n - symbol_r)]
> +                loop_vars = [(str(symbol_p), 0, embedded_degree - 1),\
> +                             (str(symbol_q), 1, float_n - symbol_p)]
>                 # Compute indices
> -                lines.append((idx0, _idx(symbol_r, symbol_s + float_1)))
> -                lines.append((idx1, _idx(symbol_r, symbol_s)))
> -                lines.append((idx2, _idx(symbol_r, symbol_s - float_1)))
> +                lines.append((idx0, _idx2D(symbol_p, symbol_q + float_1)))
> +                lines.append((idx1, _idx2D(symbol_p, symbol_q)))
> +                lines.append((idx2, _idx2D(symbol_p, symbol_q - float_1)))
>                 # Comute all helper variables
> -                jrc = _jrc(float_2*symbol_r + float_1, float_0, symbol_s)
> +                jrc = _jrc(float_2*symbol_p + float_1, float_0, symbol_q)
>                 lines.append((str(an), jrc[0]))
>                 lines.append((str(bn), jrc[1]))
>                 lines.append((str(cn), jrc[2]))
>                 # Compute value
>                 fac0 = create_product([an * symbol_y + bn, basis_idx1])
>                 fac1 = cn * basis_idx2
> -                lines.append((format_basisvalue(idx0), fac0 - fac1))
> +                lines.append((str(basis_idx0), fac0 - fac1))
>                 # Create loop (block of lines)
>                 code += generate_loop(lines, loop_vars, Indent, format)
>
> @@ -820,14 +839,14 @@
>             #    results[idx(p,1),:] = 0.5 * (1+2.0*p+(3.0+2.0*p)*y) \
>             #        * results[idx(p,0)]
>             lines = []
> -            loop_vars = [(format_r, 0, embedded_degree)]
> +            loop_vars = [(str(symbol_p), 0, embedded_degree)]
>             # Compute indices
> -            lines.append((idx0, _idx(symbol_r, float_1)))
> -            lines.append((idx1, _idx(symbol_r, float_0)))
> +            lines.append((idx0, _idx2D(symbol_p, float_1)))
> +            lines.append((idx1, _idx2D(symbol_p, float_0)))
>             # Compute value
> -            fac0 = create_product([float_3 + float_2*symbol_r, symbol_y])
> -            fac1 = create_product([float_0_5, float_1 + float_2*symbol_r + fac0, basis_idx1])
> -            lines.append((format_basisvalue(idx0), fac1.expand().reduce_ops()))
> +            fac0 = create_product([float_3 + float_2*symbol_p, symbol_y])
> +            fac1 = create_product([float_0_5, float_1 + float_2*symbol_p + fac0, basis_idx1])
> +            lines.append((str(basis_idx0), fac1.expand().reduce_ops()))
>             # Create loop (block of lines)
>             code += generate_loop(lines, loop_vars, Indent, format)
>
> @@ -836,43 +855,186 @@
>             #    for q in range(n-p+1):
>             #        results[idx(p,q),:] *= math.sqrt((p+0.5)*(p+q+1.0))
>             lines = []
> -            loop_vars = [(format_r, 0, embedded_degree + 1), \
> -                         (format_s, 0, float_n1 - symbol_r)]
> +            loop_vars = [(str(symbol_p), 0, embedded_degree + 1), \
> +                         (str(symbol_q), 0, float_n1 - symbol_p)]
>             # Compute indices
> -            lines.append((idx0, _idx(symbol_r, symbol_s)))
> +            lines.append((idx0, _idx2D(symbol_p, symbol_q)))
>             # Compute value
> -            fac0 = create_product([symbol_r + float_0_5, symbol_r + symbol_s + float_1])
> +            fac0 = create_product([symbol_p + float_0_5, symbol_p + symbol_q + float_1])
>             fac2 = create_symbol( format_sqrt(str(fac0)), CONST )
> -            lines.append((format_basisvalue(idx0), create_product([basis_idx0, fac2])))
> +            lines.append((str(basis_idx0), create_product([basis_idx0, fac2])))
>             # Create loop (block of lines)
>             code += generate_loop(lines, loop_vars, Indent, format)
>
>     # 3D
>     elif (element_cell_domain == "tetrahedron"):
> -        count = 0
> -#        for k in range(0, data.degree()+1):  # loop over degree
> -#            for i in range(0, k+1):
> -#                for j in range(0, k - i + 1):
> -#                    ii = k-i-j
> -#                    jj = j
> -#                    kk = i
> -#                    factor = math.sqrt( (ii+0.5) * (ii+jj+1.0) * (ii+jj+kk+1.5) )
> -
> -#                    symbol = format_multiply([format_psitilde_a + format_secondary_index(ii),\
> -#                             format_scalings(format_y, ii), format_psitilde_bs(ii) + format_secondary_index(jj),\
> -#                             format_scalings(format_z, (ii+jj)),\
> -#                             format_psitilde_cs(ii, jj) + format_secondary_index(kk)])
> -
> -#                    # Declare variable
> -#                    name = format_const_float + format_basisvalue(count)
> -
> -#                    # Let inner_product handle format of factor
> -#                    value = inner_product([factor],[symbol], format)
> -
> -#                    code += [(Indent.indent(name), value)]
> -#                    count += 1
> -#        if (count != poly_dim):
> -#            error("The number of basis values must be the same as the polynomium dimension of the base")
> +        # FIAT_NEW.expansions.TetrahedronExpansionSet
> +
> +        # Compute helper factors
> +        # FIAT_NEW code
> +        # factor1 = 0.5 * ( 2.0 + 2.0*x + y + z )
> +        # factor2 = (0.5*(y+z))**2
> +        # factor3 = 0.5 * ( 1 + 2.0 * y + z )
> +        # factor4 = 0.5 * ( 1 - z )
> +        # factor5 = factor4 ** 2
> +        fac1 = create_product([float_0_5, float_2 + float_2*symbol_x + symbol_y + symbol_z])
> +        fac2 = create_product([float_0_25, symbol_y + symbol_z, symbol_y + symbol_z])
> +        fac3 = create_product([float_0_5, float_1 + float_2*symbol_y + symbol_z])
> +        fac4 = create_product([float_0_5, float_1 - symbol_z])
> +        code += [(format_float_decl + str(f1), fac1)]
> +        code += [(format_float_decl + str(f2), fac2)]
> +        code += [(format_float_decl + str(f3), fac3)]
> +        code += [(format_float_decl + str(f4), fac4)]
> +        code += [(format_float_decl + str(f5), f4*f4)]
> +
> +        code += ["", Indent.indent(format["comment"]("Compute basisvalues"))]
> +        # The initial value basisvalue 0 is always 1.0
> +        # FIAT_NEW code
> +        # for ii in range( results.shape[1] ):
> +        #    results[0,ii] = 1.0 + apts[ii,0]-apts[ii,0]+apts[ii,1]-apts[ii,1]
> +        code += [(format_basisvalue(0), format_float(1.0))]
> +
> +        # Only continue if the embedded degree is larger than zero
> +        if embedded_degree > 0:
> +
> +            # The initial value of basisfunction 1 is equal to f1
> +            # FIAT_NEW code
> +            # results[idx(1,0),:] = f1
> +            code += [(format_basisvalue(1), str(f1))]
> +
> +            # Only active is embedded_degree > 1
> +            if embedded_degree > 1:
> +
> +                # FIAT_NEW code
> +                # for p in range(1,n):
> +                #    a1 = ( 2.0 * p + 1.0 ) / ( p + 1.0 )
> +                #    a2 = p / (p + 1.0)
> +                #    results[idx(p+1,0,0)] = a1 * factor1 * results[idx(p,0,0)] \
> +                #        -a2 * factor2 * results[ idx(p-1,0,0) ]
> +                lines = []
> +                loop_vars = [(str(symbol_p), 1, embedded_degree)]
> +                # Compute indices
> +                lines.append((idx0, _idx3D(symbol_p + float_1, float_0, float_0)))
> +                lines.append((idx1, _idx3D(symbol_p          , float_0, float_0)))
> +                lines.append((idx2, _idx3D(symbol_p - float_1, float_0, float_0)))
> +                # Compute value
> +                fac1 = create_fraction(float_2*symbol_p + float_1, symbol_p + float_1)
> +                fac2 = create_fraction(symbol_p, symbol_p + float_1)
> +                fac3 = create_product([fac1, f1, basis_idx1]) - create_product([fac2, f2, basis_idx2])
> +                lines.append((str(basis_idx0), fac3))
> +                # Create loop (block of lines)
> +                code += generate_loop(lines, loop_vars, Indent, format)
> +
> +                # FIAT_NEW code
> +                # for p in range(0,n-1):
> +                #    for q in range(1,n-p):
> +                #        (aq,bq,cq) = jrc(2*p+1,0,q)
> +                #        qmcoeff = aq * factor3 + bq * factor4
> +                #        qm1coeff = cq * factor5
> +                #        results[idx(p,q+1,0)] = qmcoeff * results[idx(p,q,0)] \
> +                #            - qm1coeff * results[idx(p,q-1,0)]
> +                lines = []
> +                loop_vars = [(str(symbol_p), 0, embedded_degree - 1),\
> +                             (str(symbol_q), 1, float_n - symbol_p)]
> +                # Compute indices
> +                lines.append((idx0, _idx3D(symbol_p, symbol_q + float_1, float_0)))
> +                lines.append((idx1, _idx3D(symbol_p, symbol_q          , float_0)))
> +                lines.append((idx2, _idx3D(symbol_p, symbol_q - float_1, float_0)))
> +                # Comute all helper variables
> +                jrc = _jrc(float_2*symbol_p + float_1, float_0, symbol_q)
> +                lines.append((str(an), jrc[0]))
> +                lines.append((str(bn), jrc[1]))
> +                lines.append((str(cn), jrc[2]))
> +                # Compute value
> +                fac1 = create_product([an*f3 + bn*f4, basis_idx1]) - cn*f5*basis_idx2
> +                lines.append((str(basis_idx0), fac1))
> +                # Create loop (block of lines)
> +                code += generate_loop(lines, loop_vars, Indent, format)
> +
> +                # FIAT_NEW code
> +                # general r by recurrence
> +                # for p in range(n-1):
> +                #     for q in range(0,n-p-1):
> +                #         for r in range(1,n-p-q):
> +                #             ar,br,cr = jrc(2*p+2*q+2,0,r)
> +                #             results[idx(p,q,r+1)] = \
> +                #                         (ar * z + br) * results[idx(p,q,r) ] \
> +                #                         - cr * results[idx(p,q,r-1) ]
> +                lines = []
> +                loop_vars = [(str(symbol_p), 0, embedded_degree - 1),\
> +                             (str(symbol_q), 0, float_nm1 - symbol_p),\
> +                             (str(symbol_r), 1, float_n - symbol_p - symbol_q)]
> +                # Compute indices
> +                lines.append((idx0, _idx3D(symbol_p, symbol_q, symbol_r + float_1)))
> +                lines.append((idx1, _idx3D(symbol_p, symbol_q, symbol_r)))
> +                lines.append((idx2, _idx3D(symbol_p, symbol_q, symbol_r - float_1)))
> +                # Comute all helper variables
> +                jrc = _jrc(float_2*symbol_p + float_2*symbol_q, float_0, symbol_r)
> +                lines.append((str(an), jrc[0]))
> +                lines.append((str(bn), jrc[1]))
> +                lines.append((str(cn), jrc[2]))
> +                # Compute value
> +                fac1 = create_product([an*symbol_z + bn, basis_idx1]) - cn*basis_idx2
> +                lines.append((str(basis_idx0), fac1))
> +                # Create loop (block of lines)
> +                code += generate_loop(lines, loop_vars, Indent, format)
> +
> +            # FIAT_NEW code
> +            # q = 1
> +            # for p in range(0,n):
> +            #    results[idx(p,1,0)] = results[idx(p,0,0)] \
> +            #        * ( p * (1.0 + y) + ( 2.0 + 3.0 * y + z ) / 2 )
> +            lines = []
> +            loop_vars = [(str(symbol_p), 0, embedded_degree)]
> +            # Compute indices
> +            lines.append((idx0, _idx3D(symbol_p, float_1, float_0)))
> +            lines.append((idx1, _idx3D(symbol_p, float_0, float_0)))
> +            # Compute value
> +            fac1 = create_fraction(float_2 + float_3*symbol_y + symbol_z, float_2)
> +            fac2 = create_product([symbol_p, float_1 + symbol_y])
> +            fac3 = create_product([basis_idx1, fac2 + fac1])
> +            lines.append((str(basis_idx0), fac3))
> +            # Create loop (block of lines)
> +            code += generate_loop(lines, loop_vars, Indent, format)
> +
> +            # FIAT_NEW code
> +            # now handle r=1
> +            # for p in range(n):
> +            #    for q in range(n-p):
> +            #        results[idx(p,q,1)] = results[idx(p,q,0)] \
> +            #            * ( 1.0 + p + q + ( 2.0 + q + p ) * z )
> +            lines = []
> +            loop_vars = [(str(symbol_p), 0, embedded_degree),\
> +                         (str(symbol_q), 0, float_n - symbol_p)]
> +            # Compute indices
> +            lines.append((idx0, _idx3D(symbol_p, symbol_q, float_1)))
> +            lines.append((idx1, _idx3D(symbol_p, symbol_q, float_0)))
> +            # Compute value
> +            fac1 = create_product([float_2 + symbol_p + symbol_q, symbol_z])
> +            fac2 = create_product([basis_idx1, float_1 + symbol_p + symbol_q + fac1])
> +            lines.append((str(basis_idx0), fac2))
> +            # Create loop (block of lines)
> +            code += generate_loop(lines, loop_vars, Indent, format)
> +
> +            # FIAT_NEW code
> +            # for p in range(n+1):
> +            #    for q in range(n-p+1):
> +            #        for r in range(n-p-q+1):
> +            #            results[idx(p,q,r)] *= math.sqrt((p+0.5)*(p+q+1.0)*(p+q+r+1.5))
> +            lines = []
> +            loop_vars = [(str(symbol_p), 0, float_n1),\
> +                         (str(symbol_q), 0, float_n1 - symbol_p),\
> +                         (str(symbol_r), 0, float_n1 - symbol_p - symbol_q)]
> +            # Compute indices
> +            lines.append((idx0, _idx3D(symbol_p, symbol_q, symbol_r)))
> +            # Compute value
> +            fac0 = create_product([symbol_p + float_0_5,\
> +                                   symbol_p + symbol_q + float_1,\
> +                                   symbol_p + symbol_q + symbol_r + float_1_5])
> +            fac2 = create_symbol( format_sqrt(str(fac0)), CONST )
> +            lines.append((str(basis_idx0), create_product([basis_idx0, fac2])))
> +            # Create loop (block of lines)
> +            code += generate_loop(lines, loop_vars, Indent, format)
>     else:
>         error("Cannot compute basis values for shape: %d" % elemet_cell_domain)
>
>
>
>

---

From: Marie Rognes <[email-masked]>
Date: Fri, 08 Jan 2010 17:38:21 +0100

> In FIAT.FiniteElement: set the correct mapping (self._mapping) for all 

I will fix. I'm pretty much done with my todo-list, so I'll start going 
through the demos
and checking the generated code (for the functions that should be in place).

--
Marie

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 8 Jan 2010 17:41:19 +0100 (CET)

> and checking the generated code (for the functions that should be in place).

OK, thanks.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Fri, 8 Jan 2010 18:51:06 +0100

On Fri, Jan 08, 2010 at 05:38:21PM +0100, Marie Rognes wrote:
> I will fix. I'm pretty much done with my todo-list, so I'll start
> going through the demos

Very good!

I'm also close to finished with tabulate_tensor but it needs some
further polishing and then I need to work on the code generation for
combining tabulate_tensor from tensor and quadrature representation.

--
Anders

---

From: Anders Logg <[email-masked]>
Date: Fri, 8 Jan 2010 18:54:43 +0100

> combining tabulate_tensor from tensor and quadrature representation.

And Kristian, note that the IR for integrals needs to be a tuple of
length 3 (cell, exterior, interior), where each item in the tuple is a
list of some kind of IR for each sub domain integral.

This is because some of the logic is now handled by

  generate_integrals_code

in codegeneration.py to avoid duplication of code.

Hope that works out for the quadrature representation too.

The IR should be stored as integrals_ir in ir.

It sounds confusing but it should be clear if you look at
generate_integrals_code.

--
Anders