ffc team mailing list archive
-
ffc team
-
Mailing list archive
-
Message #00572
integration over boundaries
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 <g.n.wells@xxxxxxxxxx>"
# 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):
Follow ups
