Breaks because of a bug in FIAT_NEW for 1D elements (or 2D elements with
facet integrals).
removed:
ffc/ffcquadraturerules.py
modified:
ffc/fiatinterface.py
ffc/quadratureelement.py
ffc/tensor/monomialintegration.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.
=== removed file 'ffc/ffcquadraturerules.py'
--- ffc/ffcquadraturerules.py 2010-01-08 09:49:43 +0000
+++ ffc/ffcquadraturerules.py 1970-01-01 00:00:00 +0000
@@ -1,94 +0,0 @@
-__author__ = "Anders Logg (logg@xxxxxxxxx)"
-__date__ = "2007-11-27"
-__copyright__ = "Copyright (C) 2007-2010 Anders Logg"
-__license__ = "GNU GPL version 3 or any later version"
-
-# Modified by Kristian B. Oelgaard, 2009
-# Last changed: 2010-01-08
-
-# Python modules.
-from numpy import array
-
-# FIXME: KBO: Use modules from NEW
-# FIAT modules.
-from FIAT.quadrature import make_quadrature as fiat_make_quadrature
-from FIAT.shapes import LINE
-from FIAT.shapes import TRIANGLE
-from FIAT.shapes import TETRAHEDRON
-
-# FIXME: KBO: Move to fiatinterface
-# Glue dictionary
-ufl2fiat_shape = {"interval":LINE, "triangle":TRIANGLE, "tetrahedron":TETRAHEDRON}
-
-# FFC modules.
-from log import error
-
-def make_quadrature(shape, n, quad_rule=None):
- """
- Generate quadrature rule (points, weights) for given shape with n
- points in each direction. The quadrature rule is generated using
- FIAT and then transformed and scaled to the UFC element.
- """
-
- # FIXME: Quadrature for vertices (shape == None)
- if shape is None or shape is "vertex":
- return ([()], array([1.0,]))
-
- # Set scaling and transform
- if shape == "interval":
- offset = array((1.0,))
- scaling = 0.5
- elif shape == "triangle":
- offset = array((1.0, 1.0))
- scaling = 0.25
- elif shape == "tetrahedron":
- offset = array((1.0, 1.0, 1.0))
- scaling = 0.125
- else:
- error("Unknown shape")
-
- # Get quadrature from FIAT
- q = fiat_make_quadrature(ufl2fiat_shape[shape], n)
- points = q.get_points()
- weights = q.get_weights()
-
- # Scale from FIAT reference cell to UFC reference cell
- for i in range(len(points)):
- points[i] = tuple(0.5*(array(points[i]) + offset))
- weights[i] *= scaling
-
- return (points, weights)
-
-def map_facet_points(points, facet):
- """
- Map points from the e (UFC) reference simplex of dimension d - 1
- to a given facet on the (UFC) reference simplex of dimension d.
- This may be used to transform points tabulated for example on the
- 2D reference triangle to points on a given facet of the reference
- tetrahedron.
- """
-
- # Special case, don't need to map coordinates on vertices
- dim = len(points[0]) + 1
- if dim == 1: return [(0.0,), (1.0,)][facet]
-
- # Vertex coordinates
- vertex_coordinates = \
- {1: ((0.,), (1.,)),
- 2: ((0., 0.), (1., 0.), (0., 1.)),
- 3: ((0., 0., 0.), (1., 0., 0.),(0., 1., 0.), (0., 0., 1))}
-
- # Facet vertices
- facet_vertices = \
- {2: ((1, 2), (0, 2), (0, 1)),
- 3: ((1, 2, 3), (0, 2, 3), (0, 1, 3), (0, 1, 2))}
-
- # Compute coordinates and map
- coordinates = [vertex_coordinates[dim][v] for v in facet_vertices[dim][facet]]
- new_points = []
- for point in points:
- w = (1.0 - sum(point),) + point
- x = tuple(sum([w[i]*array(coordinates[i]) for i in range(len(w))]))
- new_points += [x]
-
- return new_points
=== modified file 'ffc/fiatinterface.py'
--- ffc/fiatinterface.py 2010-01-07 13:16:49 +0000
+++ ffc/fiatinterface.py 2010-01-08 13:06:48 +0000
@@ -5,7 +5,7 @@
# Modified by Garth N. Wells, 2009.
# Modified by Marie Rognes, 2009-2010.
-# Last changed: 2010-01-07
+# Last changed: 2010-01-08
# UFL modules
from ufl import FiniteElement as UFLFiniteElement
@@ -19,6 +19,7 @@
from FIAT_NEW.lagrange import Lagrange
from FIAT_NEW.brezzi_douglas_marini import BrezziDouglasMarini
from FIAT_NEW.discontinuous_lagrange import DiscontinuousLagrange
+from FIAT_NEW.quadrature import make_quadrature
# FFC modules
from log import debug, error
@@ -40,7 +41,6 @@
"Brezzi-Douglas-Marini": BrezziDouglasMarini,
"Discontinuous Lagrange": DiscontinuousLagrange}
-
# Mapping from dimension to number of mesh sub-entities:
entities_per_dim = {1: [2, 1], 2: [3, 3, 1], 3: [4, 6, 4, 1]}
@@ -79,3 +79,45 @@
element = ElementClass(cell, ufl_element.degree())
return element
+
+def create_quadrature(shape, num_points):
+ """
+ Generate quadrature rule (points, weights) for given shape with
+ num_points points in each direction.
+ """
+ quad_rule = make_quadrature(reference_cell(shape), num_points)
+ return quad_rule.get_points(), quad_rule.get_weights()
+
+def map_facet_points(points, facet):
+ """
+ Map points from the e (UFC) reference simplex of dimension d - 1
+ to a given facet on the (UFC) reference simplex of dimension d.
+ This may be used to transform points tabulated for example on the
+ 2D reference triangle to points on a given facet of the reference
+ tetrahedron.
+ """
+
+ # Special case, don't need to map coordinates on vertices
+ dim = len(points[0]) + 1
+ if dim == 1: return [(0.0,), (1.0,)][facet]
+
+ # Vertex coordinates
+ vertex_coordinates = \
+ {1: ((0.,), (1.,)),
+ 2: ((0., 0.), (1., 0.), (0., 1.)),
+ 3: ((0., 0., 0.), (1., 0., 0.),(0., 1., 0.), (0., 0., 1))}
+
+ # Facet vertices
+ facet_vertices = \
+ {2: ((1, 2), (0, 2), (0, 1)),
+ 3: ((1, 2, 3), (0, 2, 3), (0, 1, 3), (0, 1, 2))}
+
+ # Compute coordinates and map
+ coordinates = [vertex_coordinates[dim][v] for v in facet_vertices[dim][facet]]
+ new_points = []
+ for point in points:
+ w = (1.0 - sum(point),) + point
+ x = tuple(sum([w[i]*array(coordinates[i]) for i in range(len(w))]))
+ new_points += [x]
+
+ return new_points
=== modified file 'ffc/quadratureelement.py'
--- ffc/quadratureelement.py 2009-12-17 20:25:33 +0000
+++ ffc/quadratureelement.py 2010-01-08 13:06:48 +0000
@@ -4,7 +4,7 @@
__license__ = "GNU GPL version 3 or any later version"
# Modified by Garth N. Wells 2006-2009
-# Last changed: 2009-12-16
+# Last changed: 2010-01-08
# Python modules.
import numpy
@@ -20,7 +20,7 @@
from log import error
#from finiteelement import ufl_domain2fiat_domain
#from dofrepresentation import DofRepresentation
-from ffcquadraturerules import make_quadrature
+#from ffcquadraturerules import make_quadrature
#from finiteelement import FiniteElement
#from finiteelement import AFFINE
#from finiteelement import CONTRAVARIANT_PIOLA
=== modified file 'ffc/tensor/monomialintegration.py'
--- ffc/tensor/monomialintegration.py 2010-01-08 09:49:43 +0000
+++ ffc/tensor/monomialintegration.py 2010-01-08 13:06:48 +0000
@@ -23,8 +23,7 @@
# FFC modules
from ffc.log import info, debug, error
-from ffc.ffcquadraturerules import make_quadrature, map_facet_points
-from ffc.fiatinterface import create_element
+from ffc.fiatinterface import create_element, create_quadrature
# FFC tensor representation modules.
from multiindex import build_indices
@@ -73,9 +72,9 @@
# Create quadrature rule and get points and weights
if domain_type == Measure.CELL:
- (points, weights) = make_quadrature(cell_shape, num_points)
+ (points, weights) = create_quadrature(cell_shape, num_points)
else:
- (points, weights) = make_quadrature(facet_shape, num_points)
+ (points, weights) = create_quadrature(facet_shape, num_points)
return (points, weights)
