ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jan 12, 2010 at 10:57:13AM -0000, noreply@xxxxxxxxxxxxx wrote:
> ------------------------------------------------------------
> revno: 1438
> committer: Marie E. Rognes <meg@xxxxxxxxx>
> branch nick: ffc-unstable
> timestamp: Tue 2010-01-12 11:53:45 +0100
> message:
>   Fixed interpolate vertex values for mixed and piola mapped elements. Puh.

I'm very happy I didn't have to implement interpolate_vertex_values! :-)

>   Mysterious error encountered in debug_dict -- could someone take a
>   look?

I will take a look.

Where should I look?

--
Anders


> added:
>   ffc/interpolatevertexvalues.py
> modified:
>   ffc/codegeneration.py
>   ffc/cpp.py
>   ffc/log.py
>   ffc/representation.py
>
>

> === modified file 'ffc/codegeneration.py'
> --- ffc/codegeneration.py	2010-01-11 12:47:41 +0000
> +++ ffc/codegeneration.py	2010-01-12 10:53:45 +0000
> @@ -11,7 +11,7 @@
>  __copyright__ = "Copyright (C) 2009 " + __author__
>  __license__  = "GNU GPL version 3 or any later version"
>
> -# Last changed: 2010-01-11
> +# Last changed: 2010-01-12
>
>  # FFC modules
>  from ffc.log import begin, end, debug_code
> @@ -21,6 +21,7 @@
>  from ffc.evaluatebasis import _evaluate_basis
>  from ffc.evaluatebasisderivatives import _evaluate_basis_derivatives
>  from ffc.evaluatedof import _evaluate_dof, _evaluate_dofs, affine_weights
> +from ffc.interpolatevertexvalues import _interpolate_vertex_values
>  from ffc.codesnippets import jacobian, cell_coordinates
>
>  # FFC specialized code generation modules
> @@ -78,7 +79,7 @@
>      code["value_dimension"] = _value_dimension(ir["value_dimension"])
>      code["evaluate_basis"] = _evaluate_basis(ir["evaluate_basis"])
>      code["evaluate_basis_all"] = not_implemented
> -#    code["evaluate_basis_derivatives"] = do_nothing
> +    #code["evaluate_basis_derivatives"] = do_nothing
>      code["evaluate_basis_derivatives"] = _evaluate_basis_derivatives(ir["evaluate_basis"])
>      code["evaluate_basis_derivatives_all"] = not_implemented
>      code["evaluate_dof"] = _evaluate_dof(ir["evaluate_dof"])
> @@ -133,46 +134,6 @@
>
>      return code
>
> -def _tabulate_entity_dofs(ir):
> -    "Generate code for tabulate_entity_dofs."
> -
> -    # Extract variables from ir
> -    entity_dofs, num_dofs_per_entity = ir
> -
> -    # Prefetch formats
> -    assign, component = format["assign"], format["component"]
> -
> -    # Add check that dimension and number of mesh entities is valid
> -    dim = len(num_dofs_per_entity)
> -    excpt = format["exception"]("d is larger than dimension (%d)" % (dim - 1))
> -    code = format["if"]("d > %d" % (dim-1), excpt)
> -
> -    # Generate cases for each dimension:
> -    all_cases = ["" for d in range(dim)]
> -    for d in range(dim):
> -
> -        # Ignore if no entities for this dimension
> -        if num_dofs_per_entity[d] == 0: continue
> -
> -        # Add check that given entity is valid:
> -        num_entities = len(entity_dofs[d].keys())
> -        excpt = format["exception"]("i is larger than number of entities (%d)"
> -                                    % (num_entities - 1))
> -        check = format["if"]("i > %d" % (num_entities - 1), excpt)
> -
> -        # Generate cases for each mesh entity
> -        cases = ["".join([assign(component("dofs", j), dof)
> -                          for (j, dof) in enumerate(entity_dofs[d][entity])])
> -                 for entity in range(num_entities)]
> -
> -        # Generate inner switch with preceding check
> -        all_cases[d] = "\n".join([check, format["switch"]("i", cases)])
> -
> -    # Generate outer switch
> -    code += format["switch"]("d", all_cases)
> -
> -    return code
> -
>  def generate_integrals_code(ir, prefix, options):
>      "Generate code for integrals from intermediate representation."
>
> @@ -357,45 +318,48 @@
>
>      return code
>
> -
> -def _interpolate_vertex_values(ir):
> -    "Generate code for interpolate_vertex_values."
> -
> -    code = ""
> -
> -    # Add code for Jacobian if necessary
> -    if ir["needs_jacobian"]:
> -        code += jacobian[ir["cell_dim"]]
> -
> -    # Extract formats
> -    inner_product = format["inner product"]
> -    component = format["component"]
> -    precision = format["float"]
> -    assign = format["assign"]
> -
> -    # Iterate over the elements
> -    for (k, data) in enumerate(ir["element_data"]):
> -
> -        # Extract vertex values for all basis functions
> -        vertex_values = data["basis_values"]
> -
> -        code += format["comment"]("Evaluate at vertices and map")
> -
> -        # Create code for each vertex
> -        for (j, vertex_value) in enumerate(vertex_values):
> -
> -            # FIXME: Map basis values according to mapping
> -
> -            # Contract basis values and coefficients
> -            dof_values = [component("dof_values", i)
> -                          for i in range(len(vertex_value))]
> -            value = inner_product(dof_values, vertex_value)
> -
> -            # Construct assign value to correct vertex_value
> -            code += assign(component("vertex_values", j), value)
> +def _tabulate_entity_dofs(ir):
> +    "Generate code for tabulate_entity_dofs."
> +
> +    # Extract variables from ir
> +    entity_dofs, num_dofs_per_entity = ir
> +
> +    # Prefetch formats
> +    assign, component = format["assign"], format["component"]
> +
> +    # Add check that dimension and number of mesh entities is valid
> +    dim = len(num_dofs_per_entity)
> +    excpt = format["exception"]("d is larger than dimension (%d)" % (dim - 1))
> +    code = format["if"]("d > %d" % (dim-1), excpt)
> +
> +    # Generate cases for each dimension:
> +    all_cases = ["" for d in range(dim)]
> +    for d in range(dim):
> +
> +        # Ignore if no entities for this dimension
> +        if num_dofs_per_entity[d] == 0: continue
> +
> +        # Add check that given entity is valid:
> +        num_entities = len(entity_dofs[d].keys())
> +        excpt = format["exception"]("i is larger than number of entities (%d)"
> +                                    % (num_entities - 1))
> +        check = format["if"]("i > %d" % (num_entities - 1), excpt)
> +
> +        # Generate cases for each mesh entity
> +        cases = ["".join([assign(component("dofs", j), dof)
> +                          for (j, dof) in enumerate(entity_dofs[d][entity])])
> +                 for entity in range(num_entities)]
> +
> +        # Generate inner switch with preceding check
> +        all_cases[d] = "\n".join([check, format["switch"]("i", cases)])
> +
> +    # Generate outer switch
> +    code += format["switch"]("d", all_cases)
>
>      return code
>
> +
> +
>  #--- Utility functioins ---
>
>  def _create_foo(numbers, prefix, class_name):
>
> === modified file 'ffc/cpp.py'
> --- ffc/cpp.py	2010-01-11 09:35:25 +0000
> +++ ffc/cpp.py	2010-01-12 10:53:45 +0000
> @@ -7,7 +7,7 @@
>
>  # Modified by Kristian B. Oelgaard 2009
>  # Modified by Marie E. Rognes 2010
> -# Last changed: 2010-01-11
> +# Last changed: 2010-01-12
>
>  # Python modules.
>  import re
> @@ -52,10 +52,13 @@
>                 "scale factor":    "det",
>                 "transform":       lambda t, j, k, r: _transform(t, j, k, r)})
>
> -# Mesh entity variable names
> +# Geometry related variable names
>  format.update({"entity index": "c.entity_indices",
>                 "num entities": "m.num_entities",
> -               "cell": lambda s: "ufc::%s" % s})
> +               "cell": lambda s: "ufc::%s" % s,
> +               "det(J)": "detJ",
> +               "J": lambda i, j: "J_%d%d" % (i, j),
> +               "Jinv" : lambda i, j: "Jinv_%d%d" % (i, j)})
>
>  # Misc
>  format.update({"bool":    lambda v: {True: "true", False: "false"}[v],
> @@ -82,24 +85,36 @@
>      one are ignored.
>      """
>
> -    # FIXME: Paranthesise negative numbers and fix ad-hoc solutions
> +    # FIXME: This could probably be way more robust and elegant.
>
>      cpp_str = format["str"]
>      non_zero_factors = []
>      for f in factors:
> +
> +        # Round-off if f is smaller than epsilon
>          if isinstance(f, (int, float)):
> -            if f < format["epsilon"]:
> +            if abs(f) < format["epsilon"]:
>                  return cpp_str(0)
>
>          # Convert to string
>          f = cpp_str(f)
>
> +        # Return zero if any factor is zero
>          if f == "0":
>              return cpp_str(0)
> +
> +        # Ignore 1 factors
>          if f == "1" or f == "1.0":
> -            pass
> -        else:
> -            non_zero_factors += [f]
> +            continue
> +
> +        # If sum-like, parentheseze factor
> +        if "+" in f or "-" in f:
> +            f = "(%s)" % f
> +
> +        non_zero_factors += [f]
> +
> +    if len(non_zero_factors) == 0:
> +        return cpp_str(1.0)
>
>      return "*".join(non_zero_factors)
>
> @@ -107,7 +122,10 @@
>      "Generate string summing a list of strings."
>
>      # FIXME: Subtract absolute value of negative numbers
> -    return " + ".join([str(t) for t in terms if (str(t) != "0")])
> +    result = " + ".join([str(t) for t in terms if (str(t) != "0")])
> +    if result == "":
> +        return format["str"](0)
> +    return result
>
>  def _inner_product(v, w):
>      "Generate string for v[0]*w[0] + ... + v[n]*w[n]."
>
> === added file 'ffc/interpolatevertexvalues.py'
> --- ffc/interpolatevertexvalues.py	1970-01-01 00:00:00 +0000
> +++ ffc/interpolatevertexvalues.py	2010-01-12 10:53:45 +0000
> @@ -0,0 +1,98 @@
> +"Code generation for interpolate_vertex_values."
> +
> +__author__ = "Marie E. Rognes (meg@xxxxxxxxx)"
> +__copyright__ = "Copyright (C) 2009"
> +__license__  = "GNU GPL version 3 or any later version"
> +
> +# Last changed: 2010-01-12
> +
> +from ffc.codesnippets import jacobian
> +from ffc.cpp import format
> +
> +# Extract code manipulation formats
> +inner_product = format["inner product"]
> +component = format["component"]
> +precision = format["float"]
> +assign = format["assign"]
> +multiply = format["multiply"]
> +
> +# Extract formats for the Jacobians
> +J = format["J"]
> +Jinv = format["Jinv"]
> +invdetJ = "1.0/%s" % format["det(J)"]
> +
> +def _interpolate_vertex_values(ir):
> +    "Generate code for interpolate_vertex_values."
> +
> +    # meg: Note: I think this works, but it might be more elegant to
> +    # give mixed element a tabulate function
> +
> +    code = ""
> +
> +    # Add code for Jacobian if necessary
> +    dim = ir["cell_dim"]
> +    if ir["needs_jacobian"]:
> +        code += jacobian[dim]
> +
> +    # Compute total value dimension for (mixed) element
> +    total_value_dim = sum(data["value_dim"] for data in ir["element_data"])
> +
> +    # Generate code for each element
> +    value_offset = 0
> +    space_offset = 0
> +    for data in ir["element_data"]:
> +
> +        code += format["comment"]("Evaluate at vertices and map")
> +        code += _interpolate_vertex_values_element(data, dim, total_value_dim,
> +                                                   value_offset, space_offset)
> +
> +        # Update offsets for value- and space dimension
> +        value_offset += data["value_dim"]
> +        space_offset += data["space_dim"]
> +
> +    return code
> +
> +
> +def _interpolate_vertex_values_element(data, dim, total_value_dim,
> +                                       value_offset=0, space_offset=0):
> +
> +    # Extract vertex values for all basis functions
> +    vertex_values = data["basis_values"]
> +    value_dim = data["value_dim"]
> +    space_dim = data["space_dim"]
> +    mapping = data["mapping"]
> +
> +    # Create code for each value dimension:
> +    code = ""
> +    for k in range(value_dim):
> +
> +        # Create code for each vertex x_j
> +        for (j, values_at_vertex) in enumerate(vertex_values):
> +
> +            if value_dim == 1: values_at_vertex = [values_at_vertex]
> +
> +            # Map basis values according to element mapping
> +            if mapping is "affine":
> +                components = values_at_vertex[k]
> +
> +            elif mapping == "contravariant piola":
> +                contraction = lambda i: inner_product([J(k, l) for l in range(dim)],
> +                                                      [values_at_vertex[l][i] for l in range(dim)])
> +                components = [multiply([invdetJ, contraction(i)]) for i in range(space_dim)]
> +
> +            elif mapping == "covariant piola":
> +                components = [inner_product([Jinv(k, l) for l in range(dim)],
> +                                            [values_at_vertex[l][i] for l in range(dim)])
> +                              for i in range(space_dim)]
> +            else:
> +                raise Exception, "No such mapping: %s accepted" % mapiping
> +
> +            # Contract coefficents and basis functions
> +            dof_values = [component("dof_values", i + space_offset) for i in range(space_dim)]
> +            value = inner_product(dof_values, components)
> +
> +            # Assign value to correct vertex
> +            index = j*total_value_dim + (k + value_offset)
> +            code += assign(component("vertex_values", index), value)
> +
> +    return code
>
> === modified file 'ffc/log.py'
> --- ffc/log.py	2009-12-17 10:28:31 +0000
> +++ ffc/log.py	2010-01-12 10:53:45 +0000
> @@ -11,7 +11,7 @@
>  __license__  = "GNU GPL version 3 or any later version"
>
>  # Modified by Kristian B. Oelgaard, 2009
> -# Last changed: 2009-12-17
> +# Last changed: 2010-01-12
>
>  # UFL modules.
>  from ufl.log import Logger
> @@ -43,7 +43,12 @@
>      for (key, value) in d.iteritems():
>          info(key)
>          info("-"*len(key))
> -        info(str(value))
> +
> +        try:
> +            info(str(value))
> +        except:
> +            # FIXME
> +            info("Mysterious error here")
>          info("")
>      end()
>
>
> === modified file 'ffc/representation.py'
> --- ffc/representation.py	2010-01-11 10:30:48 +0000
> +++ ffc/representation.py	2010-01-12 10:53:45 +0000
> @@ -98,9 +98,6 @@
>      element = create_element(ufl_element)
>      cell = ufl_element.cell()
>
> -    print "\nelement.entity_dofs():"
> -    print element.entity_dofs()
> -
>      # Precompute repeatedly used items
>      num_dofs_per_entity = _num_dofs_per_entity(element)
>      facet_dofs = _tabulate_facet_dofs(element, cell)
> @@ -201,6 +198,7 @@
>
>      # FIXME: Arbitrary tensor elements?
>      # FIXME: Move to FiniteElement/MixedElement
> +
>      shape = element.value_shape()
>      if shape == ():
>          return 1
> @@ -282,12 +280,11 @@
>
>      # Compute data for each constituent element
>      extract = lambda values: values[values.keys()[0]].transpose()
> -
>      ir["element_data"] = [{"value_dim": _value_dimension(e),
>                             "basis_values": extract(e.tabulate(0, vertices)),
> -                           "mapping": e.mapping()}
> +                           "mapping": e.mapping()[0],
> +                           "space_dim": e.space_dimension()}
>                            for e in all_elements(element)]
> -
>      return ir
>
>
>

Attachment: signature.asc
Description: Digital signature