ffc team mailing list archive
-
ffc team
-
Mailing list archive
-
Message #00574
Re: integration over boundaries
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
On Mon, May 01, 2006 at 07:40:00PM +0200, Garth N. Wells wrote:
> 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):
> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/ffc-dev
References
