/usr/lib/python2.7/dist-packages/ffc/interpolatevertexvalues.py is in python-ffc 1.4.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 | "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
#
# Last changed: 2010-02-09
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, str):
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 three 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)
"""
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)]
else:
raise Exception, "No such mapping: %s accepted" % mapping
return change_of_variables
|