1
"Code generation for finite element"
3
__author__ = "Anders Logg (logg@simula.no)"
4
__date__ = "2007-01-23 -- 2009-08-26"
5
__copyright__ = "Copyright (C) 2007-2009 Anders Logg"
6
__license__ = "GNU GPL version 3 or any later version"
8
# Modified by Kristian Oelgaard 2009
9
# Modified by Marie Rognes 2008
10
# Modified by Garth N. Wells 2009
13
from FIAT.shapes import LINE
16
from ffc.common.log import debug, error
19
from ffc.fem.finiteelement import *
20
from ffc.fem.vectorelement import *
21
from ffc.fem.dofmap import *
22
from ffc.fem.quadratureelement import *
24
# FFC code generation common modules
25
from evaluatebasis import evaluate_basis
26
from evaluatebasisderivatives import evaluate_basis_derivatives
27
from codeutils import inner_product
29
def generate_finite_elements(form_data, format):
30
"Generate code for finite elements, including recursively nested sub elements."
32
debug("Generating code for finite elements...")
35
# Iterate over form elements
36
for (i, element) in enumerate(form_data.ffc_elements):
38
# Extract sub elements
39
sub_elements = _extract_sub_elements(element, (i,))
41
# Generate code for each element
42
for (label, sub_element) in sub_elements:
43
code += [(label, generate_finite_element(sub_element, format))]
48
def generate_finite_element(element, format):
49
"""Generate dictionary of code for the given finite element
50
according to the given format."""
54
# Generate code for signature
55
code["signature"] = element.signature()
57
# Generate code for cell_shape
58
code["cell_shape"] = format["cell shape"](element.cell_shape())
60
# Generate code for space_dimension
61
code["space_dimension"] = "%d" % element.space_dimension()
63
# Generate code for value_rank
64
# FIXME: This is just a temporary hack to 'support' tensor elements
65
code["value_rank"] = "%d" % element.value_rank()
66
code["value_rank"] = "%d" % element._rank
68
# Generate code for value_dimension
69
code["value_dimension"] = ["%d" % element.value_dimension(i) for i in range(max(element.value_rank(), 1))]
71
# Disable code generation for unsupported functions of QuadratureElement,
72
# (or MixedElements including QuadratureElements)
73
if not True in [isinstance(e, QuadratureElement) for e in element.extract_elements()]:
74
# Generate code for evaluate_basis
75
code["evaluate_basis"] = evaluate_basis(element, format)
77
# Generate code for evaluate_basis_derivatives
78
code["evaluate_basis_derivatives"] = evaluate_basis_derivatives(element, format)
80
# Generate vectorised version of evaluate functions
81
code["evaluate_basis_all"] =\
82
format["exception"]("The vectorised version of evaluate_basis() is not yet implemented.")
83
code["evaluate_basis_derivatives_all"] =\
84
format["exception"]("The vectorised version of evaluate_basis_derivatives() is not yet implemented.")
86
# Generate code for interpolate_vertex_values
87
code["interpolate_vertex_values"] = __generate_interpolate_vertex_values(element, format)
89
code["evaluate_basis"] =\
90
format["exception"]("evaluate_basis() is not supported for QuadratureElement")
91
code["evaluate_basis_derivatives"] =\
92
format["exception"]("evaluate_basis_derivatives() is not supported for QuadratureElement")
94
# Generate vectorised version of evaluate functions
95
code["evaluate_basis_all"] =\
96
format["exception"]("evaluate_basis_all() is not supported for QuadratureElement.")
97
code["evaluate_basis_derivatives_all"] =\
98
format["exception"]("evaluate_basis_derivatives_all() is not supported for QuadratureElement.")
100
code["interpolate_vertex_values"] =\
101
format["exception"]("interpolate_vertex_values() is not supported for QuadratureElement")
103
# Generate code for evaluate_dof
104
code["evaluate_dof"] = __generate_evaluate_dof(element, format)
106
# Generate code for num_sub_elements
107
code["num_sub_elements"] = "%d" % element.num_sub_elements()
111
def _extract_sub_elements(element, parent):
112
"""Recursively extract sub elements as a list of tuples where
113
each tuple consists of a tuple labeling the sub element and
114
the sub element itself."""
116
if element.num_sub_elements() == 1:
117
return [(parent, element)]
119
for i in range(element.num_sub_elements()):
120
sub_elements += _extract_sub_elements(element.sub_element(i), parent + (i,))
121
return sub_elements + [(parent, element)]
123
# FIXME: This is C++ dependent, move relevant parts to ufc_format.py
124
def __generate_evaluate_dof(element, format):
125
"Generate code for evaluate_dof"
127
# Generate code as a list of lines
131
block = format["block"]
132
separator = format["separator"]
133
floating_point = format["floating point"]
134
comment = format["comment"]
136
# Generate dof map and get _copy_ of dof representations.
137
dof_map = DofMap(element)
138
dofs = [DofRepresentation(dof) for dof in dof_map.dual_basis()]
141
# For ease in the code generalization, pad the points, directions
142
# and weights with zeros according to the maximal number of
143
# points. (Hence, the copy of the dofs above.)
144
max_num_points = dof_map.get_max_num_of_points()
145
num_points_per_dof = [dof.pad_points_and_weights(max_num_points)
148
# Compute the value dimension of the functions
150
for i in range(element.value_rank()):
151
num_values *= element.value_dimension(i)
152
# Check that the value dimension is the same for all dofs and that
153
# it matches the dimension of the function
154
value_dim = pick_first([dof.value_dim() for dof in dofs])
155
if value_dim != num_values:
156
error("Directional component does not match vector dimension")
158
# Initialize the points, weights and directions for the dofs:
159
code += [comment("The reference points, direction and weights:")]
160
s = block(separator.join(
161
[block(separator.join([block(separator.join([floating_point(c)
163
for point in dof.points]))
165
code += ["%sX[%d][%d][%d] = %s;" % (format["table declaration"],
166
num_dofs, max_num_points,
167
element.cell_dimension(), s)]
168
s = block(separator.join(
169
[block(separator.join([floating_point(w) for w in dof.weights]))
171
code += ["%sW[%d][%d] = %s;" % (format["table declaration"],
172
num_dofs, max_num_points, s)]
174
s = block(separator.join(
175
[block(separator.join([block(separator.join([floating_point(dk)
177
for d in dof.directions]))
179
code += ["%sD[%d][%d][%d] = %s;" % (format["table declaration"],
180
num_dofs, max_num_points,
184
# Compute the declarations needed for function mapping and the
185
# code for mapping each set of function values:
186
(map_declarations, map_values_code) = \
187
__map_function_values(num_values, element, format)
188
code += [map_declarations]
190
# Loop over the number of points (if more than one) and evaluate
192
code += ["%sresult = 0.0;" % format["float declaration"]]
193
code += [comment("Iterate over the points:") ]
197
# If there is more than one point, we need to add a table of the
198
# number of points per dof, add the loop, add indentation and add
201
if max_num_points > 1:
202
num_points_per_dof_code = block(separator.join([str(n) for n in num_points_per_dof]))
203
code += ["%sns[%d] = %s;" % (format["static const uint declaration"],
204
num_dofs, num_points_per_dof_code)]
205
code += ["%s%s" % (format["loop"]("j", "0", "ns[i]"), " {")]
206
(tab, endloop, index) = (2, "\n} // End for", "j")
208
# Map the points from the reference onto the physical element
209
code += [indent(format["snippet map_onto_physical"](element.cell_dimension())
210
% {"j": index}, tab)]
212
# Evaluate the function at the physical points
213
code += [indent(comment("Evaluate function at physical points"), tab)]
214
code += [indent("double values[%d];" % num_values, tab)]
215
code += [indent("f.evaluate(values, y, c);\n", tab)]
217
# Map the function values according to the given mapping(s)
218
code += [indent(comment("Map function values using appropriate mapping"),
220
code += [indent(map_values_code, tab)]
223
# Note that we do not map the weights (yet).
224
code += [indent(comment("Note that we do not map the weights (yet)."),tab)]
227
# Take the directional components of the function values and
228
# multiply by the weights:
229
code += [indent(format["snippet calculate dof"] % {"dim": value_dim,
230
"index": index}, tab)]
234
# Return the calculated value
235
code += [format["return"]("result")]
238
# FIXME: This is C++ dependent, move relevant parts to ufc_format.py
239
def __map_function_values(num_values, element, format):
241
block = format["block"]
242
separator = format["separator"]
243
comment = format["comment"]
247
# Compute the mapping associated with each dof:
248
# Looping the extracted element should have the same effect as using the
249
# old element.space_mapping() function.
250
# mappings = [mapping_to_int[element.space_mapping(i)]
251
# for i in range(element.space_dimension())]
253
for e in element.extract_elements():
254
for d in range(e.space_dimension()):
255
mappings.append(mapping_to_int[e.mapping()])
257
whichmappings = set(mappings)
259
# If there is more than one mapping involved, we will need to
260
# keep track of them at runtime:
261
if len(whichmappings) > 1:
262
precode += ["%smappings[%d] = %s;" %
263
(format["static const uint declaration"], len(mappings),
264
block(separator.join(([str(m) for m in mappings]))))]
266
# Check whether we will need a piola
267
contrapiola_present = mapping_to_int[CONTRAVARIANT_PIOLA] in whichmappings
268
copiola_present = mapping_to_int[COVARIANT_PIOLA] in whichmappings
269
piola_present = contrapiola_present or copiola_present
271
# Add code for the jacobian if we need it for the
272
# mappings. Otherwise, just add code for the vertex coordinates
273
if contrapiola_present:
274
# If contravariant piola: Will need J, det J and J^{-1}
275
precode += [format["snippet jacobian"](element.cell_dimension())
276
% {"restriction":""}]
277
precode += ["\ndouble copyofvalues[%d];" % num_values]
278
elif copiola_present:
279
# If covariant piola: Will need J only
280
precode += [format["snippet only jacobian"](element.cell_dimension())
281
% {"restriction":""}]
282
precode += ["\ndouble copyofvalues[%d];" % num_values]
284
precode += [format["get cell vertices"]]
286
# We have to add offsets to the code if there are mixed
287
# piola-mapped elements with an offset. (Ex: DG0 + RT)
289
if element.num_sub_elements() > 1 and piola_present:
292
for i in range(element.num_sub_elements()):
293
subelement = element.sub_element(i)
294
value_offsets += [adjustment]*subelement.space_dimension()
295
adjustment += subelement.value_dimension(0)
297
# if mapping[i] != AFFINE == 0 and value_offset[i] !=
298
# 0 for all i, then need_offsets = True:
299
need_offsets = bool(max([mappings[i] and value_offsets[i]
300
for i in range(len(mappings))]))
302
precode += ["const int offsets[%d] = %s;" %
304
block(separator.join([str(o) for o in value_offsets])))]
305
offset = "offsets[i] + "
308
# Then it just remains to actually add the different mappings to the code:
309
n = element.cell_dimension()
310
mappings_code = {mapping_to_int[AFFINE]: __affine_map(),
311
mapping_to_int[CONTRAVARIANT_PIOLA]:
312
__contravariant_piola(n, offset),
313
mapping_to_int[COVARIANT_PIOLA]: __covariant_piola(n, offset)}
315
ifs = ["mappings[i] == %d" % mapping for mapping in whichmappings]
316
cases = [mappings_code[mapping] for mapping in whichmappings]
317
code += [__generate_if_block(ifs, cases,
318
comment("Other mappings not applicable."))]
319
return ("\n".join(precode), "\n".join(code))
321
# FIXME: This is C++ dependent, move to ufc_format.py
323
return "// Affine map: Do nothing"
325
# FIXME: This is C++ dependent, move to ufc_format.py
326
def __contravariant_piola(dim, offset=""):
328
code += ["// Copy old values:"]
330
code += ["copyofvalues[%s%d] = values[%s%d];" % (offset, i, offset, i)]
331
code += ["// Do the inverse of div piola "]
333
terms = ["Jinv_%d%d*copyofvalues[%s%d]" % (i, j, offset, j)
335
code += ["values[%s%d] = detJ*(%s);" % (offset, i, "+".join(terms))]
336
return "\n".join(code)
338
# FIXME: This is C++ dependent, move to ufc_format.py
339
def __covariant_piola(dim, offset=""):
341
code += ["// Copy old values:"]
343
code += ["copyofvalues[%s%d] = values[%s%d];" % (offset, i, offset, i)]
344
code += ["// Do the inverse of curl piola "]
346
terms = ["J_%d%d*copyofvalues[%s%d]" % (j, i, offset, j)
348
code += ["values[%s%d] = %s;" % (offset, i, "+".join(terms))]
349
return "\n".join(code)
351
# FIXME: This is C++ dependent, move to ufc_format.py
352
def __generate_if_block(ifs, cases, default = ""):
353
"Generate if block from given ifs and cases"
354
if len(ifs) != len(cases):
355
error("Mismatch of dimensions ifs and cases")
357
# Special case: no cases
361
# Special case: one case
366
code = "if (%s) { \n" % ifs[0]
367
code += "%s\n" % indent(cases[0], 2)
368
for i in range(len(ifs)-1):
369
code += "} else if (%s) {\n %s \n" % (ifs[i+1], indent(cases[i+1], 2))
370
code += "} else { \n %s \n}" % indent(default,2)
373
def __generate_interpolate_vertex_values(element, format):
374
"Generate code for interpolate_vertex_values"
376
# Check that we have a scalar- or vector-valued element
377
if element.value_rank() > 1:
378
return format["exception"]("interpolate_vertex_values not implemented for this type of element")
380
# Generate code as a list of declarations
383
# Set vertices (note that we need to use the FIAT reference cells)
384
if element.cell_shape() == LINE:
385
vertices = [(0,), (1,)]
386
elif element.cell_shape() == TRIANGLE:
387
vertices = [(0, 0), (1, 0), (0, 1)]
388
elif element.cell_shape() == TETRAHEDRON:
389
vertices = [(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)]
391
# Extract nested sub elements
392
sub_elements = element.extract_elements()
394
# Iterate over sub elements
395
offset_dof_values = 0
396
offset_vertex_values = 0
397
need_jacobian = False
401
for sub_element in sub_elements:
402
size += sub_element.value_dimension(0)
404
for sub_element in sub_elements:
406
# Tabulate basis functions at vertices
407
table = sub_element.tabulate(0, vertices)
409
# Check which transform we should use to map the basis functions
410
mapping = sub_element.mapping()
412
# Generate different code depending on mapping
413
if mapping == AFFINE:
415
code += [format["comment"]("Evaluate at vertices and use affine mapping")]
417
# Handle scalars and vectors
418
if sub_element.value_rank() == 0:
419
for v in range(len(vertices)):
420
coefficients = table[0][sub_element.cell_dimension()*(0,)][:, v]
421
dof_values = [format["dof values"](offset_dof_values + n) for n in range(len(coefficients))]
422
name = format["vertex values"](size*v + offset_vertex_values)
423
value = inner_product(coefficients, dof_values, format)
424
code += [(name, value)]
426
for dim in range(sub_element.value_dimension(0)):
427
for v in range(len(vertices)):
428
coefficients = table[dim][0][sub_element.cell_dimension()*(0,)][:, v]
429
dof_values = [format["dof values"](offset_dof_values + n) for n in range(len(coefficients))]
430
name = format["vertex values"](size*v + offset_vertex_values + dim)
431
value = inner_product(coefficients, dof_values, format)
432
code += [(name, value)]
434
elif (mapping == CONTRAVARIANT_PIOLA or mapping == COVARIANT_PIOLA):
436
code += [format["comment"]("Evaluate at vertices and use Piola mapping")]
438
# Remember to add code later for Jacobian
441
# Check that dimension matches for Piola transform
442
if not sub_element.value_dimension(0) == sub_element.cell_dimension():
443
error("Vector dimension of basis function does not match for Piola transform.")
445
# Get entities for the dofs
446
dof_entities = DofMap(sub_element).dof_entities()
448
for dim in range(sub_element.value_dimension(0)):
449
for v in range(len(vertices)):
451
for n in range(sub_element.space_dimension()):
452
# Get basis function values at vertices
453
coefficients = [table[j][0][sub_element.cell_dimension()*(0,)][n, v] for j in range(sub_element.value_dimension(0))]
455
if mapping == COVARIANT_PIOLA:
456
# Get row of inverse transpose Jacobian
457
jacobian_row = [format["transform"]("JINV", j, dim, None) for j in range(sub_element.cell_dimension())]
459
# mapping == CONTRAVARIANT_PIOLA:
460
# Get row of Jacobian
461
jacobian_row = [format["transform"]("J", j, dim, None) for j in range(sub_element.cell_dimension())]
463
# Multiply vector-valued basis function with Jacobian
464
basis_function = inner_product(coefficients, jacobian_row, format)
465
# Add paranthesis if necessary
466
if "+" in basis_function or "-" in basis_function: # Cheating, should use dictionary
467
basis_function = format["grouping"](basis_function)
468
# Multiply with dof value
469
factors = [format["dof values"](offset_dof_values + n), basis_function]
471
if not basis_function == format["floating point"](0):
472
terms += [format["multiply"](factors)]
474
sum = format["grouping"](format["add"](terms))
476
sum = format["add"](terms)
477
name = format["vertex values"](size*v + offset_vertex_values + dim)
478
if mapping == CONTRAVARIANT_PIOLA:
479
value = format["multiply"]([format["inverse"](format["determinant"](None)), sum])
481
value = format["multiply"]([sum])
482
code += [(name, value)]
485
error("Unknown mapping: " + str(mapping))
487
offset_dof_values += sub_element.space_dimension()
488
offset_vertex_values += sub_element.value_dimension(0)
490
# Insert code for computing quantities needed for Piola mapping
492
code.insert(0, format["snippet jacobian"](element.cell_dimension()) % {"restriction": ""})
added
removed
ffc/compiler/finiteelement.py
