This file is indexed.

/usr/lib/python2.7/dist-packages/ffc/interpolatevertexvalues.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
# -*- coding: utf-8 -*-
"Code generation for interpolate_vertex_values."

# Copyright (C) 2009 Marie E. Rognes
#
# 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 2010
# Modified by Lizao Li 2015, 2016
#
# Last changed: 2016-08-17

from six import string_types
from ffc.cpp import format, remove_unused

# Extract code manipulation formats
inner = format["inner product"]
component = format["component"]
assign = format["assign"]
multiply = format["multiply"]

# Extract formats for the Jacobians
J = format["J"]
Jinv = format["inv(J)"]
invdetJ = format["inverse"](format["det(J)"](None))

f_dof_values = format["argument dof values"]
f_vertex_values = format["argument vertex values"]


def interpolate_vertex_values(ir):
    "Generate code for interpolate_vertex_values."

    # Handle unsupported elements.
    if isinstance(ir, string_types):
        return format["exception"]("interpolate_vertex_values: %s" % ir)

    # Add code for Jacobian if necessary
    code = []
    gdim = ir["geometric_dimension"]
    tdim = ir["topological_dimension"]
    if ir["needs_jacobian"]:

        # Generate code for basic geometric quantities
        code.append(format["compute_jacobian"](tdim, gdim))
        code.append("")
        code.append(format["compute_jacobian_inverse"](tdim, gdim))
        if ir["needs_oriented"]:
            code.append("")
            code.append(format["orientation"](tdim, gdim))

    # Compute total value dimension for (mixed) element
    total_dim = ir["physical_value_size"]

    # Generate code for each element
    value_offset = 0
    space_offset = 0
    for data in ir["element_data"]:
        # Add vertex interpolation for this element
        code.append(format["comment"]("Evaluate function and change variables"))
        code.append(_interpolate_vertex_values_element(data, gdim, tdim,
                                                       total_dim,
                                                       value_offset,
                                                       space_offset))

        # Update offsets for value- and space dimension
        value_offset += data["physical_value_size"]
        space_offset += data["space_dim"]

    # Remove unused variables. (Not tracking set of used variables in
    # order to keep this code clean. Since generated code is of
    # limited size, this should be ok.)
    clean_code = remove_unused("\n".join(code))
    return clean_code


def _interpolate_vertex_values_element(data, gdim, tdim, total_value_size,
                                       value_offset=0, space_offset=0):

    # Extract vertex values for all basis functions
    vertex_values = data["basis_values"]
    value_size = data["physical_value_size"]
    space_dim = data["space_dim"]
    mapping = data["mapping"]

    # Map basis values according to element mapping. Assumes single
    # mapping for each (non-mixed) element
    change_of_variables = _change_variables(data["mapping"], gdim, tdim,
                                            space_dim)

    # Create code for each value dimension:
    code = []
    for k in range(value_size):

        # Create code for each vertex x_j
        for (j, values_at_vertex) in enumerate(vertex_values):

            if value_size == 1:
                values_at_vertex = [values_at_vertex]

            # Map basis functions using appropriate mapping
            components = change_of_variables(values_at_vertex, k)

            # Contract coefficients and basis functions
            dof_values = [component(f_dof_values, i + space_offset)
                          for i in range(space_dim)]

            value = inner(dof_values, components)

            # Assign value to correct vertex
            index = j * total_value_size + (k + value_offset)
            code.append(assign(component(f_vertex_values, index), value))

    return "\n".join(code)


def _change_variables(mapping, gdim, tdim, space_dim):
    """

    How to map a field G from the reference domain to a physical
    domain: For the converse approach -- see evaluatedof.py

    Let g be a field defined on the reference domain T_0 (of
    dimension tdim) with reference coordinates X. Let T be a a
    physical domain (of dimension gdim) with coordinates x. Assume
    that F: T_0 -> T such that

      x = F(X)

    Let J be the Jacobian of F, i.e J = dx/dX and let K denote the
    (pseudo)-inverse of the Jacobian K = J^{-1}. Note that J is gdim x
    tdim, and conversely K is tdim x gdim.  Then we (currently) have
    the following four types of mappings:

    'affine' mapping for G:

      g(x) = G(X)

    For vector fields G:

    'contravariant piola' mapping for f:

      g(x) = 1.0/det(J) J G(X)   i.e   g_i(x) = 1.0/det(J) J_ij G_j(X)

    'covariant piola' mapping for f:

      g(x) = K^T G(X)              i.e   g_i(x) = K^T_ij G_j(X) = K_ji G_j(X)

    'double covariant piola' mapping for f:

      g_il(x) = K_{ji} G_{jk} K_{kl}

    'double contravariant piola' mapping for g:

      g_il(x) = (det(J))^(-2) J_ij G_jk(X) J_lk

    """

    if mapping is "affine":
        change_of_variables = lambda G, i: G[i]
    elif mapping == "contravariant piola":
        change_of_variables = lambda G, i: [multiply([invdetJ, inner([J(i, j, gdim, tdim) for j in range(tdim)],
                                                                     [G[j][index] for j in range(tdim)])])
                                            for index in range(space_dim)]
    elif mapping == "covariant piola":
        change_of_variables = lambda G, i: [inner([Jinv(j, i, tdim, gdim) for j in range(tdim)],
                                                  [G[j][index] for j in range(tdim)])
                                            for index in range(space_dim)]
    elif mapping == "double covariant piola":
        change_of_variables = lambda G, i: [
            inner([inner([Jinv(j, i // tdim, tdim, gdim) for j in range(tdim)],
                         [G[j][k][index] for j in range(tdim)])
                   for k in range(tdim)],
                  [Jinv(k, i % tdim, tdim, gdim) for k in range(tdim)])
            for index in range(space_dim)]
    elif mapping == "double contravariant piola":
        change_of_variables = lambda G, i: [
            multiply([invdetJ, invdetJ, inner(
                [inner([J(i // tdim, j, tdim, gdim) for j in range(tdim)],
                         [G[j][k][index] for j in range(tdim)])
                   for k in range(tdim)],
                  [J(i % tdim, k, tdim, gdim) for k in range(tdim)])])
            for index in range(space_dim)]
    else:
        raise Exception("No such mapping: %s accepted" % mapping)
    return change_of_variables