Commit abb3a59472ad82278ef4428d0a42a8fba7987838 - fiat

+144

-0

.mailmap less more

	0	Anders Logg <logg@chalmers.se>
	1	Anders Logg <logg@chalmers.se> <ffc@bamse>
	2	Anders Logg <logg@chalmers.se> <logg@yavanna>
	3	Anders Logg <logg@chalmers.se> <logg@aule>
	4	Anders Logg <logg@chalmers.se> logg <logg@localhost.localdomain>
	5	Anders Logg <logg@chalmers.se> fenics <devnull@localhost>
	6	Anders Logg <logg@chalmers.se> hg <hg@numenor>
	7	Anders Logg <logg@chalmers.se> root <root@leja>
	8	Anders Logg <logg@chalmers.se> <logg@glaurung>
	9	Anders Logg <logg@chalmers.se> <logg@simula.no>
	10	Anders Logg <logg@chalmers.se> <anders.logg+bitbucket@gmail.com>
	11	Anders Logg <logg@chalmers.se> <logg@olorin>
	12	Anders Logg <logg@chalmers.se> <logg@tti-c.org>
	13	Anders Logg <logg@chalmers.se> anders <anders@localhost.localdomain>
	14	Anders Logg <logg@chalmers.se> logg <devnull@localhost>
	15	Anders Logg <logg@chalmers.se> logg <logg@bigblue.simula.no>
	16	Anders Logg <logg@chalmers.se> logg <logg@bunjil.simula.no>
	17	Anders Logg <logg@chalmers.se> logg <logg@olorin>
	18	Andy R. Terrel <andy.terrel@gmail.com>
	19	Andy R. Terrel <andy.terrel@gmail.com> <aterrel@uchicago.edu>
	20	Andy R. Terrel <andy.terrel@gmail.com> aterrel <devnull@localhost>
	21	Andy R. Terrel <andy.terrel@gmail.com> <aterrel@dhcp-11-20.cs.uchicago.edu>
	22	Benjamin Kehlet <benjamik@simula.no>
	23	Benjamin Kehlet <benjamik@simula.no> <benjamik@benjamik-hp.simula.no>
	24	Benjamin Kehlet <benjamik@simula.no> <benjamik@ifi.uio.no>
	25	Benjamin Kehlet <benjamik@simula.no> <benjamik@login.ifi.uio.no>
	26	Benjamin Kehlet <benjamik@simula.no> <bkehlet@bkehlet-laptop>
	27	Chris Richardson <chris@bpi.cam.ac.uk>
	28	Chris Richardson <chris@bpi.cam.ac.uk> chris <chris@Lenovo-Edge.(none)>
	29	Chris Richardson <chris@bpi.cam.ac.uk> Chris Richardson <bpi.cam.ac.uk>
	30	Chris Richardson <chris@bpi.cam.ac.uk> root <root@bpi.cam.ac.uk>
	31	Corrado Maurini <corrado.maurini@upmc.fr>
	32	Corrado Maurini <corrado.maurini@upmc.fr> <cmaurini@gmail.com>
	33	Dag Lindbo <dag@csc.kth.se>
	34	Dag Lindbo <dag@csc.kth.se> dag <dag@dag-laptop>
	35	Dag Lindbo <dag@csc.kth.se> dag <dag@na55.nada.kth.se>
	36	David Ham <david.ham@imperial.ac.uk>
	37	David Ham <david.ham@imperial.ac.uk> <David.Ham@imperial.ac.uk>
	38	David Ham <david.ham@imperial.ac.uk> david.ham@imperial.ac.uk <>
	39	Evan Lezar <evanlezar@gmail.com>
	40	Evan Lezar <evanlezar@gmail.com> <mail@evanlezar.com>
	41	Evan Lezar <evanlezar@gmail.com> elezar <elezar@labby>
	42	Evan Lezar <evanlezar@gmail.com> elezar <elezar@lefauve>
	43	Fredrik Valdmanis <fredrik@valdmanis.com>
	44	Fredrik Valdmanis <fredrik@valdmanis.com> Fredrik Valdmanis <fredva@ifi.uio.no>
	45	Garth N. Wells <gnw20@cam.ac.uk>
	46	Garth N. Wells <gnw20@cam.ac.uk> root <root@fornax.esc.cam.ac.uk>
	47	Garth N. Wells <gnw20@cam.ac.uk> garth <devnull@localhost>
	48	Garth N. Wells <gnw20@cam.ac.uk> <garth@debian.eng.cam.ac.uk>
	49	Garth N. Wells <gnw20@cam.ac.uk> <garth@fedora32-virtual>
	50	Garth N. Wells <gnw20@cam.ac.uk> <garth@gnw20pc>
	51	Garth N. Wells <gnw20@cam.ac.uk> <garth@home-laptop>
	52	Garth N. Wells <gnw20@cam.ac.uk> <garth@localhost.localdomain>
	53	Garth N. Wells <gnw20@cam.ac.uk> <garth@ubuntu32-virtual>
	54	Garth N. Wells <gnw20@cam.ac.uk> <garth@ubuntu.eng.cam.ac.uk>
	55	Garth N. Wells <gnw20@cam.ac.uk> <garth@garth-laptop>
	56	Garth N. Wells <gnw20@cam.ac.uk> <garth@gnw20pc>
	57	Garth N. Wells <gnw20@cam.ac.uk> <g.n.wells@tudelft.nl>
	58	Garth N. Wells <gnw20@cam.ac.uk> gnw20@cam.ac.uk <>
	59	gideonsimpson <gideonsimpson@dyn-128-59-151-14.dyn.columbia.edu>
	60	Gustav Magnus Vikström <gustavv@simula.no>
	61	Gustav Magnus Vikström <gustavv@simula.no> <gustavv@ifi.uio.no>
	62	Gustav Magnus Vikström <gustavv@simula.no> <gustavv@ifi.uio.no>
	63	Gustav Magnus Vikström <gustavv@simula.no> <gustavv@utlaan-laptop-1>
	64	Harish Narayanan <hnarayanan@gmail.com>
	65	Harish Narayanan <hnarayanan@gmail.com> Harish Narayanan <harish@simula.no>
	66	Johan Hoffman <jhoffman@csc.kth.se>
	67	Johan Hoffman <jhoffman@csc.kth.se> hoffman <devnull@localhost>
	68	Johan Hoffman <jhoffman@csc.kth.se> <jhoffman@na41.nada.kth.se>
	69	Johan Hoffman <jhoffman@csc.kth.se> <jhoffman@na42.nada.kth.se>
	70	Ilmar Wilbers <ilmarw@simula.no>
	71	Ilmar Wilbers <ilmarw@simula.no> <ilmarw@gogmagog.simula.no>
	72	Ilmar Wilbers <ilmarw@simula.no> <ilmarw@multiboot.local>
	73	Jack S. Hale <jack.hale@uni.lu>
	74	Jack S. Hale <jack.hale@uni.lu> <j.hale09@imperial.ac.uk>
	75	Johan Hake <hake.dev@gmail.com>
	76	Johan Hake <hake.dev@gmail.com> <johan.hake@gmail.com>
	77	Johan Hake <hake.dev@gmail.com> <hake@simula.no>
	78	Johan Jansson <jjan@csc.kth.se>
	79	Johan Jansson <jjan@csc.kth.se> johan <johan@localhost.localdomain>
	80	Johan Jansson <jjan@csc.kth.se> johan <johan@nova>
	81	Johan Jansson <jjan@csc.kth.se> johanjan <devnull@localhost>
	82	Johan Jansson <jjan@csc.kth.se> johanjan <johanjan@localhost.localdomain>
	83	Johan Jansson <jjan@csc.kth.se> Johan Jansson <johanjan@math.chalmers.se>
	84	Johannes Ring <johannr@simula.no>
	85	Johannes Ring <johannr@simula.no> <johannr@communalis.simula.no>
	86	Kent-Andre Mardal <kent-and@simula.no>
	87	Kent-Andre Mardal <kent-and@simula.no> <kent-and@localhost>
	88	Kristian B. Ølgaard <k.b.oelgaard@gmail.com>
	89	Kristian B. Ølgaard <k.b.oelgaard@gmail.com> <k.b.oelgaard@tudelft.nl>
	90	Kristian B. Ølgaard <k.b.oelgaard@gmail.com> <oelgaard@localhost.localdomain>
	91	Magnus Vikstrøm <gustavv@ifi.uio.no>
	92	Marco Morandini <marco.morandini@polimi.it>
	93	Marco Morandini <marco.morandini@polimi.it> <morandini@aero.polimi.it>
	94	Marie E. Rognes <meg@simula.no>
	95	Marie E. Rognes <meg@simula.no> <meg@math.uio.no>
	96	Marie E. Rognes <meg@simula.no> <meg@meg-laptop>
	97	Marie E. Rognes <meg@simula.no> Marie E. Rognes (meg@simula.no) <Marie E. Rognes (meg@simula.no)>
	98	Marie E. Rognes <meg@simula.no> meg@simula.no <>
	99	Martin Sandve Alnæs <martinal@simula.no>
	100	Martin Sandve Alnæs <martinal@simula.no> <martinal@localhost>
	101	Martin Sandve Alnæs <martinal@simula.no> <martinal@martinal-desktop>
	102	Michele Zaffalon <michele.zaffalon@gmail.com>
	103	Mikael Mortensen <mikaem@math.uio.no>
	104	Mikael Mortensen <mikaem@math.uio.no> <mikael.mortensen@gmail.com>
	105	Nate Sime <njcs4@cam.ac.uk>
	106	Nate Sime <njcs4@cam.ac.uk> <njcs4@galah.bpi.cam.ac.uk>
	107	Nuno Lopes <ndl@ptmat.fc.ul.pt>
	108	Nuno Lopes <ndl@ptmat.fc.ul.pt> N.Lopes <devnull@localhost>
	109	Patrick Farrell <patrick.farrell@maths.ox.ac.uk>
	110	Patrick Farrell <patrick.farrell@maths.ox.ac.uk> <patrick.farrell06@imperial.ac.uk>
	111	Patrick Farrell <patrick.farrell@maths.ox.ac.uk> <patrick.farrell@imperial.ac.uk>
	112	Quang Ha <qth20@cam.ac.uk>
	113	Kristoffer Selim <selim@simula.no>
	114	Kristoffer Selim <selim@simula.no> <selim@selim-laptop>
	115	Simon Funke <simon@simula.no>
	116	Simon Funke <simon@simula.no> <simon.funke@gmail.com>
	117	Simon Funke <simon@simula.no> <s.funke09@imperial.ac.uk>
	118	Solveig Bruvoll <solveio@ifi.uio.no>
	119	Solveig Masvie <smasvie@gmail.com>
	120	Steven Vandekerckhove <steven.vandekerckhove@kuleuven-kulak.be>
	121	Steven Vandekerckhove <steven.vandekerckhove@kuleuven-kulak.be> <Steven.Vandekerckhove@kuleuven-kulak.be>
	122	stockli <stockli@carleman.nada.kth.se>
	123	stockli <stockli@carleman.nada.kth.se> <stockli@localhost>
	124	Tianyi Li <tianyikillua@gmail.com>
	125	Steffen Müthing <steffen.muething@ipvs.uni-stuttgart.de>
	126	Steffen Müthing <steffen.muething@ipvs.uni-stuttgart.de> Steffen Müthing steffen.muething@ipvs.uni-stuttgart.de <>
	127	Miklós Homolya <m.homolya14@imperial.ac.uk>
	128	Åsmund Ødegård <aasmund@simula.no>
	129	Åsmund Ødegård <aasmund@simula.no> <aasmundo@manpower.local>
	130	Ola Skavhaug <skavhaug@simula.no>
	131	Andre Massing <massing@simula.no>
	132	Andrew McRae <a.mcrae12@imperial.ac.uk>
	133	Andrew McRae <a.t.t.mcrae@bath.ac.uk> <a.mcrae12@imperial.ac.uk>
	134	Lawrence Mitchell <wence@gmx.li> Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk>
	135	Robert Kirby <robert.c.kirby@gmail.com> Robert C. Kirby <kirby@uchicago.edu>
	136	Robert Kirby <robert.c.kirby@gmail.com> Robert C. Kirby <robert.c.kirby@gmail.com>
	137	Robert Kirby <robert.c.kirby@gmail.com> kirby <devnull@localhost>
	138	Robert Kirby <robert.c.kirby@gmail.com> kirby <kirby@LittleMy.local>
	139	Robert Kirby <robert.c.kirby@gmail.com> kirby <kirby@localhost>
	140	Robert Kirby <robert.c.kirby@gmail.com> kirby <kirby@ma072.dhcp.ttu.edu>
	141	Robert Kirby <robert.c.kirby@gmail.com> kirby <kirby@ma251.dhcp.ttu.edu>
	142	Robert Kirby <robert.c.kirby@gmail.com> kirby <kirby@zortzi.math.ttu.edu>
	143	Matthew Knepley <knepley@gmail.com> knepley <knepley@knepley-laptop>

+26

-0

.travis.yml less more

	0	sudo: false
	1	notifications:
	2	irc:
	3	channels: "chat.freenode.net#firedrake"
	4	skip_join: true
	5	on_success: change
	6	on_failure: always
	7	template: "%{repository}#%{build_number} (%{branch} - %{commit} : %{author}): %{message} \| %{build_url}"
	8	language: python
	9	python:
	10	- "2.7"
	11	- "3.4"
	12
	13	before_install:
	14	- pip install --upgrade pip
	15	- pip install flake8
	16	- pip install flake8-future-import
	17	- pip install pytest
	18	- pip install six
	19
	20	install:
	21	- pip install .
	22
	23	script:
	24	- flake8 .
	25	- DATA_REPO_GIT="" py.test test/

+3

-0

AUTHORS less more

66	66
67	67	Matthias Liertzer
68	68	email: matthias@liertzer.at
	69
	70	Ivan Yashchuk
	71	email: ivan.yashchuk@aalto.fi

+11

-0

ChangeLog.rst less more

0	0	Changelog
1	1	=========
	2
	3	2017.2.0 (2017-09-19)
	4	---------------------
	5
	6	- Add quadrilateral and hexahedron reference cells
	7	- Add quadrilateral and hexahedron elements (with a wrapping class for TensorProductElement)
	8
	9	2017.1.0.post1 (2017-09-12)
	10	---------------------------
	11
	12	- Change PyPI package name to fenics-fiat.
2	13
3	14	2017.1.0 (2017-05-09)
4	15	---------------------

+3

-1

FIAT/__init__.py less more

32	32	from FIAT.nodal_enriched import NodalEnrichedElement
33	33	from FIAT.discontinuous import DiscontinuousElement
34	34	from FIAT.hdiv_trace import HDivTrace
	35	from FIAT.mixed import MixedElement # noqa: F401
35	36	from FIAT.restricted import RestrictedElement # noqa: F401
	37	from FIAT.quadrature_element import QuadratureElement # noqa: F401
36	38
37	39	# Important functionality
38	40	from FIAT.quadrature import make_quadrature # noqa: F401

40	42	from FIAT.reference_element import ufc_cell, ufc_simplex # noqa: F401
41	43	from FIAT.hdivcurl import Hdiv, Hcurl # noqa: F401
42	44
43		__version__ = pkg_resources.get_distribution("FIAT").version
	45	__version__ = pkg_resources.get_distribution("fenics-fiat").version
44	46
45	47	# List of supported elements and mapping to element classes
46	48	supported_elements = {"Argyris": Argyris,

+11

-6

FIAT/argyris.py less more

17	17	from __future__ import absolute_import, print_function, division
18	18
19	19	from FIAT import finite_element, polynomial_set, dual_set, functional
	20	from FIAT.reference_element import TRIANGLE
20	21
21	22
22	23	class ArgyrisDualSet(dual_set.DualSet):

29	30	verts = ref_el.get_vertices()
30	31	sd = ref_el.get_spatial_dimension()
31	32
32		if sd != 2:
33		raise Exception("Illegal spatial dimension")
	33	if ref_el.get_shape() != TRIANGLE:
	34	raise ValueError("Argyris only defined on triangles")
34	35
35	36	pe = functional.PointEvaluation
36	37	pd = functional.PointDerivative

49	50	nodes.append(pd(ref_el, verts[v], alpha))
50	51
51	52	# second derivatives
52		alphas = [[2, 0], [0, 2], [1, 1]]
	53	alphas = [[2, 0], [1, 1], [0, 2]]
53	54	for alpha in alphas:
54	55	nodes.append(pd(ref_el, verts[v], alpha))
55	56

82	83	nodes.extend(internalnds)
83	84	entity_ids[2][0] = list(range(cur, cur + len(internalpts)))
84	85	cur += len(internalpts)
	86	else:
	87	entity_ids[2] = {0: []}
85	88
86	89	super(ArgyrisDualSet, self).__init__(nodes, ref_el, entity_ids)
87	90

97	100	top = ref_el.get_topology()
98	101	verts = ref_el.get_vertices()
99	102	sd = ref_el.get_spatial_dimension()
100		if sd != 2:
101		raise Exception("Illegal spatial dimension")
	103	if ref_el.get_shape() != TRIANGLE:
	104	raise ValueError("Argyris only defined on triangles")
102	105
103	106	pd = functional.PointDerivative
104	107

115	118	nodes.append(pd(ref_el, verts[v], alpha))
116	119
117	120	# second derivatives
118		alphas = [[2, 0], [0, 2], [1, 1]]
	121	alphas = [[2, 0], [1, 1], [0, 2]]
119	122	for alpha in alphas:
120	123	nodes.append(pd(ref_el, verts[v], alpha))
121	124

130	133	nodes.append(n)
131	134	entity_ids[1][e] = [cur]
132	135	cur += 1
	136
	137	entity_ids[2] = {0: []}
133	138
134	139	super(QuinticArgyrisDualSet, self).__init__(nodes, ref_el, entity_ids)
135	140

+36

-68

FIAT/enriched.py less more

16	16
17	17	from __future__ import absolute_import, print_function, division
18	18
19		import numpy as np
20		from copy import copy
	19	from itertools import chain
	20
	21	import numpy
21	22
22	23	from FIAT.finite_element import FiniteElement
23	24	from FIAT.dual_set import DualSet
	25	from FIAT.mixed import concatenate_entity_dofs
	26
24	27
25	28	__all__ = ['EnrichedElement']
26	29

34	37	"""
35	38
36	39	def __init__(self, *elements):
37		assert len(elements) == 2, "EnrichedElement only implemented for two subelements"
38		A, B = elements
39
40	40	# Firstly, check it makes sense to enrich. Elements must have:
41	41	# - same reference element
42	42	# - same mapping
43	43	# - same value shape
44		if not A.get_reference_element() == B.get_reference_element():
	44	if len(set(e.get_reference_element() for e in elements)) > 1:
45	45	raise ValueError("Elements must be defined on the same reference element")
46		if not A.mapping()[0] == B.mapping()[0]:
	46	if len(set(m for e in elements for m in e.mapping())) > 1:
47	47	raise ValueError("Elements must have same mapping")
48		if not A.value_shape() == B.value_shape():
	48	if len(set(e.value_shape() for e in elements)) > 1:
49	49	raise ValueError("Elements must have the same value shape")
50	50
51	51	# order is at least max, possibly more, though getting this
52	52	# right isn't important AFAIK
53		order = max(A.get_order(), B.get_order())
	53	order = max(e.get_order() for e in elements)
54	54	# form degree is essentially max (not true for Hdiv/Hcurl,
55	55	# but this will raise an error above anyway).
56	56	# E.g. an H^1 function enriched with an L^2 is now just L^2.
57		if A.get_formdegree() is None or B.get_formdegree() is None:
	57	if any(e.get_formdegree() is None for e in elements):
58	58	formdegree = None
59	59	else:
60		formdegree = max(A.get_formdegree(), B.get_formdegree())
	60	formdegree = max(e.get_formdegree() for e in elements)
61	61
62	62	# set up reference element and mapping, following checks above
63		ref_el = A.get_reference_element()
64		mapping = A.mapping()[0]
	63	ref_el, = set(e.get_reference_element() for e in elements)
	64	mapping, = set(m for e in elements for m in e.mapping())
65	65
66	66	# set up entity_ids - for each geometric entity, just concatenate
67	67	# the entities of the constituent elements
68		Adofs = A.entity_dofs()
69		Bdofs = B.entity_dofs()
70		offset = A.space_dimension() # number of entities belonging to A
71		entity_ids = {}
72
73		for ent_dim in Adofs:
74		entity_ids[ent_dim] = {}
75		for ent_dim_index in Adofs[ent_dim]:
76		entlist = copy(Adofs[ent_dim][ent_dim_index])
77		entlist += [c + offset for c in Bdofs[ent_dim][ent_dim_index]]
78		entity_ids[ent_dim][ent_dim_index] = entlist
	68	entity_ids = concatenate_entity_dofs(ref_el, elements)
79	69
80	70	# set up dual basis - just concatenation
81		nodes = A.dual_basis() + B.dual_basis()
	71	nodes = list(chain.from_iterable(e.dual_basis() for e in elements))
82	72	dual = DualSet(nodes, ref_el, entity_ids)
83	73
84	74	super(EnrichedElement, self).__init__(ref_el, dual, order, formdegree, mapping)
85	75
86		# Set up constituent elements
87		self.A = A
88		self.B = B
89
90	76	# required degree (for quadrature) is definitely max
91		self.polydegree = max(A.degree(), B.degree())
	77	self.polydegree = max(e.degree() for e in elements)
92	78
93	79	# Store subelements
94	80	self._elements = elements

115	101	"""Return tabulated values of derivatives up to given order of
116	102	basis functions at given points."""
117	103
118		# Again, simply concatenate at the basis-function level
119		# Number of array dimensions depends on whether the space
120		# is scalar- or vector-valued, so treat these separately.
	104	num_components = numpy.prod(self.value_shape())
	105	table_shape = (self.space_dimension(), num_components, len(points))
121	106
122		Asd = self.A.space_dimension()
123		Bsd = self.B.space_dimension()
124		Atab = self.A.tabulate(order, points, entity)
125		Btab = self.B.tabulate(order, points, entity)
126		npoints = len(points)
127		vs = self.A.value_shape()
128		rank = len(vs) # scalar: 0, vector: 1
	107	table = {}
	108	irange = slice(0)
	109	for element in self._elements:
129	110
130		result = {}
131		for index in Atab:
132		if rank == 0:
133		# scalar valued
134		# Atab[index] and Btab[index] look like
135		# array[basis_fn][point]
136		# We build a new array, which will be the concatenation
137		# of the two subarrays, in the first index.
	111	etable = element.tabulate(order, points, entity)
	112	irange = slice(irange.stop, irange.stop + element.space_dimension())
138	113
139		temp = np.zeros((Asd + Bsd, npoints),
140		dtype=Atab[index].dtype)
141		temp[:Asd, :] = Atab[index][:, :]
142		temp[Asd:, :] = Btab[index][:, :]
	114	# Insert element table into table
	115	for dtuple in etable.keys():
143	116
144		result[index] = temp
145		elif rank == 1:
146		# vector valued
147		# Atab[index] and Btab[index] look like
148		# array[basis_fn][x/y/z][point]
149		# We build a new array, which will be the concatenation
150		# of the two subarrays, in the first index.
	117	if dtuple not in table:
	118	if num_components == 1:
	119	table[dtuple] = numpy.zeros((self.space_dimension(), len(points)),
	120	dtype=etable[dtuple].dtype)
	121	else:
	122	table[dtuple] = numpy.zeros(table_shape,
	123	dtype=etable[dtuple].dtype)
151	124
152		temp = np.zeros((Asd + Bsd, vs[0], npoints),
153		dtype=Atab[index].dtype)
154		temp[:Asd, :, :] = Atab[index][:, :, :]
155		temp[Asd:, :, :] = Btab[index][:, :, :]
	125	table[dtuple][irange][:] = etable[dtuple]
156	126
157		result[index] = temp
158		else:
159		raise NotImplementedError("must be scalar- or vector-valued")
160		return result
	127	return table
161	128
162	129	def value_shape(self):
163	130	"""Return the value shape of the finite element functions."""
164		return self.A.value_shape()
	131	result, = set(e.value_shape() for e in self._elements)
	132	return result
165	133
166	134	def dmats(self):
167	135	"""Return dmats: expansion coefficients for basis function

+1

-1

FIAT/finite_element.py less more

55	55	return self.dual
56	56
57	57	def get_order(self):
58		"""Return the order of the element (may be different from the degree."""
	58	"""Return the order of the element (may be different from the degree)."""
59	59	return self.order
60	60
61	61	def dual_basis(self):

+18

-0

FIAT/functional.py less more

227	227	dpt_dict = {pt: [(n[i], alphas[i], tuple()) for i in range(sd)]}
228	228
229	229	Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv")
	230
	231	def to_riesz(self, poly_set):
	232	x = list(self.deriv_dict.keys())[0]
	233
	234	X = sympy.DeferredVector('x')
	235	dx = numpy.asarray([X[i] for i in range(len(x))])
	236
	237	es = poly_set.get_expansion_set()
	238	ed = poly_set.get_embedded_degree()
	239
	240	bfs = es.tabulate(ed, [dx])[:, 0]
	241
	242	# We need the gradient dotted with the normal.
	243	return numpy.asarray(
	244	[sympy.lambdify(
	245	X, sum([sympy.diff(b, dxi)*ni
	246	for dxi, ni in zip(dx, self.n)]))(x)
	247	for b in bfs])
230	248
231	249
232	250	class IntegralMoment(Functional):

+1

-1

FIAT/hdiv_trace.py less more

203	203	# If not successful, return NaNs
204	204	if not success:
205	205	for key in phivals:
206		phivals[key] = np.full(shape=(sd, len(points)), fill_value=np.nan)
	206	phivals[key] = np.full(shape=(self.space_dimension(), len(points)), fill_value=np.nan)
207	207
208	208	return phivals
209	209

+21

-13

FIAT/hermite.py less more

0	0	# Copyright (C) 2008 Robert C. Kirby (Texas Tech University)
	1	# Modified 2017 by RCK
1	2	#
2	3	# This file is part of FIAT.
3	4	#

50	51	entity_ids[0][v] = list(range(cur, cur + 1 + sd))
51	52	cur += sd + 1
52	53
	54	# now only have dofs at the barycenter, which is the
	55	# maximal dimension
53	56	# no edge dof
	57
54	58	entity_ids[1] = {}
	59	for i in top[1]:
	60	entity_ids
	61	entity_ids[1][i] = []
55	62
56		# face dof
57		# point evaluation at barycenter
58		entity_ids[2] = {}
59		for f in sorted(top[2]):
60		pt = ref_el.make_points(2, f, 3)[0]
61		n = functional.PointEvaluation(ref_el, pt)
62		nodes.append(n)
63		entity_ids[2] = list(range(cur, cur + 1))
64		cur += 1
	63	if sd > 1:
	64	# face dof
	65	# point evaluation at barycenter
	66	entity_ids[2] = {}
	67	for f in sorted(top[2]):
	68	pt = ref_el.make_points(2, f, 3)[0]
	69	n = functional.PointEvaluation(ref_el, pt)
	70	nodes.append(n)
	71	entity_ids[2][f] = list(range(cur, cur + 1))
	72	cur += 1
65	73
66		for dim in range(3, sd + 1):
67		entity_ids[dim] = {}
68		for facet in top[dim]:
69		entity_ids[dim][facet] = []
	74	for dim in range(3, sd + 1):
	75	entity_ids[dim] = {}
	76	for facet in top[dim]:
	77	entity_ids[dim][facet] = []
70	78
71	79	super(CubicHermiteDualSet, self).__init__(nodes, ref_el, entity_ids)
72	80

+124

-0

FIAT/mixed.py less more

	0	# -- coding: utf-8 --
	1	#
	2	# Copyright (C) 2005-2010 Anders Logg
	3	#
	4	# This file is part of FIAT.
	5	#
	6	# FIAT is free software: you can redistribute it and/or modify
	7	# it under the terms of the GNU Lesser General Public License as published by
	8	# the Free Software Foundation, either version 3 of the License, or
	9	# (at your option) any later version.
	10	#
	11	# FIAT is distributed in the hope that it will be useful,
	12	# but WITHOUT ANY WARRANTY; without even the implied warranty of
	13	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	14	# GNU Lesser General Public License for more details.
	15	#
	16	# You should have received a copy of the GNU Lesser General Public License
	17	# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
	18
	19	from __future__ import absolute_import, print_function, division
	20	from six import iteritems
	21	from six.moves import map
	22
	23	import numpy
	24
	25	from operator import add
	26	from functools import partial
	27
	28	from FIAT.dual_set import DualSet
	29	from FIAT.finite_element import FiniteElement
	30
	31
	32	class MixedElement(FiniteElement):
	33	"""A FIAT-like representation of a mixed element.
	34
	35	:arg elements: An iterable of FIAT elements.
	36	:arg ref_el: The reference element (optional).
	37
	38	This object offers tabulation of the concatenated basis function
	39	tables along with an entity_dofs dict."""
	40	def __init__(self, elements, ref_el=None):
	41	elements = tuple(elements)
	42
	43	cells = set(e.get_reference_element() for e in elements)
	44	if ref_el is not None:
	45	cells.add(ref_el)
	46	ref_el, = cells
	47
	48	# These functionals are absolutely wrong, they all map from
	49	# functions of the wrong shape, and potentially of different
	50	# shapes. However, they are wrong precisely as FFC hacks
	51	# expect them to be. :(
	52	nodes = [L for e in elements for L in e.dual_basis()]
	53
	54	entity_dofs = concatenate_entity_dofs(ref_el, elements)
	55
	56	dual = DualSet(nodes, ref_el, entity_dofs)
	57	super(MixedElement, self).__init__(ref_el, dual, None, mapping=None)
	58	self._elements = elements
	59
	60	def elements(self):
	61	return self._elements
	62
	63	def num_sub_elements(self):
	64	return len(self._elements)
	65
	66	def value_shape(self):
	67	return (sum(numpy.prod(e.value_shape(), dtype=int) for e in self.elements()), )
	68
	69	def mapping(self):
	70	return [m for e in self._elements for m in e.mapping()]
	71
	72	def get_nodal_basis(self):
	73	raise NotImplementedError("get_nodal_basis not implemented")
	74
	75	def tabulate(self, order, points, entity=None):
	76	"""Tabulate a mixed element by appropriately splatting
	77	together the tabulation of the individual elements.
	78	"""
	79	shape = (self.space_dimension(),) + self.value_shape() + (len(points),)
	80
	81	output = {}
	82
	83	sub_dims = [0] + list(e.space_dimension() for e in self.elements())
	84	sub_cmps = [0] + list(numpy.prod(e.value_shape(), dtype=int)
	85	for e in self.elements())
	86	irange = numpy.cumsum(sub_dims)
	87	crange = numpy.cumsum(sub_cmps)
	88
	89	for i, e in enumerate(self.elements()):
	90	table = e.tabulate(order, points, entity)
	91
	92	for d, tab in iteritems(table):
	93	try:
	94	arr = output[d]
	95	except KeyError:
	96	arr = numpy.zeros(shape, dtype=tab.dtype)
	97	output[d] = arr
	98
	99	ir = irange[i:i+2]
	100	cr = crange[i:i+2]
	101	tab = tab.reshape(ir[1] - ir[0], cr[1] - cr[0], -1)
	102	arr[slice(ir), slice(cr)] = tab
	103
	104	return output
	105
	106	def is_nodal(self):
	107	"""True if primal and dual bases are orthogonal."""
	108	return all(e.is_nodal() for e in self._elements)
	109
	110
	111	def concatenate_entity_dofs(ref_el, elements):
	112	"""Combine the entity_dofs from a list of elements into a combined
	113	entity_dof containing the information for the concatenated DoFs of
	114	all the elements."""
	115	entity_dofs = {dim: {i: [] for i in entities}
	116	for dim, entities in iteritems(ref_el.get_topology())}
	117	offsets = numpy.cumsum([0] + list(e.space_dimension()
	118	for e in elements), dtype=int)
	119	for i, d in enumerate(e.entity_dofs() for e in elements):
	120	for dim, dofs in iteritems(d):
	121	for ent, off in iteritems(dofs):
	122	entity_dofs[dim][ent] += list(map(partial(add, offsets[i]), off))
	123	return entity_dofs

+5

-3

FIAT/morley.py less more

17	17	from __future__ import absolute_import, print_function, division
18	18
19	19	from FIAT import finite_element, polynomial_set, dual_set, functional
	20	from FIAT.reference_element import TRIANGLE
20	21
21	22
22	23	class MorleyDualSet(dual_set.DualSet):

33	34	# need to do this dimension-by-dimension, facet-by-facet
34	35	top = ref_el.get_topology()
35	36	verts = ref_el.get_vertices()
36		sd = ref_el.get_spatial_dimension()
37		if sd != 2:
38		raise Exception("Illegal spatial dimension")
	37	if ref_el.get_shape() != TRIANGLE:
	38	raise ValueError("Morley only defined on triangles")
39	39
40	40	# vertex point evaluations
41	41

55	55	entity_ids[1][e] = [cur]
56	56	cur += 1
57	57
	58	entity_ids[2] = {0: []}
	59
58	60	super(MorleyDualSet, self).__init__(nodes, ref_el, entity_ids)
59	61
60	62

+82

-0

FIAT/quadrature_element.py less more

	0	# -- coding: utf-8 --
	1
	2	# Copyright (C) 2007-2016 Kristian B. Oelgaard
	3	# Copyright (C) 2017 Miklós Homolya
	4	#
	5	# This file is part of FIAT.
	6	#
	7	# FIAT is free software: you can redistribute it and/or modify
	8	# it under the terms of the GNU Lesser General Public License as published by
	9	# the Free Software Foundation, either version 3 of the License, or
	10	# (at your option) any later version.
	11	#
	12	# FIAT is distributed in the hope that it will be useful,
	13	# but WITHOUT ANY WARRANTY; without even the implied warranty of
	14	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	15	# GNU Lesser General Public License for more details.
	16	#
	17	# You should have received a copy of the GNU Lesser General Public License
	18	# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
	19	#
	20	# Modified by Garth N. Wells 2006-2009
	21
	22	from __future__ import absolute_import, print_function, division
	23	from six.moves import range
	24
	25	# Python modules.
	26	import numpy
	27	import six
	28
	29	# FIAT modules.
	30	from FIAT.dual_set import DualSet
	31	from FIAT.finite_element import FiniteElement
	32	from FIAT.functional import PointEvaluation
	33
	34
	35	class QuadratureElement(FiniteElement):
	36	"""A set of quadrature points pretending to be a finite element."""
	37
	38	def __init__(self, ref_el, points):
	39	# Create entity dofs.
	40	entity_dofs = {dim: {entity: [] for entity in entities}
	41	for dim, entities in six.iteritems(ref_el.get_topology())}
	42	entity_dofs[ref_el.get_dimension()] = {0: list(range(len(points)))}
	43
	44	# The dual nodes are PointEvaluations at the quadrature points.
	45	# FIXME: KBO: Check if this gives expected results for code like evaluate_dof.
	46	nodes = [PointEvaluation(ref_el, tuple(point)) for point in points]
	47
	48	# Construct the dual set
	49	dual = DualSet(nodes, ref_el, entity_dofs)
	50
	51	super(QuadratureElement, self).__init__(ref_el, dual, order=None)
	52	self._points = points # save the quadrature points
	53
	54	def value_shape(self):
	55	"The QuadratureElement is scalar valued"
	56	return ()
	57
	58	def tabulate(self, order, points, entity=None):
	59	"""Return the identity matrix of size (num_quad_points, num_quad_points),
	60	in a format that monomialintegration and monomialtabulation understands."""
	61
	62	if entity is not None and entity != (self.ref_el.get_dimension(), 0):
	63	raise ValueError('QuadratureElement does not "tabulate" on subentities.')
	64
	65	# Derivatives are not defined on a QuadratureElement
	66	if order:
	67	raise ValueError("Derivatives are not defined on a QuadratureElement.")
	68
	69	# Check that incoming points are equal to the quadrature points.
	70	if len(points) != len(self._points) or abs(numpy.array(points) - self._points).max() > 1e-12:
	71	raise AssertionError("Mismatch of quadrature points!")
	72
	73	# Return the identity matrix of size len(self._points).
	74	values = numpy.eye(len(self._points))
	75	dim = self.ref_el.get_spatial_dimension()
	76	return {(0,) * dim: values}
	77
	78	@staticmethod
	79	def is_nodal():
	80	# No polynomial basis, but still nodal.
	81	return True

+2

-8

FIAT/quadrature_schemes.py less more

43	43	from numpy import array, arange, float64
44	44
45	45	# FIAT
46		from FIAT.reference_element import QUADRILATERAL, TENSORPRODUCT, UFCTriangle, UFCTetrahedron
	46	from FIAT.reference_element import QUADRILATERAL, HEXAHEDRON, TENSORPRODUCT, UFCTriangle, UFCTetrahedron
47	47	from FIAT.quadrature import QuadratureRule, make_quadrature, make_tensor_product_quadrature
48	48
49	49

72	72	for c, d in zip(ref_el.cells, degree)]
73	73	return make_tensor_product_quadrature(*quad_rules)
74	74
75		if ref_el.get_shape() == QUADRILATERAL:
	75	if ref_el.get_shape() in [QUADRILATERAL, HEXAHEDRON]:
76	76	return create_quadrature(ref_el.product, degree, scheme)
77	77
78	78	if degree < 0:

98	98
99	99	# Number of points per axis for exact integration
100	100	num_points_per_axis = (degree + 1 + 1) // 2
101
102		# Check for excess
103		if num_points_per_axis > 30:
104		dim = ref_el.get_spatial_dimension()
105		raise RuntimeError("Requested a quadrature rule with %d points per direction (%d points)" %
106		(num_points_per_axis, num_points_per_axis**dim))
107	101
108	102	# Create and return FIAT quadrature rule
109	103	return make_quadrature(ref_el, num_points_per_axis)

+157

-63

FIAT/reference_element.py less more

30	30	"""
31	31	from __future__ import absolute_import, print_function, division
32	32
33		from itertools import chain, product
	33	from itertools import chain, product, count
34	34	import operator
35	35	from math import factorial
36	36
37		from six import iteritems, itervalues
	37	from six import iteritems
38	38	from six import string_types
39		from six.moves import reduce
	39	from six.moves import reduce, range
40	40	import numpy
41	41	from numpy import ravel_multi_index, transpose
	42	from collections import defaultdict
42	43
43	44
44	45	POINT = 0

46	47	TRIANGLE = 2
47	48	TETRAHEDRON = 3
48	49	QUADRILATERAL = 11
	50	HEXAHEDRON = 111
49	51	TENSORPRODUCT = 99
50	52
51	53

392	394	n = Simplex.compute_normal(self, facet_i) # skip UFC overrides
393	395	return n / numpy.linalg.norm(n, numpy.inf)
394	396
395		def get_facet_transform(self, facet_i):
396		"""Return a function f such that for a point with facet coordinates
397		x_f on facet_i, x_c = f(x_f) is the corresponding cell coordinates.
	397	def get_entity_transform(self, dim, entity):
	398	"""Returns a mapping of point coordinates from the
	399	`entity`-th subentity of dimension `dim` to the cell.
	400
	401	:arg dim: subentity dimension (integer)
	402	:arg entity: entity number (integer)
398	403	"""
399		t = self.get_topology()
	404	topology = self.get_topology()
	405	celldim = self.get_spatial_dimension()
	406	codim = celldim - dim
	407	if dim == 0:
	408	# Special case vertices.
	409	i, = topology[dim][entity]
	410	vertex = self.get_vertices()[i]
	411	return lambda point: vertex
	412	elif dim == celldim:
	413	assert entity == 0
	414	return lambda point: point
400	415
401	416	try:
402		f_el = self.get_facet_element()
	417	subcell = self.construct_subelement(dim)
403	418	except NotImplementedError:
404	419	# Special case for 1D elements.
405		x_c = self.get_vertices_of_subcomplex(t[0][facet_i])[0]
406
	420	x_c, = self.get_vertices_of_subcomplex(topology[0][entity])
407	421	return lambda x: x_c
408	422
409		sd_c = self.get_spatial_dimension()
410		sd_f = f_el.get_spatial_dimension()
411
412		# Facet vertices in facet space.
413		v_f = numpy.array(f_el.get_vertices())
414
415		A = numpy.zeros([sd_f, sd_f])
416
417		for i in range(A.shape[0]):
418		A[i, :] = (v_f[i + 1] - v_f[0])
	423	subdim = subcell.get_spatial_dimension()
	424
	425	assert subdim == celldim - codim
	426
	427	# Entity vertices in entity space.
	428	v_e = numpy.asarray(subcell.get_vertices())
	429
	430	A = numpy.zeros([subdim, subdim])
	431
	432	for i in range(subdim):
	433	A[i, :] = (v_e[i + 1] - v_e[0])
419	434	A[i, :] /= A[i, :].dot(A[i, :])
420	435
421		# Facet vertices in cell space.
422		v_c = numpy.array(self.get_vertices_of_subcomplex(t[sd_c - 1][facet_i]))
423
424		B = numpy.zeros([sd_c, sd_f])
425
426		for j in range(B.shape[1]):
	436	# Entity vertices in cell space.
	437	v_c = numpy.asarray(self.get_vertices_of_subcomplex(topology[dim][entity]))
	438
	439	B = numpy.zeros([celldim, subdim])
	440
	441	for j in range(subdim):
427	442	B[:, j] = (v_c[j + 1] - v_c[0])
428	443
429	444	C = B.dot(A)
430	445
431		offset = v_c[0] - C.dot(v_f[0])
	446	offset = v_c[0] - C.dot(v_e[0])
432	447
433	448	return lambda x: offset + C.dot(x)
434	449

436	451	"""Returns the subelement dimension of the cell. Same as the
437	452	spatial dimension."""
438	453	return self.get_spatial_dimension()
439
440		def get_entity_transform(self, dim, entity_i):
441		"""Returns a mapping of point coordinates from the
442		`entity_i`-th subentity of dimension `dim` to the cell.
443
444		:arg dim: subentity dimension (integer)
445		:arg entity_i: entity number (integer)
446		"""
447		space_dim = self.get_spatial_dimension()
448		if dim == space_dim: # cell points
449		assert entity_i == 0
450		return lambda point: point
451		elif dim == space_dim - 1: # facet points
452		return self.get_facet_transform(entity_i)
453		elif dim == 0: # vertex point
454		i, = self.get_topology()[dim][entity_i]
455		vertex = self.get_vertices()[i]
456		return lambda point: vertex
457		else:
458		raise NotImplementedError("Co-dimension >1 not implemented.")
459	454
460	455
461	456	# Backwards compatible name

772	767	True)
773	768
774	769
775		class FiredrakeQuadrilateral(Cell):
	770	class UFCQuadrilateral(Cell):
776	771	"""This is the reference quadrilateral with vertices
777	772	(0.0, 0.0), (0.0, 1.0), (1.0, 0.0) and (1.0, 1.0)."""
778	773

781	776	pt = product.get_topology()
782	777
783	778	verts = product.get_vertices()
784		topology = {0: pt[(0, 0)],
785		1: dict(enumerate(list(itervalues(pt[(0, 1)])) + list(itervalues(pt[(1, 0)])))),
786		2: pt[(1, 1)]}
787		super(FiredrakeQuadrilateral, self).__init__(QUADRILATERAL, verts, topology)
	779	topology = flatten_entities(pt)
	780
	781	super(UFCQuadrilateral, self).__init__(QUADRILATERAL, verts, topology)
	782
788	783	self.product = product
	784	self.unflattening_map = compute_unflattening_map(pt)
789	785
790	786	def get_dimension(self):
791	787	"""Returns the subelement dimension of the cell. Same as the

814	810	:arg dim: entity dimension (integer)
815	811	:arg entity_i: entity number (integer)
816	812	"""
817		d, e = {0: lambda e: ((0, 0), e),
818		1: lambda e: {0: ((0, 1), 0),
819		1: ((0, 1), 1),
820		2: ((1, 0), 0),
821		3: ((1, 0), 1)}[e],
822		2: lambda e: ((1, 1), e)}[dim](entity_i)
	813	d, e = self.unflattening_map[(dim, entity_i)]
823	814	return self.product.get_entity_transform(d, e)
824	815
825	816	def volume(self):

829	820	def compute_reference_normal(self, facet_dim, facet_i):
830	821	"""Returns the unit normal in infinity norm to facet_i."""
831	822	assert facet_dim == 1
832		d, i = {0: ((0, 1), 0),
833		1: ((0, 1), 1),
834		2: ((1, 0), 0),
835		3: ((1, 0), 1)}[facet_i]
	823	d, i = self.unflattening_map[(facet_dim, facet_i)]
	824	return self.product.compute_reference_normal(d, i)
	825
	826	def contains_point(self, point, epsilon=0):
	827	"""Checks if reference cell contains given point
	828	(with numerical tolerance)."""
	829	return self.product.contains_point(point, epsilon=epsilon)
	830
	831
	832	class UFCHexahedron(Cell):
	833	"""This is the reference hexahedron with vertices
	834	(0.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0), (0.0, 1.0, 1.0),
	835	(1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (1.0, 1.0, 0.0) and (1.0, 1.0, 1.0)."""
	836
	837	def __init__(self):
	838	product = TensorProductCell(UFCInterval(), UFCInterval(), UFCInterval())
	839	pt = product.get_topology()
	840
	841	verts = product.get_vertices()
	842	topology = flatten_entities(pt)
	843
	844	super(UFCHexahedron, self).__init__(HEXAHEDRON, verts, topology)
	845
	846	self.product = product
	847	self.unflattening_map = compute_unflattening_map(pt)
	848
	849	def get_dimension(self):
	850	"""Returns the subelement dimension of the cell. Same as the
	851	spatial dimension."""
	852	return self.get_spatial_dimension()
	853
	854	def construct_subelement(self, dimension):
	855	"""Constructs the reference element of a cell subentity
	856	specified by subelement dimension.
	857
	858	:arg dimension: subentity dimension (integer)
	859	"""
	860	if dimension == 3:
	861	return self
	862	elif dimension == 2:
	863	return UFCQuadrilateral()
	864	elif dimension == 1:
	865	return UFCInterval()
	866	elif dimension == 0:
	867	return Point()
	868	else:
	869	raise ValueError("Invalid dimension: %d" % (dimension,))
	870
	871	def get_entity_transform(self, dim, entity_i):
	872	"""Returns a mapping of point coordinates from the
	873	`entity_i`-th subentity of dimension `dim` to the cell.
	874
	875	:arg dim: entity dimension (integer)
	876	:arg entity_i: entity number (integer)
	877	"""
	878	d, e = self.unflattening_map[(dim, entity_i)]
	879	return self.product.get_entity_transform(d, e)
	880
	881	def volume(self):
	882	"""Computes the volume in the appropriate dimensional measure."""
	883	return self.product.volume()
	884
	885	def compute_reference_normal(self, facet_dim, facet_i):
	886	"""Returns the unit normal in infinity norm to facet_i."""
	887	assert facet_dim == 2
	888	d, i = self.unflattening_map[(facet_dim, facet_i)]
836	889	return self.product.compute_reference_normal(d, i)
837	890
838	891	def contains_point(self, point, epsilon=0):

918	971	# Tensor product cell
919	972	return TensorProductCell(map(ufc_cell, celltype.split(" ")))
920	973	elif celltype == "quadrilateral":
921		return FiredrakeQuadrilateral()
	974	return UFCQuadrilateral()
	975	elif celltype == "hexahedron":
	976	return UFCHexahedron()
922	977	elif celltype == "interval":
923	978	return ufc_simplex(1)
924	979	elif celltype == "triangle":

951	1006	p = numpy.prod([si for si in s if (si) > 1.e-10])
952	1007
953	1008	return p / factorial(sd)
	1009
	1010
	1011	def tuple_sum(tree):
	1012	"""
	1013	This function calculates the sum of elements in a tuple, it is needed to handle nested tuples in TensorProductCell.
	1014	Example: tuple_sum(((1, 0), 1)) returns 2
	1015	If input argument is not the tuple, returns input.
	1016	"""
	1017	if isinstance(tree, tuple):
	1018	return sum(map(tuple_sum, tree))
	1019	else:
	1020	return tree
	1021
	1022
	1023	def flatten_entities(topology_dict):
	1024	"""This function flattens topology dict of TensorProductCell and entity_dofs dict of TensorProductElement"""
	1025
	1026	flattened_entities = defaultdict(list)
	1027	for dim in sorted(topology_dict.keys()):
	1028	flat_dim = tuple_sum(dim)
	1029	flattened_entities[flat_dim] += [v for k, v in sorted(topology_dict[dim].items())]
	1030
	1031	return {dim: dict(enumerate(entities))
	1032	for dim, entities in iteritems(flattened_entities)}
	1033
	1034
	1035	def compute_unflattening_map(topology_dict):
	1036	"""This function returns unflattening map for the given tensor product topology dict."""
	1037
	1038	counter = defaultdict(count)
	1039	unflattening_map = {}
	1040
	1041	for dim, entities in sorted(iteritems(topology_dict)):
	1042	flat_dim = tuple_sum(dim)
	1043	for entity in entities:
	1044	flat_entity = next(counter[flat_dim])
	1045	unflattening_map[(flat_dim, flat_entity)] = (dim, entity)
	1046
	1047	return unflattening_map

+80

-1

FIAT/tensor_product.py less more

20	20
21	21	import numpy
22	22	from FIAT.finite_element import FiniteElement
23		from FIAT.reference_element import TensorProductCell
	23	from FIAT.reference_element import TensorProductCell, UFCQuadrilateral, UFCHexahedron, flatten_entities, compute_unflattening_map
	24	from FIAT.dual_set import DualSet
24	25	from FIAT.polynomial_set import mis
25	26	from FIAT import dual_set
26	27	from FIAT import functional

365	366	def get_num_members(self, arg):
366	367	"""Return number of members of the expansion set."""
367	368	raise NotImplementedError("get_num_members not implemented")
	369
	370	def is_nodal(self):
	371	# This element is nodal iff all factor elements are nodal.
	372	return all([self.A.is_nodal(), self.B.is_nodal()])
	373
	374
	375	class FlattenedDimensions(FiniteElement):
	376	"""A wrapper class that flattens entity dimensions of a FIAT element defined
	377	on a TensorProductCell to one with quadrilateral/hexahedron entities.
	378	TensorProductCell has dimension defined as a tuple of factor element dimensions
	379	(i, j) in 2D and (i, j, k) in 3D.
	380	Flattened dimension is a sum of the tuple elements."""
	381
	382	def __init__(self, element):
	383
	384	nodes = element.dual.nodes
	385	dim = element.ref_el.get_spatial_dimension()
	386
	387	if dim == 2:
	388	ref_el = UFCQuadrilateral()
	389	elif dim == 3:
	390	ref_el = UFCHexahedron()
	391	else:
	392	raise ValueError("Illegal element dimension %s" % dim)
	393
	394	entity_ids = element.dual.entity_ids
	395
	396	flat_entity_ids = flatten_entities(entity_ids)
	397	dual = DualSet(nodes, ref_el, flat_entity_ids)
	398	super(FlattenedDimensions, self).__init__(ref_el, dual, element.get_order(), element.get_formdegree(), element._mapping)
	399	self.element = element
	400
	401	# Construct unflattening map for passing correct values to tabulate()
	402	self.unflattening_map = compute_unflattening_map(self.element.ref_el.get_topology())
	403
	404	def degree(self):
	405	"""Return the degree of the (embedding) polynomial space."""
	406	return self.element.degree()
	407
	408	def tabulate(self, order, points, entity=None):
	409	"""Return tabulated values of derivatives up to given order of
	410	basis functions at given points."""
	411	if entity is None:
	412	entity = (self.get_reference_element().get_spatial_dimension(), 0)
	413
	414	# Entity is provided in flattened form (d, i)
	415	# Appropriate product entity is taken from the unflattening_map dict
	416	entity_dim, entity_id = entity
	417	product_entity = self.unflattening_map[(entity_dim, entity_id)]
	418
	419	return self.element.tabulate(order, points, product_entity)
	420
	421	def value_shape(self):
	422	"""Return the value shape of the finite element functions."""
	423	return self.element.value_shape()
	424
	425	def get_nodal_basis(self):
	426	"""Return the nodal basis, encoded as a PolynomialSet object,
	427	for the finite element."""
	428	raise self.element.get_nodal_basis()
	429
	430	def get_coeffs(self):
	431	"""Return the expansion coefficients for the basis of the
	432	finite element."""
	433	raise self.element.get_coeffs()
	434
	435	def dmats(self):
	436	"""Return dmats: expansion coefficients for basis function
	437	derivatives."""
	438	raise self.element.dmats()
	439
	440	def get_num_members(self, arg):
	441	"""Return number of members of the expansion set."""
	442	raise self.element.get_num_members(arg)
	443
	444	def is_nodal(self):
	445	# This element is nodal iff unflattened element is nodal.
	446	return self.element.is_nodal()

+20

-0

bitbucket-pipelines.yml less more

	0	image: ubuntu:16.04
	1
	2	pipelines:
	3	default:
	4	- step:
	5	script:
	6	- apt-get update
	7	- apt-get install -y
	8	git
	9	python-minimal python-pip python-pytest python-setuptools python-wheel
	10	python3-minimal python3-pip python3-pytest python3-setuptools python3-wheel
	11	--no-install-recommends
	12	- pip2 install flake8 flake8-future-import
	13	- pip3 install flake8 flake8-future-import
	14	- pip2 install .
	15	- pip3 install .
	16	- python2 -m flake8
	17	- python3 -m flake8
	18	- DATA_REPO_GIT="" python2 -m pytest -v test/
	19	- DATA_REPO_GIT="" python3 -m pytest -v test/

+1

-1

doc/sphinx/source/conf.py less more

55	55	project = u'FInite element Automatic Tabulator (FIAT)'
56	56	this_year = datetime.date.today().year
57	57	copyright = u'%s, FEniCS Project' % this_year
58		version = pkg_resources.get_distribution("FIAT").version
	58	version = pkg_resources.get_distribution("fenics-fiat").version
59	59	release = version
60	60
61	61	# The language for content autogenerated by Sphinx. Refer to documentation

+10

-0

doc/sphinx/source/releases/v2017.1.0.post1.rst less more

	0	=================================
	1	Changes in version 2017.1.0.post1
	2	=================================
	3
	4	FIAT 2017.1.0.post1 was released on 2017-09-12.
	5
	6	Summary of changes
	7	==================
	8
	9	- Change PyPI package name to fenics-fiat.

+11

-0

doc/sphinx/source/releases/v2017.2.0.rst less more

	0	===========================
	1	Changes in version 2017.2.0
	2	===========================
	3
	4	FIAT 2017.2.0 was released on 2017-09-19.
	5
	6	Summary of changes
	7	==================
	8
	9	- Add quadrilateral and hexahedron reference cells
	10	- Add quadrilateral and hexahedron elements (with a wrapping class for TensorProductElement)

+1

-0

doc/sphinx/source/releases.rst less more

8	8	:maxdepth: 2
9	9
10	10	releases/next
	11	releases/v2017.1.0.post1
11	12	releases/v2017.1.0
12	13	releases/v2016.2.0
13	14	releases/v2016.1.0

+8

-0

release.conf less more

	0	# Configuration file for fenics-release
	1
	2	PACKAGE="fiat"
	3	BRANCH="master"
	4	FILES="ChangeLog.rst \
	5	setup.py \
	6	doc/sphinx/source/releases/next.rst \
	7	doc/sphinx/source/releases.rst"

+3

-3

setup.py less more

12	12	print("Python 2.7 or higher required, please upgrade.")
13	13	sys.exit(1)
14	14
15		version = "2017.1.0"
	15	version = "2017.2.0"
16	16
17	17	url = "https://bitbucket.org/fenics-project/fiat/"
18	18	tarball = None
19	19	if 'dev' not in version:
20		tarball = url + "downloads/fiat-%s.tar.gz" % version
	20	tarball = url + "downloads/fenics-fiat-%s.tar.gz" % version
21	21
22		setup(name="FIAT",
	22	setup(name="fenics-fiat",
23	23	description="FInite element Automatic Tabulator",
24	24	version=version,
25	25	author="Robert C. Kirby et al.",

+15

-0

shippable.yml less more

	0	language: python
	1	python:
	2	- 2.7
	3	- 3.6
	4
	5	env:
	6	- DATA_REPO_GIT=""
	7
	8	build:
	9	ci:
	10	- pip install --upgrade pip
	11	- pip install flake8 flake8-future-import pytest
	12	- pip install .
	13	- python -m flake8 .
	14	- python -m pytest -v test/

+1

-0

test/regression/.gitignore less more

0

fiat-reference-data/

+70

-8

test/unit/test_fiat.py less more

30	30	from FIAT.raviart_thomas import RaviartThomas # noqa: F401
31	31	from FIAT.discontinuous_raviart_thomas import DiscontinuousRaviartThomas # noqa: F401
32	32	from FIAT.brezzi_douglas_marini import BrezziDouglasMarini # noqa: F401
	33	from FIAT.mixed import MixedElement
33	34	from FIAT.nedelec import Nedelec # noqa: F401
34	35	from FIAT.nedelec_second_kind import NedelecSecondKind # noqa: F401
35	36	from FIAT.regge import Regge # noqa: F401

40	41	from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre # noqa: F401
41	42	from FIAT.restricted import RestrictedElement # noqa: F401
42	43	from FIAT.tensor_product import TensorProductElement # noqa: F401
	44	from FIAT.tensor_product import FlattenedDimensions # noqa: F401
43	45	from FIAT.hdivcurl import Hdiv, Hcurl # noqa: F401
44	46	from FIAT.argyris import Argyris, QuinticArgyris # noqa: F401
45	47	from FIAT.hermite import CubicHermite # noqa: F401

194	196	" Regge(S, 1),"
195	197	" RestrictedElement(Regge(S, 2), restriction_domain='interior')"
196	198	")",
	199	"Argyris(T, 5)",
	200	"QuinticArgyris(T)",
	201	"CubicHermite(I)",
	202	"CubicHermite(T)",
	203	"CubicHermite(S)",
	204	"Morley(T)",
	205
	206	# MixedElement made of nodal elements should be nodal, but its API
	207	# is currently just broken.
	208	xfail_impl("MixedElement(["
	209	" DiscontinuousLagrange(T, 1),"
	210	" RaviartThomas(T, 2)"
	211	"])"),
197	212
198	213	# Following element do not bother implementing get_nodal_basis
199	214	# so the test would need to be rewritten using tabulate

209	224	"Hcurl(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), "
210	225	"Hcurl(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))"
211	226	")"),
212
213		# These elements have broken constructor
214		xfail_key("Argyris(T, 1)",),
215		xfail_key("QuinticArgyris(T)",),
216		xfail_key("CubicHermite(I)",),
217		xfail_key("CubicHermite(T)",),
218		xfail_key("CubicHermite(S)",),
219		xfail_key("Morley(T)",),
	227	# Following elements are checked using tabulate
	228	xfail_impl("TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))"),
	229	xfail_impl("TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))"),
	230	xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)), Lagrange(I, 1))"),
	231	xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)), Lagrange(I, 2))"),
	232	xfail_impl("FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)))"),
	233	xfail_impl("FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)))"),
	234	xfail_impl("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))), Lagrange(I, 1)))"),
	235	xfail_impl("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))), Lagrange(I, 2)))"),
220	236	]
221	237
222	238

302	318	e1.get_dual_set().to_riesz(e1.get_nodal_basis()))
303	319
304	320
	321	def test_mixed_is_nodal():
	322	element = MixedElement([DiscontinuousLagrange(T, 1), RaviartThomas(T, 2)])
	323
	324	assert element.is_nodal()
	325
	326
	327	def test_mixed_is_not_nodal():
	328	element = MixedElement([
	329	EnrichedElement(
	330	RaviartThomas(T, 1),
	331	RestrictedElement(RaviartThomas(T, 2), restriction_domain="interior")
	332	),
	333	DiscontinuousLagrange(T, 1)
	334	])
	335
	336	assert not element.is_nodal()
	337
	338
	339	@pytest.mark.parametrize('element', [
	340	"TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))",
	341	"TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))",
	342	"TensorProductElement(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)), Lagrange(I, 1))",
	343	"TensorProductElement(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)), Lagrange(I, 2))",
	344	"FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)))",
	345	"FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)))",
	346	"FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))), Lagrange(I, 1)))",
	347	"FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))), Lagrange(I, 2)))",
	348	])
	349	def test_nodality_tpe(element):
	350	"""Check that generated TensorProductElements are nodal, i.e. nodes evaluated
	351	on basis functions give Kronecker delta
	352	"""
	353	# Instantiate element
	354	element = eval(element)
	355
	356	# Get nodes coordinates
	357	nodes_coord = list(next(iter(f.get_point_dict().keys())) for f in element.dual_basis())
	358
	359	# Check nodality
	360	for j, x in enumerate(nodes_coord):
	361	table = element.tabulate(0, (x,))
	362	basis = table[list(table.keys())[0]]
	363	for i in range(len(basis)):
	364	assert np.isclose(basis[i], 1.0 if i == j else 0.0)
	365
	366
305	367	if __name__ == '__main__':
306	368	import os
307	369	pytest.main(os.path.abspath(__file__))

+15

-20

test/unit/test_quadrature.py less more

22	22	import pytest
23	23	import FIAT
24	24	from FIAT.reference_element import UFCInterval, UFCTriangle, UFCTetrahedron
25		from FIAT.reference_element import FiredrakeQuadrilateral, TensorProductCell
	25	from FIAT.reference_element import UFCQuadrilateral, UFCHexahedron, TensorProductCell
26	26
27	27
28	28	@pytest.fixture(scope='module')

42	42
43	43	@pytest.fixture(scope='module')
44	44	def quadrilateral():
45		return FiredrakeQuadrilateral()
	45	return UFCQuadrilateral()
	46
	47
	48	@pytest.fixture(scope='module')
	49	def hexahedron():
	50	return UFCHexahedron()
46	51
47	52
48	53	@pytest.fixture(scope='module')

60	65	@pytest.fixture(scope='module')
61	66	def extr_quadrilateral():
62	67	"""Extruded quadrilateral = quadrilateral x interval"""
63		return TensorProductCell(FiredrakeQuadrilateral(), UFCInterval())
	68	return TensorProductCell(UFCQuadrilateral(), UFCInterval())
64	69
65	70
66	71	@pytest.fixture(params=["canonical", "default"])

115	120	(2*(degree + 2) - 2) / ((degree + 1)(degree + 2)))
116	121
117	122
	123	@pytest.mark.parametrize("degree", range(8))
	124	def test_create_quadrature_hexahedron(hexahedron, degree, scheme):
	125	q = FIAT.create_quadrature(hexahedron, degree, scheme)
	126	assert numpy.allclose(q.integrate(lambda x: sum(x)**degree),
	127	-3 * (2(degree + 3) - 3(degree + 2) - 1) / ((degree + 1)(degree + 2)(degree + 3)))
	128
	129
118	130	@pytest.mark.parametrize("extrdeg", range(4))
119	131	@pytest.mark.parametrize("basedeg", range(5))
120	132	def test_create_quadrature_extr_quadrilateral(extr_quadrilateral, basedeg, extrdeg, scheme):

138	150	def test_invalid_quadrature_degree_tensor_prod(cell):
139	151	with pytest.raises(ValueError):
140	152	FIAT.create_quadrature(cell, (-1, -1))
141
142
143		@pytest.mark.parametrize("cell", [interval(),
144		triangle(),
145		tetrahedron(),
146		quadrilateral()])
147		def test_high_degree_runtime_error(cell):
148		with pytest.raises(RuntimeError):
149		FIAT.create_quadrature(cell, 60)
150
151
152		@pytest.mark.parametrize("cell", [extr_interval(),
153		extr_triangle(),
154		extr_quadrilateral()])
155		def test_high_degree_runtime_error_tensor_prod(cell):
156		with pytest.raises(RuntimeError):
157		FIAT.create_quadrature(cell, (60, 60))
158	153
159	154
160	155	def test_tensor_product_composition(interval, triangle, extr_triangle, scheme):

+69

-0

test/unit/test_quadrature_element.py less more

	0	# Copyright (C) 2017 Miklos Homolya
	1	#
	2	# This file is part of FIAT.
	3	#
	4	# FIAT is free software: you can redistribute it and/or modify
	5	# it under the terms of the GNU Lesser General Public License as published by
	6	# the Free Software Foundation, either version 3 of the License, or
	7	# (at your option) any later version.
	8	#
	9	# FIAT is distributed in the hope that it will be useful,
	10	# but WITHOUT ANY WARRANTY; without even the implied warranty of
	11	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	12	# GNU Lesser General Public License for more details.
	13	#
	14	# You should have received a copy of the GNU Lesser General Public License
	15	# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
	16
	17	from __future__ import absolute_import, print_function, division
	18
	19	import pytest
	20	import numpy as np
	21
	22	from FIAT import QuadratureElement, make_quadrature, ufc_simplex
	23
	24
	25	@pytest.fixture(params=[1, 2, 3])
	26	def cell(request):
	27	return ufc_simplex(request.param)
	28
	29
	30	@pytest.fixture
	31	def quadrature(cell):
	32	return make_quadrature(cell, 2)
	33
	34
	35	@pytest.fixture
	36	def element(cell, quadrature):
	37	return QuadratureElement(cell, quadrature.get_points())
	38
	39
	40	def test_order(element, quadrature):
	41	with pytest.raises(ValueError):
	42	element.tabulate(1, quadrature.get_points())
	43
	44
	45	def test_points(element, quadrature):
	46	points = quadrature.get_points()
	47	wrong_points = np.linspace(0.0, 1.0, points.size).reshape(points.shape)
	48	with pytest.raises(AssertionError):
	49	element.tabulate(0, wrong_points)
	50
	51
	52	def test_entity(element, quadrature):
	53	dim = element.get_reference_element().get_spatial_dimension()
	54	points = make_quadrature(ufc_simplex(dim - 1), 2).get_points()
	55	with pytest.raises(ValueError):
	56	element.tabulate(0, points, entity=(dim - 1, 1))
	57
	58
	59	def test_result(element, quadrature):
	60	dim = element.get_reference_element().get_spatial_dimension()
	61	points = quadrature.get_points()
	62	actual = element.tabulate(0, points)[(0,) * dim]
	63	assert np.allclose(np.eye(len(points)), actual)
	64
	65
	66	if __name__ == '__main__':
	67	import os
	68	pytest.main(os.path.abspath(__file__))

+12

-4

test/unit/test_reference_element.py less more

20	20	import numpy as np
21	21
22	22	from FIAT.reference_element import UFCInterval, UFCTriangle, UFCTetrahedron
23		from FIAT.reference_element import Point, FiredrakeQuadrilateral, TensorProductCell
	23	from FIAT.reference_element import Point, TensorProductCell, UFCQuadrilateral, UFCHexahedron
24	24
25	25	point = Point()
26	26	interval = UFCInterval()
27	27	triangle = UFCTriangle()
28		quadrilateral = FiredrakeQuadrilateral()
	28	quadrilateral = UFCQuadrilateral()
	29	hexahedron = UFCHexahedron()
29	30	tetrahedron = UFCTetrahedron()
30	31	interval_x_interval = TensorProductCell(interval, interval)
31	32	triangle_x_interval = TensorProductCell(triangle, interval)

40	41	(tetrahedron, 1/6),
41	42	(interval_x_interval, 1),
42	43	(triangle_x_interval, 1/2),
43		(quadrilateral_x_interval, 1)])
	44	(quadrilateral_x_interval, 1),
	45	(hexahedron, 1)])
44	46	def test_volume(cell, volume):
45	47	assert np.allclose(volume, cell.volume())
46	48

58	60	(tetrahedron, [[1, 1, 1],
59	61	[-1, 0, 0],
60	62	[0, -1, 0],
61		[0, 0, -1]])])
	63	[0, 0, -1]]),
	64	(hexahedron, [[-1, 0, 0],
	65	[1, 0, 0],
	66	[0, -1, 0],
	67	[0, 1, 0],
	68	[0, 0, -1],
	69	[0, 0, 1]])])
62	70	def test_reference_normal(cell, normals):
63	71	facet_dim = cell.get_spatial_dimension() - 1
64	72	for facet_number in range(len(cell.get_topology()[facet_dim])):

+41

-1

test/unit/test_tensor_product.py less more

28	28	from FIAT.discontinuous_lagrange import DiscontinuousLagrange
29	29	from FIAT.nedelec import Nedelec
30	30	from FIAT.raviart_thomas import RaviartThomas
31		from FIAT.tensor_product import TensorProductElement
	31	from FIAT.tensor_product import TensorProductElement, FlattenedDimensions
32	32	from FIAT.hdivcurl import Hdiv, Hcurl
33	33	from FIAT.enriched import EnrichedElement
34	34

521	521	assert tab[dd][7][2][0] == 0.0
522	522
523	523
	524	def test_flattened_against_tpe_quad():
	525	T = UFCInterval()
	526	P1 = Lagrange(T, 1)
	527	tpe_quad = TensorProductElement(P1, P1)
	528	flattened_quad = FlattenedDimensions(tpe_quad)
	529	assert tpe_quad.value_shape() == ()
	530	tpe_tab = tpe_quad.tabulate(1, [(0.1, 0.2)])
	531	flattened_tab = flattened_quad.tabulate(1, [(0.1, 0.2)])
	532
	533	for da, db in [[(0,), (0,)], [(1,), (0,)], [(0,), (1,)]]:
	534	dc = da + db
	535	assert np.isclose(tpe_tab[dc][0][0], flattened_tab[dc][0][0])
	536	assert np.isclose(tpe_tab[dc][1][0], flattened_tab[dc][1][0])
	537	assert np.isclose(tpe_tab[dc][2][0], flattened_tab[dc][2][0])
	538	assert np.isclose(tpe_tab[dc][3][0], flattened_tab[dc][3][0])
	539
	540
	541	def test_flattened_against_tpe_hex():
	542	T = UFCInterval()
	543	P1 = Lagrange(T, 1)
	544	tpe_quad = TensorProductElement(P1, P1)
	545	tpe_hex = TensorProductElement(tpe_quad, P1)
	546	flattened_quad = FlattenedDimensions(tpe_quad)
	547	flattened_hex = FlattenedDimensions(TensorProductElement(flattened_quad, P1))
	548	assert tpe_quad.value_shape() == ()
	549	tpe_tab = tpe_hex.tabulate(1, [(0.1, 0.2, 0.3)])
	550	flattened_tab = flattened_hex.tabulate(1, [(0.1, 0.2, 0.3)])
	551
	552	for da, db, dc in [[(0,), (0,), (0,)], [(1,), (0,), (0,)], [(0,), (1,), (0,)], [(0,), (0,), (1,)]]:
	553	dd = da + db + dc
	554	assert np.isclose(tpe_tab[dd][0][0], flattened_tab[dd][0][0])
	555	assert np.isclose(tpe_tab[dd][1][0], flattened_tab[dd][1][0])
	556	assert np.isclose(tpe_tab[dd][2][0], flattened_tab[dd][2][0])
	557	assert np.isclose(tpe_tab[dd][3][0], flattened_tab[dd][3][0])
	558	assert np.isclose(tpe_tab[dd][4][0], flattened_tab[dd][4][0])
	559	assert np.isclose(tpe_tab[dd][5][0], flattened_tab[dd][5][0])
	560	assert np.isclose(tpe_tab[dd][6][0], flattened_tab[dd][6][0])
	561	assert np.isclose(tpe_tab[dd][7][0], flattened_tab[dd][7][0])
	562
	563
524	564	if __name__ == '__main__':
525	565	import os
526	566	pytest.main(os.path.abspath(__file__))