python-ffc 1.0.0-1 / usr / share / pyshared / ffc / quadrature / quadraturerepresentation.py

This file is indexed.

/usr/share/pyshared/ffc/quadrature/quadraturerepresentation.py is in python-ffc 1.0.0-1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below. 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
"Quadrature representation class for UFL"

# Copyright (C) 2009-2010 Kristian B. Oelgaard
#
# This file is part of FFC.
#
# FFC is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# FFC is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with FFC. If not, see <http://www.gnu.org/licenses/>.
#
# Modified by Anders Logg, 2009.
#
# First added:  2009-01-07
# Last changed: 2010-05-18

# UFL modules
from ufl.classes import Form, Integral, SpatialDerivative
from ufl.algorithms import extract_unique_elements, extract_type, extract_elements, propagate_restrictions

# FFC modules
from ffc.log import ffc_assert, info, error
from ffc.fiatinterface import create_element
from ffc.fiatinterface import map_facet_points
from ffc.quadrature.quadraturetransformer import QuadratureTransformer
from ffc.quadrature.optimisedquadraturetransformer import QuadratureTransformerOpt
from ffc.quadrature_schemes import create_quadrature

def compute_integral_ir(domain_type,
                        domain_id,
                        integrals,
                        metadata,
                        form_data,
                        form_id,
                        parameters,
                        common_cell=None):
    "Compute intermediate represention of integral."

    info("Computing quadrature representation")

    # Initialise representation
    num_facets = form_data.num_facets
    ir = {"representation":       "quadrature",
          "domain_type":          domain_type,
          "domain_id":            domain_id,
          "form_id":              form_id,
          "geometric_dimension":  form_data.geometric_dimension,
          "num_facets":           num_facets,
          "geo_consts":           {}}

    # Sort integrals and tabulate basis.
    sorted_integrals = _sort_integrals(integrals, metadata, form_data)
    integrals_dict, psi_tables, quad_weights = \
        _tabulate_basis(sorted_integrals,
                        domain_type,
                        form_data.num_facets,
                        common_cell)

    # Create dimensions of primary indices, needed to reset the argument 'A'
    # given to tabulate_tensor() by the assembler.
    prim_idims = []
    for argument in form_data.arguments:
        element = create_element(argument.element())
        prim_idims.append(element.space_dimension())
    ir["prim_idims"] = prim_idims

    # Create optimise parameters.
    optimise_parameters = {"eliminate zeros":     False,
                           "ignore ones":         False,
                           "remove zero terms":   False,
                           "optimisation":        False,
                           "ignore zero tables":  False}

    if parameters["optimize"]:
        optimise_parameters["ignore ones"]        = True
        optimise_parameters["remove zero terms"]  = True
        optimise_parameters["ignore zero tables"] = True

        # Do not include this in below if/else clause since we want to be
        # able to switch on this optimisation in addition to the other
        # optimisations.
        if "eliminate_zeros" in parameters:
            optimise_parameters["eliminate zeros"] = True

        if "simplify_expressions" in parameters:
            optimise_parameters["optimisation"] = "simplify_expressions"
        elif "precompute_ip_const" in parameters:
            optimise_parameters["optimisation"] = "precompute_ip_const"
        elif "precompute_basis_const" in parameters:
            optimise_parameters["optimisation"] = "precompute_basis_const"
        # The current default optimisation (for -O) is equal to
        # '-feliminate_zeros -fsimplify_expressions'.
        else:
            # If '-O -feliminate_zeros' was given on the command line, do not
            # simplify expressions
            if not "eliminate_zeros" in parameters:
                optimise_parameters["eliminate zeros"] = True
                optimise_parameters["optimisation"]    = "simplify_expressions"

    # Save the optisation parameters.
    ir["optimise_parameters"] = optimise_parameters

    # Create transformer.
    if optimise_parameters["optimisation"]:
        transformer = QuadratureTransformerOpt(psi_tables,
                                               quad_weights,
                                               form_data.geometric_dimension,
                                               optimise_parameters)
    else:
        transformer = QuadratureTransformer(psi_tables,
                                            quad_weights,
                                            form_data.geometric_dimension,
                                            optimise_parameters)

    # Add tables for weights, name_map and basis values.
    ir["quadrature_weights"]  = quad_weights
    ir["name_map"] = transformer.name_map
    ir["unique_tables"] = transformer.unique_tables

    # Transform integrals.
    if domain_type == "cell":
        # Compute transformed integrals.
        info("Transforming cell integral")
        transformer.update_facets(None, None)
        ir["trans_integrals"] = _transform_integrals(transformer, integrals_dict, domain_type)
    elif domain_type == "exterior_facet":
        # Compute transformed integrals.
        terms = [None for i in range(num_facets)]
        for i in range(num_facets):
            info("Transforming exterior facet integral %d" % i)
            transformer.update_facets(i, None)
            terms[i] = _transform_integrals(transformer, integrals_dict, domain_type)
        ir["trans_integrals"] = terms
    elif domain_type == "interior_facet":
        # Compute transformed integrals.
        terms = [[None for j in range(num_facets)] for i in range(num_facets)]
        for i in range(num_facets):
            for j in range(num_facets):
                info("Transforming interior facet integral (%d, %d)" % (i, j))
                transformer.update_facets(i, j)
                terms[i][j] = _transform_integrals(transformer, integrals_dict, domain_type)
        ir["trans_integrals"] = terms
    else:
        error("Unhandled domain type: " + str(domain_type))

    # Save tables map, to extract table names for optimisation option -O.
    ir["psi_tables_map"] = transformer.psi_tables_map
    ir["additional_includes_set"] = transformer.additional_includes_set

    return ir

def _tabulate_basis(sorted_integrals, domain_type, num_facets, common_cell=None):
    "Tabulate the basisfunctions and derivatives."

    # Initialise return values.
    quadrature_weights = {}
    psi_tables = {}
    integrals = {}

    # Loop the quadrature points and tabulate the basis values.
    for pr, integral in sorted_integrals.iteritems():

        # Extract number of points and the rule.
        # TODO: The rule is currently unused because the fiatinterface does not
        # implement support for other rules than those defined in FIAT_NEW
        degree, rule = pr

        # Get all unique elements in integral.
        elements = extract_unique_elements(integral)

        # Create a list of equivalent FIAT elements (with same ordering of elements).
        fiat_elements = [create_element(e) for e in elements]

        # Get cell and facet domains.
        if common_cell is None:
            cell = integral.integrand().cell()
        else:
            cell = common_cell
        cell_domain = cell.domain()
        facet_domain = cell.facet_domain()

        # Make quadrature rule and get points and weights.
        # FIXME: Make create_quadrature() take a rule argument.
        if domain_type == "cell":
            (points, weights) = create_quadrature(cell_domain, degree, rule)
        elif domain_type == "exterior_facet" or domain_type == "interior_facet":
            (points, weights) = create_quadrature(facet_domain, degree, rule)
        else:
            error("Unknown integral type: " + str(domain_type))

        # Add points and rules to dictionary.
        len_weights = len(weights) # The TOTAL number of weights/points
        # TODO: This check should not be needed, remove later.
        ffc_assert(len_weights not in quadrature_weights, \
                    "This number of points is already present in the weight table: " + repr(quadrature_weights))
        quadrature_weights[len_weights] = (weights, points)

        # Add the number of points to the psi tables dictionary.
        # TODO: This check should not be needed, remove later.
        ffc_assert(len_weights not in psi_tables, \
                    "This number of points is already present in the psi table: " + repr(psi_tables))
        psi_tables[len_weights] = {}

        # Add the integral with the number of points as a key to the return integrals.
        # TODO: This check should not be needed, remove later.
        ffc_assert(len_weights not in integrals, \
                    "This number of points is already present in the integrals: " + repr(integrals))
        integrals[len_weights] = integral

        # TODO: This is most likely not the best way to get the highest
        # derivative of an element.
        # Initialise dictionary of elements and the number of derivatives.
        num_derivatives = dict([(e, 0) for e in elements])
        # Extract the derivatives from the integral.
        derivatives = set(extract_type(integral, SpatialDerivative))

        # Loop derivatives and extract multiple derivatives.
        for d in list(derivatives):
            num_deriv = len(extract_type(d, SpatialDerivative))

            # TODO: Safety check, SpatialDerivative only has one operand,
            # and there should be only one element?!
            elem = extract_elements(d.operands()[0])
            ffc_assert(len(elem) == 1, "SpatialDerivative has more than one element: " + repr(elem))
            elem = elem[0]
            # Set the number of derivatives to the highest value
            # encountered so far.
            num_derivatives[elem] = max(num_derivatives[elem], num_deriv)

        # Loop FIAT elements and tabulate basis as usual.
        for i, element in enumerate(fiat_elements):
            # Get order of derivatives.
            deriv_order = num_derivatives[elements[i]]

            # Tabulate for different integral types and insert table into
            # dictionary based on UFL elements.
            if domain_type == "cell":
                psi_tables[len_weights][elements[i]] =\
                {None: element.tabulate(deriv_order, points)}
            elif domain_type == "exterior_facet" or domain_type == "interior_facet":
                psi_tables[len_weights][elements[i]] = {}
                for facet in range(num_facets):
                    psi_tables[len_weights][elements[i]][facet] =\
                        element.tabulate(deriv_order, map_facet_points(points, facet))
            else:
                error("Unknown domain_type: %s" % domain_type)

    return (integrals, psi_tables, quadrature_weights)

def _sort_integrals(integrals, metadata, form_data):
    """Sort integrals according to the number of quadrature points needed per axis.
    Only consider those integrals defined on the given domain."""

    sorted_integrals = {}
    # TODO: We might want to take into account that a form like
    # a = f*g*h*v*u*dx(0, quadrature_order=4) + f*v*u*dx(0, quadrature_order=2),
    # although it involves two integrals of different order, will most
    # likely be integrated faster if one does
    # a = (f*g*h + f)*v*u*dx(0, quadrature_order=4)
    # It will of course only work for integrals defined on the same
    # subdomain and representation.
    for integral in integrals:
        # Get default degree and rule.
        degree = metadata["quadrature_degree"]
        rule  = metadata["quadrature_rule"]
        integral_metadata = integral.measure().metadata()
        # Override if specified in integral metadata
        if not integral_metadata is None:
            if "quadrature_degree" in integral_metadata:
                degree = integral_metadata["quadrature_degree"]
            if "quadrature_rule" in integral_metadata:
                rule = integral_metadata["quadrature_rule"]

        # Create form and add to dictionary according to degree and rule.
        form = Form([Integral(integral.integrand(), integral.measure().reconstruct(metadata={}))])
        if not (degree, rule) in sorted_integrals:
            sorted_integrals[(degree, rule)] = form
        else:
            sorted_integrals[(degree, rule)] += form
    # Extract integrals form forms.
    for key, val in sorted_integrals.items():
        if len(val.integrals()) != 1:
            error("Only expected one integral over one subdomain: %s" % repr(val))
        sorted_integrals[key] = val.integrals()[0]

    return sorted_integrals

def _transform_integrals(transformer, integrals, domain_type):
    "Transform integrals from UFL expression to quadrature representation."
    transformed_integrals = []
    for point, integral in integrals.items():
        transformer.update_points(point)
        integrand = integral.integrand()
        if domain_type == "interior_facet":
            integrand = propagate_restrictions(integrand)
        terms = transformer.generate_terms(integrand)
        transformed_integrals.append((point, terms, transformer.functions, \
                                      {}, transformer.coordinate, transformer.conditionals))
    return transformed_integrals

APT Browse - Built by Thomas Orozco.