This file is indexed.

/usr/lib/python2.7/dist-packages/ffc/cpp.py is in python-ffc 2016.2.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
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
# -*- coding: utf-8 -*-
"This module defines rules and algorithms for generating C++ code."

# Copyright (C) 2009-2016 Anders Logg
#
# 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 Kristian B. Oelgaard 2011
# Modified by Marie E. Rognes 2010
# Modified by Martin Sandve Alnæs 2013-2016

# Python modules
import re
import numpy
import platform
from six import string_types

# UFL modules
from ufl import custom_integral_types

# FFC modules
from ffc.log import debug, error
from six.moves import zip

# ufc class names

def make_classname(prefix, basename, signature):
    pre = prefix.lower() + "_" if prefix else ""
    sig = str(signature).lower()
    return "%s%s_%s" % (pre, basename, sig)


def make_integral_classname(prefix, integral_type, form_id, subdomain_id):
    basename = "%s_integral_%s" % (integral_type, str(form_id).lower())
    return make_classname(prefix, basename, subdomain_id)


# Mapping of restrictions
_fixed_map = {None: "", "+": "_0", "-": "_1"}
_choose_map = lambda r: _fixed_map[r] if r in _fixed_map else "_%s" % str(r)

# FIXME: MSA: Using a dict to collect functions in a namespace is weird
#             and makes the code harder to follow, change to a class
#             with member functions instead!
# FIXME: KBO: format is a builtin_function, i.e., we should use a different name.
# Formatting rules
format = {}

# Program flow
format.update({
    "return": lambda v: "return %s;" % str(v),
    "grouping": lambda v: "(%s)" % v,
    "block": lambda v: "{%s}" % v,
    "block begin": "{",
    "block end": "}",
    "list": lambda v: format["block"](format["list separator"].join([str(l) for l in v])),
    "switch": lambda v, cases, default=None, numbers=None: _generate_switch(v, cases, default, numbers),
    "exception": lambda v: "throw std::runtime_error(\"%s\");" % v,
    "warning": lambda v: 'std::cerr << "*** FFC warning: " << "%s" << std::endl;' % v,
    "comment": lambda v: "// %s" % v,
    "if": lambda c, v: "if (%s)\n{\n%s\n}\n" % (c, v),
    "loop": lambda i, j, k: "for (unsigned int %s = %s; %s < %s; %s++)" % (i, j, i, k, i),
    "generate loop": lambda v, w, _indent=0: _generate_loop(v, w, _indent),
    "is equal": " == ",
    "not equal": " != ",
    "less than": " < ",
    "greater than": " > ",
    "less equal": " <= ",
    "greater equal": " >= ",
    "and": " && ",
    "or": " || ",
    "not": lambda v: "!(%s)" % v,
    "do nothing": "// Do nothing"
})

# Declarations
format.update({
    "declaration": lambda t, n, v=None: _declaration(t, n, v),
    "float declaration": "double",
    "int declaration": "int",
    "uint declaration": "unsigned int",
    "static const uint declaration": "static const unsigned int",
    "static const float declaration": "static const double",
    "vector table declaration": "std::vector< std::vector<double> >",
    "double array declaration": "double*",
    "const double array declaration": "const double*",
    "const float declaration": lambda v, w: "const double %s = %s;" % (v, w),
    "const uint declaration": lambda v, w: "const unsigned int %s = %s;" % (v, w),
    "dynamic array": lambda t, n, s: "%s *%s = new %s[%s];" % (t, n, t, s),
    "static array": lambda t, n, s: "static %s %s[%d];" % (t, n, s),
    "fixed array": lambda t, n, s: "%s %s[%d];" % (t, n, s),
    "delete dynamic array": lambda n, s=None: _delete_array(n, s),
    "create foo": lambda v: "new %s()" % v,
    "create factory": lambda v: "create_%s()" % v
})

# Mathematical operators
format.update({
    "add": lambda v: " + ".join(v),
    "iadd": lambda v, w: "%s += %s;" % (str(v), str(w)),
    "sub": lambda v: " - ".join(v),
    "neg": lambda v: "-%s" % v,
    "mul": lambda v: "*".join(v),
    "imul": lambda v, w: "%s *= %s;" % (str(v), str(w)),
    "div": lambda v, w: "%s/%s" % (str(v), str(w)),
    "inverse": lambda v: "(1.0/%s)" % v,
    "std power": lambda base, exp: "std::pow(%s, %s)" % (base, exp),
    "exp": lambda v: "std::exp(%s)" % str(v),
    "ln": lambda v: "std::log(%s)" % str(v),
    "cos": lambda v: "std::cos(%s)" % str(v),
    "sin": lambda v: "std::sin(%s)" % str(v),
    "tan": lambda v: "std::tan(%s)" % str(v),
    "cosh": lambda v: "std::cosh(%s)" % str(v),
    "sinh": lambda v: "std::sinh(%s)" % str(v),
    "tanh": lambda v: "std::tanh(%s)" % str(v),
    "acos": lambda v: "std::acos(%s)" % str(v),
    "asin": lambda v: "std::asin(%s)" % str(v),
    "atan": lambda v: "std::atan(%s)" % str(v),
    "atan_2": lambda v1, v2: "std::atan2(%s,%s)" % (str(v1), str(v2)),
    "erf": lambda v: "erf(%s)" % str(v),
    "bessel_i": lambda v, n: "boost::math::cyl_bessel_i(%s, %s)" % (str(n), str(v)),
    "bessel_j": lambda v, n: "boost::math::cyl_bessel_j(%s, %s)" % (str(n), str(v)),
    "bessel_k": lambda v, n: "boost::math::cyl_bessel_k(%s, %s)" % (str(n), str(v)),
    "bessel_y": lambda v, n: "boost::math::cyl_neumann(%s, %s)" % (str(n), str(v)),
    "absolute value": lambda v: "std::abs(%s)" % str(v),
    "min value": lambda l, r: "std::min(%s, %s)" % (str(l), str(r)),
    "max value": lambda l, r: "std::max(%s, %s)" % (str(l), str(r)),
    "sqrt": lambda v: "std::sqrt(%s)" % str(v),
    "addition": lambda v: _add(v),
    "multiply": lambda v: _multiply(v),
    "power": lambda base, exp: _power(base, exp),
    "inner product": lambda v, w: _inner_product(v, w),
    "assign": lambda v, w: "%s = %s;" % (v, str(w)),
    "component": lambda v, k: _component(v, k)
})

# Formatting used in tabulate_tensor
format.update({
    "geometry tensor": lambda j, a: "G%d_%s" % (j, "_".join(["%d" % i for i in a]))
})

# Geometry related variable names (from code snippets).
format.update({
    "entity index": "entity_indices",
    "num entities": "num_global_entities",
    "cell": lambda s: "ufc::shape::%s" % s,
    "J": lambda i, j, m, n: "J[%d]" % _flatten(i, j, m, n),
    "inv(J)": lambda i, j, m, n: "K[%d]" % _flatten(i, j, m, n),
    "det(J)": lambda r=None: "detJ%s" % _choose_map(r),
    "cell volume": lambda r=None: "volume%s" % _choose_map(r),
    "circumradius": lambda r=None: "circumradius%s" % _choose_map(r),
    "facet area": "facet_area",
    "min facet edge length": lambda r: "min_facet_edge_length",
    "max facet edge length": lambda r: "max_facet_edge_length",
    "scale factor": "det",
    "transform": lambda t, i, j, m, n, r: _transform(t, i, j, m, n, r),
    "normal component": lambda r, j: "n%s%s" % (_choose_map(r), j),
    "x coordinate": "X",
    "y coordinate": "Y",
    "z coordinate": "Z",
    "ip coordinates": lambda i, j: "X%d[%d]" % (i, j),
    "affine map table": lambda i, j: "FEA%d_f%d" % (i, j),
    "coordinate_dofs": lambda r=None: "coordinate_dofs%s" % _choose_map(r)
})

# UFC function arguments and class members (names)
format.update({
    "element tensor": lambda i: "A[%s]" % i,
    "element tensor term": lambda i, j: "A%d[%s]" % (j, i),
    "coefficient": lambda j, k: format["component"]("w", [j, k]),
    "argument basis num": "i",
    "argument derivative order": "n",
    "argument values": "values",
    "argument coordinates": "dof_coordinates",
    "facet": lambda r: "facet%s" % _choose_map(r),
    "vertex": "vertex",
    "argument axis": "i",
    "argument dimension": "d",
    "argument entity": "i",
    "member global dimension": "_global_dimension",
    "argument dofs": "dofs",
    "argument dof num": "i",
    "argument dof values": "dof_values",
    "argument vertex values": "vertex_values",
    "argument sub": "i",  # sub element
    "argument subdomain": "subdomain_id",  # sub domain
})

# Formatting used in evaluatedof.
format.update({
    "dof vals": "vals",
    "dof result": "result",
    "dof X": lambda i: "X_%d" % i,
    "dof D": lambda i: "D_%d" % i,
    "dof W": lambda i: "W_%d" % i,
    "dof copy": lambda i: "copy_%d" % i,
    "dof physical coordinates": "y"
})


# Formatting used in evaluate_basis, evaluate_basis_derivatives and quadrature
# code generators.
format.update({
    # evaluate_basis and evaluate_basis_derivatives
    "tmp value": lambda i: "tmp%d" % i,
    "tmp ref value": lambda i: "tmp_ref%d" % i,
    "local dof": "dof",
    "basisvalues": "basisvalues",
    "coefficients": lambda i: "coefficients%d" % (i),
    "num derivatives": lambda t_or_g: "num_derivatives" + t_or_g,
    "derivative combinations": lambda t_or_g: "combinations" + t_or_g,
    "transform matrix": "transform",
    "transform Jinv": "Jinv",
    "dmats": lambda i: "dmats%s" % (i),
    "dmats old": "dmats_old",
    "reference derivatives": "derivatives",
    "dof values": "dof_values",
    "dof map if":                 lambda i, j: "%d <= %s && %s <= %d"\
                                  % (i, format["argument basis num"], format["argument basis num"], j),
    "dereference pointer": lambda n: "*%s" % n,
    "reference variable": lambda n: "&%s" % n,
    "call basis": lambda i, s: "_evaluate_basis(%s, %s, x, coordinate_dofs, cell_orientation);" % (i, s),
    "call basis_all": "_evaluate_basis_all(values, x, coordinate_dofs, cell_orientation);",
    "call basis_derivatives": lambda i, s: "_evaluate_basis_derivatives(%s, n, %s, x, coordinate_dofs, cell_orientation);" % (i, s),
    "call basis_derivatives_all": lambda i, s: "_evaluate_basis_derivatives_all(n, %s, x, coordinate_dofs, cell_orientation);" % s,

    # quadrature code generators
    "integration points": "ip",
    "first free index": "j",
    "second free index": "k",
    "geometry constant": lambda i: "G[%d]" % i,
    "ip constant": lambda i: "I[%d]" % i,
    "basis constant": lambda i: "B[%d]" % i,
    "conditional": lambda i: "C[%d]" % i,
    "evaluate conditional": lambda i, j, k: "(%s) ? %s : %s" % (i, j, k),
              #    "geometry constant":   lambda i: "G%d" % i,
              #    "ip constant":         lambda i: "I%d" % i,
              #    "basis constant":      lambda i: "B%d" % i,
              "function value": lambda i: "F%d" % i,
              "nonzero columns": lambda i: "nzc%d" % i,
              "weight": lambda i: "W" if i is None else "W%d" % (i),
              "psi name": lambda c, et, e, co, d, a: _generate_psi_name(c, et, e, co, d, a),
              # both
              "free indices": ["r", "s", "t", "u"],
              "matrix index": lambda i, j, range_j: _matrix_index(i, str(j), str(range_j)),
              "quadrature point": lambda i, gdim: "quadrature_points + %s*%d" % (i, gdim),
              "facet_normal_custom": lambda gdim: _generate_facet_normal_custom(gdim),
              })

# Misc
format.update({
    "bool": lambda v: {True: "true", False: "false"}[v],
    "str": lambda v: "%s" % v,
    "int": lambda v: "%d" % v,
    "list separator": ", ",
    "block separator": ",\n",
    "new line": "\\\n",
    "tabulate tensor": lambda m: _tabulate_tensor(m),
})

# Code snippets
from ffc.codesnippets import *

format.update({
    "compute_jacobian": lambda tdim, gdim, r=None:
                                compute_jacobian[tdim][gdim] % {"restriction": _choose_map(r)},
    "compute_jacobian_inverse": lambda tdim, gdim, r=None:
                                compute_jacobian_inverse[tdim][gdim] % {"restriction": _choose_map(r)},
    "orientation": lambda tdim, gdim, r=None: orientation_snippet % {"restriction": _choose_map(r)} if tdim != gdim else "",
    "facet determinant": lambda tdim, gdim, r=None: facet_determinant[tdim][gdim] % {"restriction": _choose_map(r)},
    "fiat coordinate map": lambda cell, gdim: fiat_coordinate_map[cell][gdim],
    "generate normal": lambda tdim, gdim, i: _generate_normal(tdim, gdim, i),
    "generate cell volume": lambda tdim, gdim, i, r=None: _generate_cell_volume(tdim, gdim, i, r),
    "generate circumradius": lambda tdim, gdim, i, r=None: _generate_circumradius(tdim, gdim, i, r),
    "generate facet area": lambda tdim, gdim: facet_area[tdim][gdim],
    "generate min facet edge length": lambda tdim, gdim, r=None: min_facet_edge_length[tdim][gdim] % {"restriction": _choose_map(r)},
    "generate max facet edge length": lambda tdim, gdim, r=None: max_facet_edge_length[tdim][gdim] % {"restriction": _choose_map(r)},
    "generate ip coordinates": lambda g, t, num_ip, name, ip, r=None: (ip_coordinates[t][g][0], ip_coordinates[t][g][1] %
                                                                       {"restriction": _choose_map(r), "ip": ip, "name": name, "num_ip": num_ip}),
    "scale factor snippet": scale_factor,
    "map onto physical": map_onto_physical,
    "evaluate basis snippet": eval_basis,
    "combinations": combinations_snippet,
    "transform snippet": transform_snippet,
    "evaluate function": evaluate_f,
    "ufc comment": comment_ufc,
    "dolfin comment": comment_dolfin,
    "header_h": header_h,
    "header_c": header_c,
    "footer": footer,
    "eval_basis_decl": eval_basis_decl,
    "eval_basis_init": eval_basis_init,
    "eval_basis": eval_basis,
    "eval_basis_copy": eval_basis_copy,
    "eval_derivs_decl": eval_derivs_decl,
    "eval_derivs_init": eval_derivs_init,
    "eval_derivs": eval_derivs,
    "eval_derivs_copy": eval_derivs_copy,
    "extract_cell_coordinates": lambda offset, r: "const double* coordinate_dofs_%d = coordinate_dofs + %d;" % (r, offset)
})

# Helper functions for formatting


def _declaration(type, name, value=None):
    if value is None:
        return "%s %s;" % (type, name)
    return "%s %s = %s;" % (type, name, str(value))


def _component(var, k):
    if not isinstance(k, (list, tuple)):
        k = [k]
    return "%s" % var + "".join("[%s]" % str(i) for i in k)


def _delete_array(name, size=None):
    if size is None:
        return "delete [] %s;" % name
    f_r = format["free indices"][0]
    code = format["generate loop"](["delete [] %s;" % format["component"](name, f_r)], [(f_r, 0, size)])
    code.append("delete [] %s;" % name)
    return "\n".join(code)


def _multiply(factors):
    """
    Generate string multiplying a list of numbers or strings.  If a
    factor is zero, the whole product is zero. Any factors equal to
    one are ignored.
    """

    # FIXME: This could probably be way more robust and elegant.

    cpp_str = format["str"]
    non_zero_factors = []
    for f in factors:

        # Round-off if f is smaller than epsilon
        if isinstance(f, (int, float)):
            if abs(f) < format["epsilon"]:
                return cpp_str(0)
            if abs(f - 1.0) < format["epsilon"]:
                continue

        # Convert to string
        f = cpp_str(f)

        # Return zero if any factor is zero
        if f == "0":
            return cpp_str(0)

        # If f is 1, don't add it to list of factors
        if f == "1":
            continue

        # If sum-like, parentheseze factor
        if "+" in f or "-" in f:
            f = "(%s)" % f

        non_zero_factors += [f]

    if len(non_zero_factors) == 0:
        return cpp_str(1.0)

    return "*".join(non_zero_factors)


def _add(terms):
    "Generate string summing a list of strings."

    # FIXME: Subtract absolute value of negative numbers
    result = " + ".join([str(t) for t in terms if (str(t) != "0")])
    if result == "":
        return format["str"](0)
    return result


def _power(base, exponent):
    "Generate code for base^exponent."
    if exponent >= 0:
        return _multiply(exponent * (base,))
    else:
        return "1.0 / (%s)" % _power(base, -exponent)


def _inner_product(v, w):
    "Generate string for v[0]*w[0] + ... + v[n]*w[n]."

    # Check that v and w have same length
    assert(len(v) == len(w)), "Sizes differ in inner-product!"

    # Special case, zero terms
    if len(v) == 0:
        return format["float"](0)

    # Straightforward handling when we only have strings
    if isinstance(v[0], string_types):
        return _add([_multiply([v[i], w[i]]) for i in range(len(v))])

    # Fancy handling of negative numbers etc
    result = None
    eps = format["epsilon"]
    add = format["add"]
    sub = format["sub"]
    neg = format["neg"]
    mul = format["mul"]
    fl = format["float"]
    for (c, x) in zip(v, w):
        if result:
            if abs(c - 1.0) < eps:
                result = add([result, x])
            elif abs(c + 1.0) < eps:
                result = sub([result, x])
            elif c > eps:
                result = add([result, mul([fl(c), x])])
            elif c < -eps:
                result = sub([result, mul([fl(-c), x])])
        else:
            if abs(c - 1.0) < eps:
                result = x
            elif abs(c + 1.0) < eps:
                result = neg(x)
            elif c > eps:
                result = mul([fl(c), x])
            elif c < -eps:
                result = neg(mul([fl(-c), x]))

    return result


def _transform(type, i, j, m, n, r):
    map_name = {"J": "J", "JINV": "K"}[type] + _choose_map(r)
    return (map_name + "[%d]") % _flatten(i, j, m, n)

# FIXME: Input to _generate_switch should be a list of tuples (i, case)


def _generate_switch(variable, cases, default=None, numbers=None):
    "Generate switch statement from given variable and cases"

    # Special case: no cases and no default
    if len(cases) == 0 and default is None:
        return format["do nothing"]
    elif len(cases) == 0:
        return default

    # Special case: one case and no default
    if len(cases) == 1 and default is None:
        return cases[0]

    # Create numbers for switch
    if numbers is None:
        numbers = list(range(len(cases)))

    # Create switch
    code = "switch (%s)\n{\n" % variable
    for (i, case) in enumerate(cases):
        code += "case %d:\n  {\n  %s\n    break;\n  }\n" % (numbers[i], indent(case, 2))
    code += "}\n"

    # Default value
    if default:
        code += "\n" + default

    return code


def _tabulate_tensor(vals):
    "Tabulate a multidimensional tensor. (Replace tabulate_matrix and tabulate_vector)."

    # Prefetch formats to speed up code generation
    f_block = format["block"]
    f_list_sep = format["list separator"]
    f_block_sep = format["block separator"]
    # FIXME: KBO: Change this to "float" once issue in set_float_formatting is fixed.
    f_float = format["floating point"]
    f_epsilon = format["epsilon"]

    # Create numpy array and get shape.
    tensor = numpy.array(vals)
    shape = numpy.shape(tensor)
    if len(shape) == 1:
        # Create zeros if value is smaller than tolerance.
        values = []
        for v in tensor:
            if abs(v) < f_epsilon:
                values.append(f_float(0.0))
            else:
                values.append(f_float(v))
        # Format values.
        return f_block(f_list_sep.join(values))
    elif len(shape) > 1:
        return f_block(f_block_sep.join([_tabulate_tensor(tensor[i]) for i in range(shape[0])]))
    else:
        error("Not an N-dimensional array:\n%s" % tensor)


def _generate_loop(lines, loop_vars, _indent):
    "This function generates a loop over a vector or matrix."

    # Prefetch formats to speed up code generation.
    f_loop = format["loop"]
    f_begin = format["block begin"]
    f_end = format["block end"]
    f_comment = format["comment"]

    if not loop_vars:
        return lines

    code = []
    for ls in loop_vars:
        # Get index and lower and upper bounds.
        index, lower, upper = ls
        # Loop index.
        code.append(indent(f_loop(index, lower, upper), _indent))
        code.append(indent(f_begin, _indent))

        # Increase indentation.
        _indent += 2

        # If this is the last loop, write values.
        if index == loop_vars[-1][0]:
            for l in lines:
                code.append(indent(l, _indent))

    # Decrease indentation and write end blocks.
    indices = [var[0] for var in loop_vars]
    indices.reverse()
    for index in indices:
        _indent -= 2
        code.append(indent(f_end + " " + f_comment("end loop over '%s'" % index), _indent))

    return code


def _matrix_index(i, j, range_j):
    "Map the indices in a matrix to an index in an array i.e., m[i][j] -> a[i*range(j)+j]"
    if i == 0:
        access = j
    elif i == 1:
        access = format["add"]([range_j, j])
    else:
        irj = format["mul"]([format["str"](i), range_j])
        access = format["add"]([irj, j])
    return access


def _generate_psi_name(counter, entity_type, entity, component, derivatives, avg):
    """Generate a name for the psi table of the form:
    FE#_f#_v#_C#_D###_A#, where '#' will be an integer value.

    FE  - is a simple counter to distinguish the various bases, it will be
          assigned in an arbitrary fashion.

    f   - denotes facets if applicable, range(element.num_facets()).

    v   - denotes vertices if applicable, range(num_vertices).

    C   - is the component number if any (flattened in the case of tensor valued functions)

    D   - is the number of derivatives in each spatial direction if any.
          If the element is defined in 3D, then D012 means d^3(*)/dydz^2.

    A   - denotes averaged over cell (AC) or facet (AF)
    """

    name = "FE%d" % counter

    if entity_type == "facet":
        name += "_f%d" % entity
    elif entity_type == "vertex":
        name += "_v%d" % entity

    if component != () and component != []:
        name += "_C%d" % component

    if any(derivatives):
        name += "_D" + "".join(map(str, derivatives))

    if avg == "cell":
        name += "_AC"
    elif avg == "facet":
        name += "_AF"

    return name


def _generate_normal(tdim, gdim, integral_type, reference_normal=False):
    "Generate code for computing normal"

    # Choose snippets
    direction = normal_direction[tdim][gdim]

    assert (gdim in facet_normal[tdim]),\
        "Facet normal not yet implemented for this tdim/gdim combo"
    normal = facet_normal[tdim][gdim]

    # Choose restrictions
    if integral_type == "exterior_facet":
        code = direction % {"restriction": "", "facet": "facet"}
        code += normal % {"direction": "", "restriction": ""}
    elif integral_type == "interior_facet":
        code = direction % {"restriction": _choose_map("+"), "facet": "facet_0"}
        code += normal % {"direction": "", "restriction": _choose_map("+")}
        code += normal % {"direction": "!", "restriction": _choose_map("-")}
    else:
        error("Unsupported integral_type: %s" % str(integral_type))
    return code


def _generate_facet_normal_custom(gdim):
    "Generate code for setting facet normal in custom integrals"
    code = format["comment"]("Set facet normal components for current quadrature point\n")
    for i in range(gdim):
        code += "const double n_0%d =   facet_normals[%d*ip + %d];\n" % (i, gdim, i)
        code += "const double n_1%d = - facet_normals[%d*ip + %d];\n" % (i, gdim, i)
    return code


def _generate_cell_volume(tdim, gdim, integral_type, r=None):
    "Generate code for computing cell volume."

    # Choose snippets
    volume = cell_volume[tdim][gdim]

    # Choose restrictions
    if integral_type in ("cell", "exterior_facet"):
        code = volume % {"restriction": ""}
    elif integral_type == "interior_facet":
        code = volume % {"restriction": _choose_map("+")}
        code += volume % {"restriction": _choose_map("-")}
    elif integral_type in custom_integral_types:
        code = volume % {"restriction": _choose_map(r)}
    else:
        error("Unsupported integral_type: %s" % str(integral_type))
    return code


def _generate_circumradius(tdim, gdim, integral_type, r=None):
    "Generate code for computing a cell's circumradius."

    # Choose snippets
    radius = circumradius[tdim][gdim]

    # Choose restrictions
    if integral_type in ("cell", "exterior_facet", "vertex"):
        code = radius % {"restriction": ""}
    elif integral_type == "interior_facet":
        code = radius % {"restriction": _choose_map("+")}
        code += radius % {"restriction": _choose_map("-")}
    elif integral_type in custom_integral_types:
        code = radius % {"restriction": _choose_map(r)}
    else:
        error("Unsupported integral_type: %s" % str(integral_type))
    return code


def _flatten(i, j, m, n):
    return i * n + j

# Other functions


def indent(block, num_spaces):
    "Indent each row of the given string block with n spaces."
    indentation = " " * num_spaces
    return indentation + ("\n" + indentation).join(block.split("\n"))


def count_ops(code):
    "Count the number of operations in code (multiply-add pairs)."
    num_add = code.count(" + ") + code.count(" - ")
    num_multiply = code.count("*") + code.count("/")
    return (num_add + num_multiply) // 2


def set_float_formatting(precision):
    "Set floating point formatting based on precision."

    # Options for float formatting
    # f1     = "%%.%df" % precision
    # f2     = "%%.%de" % precision
    f1 = "%%.%dg" % precision
    f2 = "%%.%dg" % precision
    f_int = "%%.%df" % 1

    eps = eval("1e-%s" % precision)

    # Regular float formatting
    def floating_point_regular(v):
        if abs(v - round(v, 1)) < eps:
            return f_int % v
        elif abs(v) < 100.0:
            return f1 % v
        else:
            return f2 % v

    # Special float formatting on Windows (remove extra leading zero)
    def floating_point_windows(v):
        return floating_point_regular(v).replace("e-0", "e-").replace("e+0", "e+")

    # Set float formatting
    if platform.system() == "Windows":
        format["float"] = floating_point_windows
    else:
        format["float"] = floating_point_regular

    # FIXME: KBO: Remove once we agree on the format of 'f1'
    format["floating point"] = format["float"]

    # Set machine precision
    format["epsilon"] = 10.0 * eval("1e-%s" % precision)

    # Hack to propagate precision to uflacs internals...
    import ffc.uflacs.language.format_value
    ffc.uflacs.language.format_value.set_float_precision(precision)


def set_exception_handling(convert_exceptions_to_warnings):
    "Set handling of exceptions."
    if convert_exceptions_to_warnings:
        format["exception"] = format["warning"]


# Declarations to examine
types = [["double"],
         ["const", "double"],
         ["const", "double", "*", "const", "*"],
         ["int"],
         ["const", "int"],
         ["unsigned", "int"],
         ["bool"],
         ["const", "bool"],
         ["static", "unsigned", "int"],
         ["const", "unsigned", "int"]]

# Special characters and delimiters
special_characters = ["+", "-", "*", "/", "=", ".", " ", ";", "(", ")", "\\", "{", "}", "[", "]", "!"]


def remove_unused(code, used_set=set()):
    """
    Remove unused variables from a given C++ code. This is useful when
    generating code that will be compiled with gcc and parameters -Wall
    -Werror, in which case gcc returns an error when seeing a variable
    declaration for a variable that is never used.

    Optionally, a set may be specified to indicate a set of variables
    names that are known to be used a priori.
    """

    # Dictionary of (declaration_line, used_lines) for variables
    variables = {}

    # List of variable names (so we can search them in order)
    variable_names = []

    lines = code.split("\n")
    for (line_number, line) in enumerate(lines):
        # Exclude commented lines.
        if line[:2] == "//" or line[:3] == "///":
            continue

        # Split words
        words = [word for word in line.split(" ") if word != ""]

        # Remember line where variable is declared
        for type in [type for type in types if " ".join(type) in " ".join(words)]:  # Fewer matches than line below.
        # for type in [type for type in types if len(words) > len(type)]:
            variable_type = words[0:len(type)]
            variable_name = words[len(type)]

            # Skip special characters
            if variable_name in special_characters:
                continue

            # Test if any of the special characters are present in the variable name
            # If this is the case, then remove these by assuming that the 'real' name
            # is the first entry in the return list. This is implemented to prevent
            # removal of e.g. 'double array[6]' if it is later used in a loop as 'array[i]'
            if variable_type == type:

                # Create correct variable name (e.g. y instead of
                # y[2]) for variables with separators
                seps_present = [sep for sep in special_characters if sep in variable_name]
                if seps_present:
                    variable_name = sorted([variable_name.split(sep)[0] for sep in seps_present])
                    variable_name = variable_name[0]

                variables[variable_name] = (line_number, [])
                if variable_name not in variable_names:
                    variable_names += [variable_name]

        # Mark line for used variables
        for variable_name in variables:
            (declaration_line, used_lines) = variables[variable_name]
            if _variable_in_line(variable_name, line) and line_number > declaration_line:
                variables[variable_name] = (declaration_line, used_lines + [line_number])

    # Reverse the order of the variable names to catch variables used
    # only by variables that are removed
    variable_names.reverse()

    # Remove declarations that are not used
    removed_lines = []
    for variable_name in variable_names:
        (declaration_line, used_lines) = variables[variable_name]
        for line in removed_lines:
            if line in used_lines:
                used_lines.remove(line)
        if not used_lines and variable_name not in used_set:
            debug("Removing unused variable: %s" % variable_name)
            lines[declaration_line] = None  # KBO: Need to completely remove line for evaluate_basis* to work
            # lines[declaration_line] = "// " + lines[declaration_line]
            removed_lines += [declaration_line]
    return "\n".join([line for line in lines if line is not None])


def _variable_in_line(variable_name, line):
    "Check if variable name is used in line"
    if variable_name not in line:
        return False
    for character in special_characters:
        line = line.replace(character, "\\" + character)
    delimiter = "[" + ",".join(["\\" + c for c in special_characters]) + "]"
    return re.search(delimiter + variable_name + delimiter, line) is not None