Commit 4e49719666854a07ea5630218eb33c42046ec8f5 - fiat

New upstream version 2019.2.0~git20200722.e16d9f3 Drew Parsons 3 years ago

14 changed file(s) with 654 addition(s) and 160 deletion(s). Raw diff Collapse all Expand all

-0

AUTHORS less more

72	72
73	73	Cyrus Cheng
74	74	email: cyruscycheng21@gmail.com
	75
	76	Reuben W. Hill
	77	email: reuben.hill10@imperial.ac.uk

-1

FIAT/P0.py less more

18	18	entity_ids = {}
19	19	nodes = []
20	20	vs = numpy.array(ref_el.get_vertices())
21		bary = tuple(numpy.average(vs, 0))
	21	if ref_el.get_dimension() == 0:
	22	bary = ()
	23	else:
	24	bary = tuple(numpy.average(vs, 0))
22	25
23	26	nodes = [functional.PointEvaluation(ref_el, bary)]
24	27	entity_ids = {}

+71

-27

FIAT/brezzi_douglas_marini.py less more

6	6
7	7	from FIAT import (finite_element, quadrature, functional, dual_set,
8	8	polynomial_set, nedelec)
	9	from FIAT.check_format_variant import check_format_variant
9	10
10	11
11	12	class BDMDualSet(dual_set.DualSet):
12		def __init__(self, ref_el, degree):
	13	def __init__(self, ref_el, degree, variant, quad_deg):
13	14
14	15	# Initialize containers for map: mesh_entity -> dof number and
15	16	# dual basis

19	20	sd = ref_el.get_spatial_dimension()
20	21	t = ref_el.get_topology()
21	22
22		# Define each functional for the dual set
23		# codimension 1 facets
24		for i in range(len(t[sd - 1])):
25		pts_cur = ref_el.make_points(sd - 1, i, sd + degree)
26		for j in range(len(pts_cur)):
27		pt_cur = pts_cur[j]
28		f = functional.PointScaledNormalEvaluation(ref_el, i, pt_cur)
29		nodes.append(f)
	23	if variant == "integral":
	24	facet = ref_el.get_facet_element()
	25	# Facet nodes are \int_F v\cdot n p ds where p \in P_{q-1}
	26	# degree is q - 1
	27	Q = quadrature.make_quadrature(facet, quad_deg)
	28	Pq = polynomial_set.ONPolynomialSet(facet, degree)
	29	Pq_at_qpts = Pq.tabulate(Q.get_points())[tuple([0]*(sd - 1))]
	30	for f in range(len(t[sd - 1])):
	31	for i in range(Pq_at_qpts.shape[0]):
	32	phi = Pq_at_qpts[i, :]
	33	nodes.append(functional.IntegralMomentOfScaledNormalEvaluation(ref_el, Q, phi, f))
30	34
31		# internal nodes
32		if degree > 1:
33		Q = quadrature.make_quadrature(ref_el, 2 * (degree + 1))
34		qpts = Q.get_points()
35		Nedel = nedelec.Nedelec(ref_el, degree - 1)
36		Nedfs = Nedel.get_nodal_basis()
37		zero_index = tuple([0 for i in range(sd)])
38		Ned_at_qpts = Nedfs.tabulate(qpts)[zero_index]
	35	# internal nodes
	36	if degree > 1:
	37	Q = quadrature.make_quadrature(ref_el, quad_deg)
	38	qpts = Q.get_points()
	39	Nedel = nedelec.Nedelec(ref_el, degree - 1, variant)
	40	Nedfs = Nedel.get_nodal_basis()
	41	zero_index = tuple([0 for i in range(sd)])
	42	Ned_at_qpts = Nedfs.tabulate(qpts)[zero_index]
39	43
40		for i in range(len(Ned_at_qpts)):
41		phi_cur = Ned_at_qpts[i, :]
42		l_cur = functional.FrobeniusIntegralMoment(ref_el, Q, phi_cur)
43		nodes.append(l_cur)
	44	for i in range(len(Ned_at_qpts)):
	45	phi_cur = Ned_at_qpts[i, :]
	46	l_cur = functional.FrobeniusIntegralMoment(ref_el, Q, phi_cur)
	47	nodes.append(l_cur)
	48
	49	elif variant == "point":
	50	# Define each functional for the dual set
	51	# codimension 1 facets
	52	for i in range(len(t[sd - 1])):
	53	pts_cur = ref_el.make_points(sd - 1, i, sd + degree)
	54	for j in range(len(pts_cur)):
	55	pt_cur = pts_cur[j]
	56	f = functional.PointScaledNormalEvaluation(ref_el, i, pt_cur)
	57	nodes.append(f)
	58
	59	# internal nodes
	60	if degree > 1:
	61	Q = quadrature.make_quadrature(ref_el, 2 * (degree + 1))
	62	qpts = Q.get_points()
	63	Nedel = nedelec.Nedelec(ref_el, degree - 1, variant)
	64	Nedfs = Nedel.get_nodal_basis()
	65	zero_index = tuple([0 for i in range(sd)])
	66	Ned_at_qpts = Nedfs.tabulate(qpts)[zero_index]
	67
	68	for i in range(len(Ned_at_qpts)):
	69	phi_cur = Ned_at_qpts[i, :]
	70	l_cur = functional.FrobeniusIntegralMoment(ref_el, Q, phi_cur)
	71	nodes.append(l_cur)
44	72
45	73	# sets vertices (and in 3d, edges) to have no nodes
46	74	for i in range(sd - 1):

70	98
71	99
72	100	class BrezziDouglasMarini(finite_element.CiarletElement):
73		"""The BDM element"""
	101	"""
	102	The BDM element
74	103
75		def __init__(self, ref_el, degree):
	104	:arg ref_el: The reference element.
	105	:arg k: The degree.
	106	:arg variant: optional variant specifying the types of nodes.
76	107
77		if degree < 1:
	108	variant can be chosen from ["point", "integral", "integral(quadrature_degree)"]
	109	"point" -> dofs are evaluated by point evaluation. Note that this variant has suboptimal
	110	convergence order in the H(div)-norm
	111	"integral" -> dofs are evaluated by quadrature rule. The quadrature degree is chosen to integrate
	112	polynomials of degree 5*k so that most expressions will be interpolated exactly. This is important
	113	when you want to have (nearly) divergence-preserving interpolation.
	114	"integral(quadrature_degree)" -> dofs are evaluated by quadrature rule of degree quadrature_degree
	115	"""
	116
	117	def __init__(self, ref_el, k, variant=None):
	118
	119	(variant, quad_deg) = check_format_variant(variant, k, "BDM")
	120
	121	if k < 1:
78	122	raise Exception("BDM_k elements only valid for k >= 1")
79	123
80	124	sd = ref_el.get_spatial_dimension()
81		poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd, ))
82		dual = BDMDualSet(ref_el, degree)
	125	poly_set = polynomial_set.ONPolynomialSet(ref_el, k, (sd, ))
	126	dual = BDMDualSet(ref_el, k, variant, quad_deg)
83	127	formdegree = sd - 1 # (n-1)-form
84		super(BrezziDouglasMarini, self).__init__(poly_set, dual, degree, formdegree,
	128	super(BrezziDouglasMarini, self).__init__(poly_set, dual, k, formdegree,
85	129	mapping="contravariant piola")

+25

-0

FIAT/check_format_variant.py less more

	0	import re
	1	import warnings
	2
	3
	4	def check_format_variant(variant, degree, element):
	5	if variant is None:
	6	variant = "point"
	7	warnings.simplefilter('always', DeprecationWarning)
	8	warnings.warn('Variant of ' + element + ' element will change from point evaluation to integral evaluation.'
	9	' You should project into variant="integral"', DeprecationWarning)
	10
	11	match = re.match(r"^integral(?:$(\d+)$)?$", variant)
	12	if match:
	13	variant = "integral"
	14	quad_degree, = match.groups()
	15	quad_degree = int(quad_degree) if quad_degree is not None else 5*(degree + 1)
	16	if quad_degree < degree + 1:
	17	raise ValueError("Warning, quadrature degree should be at least %s" % (degree + 1))
	18	elif variant == "point":
	19	quad_degree = None
	20	else:
	21	raise ValueError('Choose either variant="point" or variant="integral"'
	22	'or variant="integral(Quadrature degree)"')
	23
	24	return (variant, quad_degree)

+37

-3

FIAT/expansions.py less more

108	108	return xi1, xi2, xi3
109	109
110	110
	111	class PointExpansionSet(object):
	112	"""Evaluates the point basis on a point reference element."""
	113
	114	def __init__(self, ref_el):
	115	if ref_el.get_spatial_dimension() != 0:
	116	raise ValueError("Must have a point")
	117	self.ref_el = ref_el
	118	self.base_ref_el = reference_element.Point()
	119
	120	def get_num_members(self, n):
	121	return 1
	122
	123	def tabulate(self, n, pts):
	124	"""Returns a numpy array A[i,j] = phi_i(pts[j]) = 1.0."""
	125	assert n == 0
	126	return numpy.ones((1, len(pts)))
	127
	128	def tabulate_derivatives(self, n, pts):
	129	"""Returns a numpy array of size A where A[i,j] = phi_i(pts[j])
	130	but where each element is an empty tuple (). This maintains
	131	compatibility with the interfaces of the interval, triangle and
	132	tetrahedron expansions."""
	133	deriv_vals = numpy.empty_like(self.tabulate(n, pts), dtype=tuple)
	134	deriv_vals.fill(())
	135
	136	return deriv_vals
	137
	138
111	139	class LineExpansionSet(object):
112	140	"""Evaluates the Legendre basis on a line reference element."""
113	141

376	404	def get_expansion_set(ref_el):
377	405	"""Returns an ExpansionSet instance appopriate for the given
378	406	reference element."""
379		if ref_el.get_shape() == reference_element.LINE:
	407	if ref_el.get_shape() == reference_element.POINT:
	408	return PointExpansionSet(ref_el)
	409	elif ref_el.get_shape() == reference_element.LINE:
380	410	return LineExpansionSet(ref_el)
381	411	elif ref_el.get_shape() == reference_element.TRIANGLE:
382	412	return TriangleExpansionSet(ref_el)

389	419	def polynomial_dimension(ref_el, degree):
390	420	"""Returns the dimension of the space of polynomials of degree no
391	421	greater than degree on the reference element."""
392		if ref_el.get_shape() == reference_element.LINE:
	422	if ref_el.get_shape() == reference_element.POINT:
	423	if degree > 0:
	424	raise ValueError("Only degree zero polynomials supported on point elements.")
	425	return 1
	426	elif ref_el.get_shape() == reference_element.LINE:
393	427	return max(0, degree + 1)
394	428	elif ref_el.get_shape() == reference_element.TRIANGLE:
395	429	return max((degree + 1) * (degree + 2) // 2, 0)
396	430	elif ref_el.get_shape() == reference_element.TETRAHEDRON:
397	431	return max(0, (degree + 1) * (degree + 2) * (degree + 3) // 6)
398	432	else:
399		raise Exception("Unknown reference element type.")
	433	raise ValueError("Unknown reference element type.")

+99

-0

FIAT/functional.py less more

387	387	return "(u.t)(%s)" % (','.join(x),)
388	388
389	389
	390	class IntegralMomentOfEdgeTangentEvaluation(Functional):
	391	r"""
	392	\int_e v\cdot t p ds
	393
	394	p \in Polynomials
	395
	396	:arg ref_el: reference element for which e is a dim-1 entity
	397	:arg Q: quadrature rule on the face
	398	:arg P_at_qpts: polynomials evaluated at quad points
	399	:arg edge: which edge.
	400	"""
	401	def __init__(self, ref_el, Q, P_at_qpts, edge):
	402	t = ref_el.compute_edge_tangent(edge)
	403	sd = ref_el.get_spatial_dimension()
	404	transform = ref_el.get_entity_transform(1, edge)
	405	pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points()))
	406	weights = Q.get_weights()
	407	pt_dict = OrderedDict()
	408	for pt, wgt, phi in zip(pts, weights, P_at_qpts):
	409	pt_dict[pt] = [(wgtphit[i], (i, )) for i in range(sd)]
	410	super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfEdgeTangentEvaluation")
	411
	412
390	413	class PointFaceTangentEvaluation(Functional):
391	414	"""Implements the evaluation of a tangential component of a
392	415	vector at a point on a facet of codimension 1."""

405	428	return "(u.t%d)(%s)" % (self.tno, ','.join(x),)
406	429
407	430
	431	class IntegralMomentOfFaceTangentEvaluation(Functional):
	432	r"""
	433	\int_F v \times n \cdot p ds
	434
	435	p \in Polynomials
	436
	437	:arg ref_el: reference element for which F is a codim-1 entity
	438	:arg Q: quadrature rule on the face
	439	:arg P_at_qpts: polynomials evaluated at quad points
	440	:arg facet: which facet.
	441	"""
	442	def __init__(self, ref_el, Q, P_at_qpts, facet):
	443	P_at_qpts = [[P_at_qpts[0][i], P_at_qpts[1][i], P_at_qpts[2][i]]
	444	for i in range(P_at_qpts.shape[1])]
	445	n = ref_el.compute_scaled_normal(facet)
	446	sd = ref_el.get_spatial_dimension()
	447	transform = ref_el.get_entity_transform(sd-1, facet)
	448	pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points()))
	449	weights = Q.get_weights()
	450	pt_dict = OrderedDict()
	451	for pt, wgt, phi in zip(pts, weights, P_at_qpts):
	452	phixn = [phi[1]n[2] - phi[2]n[1],
	453	phi[2]n[0] - phi[0]n[2],
	454	phi[0]n[1] - phi[1]n[0]]
	455	pt_dict[pt] = [(wgt(-n[2]phixn[1]+n[1]*phixn[2]), (0, )),
	456	(wgt(n[2]phixn[0]-n[0]*phixn[2]), (1, )),
	457	(wgt(-n[1]phixn[0]+n[0]*phixn[1]), (2, ))]
	458	super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfFaceTangentEvaluation")
	459
	460
	461	class MonkIntegralMoment(Functional):
	462	r"""
	463	face nodes are \int_F v\cdot p dA where p \in P_{q-2}(f)^3 with p \cdot n = 0
	464	(cmp. Peter Monk - Finite Element Methods for Maxwell's equations p. 129)
	465	Note that we don't scale by the area of the facet
	466
	467	:arg ref_el: reference element for which F is a codim-1 entity
	468	:arg Q: quadrature rule on the face
	469	:arg P_at_qpts: polynomials evaluated at quad points
	470	:arg facet: which facet.
	471	"""
	472
	473	def __init__(self, ref_el, Q, P_at_qpts, facet):
	474	sd = ref_el.get_spatial_dimension()
	475	weights = Q.get_weights()
	476	pt_dict = OrderedDict()
	477	transform = ref_el.get_entity_transform(sd-1, facet)
	478	pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points()))
	479	for pt, wgt, phi in zip(pts, weights, P_at_qpts):
	480	pt_dict[pt] = [(wgt*phi[i], (i, )) for i in range(sd)]
	481	super().__init__(ref_el, (sd, ), pt_dict, {}, "MonkIntegralMoment")
	482
	483
408	484	class PointScaledNormalEvaluation(Functional):
409	485	"""Implements the evaluation of the normal component of a vector at a
410	486	point on a facet of codimension 1, where the normal is scaled by

421	497	def tostr(self):
422	498	x = list(map(str, list(self.pt_dict.keys())[0]))
423	499	return "(u.n)(%s)" % (','.join(x),)
	500
	501
	502	class IntegralMomentOfScaledNormalEvaluation(Functional):
	503	r"""
	504	\int_F v\cdot n p ds
	505
	506	p \in Polynomials
	507
	508	:arg ref_el: reference element for which F is a codim-1 entity
	509	:arg Q: quadrature rule on the face
	510	:arg P_at_qpts: polynomials evaluated at quad points
	511	:arg facet: which facet.
	512	"""
	513	def __init__(self, ref_el, Q, P_at_qpts, facet):
	514	n = ref_el.compute_scaled_normal(facet)
	515	sd = ref_el.get_spatial_dimension()
	516	transform = ref_el.get_entity_transform(sd - 1, facet)
	517	pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points()))
	518	weights = Q.get_weights()
	519	pt_dict = OrderedDict()
	520	for pt, wgt, phi in zip(pts, weights, P_at_qpts):
	521	pt_dict[pt] = [(wgtphin[i], (i, )) for i in range(sd)]
	522	super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfScaledNormalEvaluation")
424	523
425	524
426	525	class PointwiseInnerProductEvaluation(Functional):

+152

-63

FIAT/nedelec.py less more

8	8	finite_element, functional)
9	9	from itertools import chain
10	10	import numpy
	11	from FIAT.check_format_variant import check_format_variant
11	12
12	13
13	14	def NedelecSpace2D(ref_el, k):

142	143	class NedelecDual2D(dual_set.DualSet):
143	144	"""Dual basis for first-kind Nedelec in 2D."""
144	145
145		def __init__(self, ref_el, degree):
	146	def __init__(self, ref_el, degree, variant, quad_deg):
146	147	sd = ref_el.get_spatial_dimension()
147	148	if sd != 2:
148	149	raise Exception("Nedelec2D only works on triangles")

151	152
152	153	t = ref_el.get_topology()
153	154
154		num_edges = len(t[1])
155
156		# edge tangents
157		for i in range(num_edges):
158		pts_cur = ref_el.make_points(1, i, degree + 2)
159		for j in range(len(pts_cur)):
160		pt_cur = pts_cur[j]
161		f = functional.PointEdgeTangentEvaluation(ref_el, i, pt_cur)
162		nodes.append(f)
163
164		# internal moments
165		if degree > 0:
166		Q = quadrature.make_quadrature(ref_el, 2 * (degree + 1))
167		qpts = Q.get_points()
168		Pkm1 = polynomial_set.ONPolynomialSet(ref_el, degree - 1)
169		zero_index = tuple([0 for i in range(sd)])
170		Pkm1_at_qpts = Pkm1.tabulate(qpts)[zero_index]
171
172		for d in range(sd):
173		for i in range(Pkm1_at_qpts.shape[0]):
174		phi_cur = Pkm1_at_qpts[i, :]
175		l_cur = functional.IntegralMoment(ref_el, Q, phi_cur, (d,))
176		nodes.append(l_cur)
	155	if variant == "integral":
	156	# edge nodes are \int_F v\cdot t p ds where p \in P_{q-1}(edge)
	157	# degree is q - 1
	158	edge = ref_el.get_facet_element()
	159	Q = quadrature.make_quadrature(edge, quad_deg)
	160	Pq = polynomial_set.ONPolynomialSet(edge, degree)
	161	Pq_at_qpts = Pq.tabulate(Q.get_points())[tuple([0]*(sd - 1))]
	162	for e in range(len(t[sd - 1])):
	163	for i in range(Pq_at_qpts.shape[0]):
	164	phi = Pq_at_qpts[i, :]
	165	nodes.append(functional.IntegralMomentOfEdgeTangentEvaluation(ref_el, Q, phi, e))
	166
	167	# internal nodes. These are \int_T v \cdot p dx where p \in P_{q-2}^2
	168	if degree > 0:
	169	Q = quadrature.make_quadrature(ref_el, quad_deg)
	170	qpts = Q.get_points()
	171	Pkm1 = polynomial_set.ONPolynomialSet(ref_el, degree - 1)
	172	zero_index = tuple([0 for i in range(sd)])
	173	Pkm1_at_qpts = Pkm1.tabulate(qpts)[zero_index]
	174
	175	for d in range(sd):
	176	for i in range(Pkm1_at_qpts.shape[0]):
	177	phi_cur = Pkm1_at_qpts[i, :]
	178	l_cur = functional.IntegralMoment(ref_el, Q, phi_cur, (d,), (sd,))
	179	nodes.append(l_cur)
	180
	181	elif variant == "point":
	182	num_edges = len(t[1])
	183
	184	# edge tangents
	185	for i in range(num_edges):
	186	pts_cur = ref_el.make_points(1, i, degree + 2)
	187	for j in range(len(pts_cur)):
	188	pt_cur = pts_cur[j]
	189	f = functional.PointEdgeTangentEvaluation(ref_el, i, pt_cur)
	190	nodes.append(f)
	191
	192	# internal moments
	193	if degree > 0:
	194	Q = quadrature.make_quadrature(ref_el, 2 * (degree + 1))
	195	qpts = Q.get_points()
	196	Pkm1 = polynomial_set.ONPolynomialSet(ref_el, degree - 1)
	197	zero_index = tuple([0 for i in range(sd)])
	198	Pkm1_at_qpts = Pkm1.tabulate(qpts)[zero_index]
	199
	200	for d in range(sd):
	201	for i in range(Pkm1_at_qpts.shape[0]):
	202	phi_cur = Pkm1_at_qpts[i, :]
	203	l_cur = functional.IntegralMoment(ref_el, Q, phi_cur, (d,))
	204	nodes.append(l_cur)
177	205
178	206	entity_ids = {}
179	207

203	231	class NedelecDual3D(dual_set.DualSet):
204	232	"""Dual basis for first-kind Nedelec in 3D."""
205	233
206		def __init__(self, ref_el, degree):
	234	def __init__(self, ref_el, degree, variant, quad_deg):
207	235	sd = ref_el.get_spatial_dimension()
208	236	if sd != 3:
209	237	raise Exception("NedelecDual3D only works on tetrahedra")

212	240
213	241	t = ref_el.get_topology()
214	242
215		# how many edges
216		num_edges = len(t[1])
217
218		for i in range(num_edges):
219		# points to specify P_k on each edge
220		pts_cur = ref_el.make_points(1, i, degree + 2)
221		for j in range(len(pts_cur)):
222		pt_cur = pts_cur[j]
223		f = functional.PointEdgeTangentEvaluation(ref_el, i, pt_cur)
224		nodes.append(f)
225
226		if degree > 0: # face tangents
227		num_faces = len(t[2])
228		for i in range(num_faces): # loop over faces
229		pts_cur = ref_el.make_points(2, i, degree + 2)
230		for j in range(len(pts_cur)): # loop over points
	243	if variant == "integral":
	244	# edge nodes are \int_F v\cdot t p ds where p \in P_{q-1}(edge)
	245	# degree is q - 1
	246	edge = ref_el.get_facet_element().get_facet_element()
	247	Q = quadrature.make_quadrature(edge, quad_deg)
	248	Pq = polynomial_set.ONPolynomialSet(edge, degree)
	249	Pq_at_qpts = Pq.tabulate(Q.get_points())[tuple([0]*(1))]
	250	for e in range(len(t[1])):
	251	for i in range(Pq_at_qpts.shape[0]):
	252	phi = Pq_at_qpts[i, :]
	253	nodes.append(functional.IntegralMomentOfEdgeTangentEvaluation(ref_el, Q, phi, e))
	254
	255	# face nodes are \int_F v\cdot p dA where p \in P_{q-2}(f)^3 with p \cdot n = 0 (cmp. Monk)
	256	# these are equivalent to dofs from Fenics book defined by
	257	# \int_F v\times n \cdot p ds where p \in P_{q-2}(f)^2
	258	if degree > 0:
	259	facet = ref_el.get_facet_element()
	260	Q = quadrature.make_quadrature(facet, quad_deg)
	261	Pq = polynomial_set.ONPolynomialSet(facet, degree-1, (sd,))
	262	Pq_at_qpts = Pq.tabulate(Q.get_points())[(0, 0)]
	263
	264	for f in range(len(t[2])):
	265	# R is used to map [1,0,0] to tangent1 and [0,1,0] to tangent2
	266	R = ref_el.compute_face_tangents(f)
	267
	268	# Skip last functionals because we only want p with p \cdot n = 0
	269	for i in range(2 * Pq.get_num_members() // 3):
	270	phi = Pq_at_qpts[i, ...]
	271	phi = numpy.matmul(phi[:-1, ...].T, R)
	272	nodes.append(functional.MonkIntegralMoment(ref_el, Q, phi, f))
	273
	274	# internal nodes. These are \int_T v \cdot p dx where p \in P_{q-3}^3(T)
	275	if degree > 1:
	276	Q = quadrature.make_quadrature(ref_el, quad_deg)
	277	qpts = Q.get_points()
	278	Pkm2 = polynomial_set.ONPolynomialSet(ref_el, degree - 2)
	279	zero_index = tuple([0 for i in range(sd)])
	280	Pkm2_at_qpts = Pkm2.tabulate(qpts)[zero_index]
	281
	282	for d in range(sd):
	283	for i in range(Pkm2_at_qpts.shape[0]):
	284	phi_cur = Pkm2_at_qpts[i, :]
	285	l_cur = functional.IntegralMoment(ref_el, Q, phi_cur, (d,), (sd,))
	286	nodes.append(l_cur)
	287
	288	elif variant == "point":
	289	num_edges = len(t[1])
	290
	291	for i in range(num_edges):
	292	# points to specify P_k on each edge
	293	pts_cur = ref_el.make_points(1, i, degree + 2)
	294	for j in range(len(pts_cur)):
231	295	pt_cur = pts_cur[j]
232		for k in range(2): # loop over tangents
233		f = functional.PointFaceTangentEvaluation(ref_el, i, k, pt_cur)
	296	f = functional.PointEdgeTangentEvaluation(ref_el, i, pt_cur)
	297	nodes.append(f)
	298
	299	if degree > 0: # face tangents
	300	num_faces = len(t[2])
	301	for i in range(num_faces): # loop over faces
	302	pts_cur = ref_el.make_points(2, i, degree + 2)
	303	for j in range(len(pts_cur)): # loop over points
	304	pt_cur = pts_cur[j]
	305	for k in range(2): # loop over tangents
	306	f = functional.PointFaceTangentEvaluation(ref_el, i, k, pt_cur)
	307	nodes.append(f)
	308
	309	if degree > 1: # internal moments
	310	Q = quadrature.make_quadrature(ref_el, 2 * (degree + 1))
	311	qpts = Q.get_points()
	312	Pkm2 = polynomial_set.ONPolynomialSet(ref_el, degree - 2)
	313	zero_index = tuple([0 for i in range(sd)])
	314	Pkm2_at_qpts = Pkm2.tabulate(qpts)[zero_index]
	315
	316	for d in range(sd):
	317	for i in range(Pkm2_at_qpts.shape[0]):
	318	phi_cur = Pkm2_at_qpts[i, :]
	319	f = functional.IntegralMoment(ref_el, Q, phi_cur, (d,))
234	320	nodes.append(f)
235
236		if degree > 1: # internal moments
237		Q = quadrature.make_quadrature(ref_el, 2 * (degree + 1))
238		qpts = Q.get_points()
239		Pkm2 = polynomial_set.ONPolynomialSet(ref_el, degree - 2)
240		zero_index = tuple([0 for i in range(sd)])
241		Pkm2_at_qpts = Pkm2.tabulate(qpts)[zero_index]
242
243		for d in range(sd):
244		for i in range(Pkm2_at_qpts.shape[0]):
245		phi_cur = Pkm2_at_qpts[i, :]
246		f = functional.IntegralMoment(ref_el, Q, phi_cur, (d,))
247		nodes.append(f)
248	321
249	322	entity_ids = {}
250	323	# set to empty

276	349
277	350
278	351	class Nedelec(finite_element.CiarletElement):
279		"""Nedelec finite element"""
280
281		def __init__(self, ref_el, q):
282
283		degree = q - 1
	352	"""
	353	Nedelec finite element
	354
	355	:arg ref_el: The reference element.
	356	:arg k: The degree.
	357	:arg variant: optional variant specifying the types of nodes.
	358
	359	variant can be chosen from ["point", "integral", "integral(quadrature_degree)"]
	360	"point" -> dofs are evaluated by point evaluation. Note that this variant has suboptimal
	361	convergence order in the H(curl)-norm
	362	"integral" -> dofs are evaluated by quadrature rule. The quadrature degree is chosen to integrate
	363	polynomials of degree 5*k so that most expressions will be interpolated exactly. This is important
	364	when you want to have (nearly) curl-preserving interpolation.
	365	"integral(quadrature_degree)" -> dofs are evaluated by quadrature rule of degree quadrature_degree
	366	"""
	367
	368	def __init__(self, ref_el, k, variant=None):
	369
	370	degree = k - 1
	371
	372	(variant, quad_deg) = check_format_variant(variant, degree, "Nedelec")
284	373
285	374	if ref_el.get_spatial_dimension() == 3:
286	375	poly_set = NedelecSpace3D(ref_el, degree)
287		dual = NedelecDual3D(ref_el, degree)
	376	dual = NedelecDual3D(ref_el, degree, variant, quad_deg)
288	377	elif ref_el.get_spatial_dimension() == 2:
289	378	poly_set = NedelecSpace2D(ref_el, degree)
290		dual = NedelecDual2D(ref_el, degree)
	379	dual = NedelecDual2D(ref_el, degree, variant, quad_deg)
291	380	else:
292	381	raise Exception("Not implemented")
293	382	formdegree = 1 # 1-form

+61

-31

FIAT/nedelec_second_kind.py less more

14	14	from FIAT.raviart_thomas import RaviartThomas
15	15	from FIAT.quadrature import make_quadrature, UFCTetrahedronFaceQuadratureRule
16	16	from FIAT.reference_element import UFCTetrahedron
	17	from FIAT.check_format_variant import check_format_variant
	18
	19	from FIAT import polynomial_set, quadrature, functional
17	20
18	21
19	22	class NedelecSecondKindDual(DualSet):

45	48	these elements coincide with the CG_k elements.)
46	49	"""
47	50
48		def __init__(self, cell, degree):
	51	def __init__(self, cell, degree, variant, quad_deg):
49	52
50	53	# Define degrees of freedom
51		(dofs, ids) = self.generate_degrees_of_freedom(cell, degree)
52
	54	(dofs, ids) = self.generate_degrees_of_freedom(cell, degree, variant, quad_deg)
53	55	# Call init of super-class
54	56	super(NedelecSecondKindDual, self).__init__(dofs, cell, ids)
55	57
56		def generate_degrees_of_freedom(self, cell, degree):
	58	def generate_degrees_of_freedom(self, cell, degree, variant, quad_deg):
57	59	"Generate dofs and geometry-to-dof maps (ids)."
58	60
59	61	dofs = []

67	69	ids[0] = dict(list(zip(list(range(d + 1)), ([] for i in range(d + 1)))))
68	70
69	71	# (d+1) degrees of freedom per entity of codimension 1 (edges)
70		(edge_dofs, edge_ids) = self._generate_edge_dofs(cell, degree, 0)
	72	(edge_dofs, edge_ids) = self._generate_edge_dofs(cell, degree, 0, variant, quad_deg)
71	73	dofs.extend(edge_dofs)
72	74	ids[1] = edge_ids
73	75
74	76	# Include face degrees of freedom if 3D
75	77	if d == 3:
76	78	(face_dofs, face_ids) = self._generate_face_dofs(cell, degree,
77		len(dofs))
	79	len(dofs), variant, quad_deg)
78	80	dofs.extend(face_dofs)
79	81	ids[2] = face_ids
80	82
81	83	# Varying degrees of freedom (possibly zero) per cell
82		(cell_dofs, cell_ids) = self._generate_cell_dofs(cell, degree, len(dofs))
	84	(cell_dofs, cell_ids) = self._generate_cell_dofs(cell, degree, len(dofs), variant, quad_deg)
83	85	dofs.extend(cell_dofs)
84	86	ids[d] = cell_ids
85	87
86	88	return (dofs, ids)
87	89
88		def _generate_edge_dofs(self, cell, degree, offset):
	90	def _generate_edge_dofs(self, cell, degree, offset, variant, quad_deg):
89	91	"""Generate degrees of freedoms (dofs) for entities of
90	92	codimension 1 (edges)."""
91	93

93	95	# freedom per entity of codimension 1 (edges)
94	96	dofs = []
95	97	ids = {}
96		for edge in range(len(cell.get_topology()[1])):
97
98		# Create points for evaluation of tangential components
99		points = cell.make_points(1, edge, degree + 2)
100
101		# A tangential component evaluation for each point
102		dofs += [Tangent(cell, edge, point) for point in points]
103
104		# Associate these dofs with this edge
105		i = len(points) * edge
106		ids[edge] = list(range(offset + i, offset + i + len(points)))
107
108		return (dofs, ids)
109
110		def _generate_face_dofs(self, cell, degree, offset):
	98
	99	if variant == "integral":
	100	edge = cell.construct_subelement(1)
	101	Q = quadrature.make_quadrature(edge, quad_deg)
	102	Pq = polynomial_set.ONPolynomialSet(edge, degree)
	103	Pq_at_qpts = Pq.tabulate(Q.get_points())[tuple([0]*(1))]
	104	for e in range(len(cell.get_topology()[1])):
	105	for i in range(Pq_at_qpts.shape[0]):
	106	phi = Pq_at_qpts[i, :]
	107	dofs.append(functional.IntegralMomentOfEdgeTangentEvaluation(cell, Q, phi, e))
	108	jj = Pq_at_qpts.shape[0] * e
	109	ids[e] = list(range(offset + jj, offset + jj + Pq_at_qpts.shape[0]))
	110
	111	elif variant == "point":
	112	for edge in range(len(cell.get_topology()[1])):
	113
	114	# Create points for evaluation of tangential components
	115	points = cell.make_points(1, edge, degree + 2)
	116
	117	# A tangential component evaluation for each point
	118	dofs += [Tangent(cell, edge, point) for point in points]
	119
	120	# Associate these dofs with this edge
	121	i = len(points) * edge
	122	ids[edge] = list(range(offset + i, offset + i + len(points)))
	123
	124	return (dofs, ids)
	125
	126	def _generate_face_dofs(self, cell, degree, offset, variant, quad_deg):
111	127	"""Generate degrees of freedoms (dofs) for faces."""
112	128
113	129	# Initialize empty dofs and identifiers (ids)

133	149	# Construct Raviart-Thomas of (degree - 1) on the
134	150	# reference face
135	151	reference_face = Q_face.reference_rule().ref_el
136		RT = RaviartThomas(reference_face, degree - 1)
	152	RT = RaviartThomas(reference_face, degree - 1, variant)
137	153	num_rts = RT.space_dimension()
138	154
139	155	# Evaluate RT basis functions at reference quadrature

167	183
168	184	return (dofs, ids)
169	185
170		def _generate_cell_dofs(self, cell, degree, offset):
	186	def _generate_cell_dofs(self, cell, degree, offset, variant, quad_deg):
171	187	"""Generate degrees of freedoms (dofs) for entities of
172	188	codimension d (cells)."""
173	189

181	197	qs = Q.get_points()
182	198
183	199	# Create Raviart-Thomas nodal basis
184		RT = RaviartThomas(cell, degree + 1 - d)
	200	RT = RaviartThomas(cell, degree + 1 - d, variant)
185	201	phi = RT.get_nodal_basis()
186	202
187	203	# Evaluate Raviart-Thomas basis at quadrature points

202	218	tetrahedra: the polynomial space described by the full polynomials
203	219	of degree k, with a suitable set of degrees of freedom to ensure
204	220	H(curl) conformity.
	221
	222	:arg ref_el: The reference element.
	223	:arg k: The degree.
	224	:arg variant: optional variant specifying the types of nodes.
	225
	226	variant can be chosen from ["point", "integral", "integral(quadrature_degree)"]
	227	"point" -> dofs are evaluated by point evaluation. Note that this variant has suboptimal
	228	convergence order in the H(curl)-norm
	229	"integral" -> dofs are evaluated by quadrature rule. The quadrature degree is chosen to integrate
	230	polynomials of degree 5*k so that most expressions will be interpolated exactly. This is important
	231	when you want to have (nearly) curl-preserving interpolation.
	232	"integral(quadrature_degree)" -> dofs are evaluated by quadrature rule of degree quadrature_degree
205	233	"""
206	234
207		def __init__(self, cell, degree):
	235	def __init__(self, cell, k, variant=None):
	236
	237	(variant, quad_deg) = check_format_variant(variant, k, "Nedelec Second Kind")
208	238
209	239	# Check degree
210		assert degree >= 1, "Second kind Nedelecs start at 1!"
	240	assert k >= 1, "Second kind Nedelecs start at 1!"
211	241
212	242	# Get dimension
213	243	d = cell.get_spatial_dimension()
214	244
215	245	# Construct polynomial basis for d-vector fields
216		Ps = ONPolynomialSet(cell, degree, (d, ))
	246	Ps = ONPolynomialSet(cell, k, (d, ))
217	247
218	248	# Construct dual space
219		Ls = NedelecSecondKindDual(cell, degree)
	249	Ls = NedelecSecondKindDual(cell, k, variant, quad_deg)
220	250
221	251	# Set form degree
222	252	formdegree = 1 # 1-form

225	255	mapping = "covariant piola"
226	256
227	257	# Call init of super-class
228		super(NedelecSecondKind, self).__init__(Ps, Ls, degree, formdegree, mapping=mapping)
	258	super(NedelecSecondKind, self).__init__(Ps, Ls, k, formdegree, mapping=mapping)

-1

FIAT/polynomial_set.py less more

74	74	for i in range(jet_order + 1):
75	75	alphas = mis(self.ref_el.get_spatial_dimension(), i)
76	76	for alpha in alphas:
77		D = form_matrix_product(self.dmats, alpha)
	77	if len(self.dmats) > 0:
	78	D = form_matrix_product(self.dmats, alpha)
	79	else:
	80	# special for vertex without defined point location
	81	assert pts == [()]
	82	D = numpy.eye(1)
78	83	result[alpha] = numpy.dot(self.coeffs,
79	84	numpy.dot(numpy.transpose(D),
80	85	base_vals))

+63

-30

FIAT/raviart_thomas.py less more

8	8	finite_element, functional)
9	9	import numpy
10	10	from itertools import chain
	11	from FIAT.check_format_variant import check_format_variant
11	12
12	13
13	14	def RTSpace(ref_el, deg):

64	65	evaluation of normals on facets of codimension 1 and internal
65	66	moments against polynomials"""
66	67
67		def __init__(self, ref_el, degree):
	68	def __init__(self, ref_el, degree, variant, quad_deg):
68	69	entity_ids = {}
69	70	nodes = []
70	71
71	72	sd = ref_el.get_spatial_dimension()
72	73	t = ref_el.get_topology()
73	74
74		# codimension 1 facets
75		for i in range(len(t[sd - 1])):
76		pts_cur = ref_el.make_points(sd - 1, i, sd + degree)
77		for j in range(len(pts_cur)):
78		pt_cur = pts_cur[j]
79		f = functional.PointScaledNormalEvaluation(ref_el, i, pt_cur)
80		nodes.append(f)
	75	if variant == "integral":
	76	facet = ref_el.get_facet_element()
	77	# Facet nodes are \int_F v\cdot n p ds where p \in P_{q-1}
	78	# degree is q - 1
	79	Q = quadrature.make_quadrature(facet, quad_deg)
	80	Pq = polynomial_set.ONPolynomialSet(facet, degree)
	81	Pq_at_qpts = Pq.tabulate(Q.get_points())[tuple([0]*(sd - 1))]
	82	for f in range(len(t[sd - 1])):
	83	for i in range(Pq_at_qpts.shape[0]):
	84	phi = Pq_at_qpts[i, :]
	85	nodes.append(functional.IntegralMomentOfScaledNormalEvaluation(ref_el, Q, phi, f))
81	86
82		# internal nodes. Let's just use points at a lattice
83		if degree > 0:
84		cpe = functional.ComponentPointEvaluation
85		pts = ref_el.make_points(sd, 0, degree + sd)
86		for d in range(sd):
87		for i in range(len(pts)):
88		l_cur = cpe(ref_el, d, (sd,), pts[i])
89		nodes.append(l_cur)
	87	# internal nodes. These are \int_T v \cdot p dx where p \in P_{q-2}^d
	88	if degree > 0:
	89	Q = quadrature.make_quadrature(ref_el, quad_deg)
	90	qpts = Q.get_points()
	91	Pkm1 = polynomial_set.ONPolynomialSet(ref_el, degree - 1)
	92	zero_index = tuple([0 for i in range(sd)])
	93	Pkm1_at_qpts = Pkm1.tabulate(qpts)[zero_index]
90	94
91		# Q = quadrature.make_quadrature(ref_el, 2 * ( degree + 1 ))
92		# qpts = Q.get_points()
93		# Pkm1 = polynomial_set.ONPolynomialSet(ref_el, degree - 1)
94		# zero_index = tuple([0 for i in range(sd)])
95		# Pkm1_at_qpts = Pkm1.tabulate(qpts)[zero_index]
	95	for d in range(sd):
	96	for i in range(Pkm1_at_qpts.shape[0]):
	97	phi_cur = Pkm1_at_qpts[i, :]
	98	l_cur = functional.IntegralMoment(ref_el, Q, phi_cur, (d,), (sd,))
	99	nodes.append(l_cur)
96	100
97		# for d in range(sd):
98		# for i in range(Pkm1_at_qpts.shape[0]):
99		# phi_cur = Pkm1_at_qpts[i, :]
100		# l_cur = functional.IntegralMoment(ref_el, Q, phi_cur, (d,), (sd,))
101		# nodes.append(l_cur)
	101	elif variant == "point":
	102	# codimension 1 facets
	103	for i in range(len(t[sd - 1])):
	104	pts_cur = ref_el.make_points(sd - 1, i, sd + degree)
	105	for j in range(len(pts_cur)):
	106	pt_cur = pts_cur[j]
	107	f = functional.PointScaledNormalEvaluation(ref_el, i, pt_cur)
	108	nodes.append(f)
	109
	110	# internal nodes. Let's just use points at a lattice
	111	if degree > 0:
	112	cpe = functional.ComponentPointEvaluation
	113	pts = ref_el.make_points(sd, 0, degree + sd)
	114	for d in range(sd):
	115	for i in range(len(pts)):
	116	l_cur = cpe(ref_el, d, (sd,), pts[i])
	117	nodes.append(l_cur)
102	118
103	119	# sets vertices (and in 3d, edges) to have no nodes
104	120	for i in range(sd - 1):

127	143
128	144
129	145	class RaviartThomas(finite_element.CiarletElement):
130		"""The Raviart-Thomas finite element"""
	146	"""
	147	The Raviart Thomas element
131	148
132		def __init__(self, ref_el, q):
	149	:arg ref_el: The reference element.
	150	:arg k: The degree.
	151	:arg variant: optional variant specifying the types of nodes.
133	152
134		degree = q - 1
	153	variant can be chosen from ["point", "integral", "integral(quadrature_degree)"]
	154	"point" -> dofs are evaluated by point evaluation. Note that this variant has suboptimal
	155	convergence order in the H(div)-norm
	156	"integral" -> dofs are evaluated by quadrature rule. The quadrature degree is chosen to integrate
	157	polynomials of degree 5*k so that most expressions will be interpolated exactly. This is important
	158	when you want to have (nearly) divergence-preserving interpolation.
	159	"integral(quadrature_degree)" -> dofs are evaluated by quadrature rule of degree quadrature_degree
	160	"""
	161
	162	def __init__(self, ref_el, k, variant=None):
	163
	164	degree = k - 1
	165
	166	(variant, quad_deg) = check_format_variant(variant, degree, "Raviart Thomas")
	167
135	168	poly_set = RTSpace(ref_el, degree)
136		dual = RTDualSet(ref_el, degree)
	169	dual = RTDualSet(ref_el, degree, variant, quad_deg)
137	170	formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form
138	171	super(RaviartThomas, self).__init__(poly_set, dual, degree, formdegree,
139	172	mapping="contravariant piola")

+12

-0

FIAT/reference_element.py less more

5	5	#
6	6	# Modified by David A. Ham (david.ham@imperial.ac.uk), 2014
7	7	# Modified by Lizao Li (lzlarryli@gmail.com), 2016
	8
8	9	"""
9	10	Abstract class and particular implementations of finite element
10	11	reference simplex geometry/topology.

480	481	verts = ((),)
481	482	topology = {0: {0: (0,)}}
482	483	super(Point, self).__init__(POINT, verts, topology)
	484
	485	def construct_subelement(self, dimension):
	486	"""Constructs the reference element of a cell subentity
	487	specified by subelement dimension.
	488
	489	:arg dimension: subentity dimension (integer). Must be zero.
	490	"""
	491	assert dimension == 0
	492	return self
483	493
484	494
485	495	class DefaultLine(Simplex):

967	977	return UFCQuadrilateral()
968	978	elif celltype == "hexahedron":
969	979	return UFCHexahedron()
	980	elif celltype == "vertex":
	981	return ufc_simplex(0)
970	982	elif celltype == "interval":
971	983	return ufc_simplex(1)
972	984	elif celltype == "triangle":

-2

FIAT/serendipity.py less more

149	149	raise NotImplementedError('no tabulate method for serendipity elements of dimension 1 or less.')
150	150	if dim >= 4:
151	151	raise NotImplementedError('tabulate does not support higher dimensions than 3.')
	152	points = np.asarray(points)
	153	npoints, pointdim = points.shape
152	154	for o in range(order + 1):
153	155	alphas = mis(dim, o)
154	156	for alpha in alphas:

159	161	callable = lambdify(variables[:dim], polynomials, modules="numpy", dummify=True)
160	162	self.basis[alpha] = polynomials
161	163	self.basis_callable[alpha] = callable
162		points = np.asarray(points)
163		T = np.asarray(callable(*(points[:, i] for i in range(points.shape[1]))))
	164	tabulation = callable(*(points[:, i] for i in range(pointdim)))
	165	T = np.asarray([np.broadcast_to(tab, (npoints, ))
	166	for tab in tabulation])
164	167	phivals[alpha] = T
165	168	return phivals
166	169

+84

-2

test/unit/test_fiat.py less more

19	19	import pytest
20	20
21	21	from FIAT.reference_element import LINE, ReferenceElement
22		from FIAT.reference_element import UFCInterval, UFCTriangle, UFCTetrahedron
	22	from FIAT.reference_element import Point, UFCInterval, UFCTriangle, UFCTetrahedron
23	23	from FIAT.lagrange import Lagrange
24	24	from FIAT.discontinuous_lagrange import DiscontinuousLagrange # noqa: F401
25	25	from FIAT.discontinuous_taylor import DiscontinuousTaylor # noqa: F401

48	48	from FIAT.enriched import EnrichedElement # noqa: F401
49	49	from FIAT.nodal_enriched import NodalEnrichedElement
50	50
51
	51	P = Point()
52	52	I = UFCInterval() # noqa: E741
53	53	T = UFCTriangle()
54	54	S = UFCTetrahedron()

109	109	"P0(I)",
110	110	"P0(T)",
111	111	"P0(S)",
	112	"DiscontinuousLagrange(P, 0)",
112	113	"DiscontinuousLagrange(I, 0)",
113	114	"DiscontinuousLagrange(I, 1)",
114	115	"DiscontinuousLagrange(I, 2)",

136	137	"RaviartThomas(S, 1)",
137	138	"RaviartThomas(S, 2)",
138	139	"RaviartThomas(S, 3)",
	140	'RaviartThomas(T, 1, variant="integral")',
	141	'RaviartThomas(T, 2, variant="integral")',
	142	'RaviartThomas(T, 3, variant="integral")',
	143	'RaviartThomas(S, 1, variant="integral")',
	144	'RaviartThomas(S, 2, variant="integral")',
	145	'RaviartThomas(S, 3, variant="integral")',
	146	'RaviartThomas(T, 1, variant="integral(2)")',
	147	'RaviartThomas(T, 2, variant="integral(3)")',
	148	'RaviartThomas(T, 3, variant="integral(4)")',
	149	'RaviartThomas(S, 1, variant="integral(2)")',
	150	'RaviartThomas(S, 2, variant="integral(3)")',
	151	'RaviartThomas(S, 3, variant="integral(4)")',
	152	'RaviartThomas(T, 1, variant="point")',
	153	'RaviartThomas(T, 2, variant="point")',
	154	'RaviartThomas(T, 3, variant="point")',
	155	'RaviartThomas(S, 1, variant="point")',
	156	'RaviartThomas(S, 2, variant="point")',
	157	'RaviartThomas(S, 3, variant="point")',
139	158	"DiscontinuousRaviartThomas(T, 1)",
140	159	"DiscontinuousRaviartThomas(T, 2)",
141	160	"DiscontinuousRaviartThomas(T, 3)",

148	167	"BrezziDouglasMarini(S, 1)",
149	168	"BrezziDouglasMarini(S, 2)",
150	169	"BrezziDouglasMarini(S, 3)",
	170	'BrezziDouglasMarini(T, 1, variant="integral")',
	171	'BrezziDouglasMarini(T, 2, variant="integral")',
	172	'BrezziDouglasMarini(T, 3, variant="integral")',
	173	'BrezziDouglasMarini(S, 1, variant="integral")',
	174	'BrezziDouglasMarini(S, 2, variant="integral")',
	175	'BrezziDouglasMarini(S, 3, variant="integral")',
	176	'BrezziDouglasMarini(T, 1, variant="integral(2)")',
	177	'BrezziDouglasMarini(T, 2, variant="integral(3)")',
	178	'BrezziDouglasMarini(T, 3, variant="integral(4)")',
	179	'BrezziDouglasMarini(S, 1, variant="integral(2)")',
	180	'BrezziDouglasMarini(S, 2, variant="integral(3)")',
	181	'BrezziDouglasMarini(S, 3, variant="integral(4)")',
	182	'BrezziDouglasMarini(T, 1, variant="point")',
	183	'BrezziDouglasMarini(T, 2, variant="point")',
	184	'BrezziDouglasMarini(T, 3, variant="point")',
	185	'BrezziDouglasMarini(S, 1, variant="point")',
	186	'BrezziDouglasMarini(S, 2, variant="point")',
	187	'BrezziDouglasMarini(S, 3, variant="point")',
151	188	"Nedelec(T, 1)",
152	189	"Nedelec(T, 2)",
153	190	"Nedelec(T, 3)",
154	191	"Nedelec(S, 1)",
155	192	"Nedelec(S, 2)",
156	193	"Nedelec(S, 3)",
	194	'Nedelec(T, 1, variant="integral")',
	195	'Nedelec(T, 2, variant="integral")',
	196	'Nedelec(T, 3, variant="integral")',
	197	'Nedelec(S, 1, variant="integral")',
	198	'Nedelec(S, 2, variant="integral")',
	199	'Nedelec(S, 3, variant="integral")',
	200	'Nedelec(T, 1, variant="integral(2)")',
	201	'Nedelec(T, 2, variant="integral(3)")',
	202	'Nedelec(T, 3, variant="integral(4)")',
	203	'Nedelec(S, 1, variant="integral(2)")',
	204	'Nedelec(S, 2, variant="integral(3)")',
	205	'Nedelec(S, 3, variant="integral(4)")',
	206	'Nedelec(T, 1, variant="point")',
	207	'Nedelec(T, 2, variant="point")',
	208	'Nedelec(T, 3, variant="point")',
	209	'Nedelec(S, 1, variant="point")',
	210	'Nedelec(S, 2, variant="point")',
	211	'Nedelec(S, 3, variant="point")',
157	212	"NedelecSecondKind(T, 1)",
158	213	"NedelecSecondKind(T, 2)",
159	214	"NedelecSecondKind(T, 3)",
160	215	"NedelecSecondKind(S, 1)",
161	216	"NedelecSecondKind(S, 2)",
162	217	"NedelecSecondKind(S, 3)",
	218	'NedelecSecondKind(T, 1, variant="integral")',
	219	'NedelecSecondKind(T, 2, variant="integral")',
	220	'NedelecSecondKind(T, 3, variant="integral")',
	221	'NedelecSecondKind(S, 1, variant="integral")',
	222	'NedelecSecondKind(S, 2, variant="integral")',
	223	'NedelecSecondKind(S, 3, variant="integral")',
	224	'NedelecSecondKind(T, 1, variant="integral(2)")',
	225	'NedelecSecondKind(T, 2, variant="integral(3)")',
	226	'NedelecSecondKind(T, 3, variant="integral(4)")',
	227	'NedelecSecondKind(S, 1, variant="integral(2)")',
	228	'NedelecSecondKind(S, 2, variant="integral(3)")',
	229	'NedelecSecondKind(S, 3, variant="integral(4)")',
	230	'NedelecSecondKind(T, 1, variant="point")',
	231	'NedelecSecondKind(T, 2, variant="point")',
	232	'NedelecSecondKind(T, 3, variant="point")',
	233	'NedelecSecondKind(S, 1, variant="point")',
	234	'NedelecSecondKind(S, 2, variant="point")',
	235	'NedelecSecondKind(S, 3, variant="point")',
163	236	"Regge(T, 0)",
164	237	"Regge(T, 1)",
165	238	"Regge(T, 2)",

422	495	assert np.isclose(basis[i], 1.0 if i == j else 0.0)
423	496
424	497
	498	@pytest.mark.parametrize('element', [
	499	'Nedelec(S, 3, variant="integral(2)")',
	500	'NedelecSecondKind(S, 3, variant="integral(3)")'
	501	])
	502	def test_error_quadrature_degree(element):
	503	with pytest.raises(ValueError):
	504	eval(element)
	505
	506
425	507	if __name__ == '__main__':
426	508	import os
427	509	pytest.main(os.path.abspath(__file__))

+32

-0

test/unit/test_serendipity.py less more

	0	from FIAT.reference_element import UFCQuadrilateral
	1	from FIAT import Serendipity
	2	import numpy as np
	3	import sympy
	4
	5
	6	def test_serendipity_derivatives():
	7	cell = UFCQuadrilateral()
	8	S = Serendipity(cell, 2)
	9
	10	x = sympy.DeferredVector("X")
	11	X, Y = x[0], x[1]
	12	basis_functions = [
	13	(1 - X)*(1 - Y),
	14	Y*(1 - X),
	15	X*(1 - Y),
	16	X*Y,
	17	Y(1 - X)(Y - 1),
	18	XY(Y - 1),
	19	X(1 - Y)(X - 1),
	20	XY(X - 1),
	21	]
	22	points = [[0.5, 0.5], [0.25, 0.75]]
	23	for alpha, actual in S.tabulate(2, points).items():
	24	expect = list(sympy.diff(basis, *zip([X, Y], alpha))
	25	for basis in basis_functions)
	26	expect = list([basis.subs(dict(zip([X, Y], point)))
	27	for point in points]
	28	for basis in expect)
	29	assert actual.shape == (8, 2)
	30	assert np.allclose(np.asarray(expect, dtype=float),
	31	actual.reshape(8, 2))