Subject: [FFC-dev] integration over boundaries

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 01 May 2006 19:40:00 +0200

Attached is a changeset for integrating terms over boundaries. Output is
only in the default DOLFIN format at this stage. 

The CVS version of FIAT is required.

The output looks OK, but I haven't tested it yet with DOLFIN. On the
DOLFIN side, the appropriate "det", which is dependent on the facet,
needs to be computed in AffineMap.

Garth
# HG changeset patch
# User "Garth N. Wells <[email-masked]>"
# Node ID 436ff708086d3e065c65f9d2cdd29d9c4d6b4cfc
# Parent  55828e695e1461e902bfa129e70d358cbc56b4a1
Support integrals over boundaries.

Output only in standard DOLFIN output.

diff -r 55828e695e14 -r 436ff708086d src/demo/NeumannProblem.form
--- a/src/demo/NeumannProblem.form	Thu Apr  6 11:56:57 2006 -0500
+++ b/src/demo/NeumannProblem.form	Mon May  1 19:31:20 2006 +0200
@@ -6,7 +6,7 @@
 #
 # Compile this form with FFC: ffc NeumannProblem.form
 
-element = FiniteElement("Lagrange", "triangle", 1)
+element = FiniteElement("Vector Lagrange", "triangle", 1)
 
 v = TestFunction(element)
 U = TrialFunction(element)
@@ -14,4 +14,5 @@ g = Function(element)
 g = Function(element)
 
 a = dot(grad(v), grad(U))*dx
-L = v*f*dx + v*g*ds
+#L = v*f*dx + v*g*ds
+L = dot(v, f)*dx + dot(v,g)*ds
diff -r 55828e695e14 -r 436ff708086d src/ffc/compiler/compiler.py
--- a/src/ffc/compiler/compiler.py	Thu Apr  6 11:56:57 2006 -0500
+++ b/src/ffc/compiler/compiler.py	Mon May  1 19:31:20 2006 +0200
@@ -116,11 +116,31 @@ def build(sums, name = "Form", language 
 
         # Create element tensors
         debug("Compiling tensor representation for interior")
-        form.AKi = ElementTensor(form.sum, "interior", format, cK_used, options)
+        form.AKi = ElementTensor(form.sum, "interior", format, cK_used, options, None)
+
         debug("Compiling tensor representation for boundary")
-        form.AKb = ElementTensor(form.sum, "boundary", format, cK_used, options)
+        # Determine number of facets (FIXME: this could be provided by finiteelement)
+        if len(form.elements) != 0:
+            shape = form.elements[0].shape()
+        elif form.test != None:
+            shape = form.test.shape() 
+        else:
+            print "Unable to determine element shape."
+
+        if shape == TRIANGLE:
+            facets = 3
+        elif shape == TETRAHEDRON:
+            facets = 4
+
+         # Compute reference tensor for each facet
+        form.AKb = []
+        debug("Compiling tensor representation for boundaries")
+        for i in range(facets):
+            form.AKb.append(ElementTensor(form.sum, "boundary", format, cK_used, options, i))
+            form.AKb[i].facet = i
 
         # Compute coefficient declarations
+        # FIXME: projections not programmed for boundaries (are changes is anything needed?)
         form.cK = __compute_coefficients(form.projections, format, cK_used)
 
         # Check primary ranks
@@ -181,7 +201,10 @@ def writeFiniteElement(element, name = "
 
 def __check_primary_ranks(form):
     "Check that all primary ranks are equal."
-    terms = form.AKi.terms + form.AKb.terms
+    terms = form.AKi.terms
+    for AKb in form.AKb:
+        terms += AKb.terms
+
     ranks = [term.A0.i.rank for term in terms]
     if not ranks[1:] == ranks[:-1]:
         "Form must be linear in each of its arguments."
diff -r 55828e695e14 -r 436ff708086d src/ffc/compiler/elementtensor.py
--- a/src/ffc/compiler/elementtensor.py	Thu Apr  6 11:56:57 2006 -0500
+++ b/src/ffc/compiler/elementtensor.py	Mon May  1 19:31:20 2006 +0200
@@ -30,7 +30,7 @@ class ElementTensor:
         aK    - a list of precomputed element tensor declarations
     """
 
-    def __init__(self, sum, type, format, cK_used, options):
+    def __init__(self, sum, type, format, cK_used, options, facet):
         "Create ElementTensor."
 
         # Check that all Products have integrals
@@ -50,7 +50,7 @@ class ElementTensor:
                 # Check if reference tensor should be computed
                 if factorization[i] == None:
                     # Compute reference tensor and add term
-                    A0 = ReferenceTensor(p)
+                    A0 = ReferenceTensor(p, facet)
                     self.terms[i] = Term(p, A0, [GK])
                 else:
                     # Add geometry tensor to previous term
diff -r 55828e695e14 -r 436ff708086d src/ffc/compiler/finiteelement.py
--- a/src/ffc/compiler/finiteelement.py	Thu Apr  6 11:56:57 2006 -0500
+++ b/src/ffc/compiler/finiteelement.py	Mon May  1 19:31:20 2006 +0200
@@ -148,10 +148,15 @@ http://www.fenics.org/pipermail/fiat-dev
         else:
             raise RuntimeError, "Can only handle scalar or vector-valued elements."
 
-    def tabulate(self, order, points):
+    def tabulate(self, order, points, facet):
         """Return tabulated values of derivatives up to given order of
         basis functions at given points."""
-        return self.element.function_space().tabulate_jet(order, points)
+        if facet == None:
+            return self.element.function_space().tabulate_jet(order, points)
+        else:
+            # FIXME: compute facet shape elsewhere
+            shape = self.shape() - 1
+            return self.element.function_space().trace_tabulate_jet(shape, facet, order, points)
 
     def __add__(self, other):
         "Create mixed element."
diff -r 55828e695e14 -r 436ff708086d src/ffc/compiler/form.py
--- a/src/ffc/compiler/form.py	Thu Apr  6 11:56:57 2006 -0500
+++ b/src/ffc/compiler/form.py	Mon May  1 19:31:20 2006 +0200
@@ -61,6 +61,7 @@ class Form:
         self.name        = name
         self.AKi         = None
         self.AKb         = None
+        self.facet       = None
         self.cK          = None
         self.rank        = None
         self.dims        = None
diff -r 55828e695e14 -r 436ff708086d src/ffc/compiler/monomialintegration.py
--- a/src/ffc/compiler/monomialintegration.py	Thu Apr  6 11:56:57 2006 -0500
+++ b/src/ffc/compiler/monomialintegration.py	Mon May  1 19:31:20 2006 +0200
@@ -23,17 +23,20 @@ from algebra import *
 from algebra import *
 from multiindex import *
 
-def integrate(product):
+def integrate(product, facet):
     """Compute the reference tensor for a given monomial term of a
     multilinear form, given as a Product."""
 
     debug("Pretabulating basis functions at quadrature points")
 
+    # Check for integration type (interior or boundary)
+    type = product.integral.type
+
     # Initialize quadrature points and weights
-    (points, weights, vscaling, dscaling) = __init_quadrature(product.basisfunctions)
+    (points, weights, vscaling, dscaling) = __init_quadrature(product.basisfunctions, type)
 
     # Initialize quadrature table for basis functions
-    table = __init_table(product.basisfunctions, points)
+    table = __init_table(product.basisfunctions, points, facet)
 
     # Compute table Psi for each factor
     psis = [__compute_psi(v, table, len(points), dscaling) for v in product.basisfunctions]
@@ -43,7 +46,7 @@ def integrate(product):
 
     return A0
 
-def __init_quadrature(basisfunctions):
+def __init_quadrature(basisfunctions, type):
     "Initialize quadrature for given monomial term."
 
     debug("Initializing quadrature.", 1)
@@ -57,22 +60,34 @@ def __init_quadrature(basisfunctions):
     debug("Total degree is %d, using %d quadrature point(s) in each dimension" % (q, m), 1)
 
     # Create quadrature rule and get points and weights
-    quadrature = make_quadrature(shape, m)
+    # FIXME: FIAT ot finiteelement should return shape of facet
+    if type == 'interior':
+        quadrature = make_quadrature(shape, m)
+    elif type == 'boundary':
+        quadrature = make_quadrature(shape-1, m)
     points = quadrature.get_points()
     weights = quadrature.get_weights()
 
     # Compensate for different choice of reference cells in FIAT
     # FIXME: Convince Rob to change his reference elements
     if shape == TRIANGLE:
-        vscaling = 0.25  # Area 1/2 instead of 2
-        dscaling = 2.0   # Scaling of derivative
+        if type == 'interior':
+            vscaling = 0.25  # Area 1/2 instead of 2
+            dscaling = 2.0   # Scaling of derivative
+        elif type == 'boundary':
+            vscaling = 0.5   # Length 1 instead of 2
+            dscaling = 2.0   # Scaling of derivative        
     elif shape == TETRAHEDRON:
-        vscaling = 0.125 # Volume 1/6 instead of 4/3
-        dscaling = 2.0   # Scaling of derivative
-        
+        if type == 'interior':
+            vscaling = 0.125 # Volume 1/6 instead of 4/3
+            dscaling = 2.0   # Scaling of derivative
+        elif type == 'boundary':
+            vscaling = 0.25  # Area 1/2 instead of 2
+            dscaling = 2.0   # Scaling of derivative
+
     return (points, weights, vscaling, dscaling)
 
-def __init_table(basisfunctions, points):
+def __init_table(basisfunctions, points, facet):
     """Initialize table of basis functions and their derivatives at
     the given quadrature points for each element."""
 
@@ -92,7 +107,7 @@ def __init_table(basisfunctions, points)
     table = {}
     for element in num_derivatives:
         order = num_derivatives[element]
-        table[element] = element.tabulate(order, points)
+        table[element] = element.tabulate(order, points, facet)
 
     return table
 
diff -r 55828e695e14 -r 436ff708086d src/ffc/compiler/referencetensor.py
--- a/src/ffc/compiler/referencetensor.py	Thu Apr  6 11:56:57 2006 -0500
+++ b/src/ffc/compiler/referencetensor.py	Mon May  1 19:31:20 2006 +0200
@@ -37,7 +37,7 @@ class ReferenceTensor:
         integral       - an Integral
         cputime        - time to compute the reference tensor"""
 
-    def __init__(self, product):
+    def __init__(self, product, facet):
         "Create ReferenceTensor."
 
         # Check that we get a Product
@@ -60,7 +60,7 @@ class ReferenceTensor:
 
         # Compute reference tensor (new version)
         t = time.time()
-        self.A0 = integrate(product)
+        self.A0 = integrate(product, facet)
 
         # Report time to compute the reference tensor
         self.cputime = time.time() - t
diff -r 55828e695e14 -r 436ff708086d src/ffc/format/dolfin.py
--- a/src/ffc/format/dolfin.py	Thu Apr  6 11:56:57 2006 -0500
+++ b/src/ffc/format/dolfin.py	Mon May  1 19:31:20 2006 +0200
@@ -462,11 +462,11 @@ public:
 """
 
     # Boundary contribution (if any)
-    if form.AKb.terms:
+    if form.AKb[0].terms:
         eval = __eval_boundary(form, options)
         output += """\
 
-  void eval(real block[], const AffineMap& map, unsigned int boundary) const
+  void eval(real block[], const AffineMap& map, unsigned int facet) const
   {
 %s  }
 """ % eval
@@ -474,7 +474,7 @@ public:
         output += """\
 
   // No contribution from the boundary
-  void eval(real block[], const AffineMap& map, unsigned int boundary) const {}   
+  void eval(real block[], const AffineMap& map, unsigned int facet) const {}   
 """
 
     # Declare class members (if any)
@@ -582,18 +582,26 @@ def __eval_boundary_default(form, option
 """ % "".join(["    const real %s = %s;\n" % (cK.name, cK.value) for cK in form.cK if cK.used])
         output += """\
     // Compute geometry tensors
-%s""" % "".join(["    const real %s = %s;\n" % (gK.name, gK.value) for gK in form.AKb.gK if gK.used])
+%s""" % "".join(["    const real %s = %s;\n" % (gK.name, gK.value) for gK in form.AKb[0].gK if gK.used])
     else:
         output += """\
     // Compute geometry tensors
-%s""" % "".join(["    const real %s = 0.0;\n" % gK.name for gK in form.AKi.gK if gK.used])
+%s""" % "".join(["    const real %s = 0.0;\n" % gK.name for gK in form.AKb[0].gK if gK.used])
 
     if not options["debug-no-element-tensor"]:
         output += """\
 
     // Compute element tensor
-%s""" % "".join(["    %s = %s;\n" % (aK.name, aK.value) for aK in form.AKb.aK])
-
+    switch ( facet )
+    { """
+        for akb in form.AKb:
+          output += """ 
+    case %s:"""  % akb.facet   
+          output += """ 
+%s""" % "".join(["      %s = %s;\n" % (aK.name, aK.value) for aK in akb.aK])
+
+        output += """\
+    } \n"""
     return output
 
 def __eval_boundary_blas(form, options):

---

From: Anders Logg <[email-masked]>
Date: Tue, 2 May 2006 13:39:07 +0200

Looks very nice! I read through the patch and it looks very similar to
what I had in mind. Impressive that you were able to figure out what to
add where.

I just arrived in Sweden yesterday and I won't have access to a computer
(other than for reading email) until Friday when I start at Simula.

I'll take a closer look at the patch then, maybe make some minor changes
and commit it. (Also assuming I have internet access on Svalbard where
I'm going to a meeting for a week starting Saturday.)

When this is added to both FFC and DOLFIN, we can make a new release
of both.

Maybe you can put your FFC tree on fenics.org (in say ~garth/hg/ffc/)
and I can pull from there?

/Anders

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