~ubuntu-branches/ubuntu/trusty/ffc/trusty

« back to all changes in this revision

Viewing changes to ffc/evaluatebasisderivatives.py

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2010-02-03 20:22:35 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20100203202235-fe8d0kajuvgy2sqn
Tags: 0.9.0-1
* New upstream release.
* debian/control: Bump Standards-Version (no changes needed).
* Update debian/copyright and debian/copyright_hints.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
"""Code generation for evaluation of derivatives of finite element basis values.
 
2
This module generates code which is more or less a C++ representation of the code
 
3
found in FIAT_NEW."""
 
4
 
 
5
__author__ = "Kristian B. Oelgaard (k.b.oelgaard@gmail.com)"
 
6
__date__ = "2007-04-16"
 
7
__copyright__ = "Copyright (C) 2007-2010 Kristian B. Oelgaard"
 
8
__license__  = "GNU GPL version 3 or any later version"
 
9
 
 
10
# Last changed: 2010-01-30
 
11
 
 
12
# Python modules
 
13
import math
 
14
import numpy
 
15
 
 
16
# FFC modules
 
17
from ffc.log import error, ffc_assert
 
18
from ffc.evaluatebasis import _map_dof
 
19
from ffc.evaluatebasis import _compute_basisvalues
 
20
from ffc.evaluatebasis import _tabulate_coefficients
 
21
from ffc.cpp import tabulate_matrix, IndentControl, remove_unused, tabulate_vector
 
22
from ffc.quadrature.quadraturegenerator_utils import generate_loop
 
23
#from ffc.quadrature.symbolics import set_format
 
24
 
 
25
# Temporary import
 
26
#from cpp import format_old as format
 
27
from ffc.cpp import format
 
28
 
 
29
def _evaluate_basis_derivatives_all(data_list):
 
30
    """Like evaluate_basis, but return the values of all basis functions (dofs)."""
 
31
 
 
32
    if isinstance(data_list, str):
 
33
        return format["exception"]("evaluate_basis_derivatives_all: %s" % data_list)
 
34
 
 
35
    f_r, f_s = format["free indices"][:2]
 
36
    f_assign = format["assign"]
 
37
 
 
38
    # Initialise objects
 
39
    Indent = IndentControl()
 
40
    code = []
 
41
 
 
42
    # FIXME: KBO: Figure out which return format to use, either:
 
43
    # [dN0[0]/dx, dN0[0]/dy, dN0[1]/dx, dN0[1]/dy, dN1[0]/dx, dN1[0]/dy, dN1[1]/dx, dN1[1]/dy, ...]
 
44
    # or
 
45
    # [dN0[0]/dx, dN1[0]/dx, ..., dN0[1]/dx, dN1[1]/dx, ..., dN0[0]/dy, dN1[0]/dy, ..., dN0[1]/dy, dN1[1]/dy, ...]
 
46
    # or
 
47
    # [dN0[0]/dx, dN0[1]/dx, ..., dN1[0]/dx, dN1[1]/dx, ..., dN0[0]/dy, dN0[1]/dy, ..., dN1[0]/dy, dN1[1]/dy, ...]
 
48
    # for vector (tensor elements), currently returning option 1.
 
49
 
 
50
    # FIXME: KBO: For now, just call evaluate_basis_derivatives and map values
 
51
    # accordingly, this will keep the amount of code at a minimum. If it turns
 
52
    # out that speed is an issue (overhead from calling evaluate_basis), we can
 
53
    # easily generate all the code.
 
54
 
 
55
    # Get total value shape and space dimension for entire element (possibly mixed).
 
56
    value_shape = sum(sum(data["value_shape"] or (1,)) for data in data_list)
 
57
    space_dimension = sum(data["space_dimension"] for data in data_list)
 
58
 
 
59
    # Special case where space dimension is one (constant elements)
 
60
    if space_dimension == 1:
 
61
        code += [format["comment"]("Element is constant, calling evaluate_basis_derivatives.")]
 
62
        code += ["evaluate_basis_derivatives(0, n, %s, coordinates, c);" % format["argument values"]]
 
63
        return "\n".join(code)
 
64
 
 
65
    # Compute number of derivatives
 
66
    # Get the topological dimension.
 
67
    topological_dimension = data_list[0]["topological_dimension"]
 
68
    code += _compute_num_derivatives(topological_dimension, Indent, format)
 
69
 
 
70
    # Declare helper value to hold single dof values and reset
 
71
    code += ["", format["comment"]("Helper variable to hold values of a single dof.")]
 
72
    if (value_shape == 1):
 
73
        num_vals = format["num derivatives"]
 
74
    else:
 
75
        # FIXME: KBO: Should the str(int()) be in format?
 
76
        num_vals = format["multiply"]([str(int(value_shape)), format["num derivatives"]])
 
77
    code += [f_assign(Indent.indent(format["float declaration"] + format["pointer"] + "dof_values"),\
 
78
              format["component"](format["new"] + format["float declaration"], num_vals))]
 
79
    loop_vars = [(f_r, 0, num_vals)]
 
80
    line = [f_assign(format["component"]("dof_values", f_r), format["floating point"](0.0))]
 
81
    code += generate_loop(line, loop_vars, Indent, format)
 
82
 
 
83
    # Create loop over dofs that calls evaluate_basis_derivatives for a single dof and
 
84
    # inserts the values into the global array.
 
85
    code += ["", format["comment"]("Loop dofs and call evaluate_basis_derivatives.")]
 
86
    lines_r = []
 
87
    loop_vars_r = [(f_r, 0, space_dimension)]
 
88
 
 
89
    # FIXME: KBO: Move evaluate_basis_derivatives string to cpp.py
 
90
    lines_r += ["evaluate_basis_derivatives(%s, n, dof_values, coordinates, c);" % f_r]
 
91
 
 
92
    loop_vars_s = [(f_s, 0, num_vals)]
 
93
    index = format["add"]([format["multiply"]([f_r, num_vals]), f_s])
 
94
    name = format["component"](format["argument values"], index)
 
95
    value = format["component"]("dof_values", f_s)
 
96
    lines_s = [f_assign(name, value)]
 
97
    lines_r += generate_loop(lines_s, loop_vars_s, Indent, format)
 
98
    code += generate_loop(lines_r, loop_vars_r, Indent, format)
 
99
 
 
100
    code += ["", format["comment"]("Delete pointer.")]
 
101
    code += [Indent.indent(format["delete pointer"]("dof_values", ""))]
 
102
 
 
103
    # Generate bode (no need to remove unused)
 
104
    return "\n".join(code)
 
105
 
 
106
def _evaluate_basis_derivatives(data_list):
 
107
    """Evaluate the derivatives of an element basisfunction at a point. The values are
 
108
    computed as in FIAT as the dot product of the coefficients (computed at compile time)
 
109
    and basisvalues which are dependent on the coordinate and thus have to be computed at
 
110
    run time.
 
111
 
 
112
    Currently the following elements are supported in 1D:
 
113
 
 
114
    Lagrange                + mixed/vector valued
 
115
    Discontinuous Lagrange  + mixed/vector valued
 
116
 
 
117
    Currently the following elements are supported in 2D and 3D:
 
118
 
 
119
    Lagrange                + mixed/vector valued
 
120
    Discontinuous Lagrange  + mixed/vector valued
 
121
    Crouzeix-Raviart        + mixed/vector valued
 
122
    Brezzi-Douglas-Marini   + mixed/vector valued
 
123
 
 
124
    Not supported in 2D or 3D:
 
125
    Raviart-Thomas ? (not tested since it is broken in FFC, but should work)
 
126
    Nedelec (broken?)"""
 
127
 
 
128
    if isinstance(data_list, str):
 
129
        return format["exception"]("evaluate_basis_derivatives: %s" % data_list)
 
130
 
 
131
    # Init return code and indent object
 
132
    code = []
 
133
    Indent = IndentControl()
 
134
 
 
135
    # Get the element cell domain, geometric and topological dimension.
 
136
    element_cell_domain = data_list[0]["cell_domain"]
 
137
    geometric_dimension = data_list[0]["geometric_dimension"]
 
138
    topological_dimension = data_list[0]["topological_dimension"]
 
139
 
 
140
    # Get code snippets for Jacobian, Inverse of Jacobian and mapping of
 
141
    # coordinates from physical element to the FIAT reference element.
 
142
    # FIXME: KBO: Change this when supporting R^2 in R^3 elements.
 
143
    code += [Indent.indent(format["jacobian and inverse"](geometric_dimension))]
 
144
    code += ["", Indent.indent(format["fiat coordinate map"](element_cell_domain))]
 
145
 
 
146
    # Compute number of derivatives that has to be computed, and declare an array to hold
 
147
    # the values of the derivatives on the reference element
 
148
    code += [""]
 
149
    code += _compute_num_derivatives(topological_dimension, Indent, format)
 
150
 
 
151
    # Generate all possible combinations of derivatives
 
152
    code += _generate_combinations(topological_dimension, Indent, format)
 
153
 
 
154
    # Generate the transformation matrix
 
155
    code += _generate_transform(element_cell_domain, Indent, format)
 
156
 
 
157
    # Reset all values
 
158
    code += _reset_values(data_list, Indent, format)
 
159
 
 
160
    if len(data_list) == 1:
 
161
        data = data_list[0]
 
162
 
 
163
        # Map degree of freedom to local degree.
 
164
        code += [_map_dof(0, Indent, format)]
 
165
 
 
166
        # Generate element code.
 
167
        code += _generate_element_code(data, 0, Indent, format)
 
168
 
 
169
    # If the element is of type MixedElement (including Vector- and TensorElement).
 
170
    else:
 
171
        # Generate element code, for all sub-elements.
 
172
        code += _mixed_elements(data_list, Indent, format)
 
173
    code = remove_unused("\n".join(code))
 
174
    return code
 
175
 
 
176
def _compute_num_derivatives(topological_dimension, Indent, format):
 
177
    "Computes the number of derivatives of order 'n' as: element.cell_shape()^n."
 
178
 
 
179
    code = [format["comment"]("Compute number of derivatives.")]
 
180
    # FIXME: KBO: Should the str(int()) be in format?
 
181
    code.append(format["declaration"](format["uint declaration"], format["num derivatives"], str(int(1))))
 
182
 
 
183
    loop_vars = [(format["free indices"][0], 0, format["argument derivative order"])]
 
184
    # FIXME: KBO: Should the str(int()) be in format?
 
185
    lines = [format["imul"](format["num derivatives"], str(int(topological_dimension)))]
 
186
 
 
187
    code += generate_loop(lines, loop_vars, Indent, format)
 
188
 
 
189
    return code
 
190
 
 
191
def _generate_combinations(topological_dimension, Indent, format):
 
192
    "Generate all possible combinations of derivatives of order 'n'"
 
193
 
 
194
    # Use code from format
 
195
    code = ["", Indent.indent(format["combinations"]\
 
196
            % {"combinations": format["derivative combinations"],\
 
197
               "topological_dimension-1": topological_dimension-1,\
 
198
               "num_derivatives" : format["num derivatives"],\
 
199
               "n": format["argument derivative order"]})]
 
200
    return code
 
201
 
 
202
def _generate_transform(element_cell_domain, Indent, format):
 
203
    """Generate the transformation matrix, whic is used to transform derivatives from reference
 
204
    element back to the physical element."""
 
205
 
 
206
    # Generate code to construct the inverse of the Jacobian
 
207
    if (element_cell_domain in ["interval", "triangle", "tetrahedron"]):
 
208
        code = ["", Indent.indent(format["transform snippet"][element_cell_domain]\
 
209
        % {"transform": format["transform matrix"],\
 
210
           "num_derivatives" : format["num derivatives"],\
 
211
           "n": format["argument derivative order"],\
 
212
           "combinations": format["derivative combinations"],\
 
213
           "K":format["transform Jinv"]})]
 
214
    else:
 
215
        error("Cannot generate transform for shape: %s" % element_cell_domain)
 
216
 
 
217
    return code
 
218
 
 
219
def _reset_values(data_list, Indent, format):
 
220
    "Reset all components of the 'values' array as it is a pointer to an array."
 
221
    f_assign = format["assign"]
 
222
    code = ["", Indent.indent(format["comment"]("Reset values. Assuming that values is always an array."))]
 
223
 
 
224
    # Get value shape and reset values. This should also work for TensorElement,
 
225
    # scalar are empty tuples, therefore (1,) in which case value_shape = 1.
 
226
    value_shape = sum(sum(data["value_shape"] or (1,)) for data in data_list)
 
227
 
 
228
    # Only multiply by value shape if different from 1.
 
229
    if value_shape == 1:
 
230
        num_vals = format["num derivatives"]
 
231
    else:
 
232
        # FIXME: KBO: Should the str(int()) be in format?
 
233
        num_vals = format["multiply"]([str(int(value_shape)), format["num derivatives"]])
 
234
    name = format["component"](format["argument values"], format["free indices"][0])
 
235
    loop_vars = [(format["free indices"][0], 0, num_vals)]
 
236
    lines = [f_assign(name, format["floating point"](0))]
 
237
    code += generate_loop(lines, loop_vars, Indent, format)
 
238
 
 
239
    return code + [""]
 
240
 
 
241
def _generate_element_code(data, sum_value_dim, Indent, format):
 
242
    "Generate code for each basis element"
 
243
 
 
244
    code = []
 
245
 
 
246
    # Compute basisvalues, from evaluatebasis.py
 
247
    code += _compute_basisvalues(data, Indent, format)
 
248
 
 
249
    # Tabulate coefficients
 
250
    code += _tabulate_coefficients(data, Indent, format)
 
251
 
 
252
    # Tabulate coefficients for derivatives
 
253
    code += _tabulate_dmats(data, Indent, format)
 
254
 
 
255
    # Compute the derivatives of the basisfunctions on the reference (FIAT) element,
 
256
    # as the dot product of the new coefficients and basisvalues
 
257
    code += _compute_reference_derivatives(data, Indent, format)
 
258
 
 
259
    # Transform derivatives to physical element by multiplication with the transformation matrix
 
260
    code += _transform_derivatives(data, sum_value_dim, Indent, format)
 
261
 
 
262
    # Delete pointers
 
263
    code += _delete_pointers(data, Indent, format)
 
264
 
 
265
    return code
 
266
 
 
267
def _mixed_elements(data_list, Indent, format):
 
268
    "Generate code for each sub-element in the event of vector valued elements or mixed elements"
 
269
 
 
270
    # Prefetch formats to speed up code generation
 
271
    f_dof_map_if = format["dof map if"]
 
272
    f_if         = format["if"]
 
273
 
 
274
    sum_value_dim = 0
 
275
    sum_space_dim = 0
 
276
 
 
277
    # Init return code.
 
278
    code = []
 
279
 
 
280
    # Generate code for each element
 
281
    for data in data_list:
 
282
 
 
283
        # Get value and space dimension (should be tensor ready).
 
284
        value_dim = sum(data["value_shape"] or (1,))
 
285
        space_dim = data["space_dimension"]
 
286
 
 
287
        # Generate map from global to local dof
 
288
        element_code = [_map_dof(sum_space_dim, Indent, format)]
 
289
 
 
290
        # Generate code for basis element
 
291
        element_code += _generate_element_code(data, sum_value_dim, Indent, format)
 
292
 
 
293
        # Increase indentation, indent code and decrease indentation.
 
294
        Indent.increase()
 
295
        if_code = remove_unused(Indent.indent("\n".join(element_code)))
 
296
        # if_code = Indent.indent("\n".join(element_code))
 
297
        Indent.decrease()
 
298
 
 
299
        # Create if statement and add to code.
 
300
        code += [f_if(f_dof_map_if(sum_space_dim, sum_space_dim + space_dim -1), if_code)]
 
301
 
 
302
        # Increase sum of value dimension, and space dimension
 
303
        sum_value_dim += value_dim
 
304
        sum_space_dim += space_dim
 
305
 
 
306
    return code
 
307
 
 
308
def _tabulate_dmats(data, Indent, format):
 
309
    "Tabulate the derivatives of the polynomial base"
 
310
 
 
311
    code = []
 
312
 
 
313
    # Prefetch formats to speed up code generation
 
314
    f_table          = format["static const float declaration"]
 
315
    f_dmats          = format["dmats"]
 
316
    f_component  = format["component"]
 
317
    f_assign = format["assign"]
 
318
 
 
319
    # Get derivative matrices (coefficients) of basis functions, computed by FIAT at compile time
 
320
    derivative_matrices = data["dmats"]
 
321
 
 
322
    code += [Indent.indent(format["comment"]("Tables of derivatives of the polynomial base (transpose)."))]
 
323
 
 
324
    # Generate tables for each spatial direction
 
325
    for i, dmat in enumerate(derivative_matrices):
 
326
 
 
327
        # Extract derivatives for current direction (take transpose, FIAT_NEW PolynomialSet.tabulate())
 
328
        matrix = numpy.transpose(dmat)
 
329
 
 
330
        # Get shape and check dimension (This is probably not needed)
 
331
        shape = numpy.shape(matrix)
 
332
        ffc_assert(shape[0] == shape[1] == data["num_expansion_members"], "Something is wrong with the shape of dmats.")
 
333
 
 
334
        # Declare varable name for coefficients
 
335
        name = f_component(f_table + f_dmats(i), [shape[0], shape[1]])
 
336
        value = tabulate_matrix(matrix, format)
 
337
        code += [f_assign(Indent.indent(name), Indent.indent(value)), ""]
 
338
 
 
339
    return code
 
340
 
 
341
def _reset_dmats(shape_dmats, indices, Indent, format):
 
342
    f_assign = format["assign"]
 
343
    code = [format["comment"]("Resetting dmats values to compute next derivative.")]
 
344
 
 
345
    loop_vars = [(indices[0], 0, shape_dmats[0]), (indices[1], 0, shape_dmats[1])]
 
346
    dmats_old = format["component"](format["dmats"](""), [indices[0], indices[1]])
 
347
    lines = [f_assign(dmats_old, format["floating point"](0.0))]
 
348
    lines += [format["if"](indices[0] + format["is equal"] + indices[1],\
 
349
              format["assign"](Indent.indent(dmats_old), format["floating point"](1.0)))]
 
350
    code += generate_loop(lines, loop_vars, Indent, format)
 
351
    return code
 
352
 
 
353
def _update_dmats(shape_dmats, indices, Indent, format):
 
354
    f_assign = format["assign"]
 
355
    code = [format["comment"]("Updating dmats_old with new values and resetting dmats.")]
 
356
    dmats = format["component"](format["dmats"](""), [indices[0], indices[1]])
 
357
    dmats_old = format["component"](format["dmats old"], [indices[0], indices[1]])
 
358
    loop_vars = [(indices[0], 0, shape_dmats[0]), (indices[1], 0, shape_dmats[1])]
 
359
    lines = [f_assign(dmats_old, dmats), f_assign(dmats, format["floating point"](0.0))]
 
360
    code += generate_loop(lines, loop_vars, Indent, format)
 
361
    return code
 
362
 
 
363
def _compute_dmats(num_dmats, shape_dmats, available_indices, deriv_index, Indent, format):
 
364
 
 
365
    s, t, u = available_indices
 
366
 
 
367
    # Reset dmats_old
 
368
    code = _reset_dmats(shape_dmats, [t, u], Indent, format)
 
369
    code += ["", format["comment"]("Looping derivative order to generate dmats.")]
 
370
 
 
371
    # Set dmats matrix equal to dmats_old
 
372
    lines = _update_dmats(shape_dmats, [t, u], Indent, format)
 
373
    loop_vars = [(s, 0, format["argument derivative order"])]
 
374
 
 
375
    lines += ["", format["comment"]("Update dmats using an inner product.")]
 
376
    # Create dmats matrix by multiplication
 
377
    comb = format["component"](format["derivative combinations"], [deriv_index, s])
 
378
    for i in range(num_dmats):
 
379
        lines += _dmats_product(shape_dmats, comb, i, [t, u], Indent, format)
 
380
 
 
381
    code += generate_loop(lines, loop_vars, Indent, format)
 
382
 
 
383
    return code
 
384
 
 
385
def _dmats_product(shape_dmats, index, i, indices, Indent, format):
 
386
 
 
387
    t, u = indices
 
388
    tu = t + u
 
389
    loop_vars = [(t, 0, shape_dmats[0]), (u, 0, shape_dmats[1])]
 
390
    dmats = format["component"](format["dmats"](""), [t, u])
 
391
    dmats_old = format["component"](format["dmats old"], [tu, u])
 
392
    value = format["multiply"]([format["component"](format["dmats"](i), [t, tu]), dmats_old])
 
393
    name = Indent.indent(format["iadd"](dmats, value))
 
394
    lines = generate_loop([name], [(tu, 0, shape_dmats[0])], Indent, format)
 
395
    code = [format["if"](index + format["is equal"] + str(i),\
 
396
            "\n".join(generate_loop(lines, loop_vars, Indent, format)))]
 
397
 
 
398
    return code
 
399
 
 
400
def _matrix_index(i, j, range_j):
 
401
    "Compute the index in an array from indices in a matrix i.e., m[i][j] -> a[i*range(j)+j]"
 
402
    if i == 0:
 
403
        access = j
 
404
    elif i == 1:
 
405
        access = format["add"]([range_j, j])
 
406
    else:
 
407
        # FIXME: KBO: Should the str(int()) be in format?
 
408
        irj = format["mul"]([str(i), range_j])
 
409
        access = format["add"]([irj, j])
 
410
    return access
 
411
 
 
412
def _compute_reference_derivatives(data, Indent, format):
 
413
    """Compute derivatives on the reference element by recursively multiply coefficients with
 
414
    the relevant derivatives of the polynomial base until the requested order of derivatives
 
415
    has been reached. After this take the dot product with the basisvalues."""
 
416
 
 
417
    code = []
 
418
 
 
419
    # Prefetch formats to speed up code generation
 
420
    f_comment    = format["comment"]
 
421
    f_float      = format["float declaration"]
 
422
    f_component  = format["component"]
 
423
    f_tmp        = format["tmp ref value"]
 
424
    f_dmats      = format["dmats"]
 
425
    f_assign     = format["assign"]
 
426
    f_iadd       = format["iadd"]
 
427
 
 
428
    f_r, f_s, f_t, f_u = format["free indices"]
 
429
 
 
430
    # Get number of components, change for tensor valued elements.
 
431
    shape = data["value_shape"]
 
432
    if shape == ():
 
433
        num_components = 1
 
434
    elif len(shape) == 1:
 
435
        num_components = shape[0]
 
436
    else:
 
437
        error("Tensor valued elements are not supported yet: %s" % data["family"])
 
438
 
 
439
    # Get shape of derivative matrix (they should all have the same shape).
 
440
    shape_dmats = numpy.shape(data["dmats"][0])
 
441
 
 
442
    code += [Indent.indent(f_comment("Compute reference derivatives"))]
 
443
 
 
444
    # Declare pointer to array that holds derivatives on the FIAT element
 
445
    code += [Indent.indent(f_comment("Declare pointer to array of derivatives on FIAT element"))]
 
446
    # The size of the array of reference derivatives is equal to the number of derivatives
 
447
    # times the number of components of the basis element
 
448
    if (num_components == 1):
 
449
        num_vals = format["num derivatives"]
 
450
    else:
 
451
        # FIXME: KBO: Should the str(int()) be in format?
 
452
        num_vals = format["multiply"]([str(int(num_components)), format["num derivatives"]])
 
453
    code += [f_assign(Indent.indent(f_float + format["pointer"] + format["reference derivatives"]),\
 
454
              format["component"](format["new"] + f_float, num_vals))]
 
455
    # Reset values of reference derivatives.
 
456
    name = format["component"](format["reference derivatives"], f_r)
 
457
    lines = [f_assign(name, format["floating point"](0))]
 
458
    code += generate_loop(lines, [(f_r, 0, num_vals)], Indent, format)
 
459
 
 
460
    code += [""]
 
461
 
 
462
    # Declare matrix of dmats (which will hold the matrix product of all combinations)
 
463
    # and dmats_old which is needed in order to perform the matrix product.
 
464
    code += [Indent.indent(f_comment("Declare derivative matrix (of polynomial basis)."))]
 
465
    matrix = numpy.eye(shape_dmats[0])
 
466
    name = f_component(f_float + f_dmats(""), [shape_dmats[0], shape_dmats[1]])
 
467
    value = tabulate_matrix(matrix, format)
 
468
    code += [f_assign(Indent.indent(name), Indent.indent(value)), ""]
 
469
    code += [Indent.indent(f_comment("Declare (auxiliary) derivative matrix (of polynomial basis)."))]
 
470
    name = f_component(f_float + format["dmats old"], [shape_dmats[0], shape_dmats[1]])
 
471
    code += [f_assign(Indent.indent(name), Indent.indent(value)), ""]
 
472
 
 
473
    # Loop all derivatives and compute value of the derivative as:
 
474
    # deriv_on_ref[r] = coeff[dof][s]*dmat[s][t]*basis[t]
 
475
    code += [Indent.indent(f_comment("Loop possible derivatives."))]
 
476
    loop_vars = [(f_r, 0, format["num derivatives"])]
 
477
 
 
478
    # Compute dmats as a recursive matrix product
 
479
    lines = _compute_dmats(len(data["dmats"]), shape_dmats, [f_s, f_t, f_u], f_r, Indent, format)
 
480
 
 
481
    # Compute derivatives for all components
 
482
    loop_vars_c = [(f_s, 0, shape_dmats[0]),(f_t, 0, shape_dmats[1])]
 
483
    lines_c = []
 
484
    for i in range(num_components):
 
485
        comp_access = _matrix_index(i, f_r, format["num derivatives"])
 
486
        name = format["component"](format["reference derivatives"], comp_access)
 
487
        coeffs = format["component"](format["coefficients"](i), [format["local dof"], f_s])
 
488
        dmats = format["component"](format["dmats"](""), [f_s, f_t])
 
489
        basis = format["component"](format["basisvalues"], f_t)
 
490
        value = format["multiply"]([coeffs, dmats, basis])
 
491
        lines_c.append(f_iadd(name, value))
 
492
    lines += generate_loop(lines_c, loop_vars_c, Indent, format)
 
493
 
 
494
    # Apply transformation if applicable.
 
495
    mapping = data["mapping"]
 
496
    if mapping == "affine":
 
497
        pass
 
498
    elif mapping == "contravariant piola":
 
499
        lines += ["", Indent.indent(format["comment"]\
 
500
                ("Using contravariant Piola transform to map values back to the physical element"))]
 
501
        # Get temporary values before mapping.
 
502
        lines += [format["const float declaration"](Indent.indent(f_tmp(i)),\
 
503
                  f_component(format["reference derivatives"], _matrix_index(i, f_r, format["num derivatives"]))) for i in range(num_components)]
 
504
 
 
505
        # Create names for inner product.
 
506
        topological_dimension = data["topological_dimension"]
 
507
        basis_col = [f_tmp(j) for j in range(topological_dimension)]
 
508
        for i in range(num_components):
 
509
            # Create Jacobian.
 
510
            jacobian_row = [format["transform"]("J", i, j, None) for j in range(topological_dimension)]
 
511
 
 
512
            # Create inner product and multiply by inverse of Jacobian.
 
513
            inner = [format["multiply"]([jacobian_row[j], basis_col[j]]) for j in range(topological_dimension)]
 
514
            sum_ = format["grouping"](format["add"](inner))
 
515
            value = format["multiply"]([format["inverse"](format["det(J)"]("")), sum_])
 
516
            name = f_component(format["reference derivatives"], _matrix_index(i, f_r, format["num derivatives"]))
 
517
            lines += [f_assign(name, value)]
 
518
    elif mapping == "covariant piola":
 
519
        lines += ["", Indent.indent(format["comment"]\
 
520
                ("Using covariant Piola transform to map values back to the physical element"))]
 
521
        # Get temporary values before mapping.
 
522
        lines += [format["const float declaration"](Indent.indent(f_tmp(i)),\
 
523
                  f_component(format["reference derivatives"], _matrix_index(i, f_r, format["num derivatives"]))) for i in range(num_components)]
 
524
        # Create names for inner product.
 
525
        topological_dimension = data["topological_dimension"]
 
526
        basis_col = [f_tmp(j) for j in range(topological_dimension)]
 
527
        for i in range(num_components):
 
528
            # Create inverse of Jacobian.
 
529
            inv_jacobian_column = [format["transform"]("JINV", j, i, None) for j in range(topological_dimension)]
 
530
 
 
531
            # Create inner product of basis values and inverse of Jacobian.
 
532
            inner = [format["multiply"]([inv_jacobian_column[j], basis_col[j]]) for j in range(topological_dimension)]
 
533
            value = format["grouping"](format["add"](inner))
 
534
            name = f_component(format["reference derivatives"], _matrix_index(i, f_r, format["num derivatives"]))
 
535
            lines += [f_assign(name, value)]
 
536
    else:
 
537
        error("Unknown mapping: %s" % mapping)
 
538
 
 
539
    # Generate loop over number of derivatives.
 
540
    code += generate_loop(lines, loop_vars, Indent, format)
 
541
 
 
542
    return code + [""]
 
543
 
 
544
def _transform_derivatives(data, sum_value_dim, Indent, format):
 
545
    """"""
 
546
 
 
547
    code = []
 
548
 
 
549
    # Prefetch formats to speed up code generation
 
550
    f_loop            = format["loop"]
 
551
    f_num_derivatives = format["num derivatives"]
 
552
    f_derivatives     = format["reference derivatives"]
 
553
    f_values          = format["argument values"]
 
554
    f_mul             = format["mul"]
 
555
    f_add             = format["add"]
 
556
    f_component       = format["component"]
 
557
    f_transform       = format["transform matrix"]
 
558
 
 
559
    code += [Indent.indent(format["comment"]("Transform derivatives back to physical element"))]
 
560
 
 
561
    code += [Indent.indent(f_loop("row", 0, f_num_derivatives))]
 
562
    code += [Indent.indent(format["block begin"])]
 
563
    # Increase indentation
 
564
    Indent.increase()
 
565
    code += [Indent.indent(f_loop("col", 0, f_num_derivatives))]
 
566
    code += [Indent.indent(format["block begin"])]
 
567
    # Increase indentation
 
568
    Indent.increase()
 
569
 
 
570
    # Get number of components, change for tensor valued elements.
 
571
    shape = data["value_shape"]
 
572
    if shape == ():
 
573
        num_components = 1
 
574
    elif len(shape) == 1:
 
575
        num_components = shape[0]
 
576
    else:
 
577
        error("Tensor valued elements are not supported yet: %s" % data["family"])
 
578
 
 
579
    # Compute offset in array values if any
 
580
    for i in range(num_components):
 
581
        access_name = _matrix_index(sum_value_dim + i, "row", f_num_derivatives)
 
582
        name = f_component(f_values, access_name)
 
583
#        if (sum_value_dim + i == 0):
 
584
#            name = f_component(f_values, "row")
 
585
#        elif (sum_value_dim + i == 1):
 
586
#            name = f_component(f_values, f_add([f_num_derivatives, "row"]))
 
587
#        else:
 
588
#            offset_name = f_mul(["%d" %(sum_value_dim + i), f_num_derivatives])
 
589
#            name = f_component(f_values, f_add([offset_name, "row"]))
 
590
        access_val = _matrix_index(i, "col", f_num_derivatives)
 
591
        value = f_mul([f_component(f_transform, ["row","col"]),\
 
592
                    f_component(f_derivatives, access_val)])
 
593
#        if (i == 0):
 
594
#            value = f_mul([f_component(f_transform, ["row","col"]),\
 
595
#                    f_component(f_derivatives, "col")])
 
596
#        elif (i == 1):
 
597
#            value = f_mul([f_component(f_transform, ["row","col"]),\
 
598
#                    f_component(f_derivatives, f_add([f_num_derivatives, "col"]))])
 
599
#        else:
 
600
#            offset_value = f_mul(["%d" %i, f_num_derivatives])
 
601
 
 
602
#            value = f_mul([f_component(f_transform, ["row","col"]),\
 
603
#                    f_component(f_derivatives, f_add([offset_value, "col"]))])
 
604
 
 
605
        # Compute values
 
606
        code += [Indent.indent(format["iadd"](name, value))]
 
607
 
 
608
    # Decrease indentation
 
609
    Indent.decrease()
 
610
    code += [Indent.indent(format["block end"])]
 
611
 
 
612
    # Decrease indentation
 
613
    Indent.decrease()
 
614
 
 
615
    code += [Indent.indent(format["block end"])]
 
616
 
 
617
    return code
 
618
 
 
619
def _delete_pointers(data, Indent, format):
 
620
    "Delete the pointers to arrays."
 
621
 
 
622
    code = []
 
623
    f_r = format["free indices"][0]
 
624
 
 
625
    code += ["", Indent.indent(format["comment"]("Delete pointer to array of derivatives on FIAT element"))]
 
626
    code += [Indent.indent(format["delete pointer"](format["reference derivatives"], "")), ""]
 
627
 
 
628
    code += [Indent.indent(format["comment"]("Delete pointer to array of combinations of derivatives and transform"))]
 
629
    loop_vars = [(f_r, 0, format["num derivatives"])]
 
630
    lines =  [format["delete pointer"](format["derivative combinations"], format["component"]("", f_r))]
 
631
    lines += [format["delete pointer"](format["transform matrix"], format["component"]("", f_r))]
 
632
    code += generate_loop(lines, loop_vars, Indent, format)
 
633
 
 
634
    code += [Indent.indent(format["delete pointer"](format["derivative combinations"], ""))]
 
635
    code += [Indent.indent(format["delete pointer"](format["transform matrix"], ""))]
 
636
 
 
637
    return code
 
638