Subject: [Ffc] [[email-masked]: [Branch ~ffc-core/ffc/dev] Rev 1412: Implemented interpolate_vertex_values for scalar CG_1.]

From: Anders Logg <[email-masked]>
Date: Thu, 7 Jan 2010 14:45:10 +0100

Which ones need to be cleaned up?

I'll be working on getting tabulate_tensor working. Let me know if I
should help with the cleaning up.

--
Anders

Subject: [Branch ~ffc-core/ffc/dev] Rev 1412: Implemented interpolate_vertex_values for scalar CG_1.
From: [email-masked]
To: Anders Logg <[email-masked]>
Date: Thu, 07 Jan 2010 13:19:14 -0000

------------------------------------------------------------
revno: 1412
committer: Marie E. Rognes <[email-masked]>
branch nick: ffc-unstable
timestamp: Thu 2010-01-07 14:16:49 +0100
message:
  Implemented interpolate_vertex_values for scalar CG_1.
  
  Time to do some clean-ups with regard to code generation utilities
  before continuing.
modified:
  ffc/codegeneration.py
  ffc/fiatinterface.py
  ffc/representation.py

--
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/codegeneration.py'
--- ffc/codegeneration.py	2010-01-06 18:32:59 +0000
+++ ffc/codegeneration.py	2010-01-07 13:16:49 +0000
@@ -11,7 +11,7 @@
 __copyright__ = "Copyright (C) 2009 " + __author__
 __license__  = "GNU GPL version 3 or any later version"
 
-# Last changed: 2010-01-06
+# Last changed: 2010-01-07
 
 # FFC modules
 from ffc.log import begin, end, debug_code
@@ -20,6 +20,7 @@
 # FFC code generation modules
 from ffc.evaluatebasis import _evaluate_basis
 from ffc.evaluatedof import _evaluate_dof, _evaluate_dofs, affine_weights
+from ffc.codesnippets import jacobian
 
 # FFC specialized code generation modules
 #from ffc.quadrature import generate_quadrature_integrals
@@ -98,7 +99,7 @@
     code["evaluate_basis_derivatives_all"] = not_implemented
     code["evaluate_dof"] = _evaluate_dof(ir["evaluate_dof"])
     code["evaluate_dofs"] = _evaluate_dofs(ir["evaluate_dofs"])
-    code["interpolate_vertex_values"] = not_implemented
+    code["interpolate_vertex_values"] = _interpolate_vertex_values(ir["interpolate_vertex_values"])
     code["num_sub_elements"] = ret(ir["num_sub_elements"])
     code["create_sub_element"] = _create_foo(ir["create_sub_element"], prefix, "finite_element")
 
@@ -264,31 +265,64 @@
 
 def _tabulate_coordinates(ir):
 
+    # Raise error if tabulate_coordinates is ill-defined
     if ir is None:
         return "// Raise error here"
 
     # Extract formats:
-    add = format["add"]
-    multiply = format["multiply"]
+    add, multiply = format["add"], format["multiply"]
     precision = format["float"]
 
     # Extract coordinates and cell dimension
-    coordinates = ir
-    cell_dim = len(coordinates[0])
+    cell_dim = len(ir[0])
 
     # Aid mapping points from reference to physical element
     coefficients = affine_weights(cell_dim)
 
+    # Generate code for each point and each component
     code = ["const double * const * x = c.coordinates;"]
-
-    for (i, c) in enumerate(coordinates):
-        w = coefficients(c)
+    for (i, coordinate) in enumerate(ir):
+        w = coefficients(coordinate)
         for j in range(cell_dim):
             value = add([multiply([precision(w[k]), "x[%d][%d]" % (k, j)])
                          for k in range(cell_dim + 1)])
             code += ["coordinates[%d][%d] = %s;" % (i, j, value)]
     return "\n".join(code)
 
+def _interpolate_vertex_values(ir):
+
+    code = []
+
+    # Add code for Jacobian if necessary
+    if ir["needs_jacobian"]:
+        code += [jacobian(ir["cell_dim"])]
+
+    # Extract formats
+    add, multiply = format["add"], format["multiply"]
+    precision = format["float"]
+
+    # Iterate over the elements
+    for (k, data) in enumerate(ir["element_data"]):
+
+        # Extract vertex values for all basis functions
+        vertex_values = data["basis_values"]
+
+        # Create code for each vertex
+        for j in range(len(vertex_values)):
+
+            # Map basis values according to mapping
+            vertex_value = [precision(v) for v in vertex_values[j]]
+
+            # Contract basis values and coefficients
+            value = add([multiply(["dof_values[%d]" % i, v])
+                         for (i, v) in enumerate(vertex_value)])
+
+            # Construct correct vertex value label
+            name = "vertex_values[%d]" % j
+            code += [name + " = " +  value + ";"]
+
+    return "\n".join(code)
+
 #--- Utility functioins ---
 
 def _create_foo(numbers, prefix, class_name):

=== modified file 'ffc/fiatinterface.py'
--- ffc/fiatinterface.py	2010-01-04 12:48:48 +0000
+++ ffc/fiatinterface.py	2010-01-07 13:16:49 +0000
@@ -5,7 +5,7 @@
 
 # Modified by Garth N. Wells, 2009.
 # Modified by Marie Rognes, 2009-2010.
-# Last changed: 2010-01-04
+# Last changed: 2010-01-07
 
 # UFL modules
 from ufl import FiniteElement as UFLFiniteElement
@@ -44,6 +44,11 @@
 # Mapping from dimension to number of mesh sub-entities:
 entities_per_dim = {1: [2, 1], 2: [3, 3, 1], 3: [4, 6, 4, 1]}
 
+def reference_cell(dim):
+    if isinstance(dim, int):
+        return ufc_simplex(dim)
+    else:
+        return ufc_simplex(domain2dim[dim])
 
 def create_element(ufl_element):
 
@@ -70,7 +75,7 @@
     "Create FIAT element corresponding to given UFL finite element."
 
     ElementClass = family2class[ufl_element.family()]
-    reference_cell = ufc_simplex(domain2dim[ufl_element.cell().domain()])
-    element = ElementClass(reference_cell, ufl_element.degree())
+    cell = reference_cell(ufl_element.cell().domain())
+    element = ElementClass(cell, ufl_element.degree())
 
     return element

=== modified file 'ffc/representation.py'
--- ffc/representation.py	2010-01-06 18:32:59 +0000
+++ ffc/representation.py	2010-01-07 13:16:49 +0000
@@ -28,7 +28,7 @@
 # FFC modules
 from ffc.utils import compute_permutations
 from ffc.log import info, error, begin, end, debug_ir, ffc_assert
-from ffc.fiatinterface import create_element, entities_per_dim
+from ffc.fiatinterface import create_element, entities_per_dim, reference_cell
 from ffc.mixedelement import MixedElement

@@ -63,11 +63,12 @@
 
     # Create FIAT element
     element = create_element(ufl_element)
+    cell = ufl_element.cell()
 
     # Compute data for each function
     ir = {}
     ir["signature"] = repr(ufl_element)
-    ir["cell_shape"] = ufl_element.cell().domain()
+    ir["cell_shape"] = cell.domain()
     ir["space_dimension"] = element.space_dimension()
     ir["value_rank"] = len(ufl_element.value_shape())
     ir["value_dimension"] = ufl_element.value_shape()
@@ -75,11 +76,12 @@
     ir["evaluate_basis_all"] = not_implemented #element.get_coeffs()
     ir["evaluate_basis_derivatives"] = element
     ir["evaluate_basis_derivatives_all"] = not_implemented #element.get_coeffs()
-    ir["evaluate_dof"] = _evaluate_dof(element, ufl_element.cell())
-    ir["evaluate_dofs"] = None
-    ir["interpolate_vertex_values"] = None
+    ir["evaluate_dof"] = _evaluate_dof(element, cell)
+    ir["evaluate_dofs"] = not_implemented
+    ir["interpolate_vertex_values"] = _interpolate_vertex_values(element, cell)
     ir["num_sub_elements"] = ufl_element.num_sub_elements()
-    ir["create_sub_element"] = [form_data.element_map[e] for e in ufl_element.sub_elements()]
+    ir["create_sub_element"] = [form_data.element_map[e]
+                                for e in ufl_element.sub_elements()]
 
     debug_ir(ir, "finite_element")
 
@@ -116,7 +118,8 @@
     ir["tabulate_entity_dofs"] = not_implemented
     ir["tabulate_coordinates"] = _tabulate_coordinates(element)
     ir["num_sub_dof_maps"] = ufl_element.num_sub_elements()
-    ir["create_sub_dof_map"] = [form_data.element_map[e] for e in ufl_element.sub_elements()]
+    ir["create_sub_dof_map"] = [form_data.element_map[e]
+                                for e in ufl_element.sub_elements()]
 
     debug_ir(ir, "dofmap")
 
@@ -202,24 +205,18 @@
 def _evaluate_dof(element, cell):
     "Compute intermediate representation of evaluate_dof."
 
-    # Setup ir
-    ir = {"mappings": element.mapping(),
-          "value_dim": _value_dimension(element),
-          "cell_dimension": cell.geometric_dimension(),
-          "dofs": [L.pt_dict for L in element.dual_basis()]}
-
     # Generate offsets: i.e value offset for each basis function
-    if not isinstance(element, MixedElement):
-        ir["offsets"] = [0]*element.space_dimension()
-    else:
-        offsets = []
-        offset = 0
-        for e in element.elements():
-            offsets += [offset]*e.space_dimension()
-            offset += _value_dimension(e)
-        ir["offsets"] = offsets
+    offsets = []
+    offset = 0
+    for e in all_elements(element):
+        offsets += [offset]*e.space_dimension()
+        offset += _value_dimension(e)
 
-    return ir
+    return {"mappings": element.mapping(),
+            "value_dim": _value_dimension(element),
+            "cell_dimension": cell.geometric_dimension(),
+            "dofs": [L.pt_dict for L in element.dual_basis()],
+            "offsets": offsets}
 
 def _tabulate_coordinates(element):
     "Compute intermediate representation of tabulate_coordinates."
@@ -230,16 +227,10 @@
 def _tabulate_dofs(element, cell):
     "Compute intermediate representation of tabulate_dofs."
 
-    if isinstance(element, MixedElement):
-        ir = [{"entites_per_dim": entities_per_dim[cell.geometric_dimension()],
+    elements = all_elements(element)
+    return [{"entites_per_dim": entities_per_dim[cell.geometric_dimension()],
                "num_dofs_per_entity": _num_dofs_per_entity(e)}
-              for e in element.elements()]
-    else:
-        ir = [{"entites_per_dim": entities_per_dim[cell.geometric_dimension()],
-               "num_dofs_per_entity": _num_dofs_per_entity(element)}]
-
-    return ir
-
+              for e in elements]
 
 def _tabulate_facet_dofs(element, cell):
     "Compute intermediate representation of tabulate_facet_dofs."
@@ -271,8 +262,39 @@
 
     return facet_dofs
 
+def _interpolate_vertex_values(element, cell):
+
+    ir = {}
+    ir["cell_dim"] = cell.geometric_dimension()
+
+    # Check whether computing the Jacobian is necessary
+    mappings = element.mapping()
+    ir["needs_jacobian"] = any("piola" in m for m in mappings)
+
+    # Get vertices of reference cell
+    vertices = reference_cell(cell.domain()).get_vertices()
+
+    # Compute data for each constituent element
+    extract = lambda values: [a for a in values[(0, 0)]]
+
+
+    ir["element_data"] = [{"space_dim": e.space_dimension(),
+                           "value_dim": _value_dimension(e),
+                           "basis_values": extract(e.tabulate(0, vertices)),
+                           "mapping": e.mapping()}
+                          for e in all_elements(element)]
+
+    return ir
+
+
 #--- Utility functions ---
 
+def all_elements(element):
+    try:
+        return element.elements()
+    except:
+        return [element]
+
 def _num_dofs_per_entity(element):
     """
     Compute list of integers representing the number of dofs

---

From: Marie Rognes <[email-masked]>
Date: Thu, 07 Jan 2010 15:07:28 +0100

> should help with the cleaning up.

I'm boldly throwing away anything that seems obsolete....

--
Marie

> --
> Anders
>   
>
> ------------------------------------------------------------------------
>
> Subject:
> [Branch ~ffc-core/ffc/dev] Rev 1412: Implemented 
> interpolate_vertex_values for scalar CG_1.
> From:
> [email-masked]
> Date:
> ------------------------------------------------------------------------
>
> _______________________________________________
> Mailing list: https://launchpad.net/~ffc
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~ffc
> More help   : https://help.launchpad.net/ListHelp

---

From: Anders Logg <[email-masked]>
Date: Thu, 7 Jan 2010 15:17:12 +0100

> I'm boldly throwing away anything that seems obsolete....

Good!

--
Anders