2010/1/8 <nore...@launchpad.net>:
------------------------------------------------------------ revno: 1426 committer: Anders Logg <l...@simula.no> branch nick: ffc-dev timestamp: Fri 2010-01-08 14:06:48 +0100 message: Use new FIAT interface for quadrature rules and remove ffcquadraturerules.py. Old make_quadrature replaced by create_quadrature now in fiatinterface.py.
If we want to support other quadrature rules should we then implement those in FIAT instead of FFC? The design of ffcfquadraturerules.py was intended to support easy addition of other rules. Kristian
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 (l...@simula.no)" -__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)
signature.asc
Description: OpenPGP digital signature
_______________________________________________ Mailing list: https://launchpad.net/~ffc Post to : ffc@lists.launchpad.net Unsubscribe : https://launchpad.net/~ffc More help : https://help.launchpad.net/ListHelp
