/usr/lib/python2.7/dist-packages/ffc/representation.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 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 | """
Compiler stage 2: Code representation
-------------------------------------
This module computes intermediate representations of forms,
elements and dofmaps. For each UFC function, we extract the
data needed for code generation at a later stage.
The representation should conform strictly to the naming and
order of functions in UFC. Thus, for code generation of the
function "foo", one should only need to use the data stored
in the intermediate representation under the key "foo".
"""
# Copyright (C) 2009-2013 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 Marie E. Rognes 2010
# Modified by Kristian B. Oelgaard 2010
# Modified by Martin Alnaes, 2013
#
# First added: 2009-12-16
# Last changed: 2014-04-24
# Python modules
from itertools import chain
# Import UFL
import ufl
from ufl.classes import Measure
from ufl.geometry import cellname2dim
# FFC modules
from ffc.utils import compute_permutations, product
from ffc.log import info, error, begin, end, debug_ir, ffc_assert, warning
from ffc.fiatinterface import create_element, cellname_to_num_entities, reference_cell
from ffc.mixedelement import MixedElement
from ffc.enrichedelement import EnrichedElement, SpaceOfReals
from ffc.quadratureelement import QuadratureElement
from ffc.cpp import set_float_formatting
def pick_representation(representation):
"Return one of the specialized code generation modules from a representation string."
if representation == "quadrature":
from ffc import quadrature
r = quadrature
elif representation == "tensor":
from ffc import tensor
r = tensor
elif representation == "uflacs":
from ffc import uflacsrepr
r = uflacsrepr
else:
error("Unknown representation: %s" % str(representation))
return r
not_implemented = None
def compute_ir(analysis, parameters):
"Compute intermediate representation."
begin("Compiler stage 2: Computing intermediate representation")
# Set code generation parameters
set_float_formatting(int(parameters["precision"]))
# Extract data from analysis
form_datas, elements, element_numbers = analysis
# Compute representation of elements
info("Computing representation of %d elements" % len(elements))
ir_elements = [_compute_element_ir(e, i, element_numbers) \
for (i, e) in enumerate(elements)]
# Compute representation of dofmaps
info("Computing representation of %d dofmaps" % len(elements))
ir_dofmaps = [_compute_dofmap_ir(e, i, element_numbers)
for (i, e) in enumerate(elements)]
# Compute and flatten representation of integrals
info("Computing representation of integrals")
irs = [_compute_integral_ir(fd, i, parameters) \
for (i, fd) in enumerate(form_datas)]
ir_integrals = [ir for ir in chain(*irs) if not ir is None]
# Compute representation of forms
info("Computing representation of forms")
ir_forms = [_compute_form_ir(fd, i, element_numbers) \
for (i, fd) in enumerate(form_datas)]
end()
return ir_elements, ir_dofmaps, ir_integrals, ir_forms
def _compute_element_ir(ufl_element, element_id, element_numbers):
"Compute intermediate representation of element."
# Create FIAT element
element = create_element(ufl_element)
domain, = ufl_element.domains() # Assuming single domain
cellname = domain.cell().cellname()
# Store id
ir = {"id": element_id}
# Compute data for each function
ir["signature"] = ufl_element.reconstruction_signature()
ir["cell_shape"] = cellname
ir["topological_dimension"] = domain.topological_dimension()
ir["geometric_dimension"] = domain.geometric_dimension()
ir["space_dimension"] = element.space_dimension()
ir["value_rank"] = len(ufl_element.value_shape())
ir["value_dimension"] = ufl_element.value_shape()
ir["evaluate_basis"] = _evaluate_basis(ufl_element, element)
ir["evaluate_dof"] = _evaluate_dof(ufl_element, element)
ir["interpolate_vertex_values"] = _interpolate_vertex_values(ufl_element, element)
ir["num_sub_elements"] = ufl_element.num_sub_elements()
ir["create_sub_element"] = _create_sub_foo(ufl_element, element_numbers)
#debug_ir(ir, "finite_element")
return ir
def _compute_dofmap_ir(ufl_element, element_id, element_numbers):
"Compute intermediate representation of dofmap."
# Create FIAT element
element = create_element(ufl_element)
domain, = ufl_element.domains() # Assuming single domain
cellname = domain.cell().cellname()
# Precompute repeatedly used items
num_dofs_per_entity = _num_dofs_per_entity(element)
facet_dofs = _tabulate_facet_dofs(element, cellname)
# Store id
ir = {"id": element_id}
# Compute data for each function
ir["signature"] = "FFC dofmap for " + ufl_element.reconstruction_signature()
ir["needs_mesh_entities"] = _needs_mesh_entities(element)
ir["topological_dimension"] = domain.topological_dimension()
ir["geometric_dimension"] = domain.geometric_dimension()
ir["global_dimension"] = _global_dimension(element)
ir["local_dimension"] = element.space_dimension()
ir["num_facet_dofs"] = len(facet_dofs[0])
ir["num_entity_dofs"] = num_dofs_per_entity
ir["tabulate_dofs"] = _tabulate_dofs(element, cellname)
ir["tabulate_facet_dofs"] = facet_dofs
ir["tabulate_entity_dofs"] = (element.entity_dofs(), num_dofs_per_entity)
ir["tabulate_coordinates"] = _tabulate_coordinates(ufl_element, element)
ir["num_sub_dofmaps"] = ufl_element.num_sub_elements()
ir["create_sub_dofmap"] = _create_sub_foo(ufl_element, element_numbers)
#debug_ir(ir, "dofmap")
return ir
def _global_dimension(element):
"Compute intermediate representation for global_dimension."
if not isinstance(element, MixedElement):
if isinstance(element, SpaceOfReals):
return ([], 1)
return (_num_dofs_per_entity(element), 0)
elements = []
reals = []
num_reals = 0
for (i, e) in enumerate(element.elements()):
if not isinstance(e, SpaceOfReals):
elements += [e]
else:
num_reals += 1
element = MixedElement(elements)
return (_num_dofs_per_entity(element), num_reals)
def _needs_mesh_entities(element):
"Compute intermediate representation for needs_mesh_entities."
# Note: The dof map for Real elements does not depend on the mesh
num_dofs_per_entity = _num_dofs_per_entity(element)
if isinstance(element, SpaceOfReals):
return [False for d in num_dofs_per_entity]
else:
return [d > 0 for d in num_dofs_per_entity]
def _compute_integral_ir(form_data, form_id, parameters):
"Compute intermediate represention for form integrals."
irs = []
# Iterate over integrals
for itg_data in form_data.integral_data:
# Select representation
# TODO: Is it possible to detach this metadata from IntegralData? It's a bit strange from the ufl side.
r = pick_representation(itg_data.metadata["representation"])
# Compute representation
ir = r.compute_integral_ir(itg_data,
form_data,
form_id,
parameters)
# Append representation
irs.append(ir)
return irs
def _compute_form_ir(form_data, form_id, element_numbers):
"Compute intermediate representation of form."
# Store id
ir = {"id": form_id}
# Compute common data
ir["classname"] = "FooForm"
ir["members"] = not_implemented
ir["constructor"] = not_implemented
ir["destructor"] = not_implemented
ir["signature"] = form_data.signature
ir["rank"] = form_data.rank
ir["num_coefficients"] = form_data.num_coefficients
ir["num_cell_domains"] = form_data.num_sub_domains.get("cell",0)
ir["num_exterior_facet_domains"] = form_data.num_sub_domains.get("exterior_facet",0)
ir["num_interior_facet_domains"] = form_data.num_sub_domains.get("interior_facet",0)
ir["num_point_domains"] = form_data.num_sub_domains.get("point",0)
ir["num_custom_domains"] = form_data.num_sub_domains.get("custom",0)
ir["has_cell_integrals"] = _has_foo_integrals("cell", form_data)
ir["has_exterior_facet_integrals"] = _has_foo_integrals("exterior_facet", form_data)
ir["has_interior_facet_integrals"] = _has_foo_integrals("interior_facet", form_data)
ir["has_point_integrals"] = _has_foo_integrals("point", form_data)
ir["has_custom_integrals"] = _has_foo_integrals("custom", form_data)
ir["create_finite_element"] = [element_numbers[e] for e in form_data.elements]
ir["create_dofmap"] = [element_numbers[e] for e in form_data.elements]
ir["create_cell_integral"] = _create_foo_integral("cell", form_data)
ir["create_exterior_facet_integral"] = _create_foo_integral("exterior_facet", form_data)
ir["create_interior_facet_integral"] = _create_foo_integral("interior_facet", form_data)
ir["create_point_integral"] = _create_foo_integral("point", form_data)
ir["create_custom_integral"] = _create_foo_integral("custom", form_data)
ir["create_default_cell_integral"] = _create_default_foo_integral("cell", form_data)
ir["create_default_exterior_facet_integral"] = _create_default_foo_integral("exterior_facet", form_data)
ir["create_default_interior_facet_integral"] = _create_default_foo_integral("interior_facet", form_data)
ir["create_default_point_integral"] = _create_default_foo_integral("point", form_data)
ir["create_default_custom_integral"] = _create_default_foo_integral("custom", form_data)
return ir
#--- Computation of intermediate representation for non-trivial functions ---
# FIXME: Move to FiniteElement/MixedElement
def _value_size(element):
"""Compute value size of element, aka the number of components.
The value size of a scalar field is 1, the value size of a vector
field (is the number of components), the value size of a higher
dimensional tensor field is the product of the value_shape of the
field. Recall that all mixed elements are flattened.
"""
shape = element.value_shape()
if shape == ():
return 1
return product(shape)
def _generate_reference_offsets(element, offset=0):
"""Generate offsets: i.e value offset for each basis function
relative to a reference element representation."""
offsets = []
if isinstance(element, MixedElement):
for e in element.elements():
offsets += _generate_reference_offsets(e, offset)
offset += _value_size(e)
elif isinstance(element, EnrichedElement):
for e in element.elements():
offsets += _generate_reference_offsets(e, offset)
else:
offsets = [offset]*element.space_dimension()
return offsets
def _generate_physical_offsets(ufl_element, offset=0):
"""Generate offsets: i.e value offset for each basis function
relative to a physical element representation."""
offsets = []
# Refer to reference if gdim == tdim. This is a hack to support
# more stuff (in particular restricted elements)
domain, = ufl_element.domains() # Assuming single domain
gdim = domain.geometric_dimension()
tdim = domain.topological_dimension()
if (gdim == tdim):
return _generate_reference_offsets(create_element(ufl_element))
if isinstance(ufl_element, ufl.MixedElement):
for e in ufl_element.sub_elements():
offsets += _generate_physical_offsets(e, offset)
offset += _value_size(e)
elif isinstance(ufl_element, ufl.EnrichedElement):
for e in ufl_element._elements:
offsets += _generate_physical_offsets(e, offset)
elif isinstance(ufl_element, ufl.FiniteElement):
element = create_element(ufl_element)
offsets = [offset]*element.space_dimension()
else:
raise NotImplementedError, \
"This element combination is not implemented"
return offsets
def _evaluate_dof(ufl_element, element):
"Compute intermediate representation of evaluate_dof."
# With regard to reference_value_size vs physical_value_size: Note
# that 'element' is the FFC/FIAT representation of the finite
# element, while 'ufl_element' is the UFL representation. In
# particular, UFL only knows about physical dimensions, so the
# value shape of the 'ufl_element' (which is used to compute the
# _value_size) will be correspond to the value size in physical
# space. FIAT however only knows about the reference element, and
# so the FIAT value shape of the 'element' will be the reference
# value size. This of course only matters for elements that have
# different physical and reference value shapes and sizes.
domain, = ufl_element.domains() # Assuming single domain
return {"mappings": element.mapping(),
"reference_value_size": _value_size(element),
"physical_value_size": _value_size(ufl_element),
"geometric_dimension": domain.geometric_dimension(),
"topological_dimension": domain.topological_dimension(),
"dofs": [L.pt_dict for L in element.dual_basis()],
"physical_offsets": _generate_physical_offsets(ufl_element)}
def _extract_elements(element):
new_elements = []
if isinstance(element, (MixedElement, EnrichedElement)):
for e in element.elements():
new_elements += _extract_elements(e)
else:
new_elements.append(element)
return new_elements
# def _num_components(element):
# """Compute the number of components of element, like _value_size, but
# does not support tensor elements."""
# shape = element.value_shape()
# if shape == ():
# return 1
# elif len(shape) == 1:
# return shape[0]
# else:
# error("Tensor valued elements are not supported yet: %d " % shape)
def _evaluate_basis(ufl_element, element):
"Compute intermediate representation for evaluate_basis."
domain, = ufl_element.domains() # Assuming single domain
cellname = domain.cell().cellname()
# Handle Mixed and EnrichedElements by extracting 'sub' elements.
elements = _extract_elements(element)
offsets = _generate_reference_offsets(element) # Must check?
mappings = element.mapping()
# This function is evidently not implemented for TensorElements
for e in elements:
if len(e.value_shape()) > 1:
return "Function not supported/implemented for TensorElements."
# Handle QuadratureElement, not supported because the basis is only defined
# at the dof coordinates where the value is 1, so not very interesting.
for e in elements:
if isinstance(e, QuadratureElement):
return "Function not supported/implemented for QuadratureElement."
# Initialise data with 'global' values.
data = {"reference_value_size": _value_size(element),
"physical_value_size": _value_size(ufl_element),
"cellname" : cellname,
"topological_dimension" : domain.topological_dimension(),
"geometric_dimension" : domain.geometric_dimension(),
"space_dimension" : element.space_dimension(),
"needs_oriented": needs_oriented_jacobian(element),
"max_degree": max([e.degree() for e in elements])
}
# Loop element and space dimensions to generate dof data.
dof = 0
dof_data = []
for e in elements:
for i in range(e.space_dimension()):
num_components = _value_size(e)
coefficients = []
coeffs = e.get_coeffs()
# Handle coefficients for vector valued basis elements
# [Raviart-Thomas, Brezzi-Douglas-Marini (BDM)].
if num_components > 1:
for c in range(num_components):
coefficients.append(coeffs[i][c])
else:
coefficients.append(coeffs[i])
dof_data.append(
{
"embedded_degree" : e.degree(),
"coeffs" : coefficients,
"num_components" : num_components,
"dmats" : e.dmats(),
"mapping" : mappings[dof],
"offset" : offsets[dof],
"num_expansion_members": e.get_num_members(e.degree())
})
dof += 1
data["dof_data"] = dof_data
return data
def _tabulate_coordinates(ufl_element, element):
"Compute intermediate representation of tabulate_coordinates."
if uses_integral_moments(element):
return {}
domain, = ufl_element.domains() # Assuming single domain
data = {}
data["tdim"] = domain.topological_dimension()
data["gdim"] = domain.geometric_dimension()
data["points"] = [L.pt_dict.keys()[0] for L in element.dual_basis()]
return data
def _tabulate_dofs(element, cellname):
"Compute intermediate representation of tabulate_dofs."
if isinstance(element, SpaceOfReals):
return None
# Extract number of entities for each dimension for this cell
num_entities = cellname_to_num_entities[cellname]
# Extract number of dofs per entity for each element
elements = all_elements(element)
num_dofs_per_element = [_num_dofs_per_entity(e) for e in elements]
# Extract local dof numbers per entity for each element
all_entity_dofs = [e.entity_dofs() for e in elements]
dofs_per_element = [[[list(dofs[dim][entity])
for entity in sorted(dofs[dim].keys())]
for dim in sorted(dofs.keys())]
for dofs in all_entity_dofs]
# Check whether we need offset
multiple_entities = any([sum(n > 0 for n in num_dofs) - 1
for num_dofs in num_dofs_per_element])
need_offset = len(elements) > 1 or multiple_entities
num_dofs_per_element = [e.space_dimension() for e in elements]
# Handle global "elements"
fakes = [isinstance(e, SpaceOfReals) for e in elements]
return (dofs_per_element, num_dofs_per_element, num_entities, need_offset, fakes)
def _tabulate_facet_dofs(element, cellname):
"Compute intermediate representation of tabulate_facet_dofs."
# Compute incidences
incidence = __compute_incidence(cellname2dim[cellname])
# Get topological dimension
D = max([pair[0][0] for pair in incidence])
# Get the number of facets
num_facets = cellname_to_num_entities[cellname][-2]
# Find out which entities are incident to each facet
incident = num_facets*[None]
for facet in range(num_facets):
incident[facet] = [pair[1] for pair in incidence if incidence[pair] == True and pair[0] == (D - 1, facet)]
# Make list of dofs
facet_dofs = []
entity_dofs = element.entity_dofs()
for facet in range(num_facets):
facet_dofs += [[]]
for dim in entity_dofs:
for entity in entity_dofs[dim]:
if (dim, entity) in incident[facet]:
facet_dofs[facet] += entity_dofs[dim][entity]
facet_dofs[facet].sort()
return facet_dofs
def _interpolate_vertex_values(ufl_element, element):
"Compute intermediate representation of interpolate_vertex_values."
# Check for QuadratureElement
for e in all_elements(element):
if isinstance(e, QuadratureElement):
return "Function is not supported/implemented for QuadratureElement."
domain, = ufl_element.domains() # Assuming single domain
cellname = domain.cell().cellname()
ir = {}
ir["geometric_dimension"] = domain.geometric_dimension()
ir["topological_dimension"] = domain.topological_dimension()
# Check whether computing the Jacobian is necessary
mappings = element.mapping()
ir["needs_jacobian"] = any("piola" in m for m in mappings)
ir["needs_oriented"] = needs_oriented_jacobian(element)
# See note in _evaluate_dofs
ir["reference_value_size"] = _value_size(element)
ir["physical_value_size"] = _value_size(ufl_element)
# Get vertices of reference cell
fiat_cell = reference_cell(cellname)
vertices = fiat_cell.get_vertices()
# Compute data for each constituent element
extract = lambda values: values[values.keys()[0]].transpose()
all_fiat_elm = all_elements(element)
ir["element_data"] = [{
# See note in _evaluate_dofs
"reference_value_size": _value_size(e),
"physical_value_size": _value_size(e), # FIXME: Get from corresponding ufl element
"basis_values": extract(e.tabulate(0, vertices)),
"mapping": e.mapping()[0],
"space_dim": e.space_dimension()}
for e in all_fiat_elm]
# FIXME: Temporary hack!
if len(ir["element_data"]) == 1:
ir["element_data"][0]["physical_value_size"] = ir["physical_value_size"]
# Consistency check, related to note in _evaluate_dofs
# This will fail for e.g. (RT1 x DG0) on a manifold
if sum(data["physical_value_size"] for data in ir["element_data"]) != ir["physical_value_size"]:
ir = "Failed to set physical value size correctly for subelements."
elif sum(data["reference_value_size"] for data in ir["element_data"]) != ir["reference_value_size"]:
ir = "Failed to set reference value size correctly for subelements."
return ir
def _create_sub_foo(ufl_element, element_numbers):
"Compute intermediate representation of create_sub_element/dofmap."
return [element_numbers[e] for e in ufl_element.sub_elements()]
def _create_foo_integral(integral_type, form_data):
"Compute intermediate representation of create_foo_integral."
return [itg_data.subdomain_id for itg_data in form_data.integral_data
if itg_data.integral_type == integral_type and isinstance(itg_data.subdomain_id, int)]
def _has_foo_integrals(integral_type, form_data):
"Compute intermediate representation of has_foo_integrals."
v = (form_data.num_sub_domains.get(integral_type,0) > 0
or _create_default_foo_integral(integral_type, form_data) is not None)
return bool(v)
def _create_default_foo_integral(integral_type, form_data):
"Compute intermediate representation of create_default_foo_integral."
itg_data = [itg_data for itg_data in form_data.integral_data
if (itg_data.subdomain_id == "otherwise" and
itg_data.integral_type == integral_type)]
ffc_assert(len(itg_data) in (0,1), "Expecting at most one default integral of each type.")
return "otherwise" if itg_data else None
#--- Utility functions ---
# FIXME: KBO: This could go somewhere else, like in UFL?
# MSA: There is probably something related in ufl somewhere,
# but I don't understand quite what this does.
# In particular it does not cover sub-sub-elements? Is that a bug?
# Also look at function naming, use single '_' for utility functions.
def all_elements(element):
if isinstance(element, MixedElement):
return element.elements()
return [element]
def _num_dofs_per_entity(element):
"""
Compute list of integers representing the number of dofs
associated with a single mesh entity.
Example: Lagrange of degree 3 on triangle: [1, 2, 1]
"""
entity_dofs = element.entity_dofs()
return [len(entity_dofs[e][0]) for e in range(len(entity_dofs.keys()))]
# These two are copied from old ffc
def __compute_incidence(D):
"Compute which entities are incident with which"
# Compute the incident vertices for each entity
sub_simplices = []
for dim in range(D + 1):
sub_simplices += [__compute_sub_simplices(D, dim)]
# Check which entities are incident, d0 --> d1 for d0 >= d1
incidence = {}
for d0 in range(0, D + 1):
for i0 in range(len(sub_simplices[d0])):
for d1 in range(d0 + 1):
for i1 in range(len(sub_simplices[d1])):
if min([v in sub_simplices[d0][i0] for v in sub_simplices[d1][i1]]) == True:
incidence[((d0, i0), (d1, i1))] = True
else:
incidence[((d0, i0), (d1, i1))] = False
return incidence
def __compute_sub_simplices(D, d):
"""Compute vertices for all sub simplices of dimension d (code
taken from Exterior)."""
# Number of vertices
num_vertices = D + 1
# Special cases: d = 0 and d = D
if d == 0:
return [[i] for i in range(num_vertices)]
elif d == D:
return [range(num_vertices)]
# Compute all permutations of num_vertices - (d + 1)
permutations = compute_permutations(num_vertices - d - 1, num_vertices)
# Iterate over sub simplices
sub_simplices = []
for i in range(len(permutations)):
# Pick tuple i among permutations (non-incident vertices)
remove = permutations[i]
# Remove vertices, keeping d + 1 vertices
vertices = [v for v in range(num_vertices) if not v in remove]
sub_simplices += [vertices]
return sub_simplices
def uses_integral_moments(element):
"True if element uses integral moments for its degrees of freedom."
integrals = set(["IntegralMoment", "FrobeniusIntegralMoment"])
tags = set([L.get_type_tag() for L in element.dual_basis()])
return len(integrals & tags) > 0
def needs_oriented_jacobian(element):
# Check whether this element needs an oriented jacobian (only
# contravariant piolas seem to need it)
return ("contravariant piola" in element.mapping())
|