2010/1/8 <nore...@launchpad.net>:
------------------------------------------------------------ revno: 1429 committer: Kristian B. Ølgaard <k.b.oelga...@gmail.com> 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)
signature.asc
Description: OpenPGP digital signature
_______________________________________________ Mailing list: https://launchpad.net/~ffc Post to : ffc@lists.launchpad.net Unsubscribe : https://launchpad.net/~ffc More help : https://help.launchpad.net/ListHelp
