shapefuncs Module

Lagrange shape functions and their first and second derivatives on the reference elements of the finite element method, including mapped infinite elements.

A shapefunc evaluates the $n$ shape functions $N_i(\xi)$ of an element at the points of a cubature, together with $\partial N_i/\partial\xi_j$ and $\partial^2 N_i/\partial\xi_j\partial\xi_k$. The reference elements and coordinates are those of the cubatures package.

Elements and orders

Constant	Element	Reference domain	Nodes, order $p$	Orders
CUB_LIN	line	$[-1,1]$	$p+1$	$\ge 0$
CUB_QUA	quadrilateral	$[-1,1]^2$	$(p+1)^2$	$\ge 0$
CUB_HEX	hexahedron	$[-1,1]^3$	$(p+1)^3$	$\ge 0$
CUB_TRI	triangle	$x,y \ge 0,\ x+y \le 1$	$(p+1)(p+2)/2$	$\ge 1$
CUB_TET	tetrahedron	$x,y,z \ge 0,\ x+y+z \le 1$	$(p+1)(p+2)(p+3)/6$	$\ge 1$
CUB_WED	wedge / prism	triangle $\times\ [-1,1]$	$(p+1)^2(p+2)/2$	$\ge 1$

Lines, quadrilaterals and hexahedra use tensor products of Lagrange polynomials on equispaced nodes; order 0 is the constant function with one node at the centre. Triangles and tetrahedra use the complete polynomials of total degree $p$ in the barycentric coordinates (Silvester's polynomials [5]); wedges use the triangle functions times Lagrange polynomials in the axial coordinate. There is no upper limit on the order, but equispaced interpolation becomes ill-conditioned at high orders.

Serendipity elements

With family = SHP_SERENDIPITY, quadrilaterals and hexahedra of order 1 to 3 have nodes on their vertices and edges only: 4, 8 and 12 nodes for quadrilaterals, 8, 20 and 32 for hexahedra. They span the serendipity space: the monomials of superlinear degree $\le p$ [7]. Each serendipity function is a fixed combination of the tensor-product Lagrange functions of the same order, computed once by set. Serendipity elements must be finite.

Factor coordinates

Every shape function is a product of one-dimensional functions of factor coordinates $s_k$, which are affine in $\xi$:

• line, quadrilateral, hexahedron: $s = \xi$;
• triangle: $s = (\xi_1, \xi_2, 1-\xi_1-\xi_2)$, the barycentric coordinates;
• tetrahedron: $s = (\xi_1, \xi_2, \xi_3, 1-\xi_1-\xi_2-\xi_3)$;
• wedge: $s = (\xi_1, \xi_2, 1-\xi_1-\xi_2, \xi_3)$.

Node $i$ has one integer lattice index $a_{ki}$ per factor coordinate, and $N_i(\xi) = \prod_k \phi_{a_{ki}}(s_k)$. Derivatives follow from the product rule and the constant $\partial s/\partial\xi$.

Node numbering

Nodes are numbered hierarchically: vertices, then edge nodes, then face nodes, then interior nodes.

• Vertices of a quadrilateral run counter-clockwise from $(-1,-1)$. A hexahedron has the bottom face $\xi_3=-1$ counter-clockwise from $(-1,-1,-1)$, then the top face in the same order. Triangle and tetrahedron vertex $k$ is where barycentric coordinate $k$ equals 1, so vertex 1 is $\xi = (1,0,\dots)$ and the last vertex is the origin. Wedge vertices are the triangle vertices at $\xi_3=-1$, then at $\xi_3=1$.
• Edges are listed below. Edge nodes run from the first vertex of the edge to the second.

Element	Edges
quadrilateral	1-2, 2-3, 3-4, 4-1
hexahedron	1-2, 2-3, 3-4, 4-1, 5-6, 6-7, 7-8, 8-5, 1-5, 2-6, 3-7, 4-8
triangle	1-2, 2-3, 3-1
tetrahedron	1-2, 2-3, 3-1, 1-4, 2-4, 3-4
wedge	1-2, 2-3, 3-1, 4-5, 5-6, 6-4, 1-4, 2-5, 3-6

• Faces of a hexahedron are $\xi_3=-1$, $\xi_3=1$, $\xi_1=-1$, $\xi_1=1$, $\xi_2=-1$, $\xi_2=1$. Tetrahedron faces are 1-2-3, 1-2-4, 2-3-4, 1-3-4. Wedge faces are the bottom and top triangles, then the quadrilaterals on edges 1-2, 2-3, 3-1.
• Interiors of quadrilateral faces and elements, and of triangles and tetrahedra, are numbered recursively as an element of the same type of order $p-2$ (quadrilateral, hexahedron), $p-3$ (triangle) or $p-4$ (tetrahedron). A quadrilateral face in $(\xi_j, \xi_k)$, (j

Variables

Type	Visibility	Attributes		Name		Initial
integer,	public,	parameter	::	SHP_FIN	=	1	Finite direction: Lagrange polynomials on $[-1,1]$
integer,	public,	parameter	::	SHP_INF	=	2	Infinite direction, mapping functions for coordinates
integer,	public,	parameter	::	SHP_CHP	=	3	Infinite direction, chopped functions for the solution
integer,	public,	parameter	::	SHP_LAGRANGE	=	1	Full Lagrange (tensor-product or complete) elements
integer,	public,	parameter	::	SHP_SERENDIPITY	=	2	Serendipity quadrilaterals and hexahedra, orders 1 to 3

Interfaces

public interface shapefunc

Construct shape functions, e.g. shapefunc(q, 2), shapefunc(q, 1, [SHP_FIN, SHP_INF]) or shapefunc(q, 2, family=SHP_SERENDIPITY, nodes=perm)

• private pure function new_shapefunc(q, order, infin, family, nodes) result(self)

  Construct shape functions of order order at the points of q; see set

  Arguments

  Type	Intent	Optional	Attributes		Name
  type(cubature),	intent(in)			::	q	Cubature on the reference element
  integer,	intent(in)			::	order	Polynomial order
  integer,	intent(in),	optional		::	infin(:)	Infinitude, SHP_FIN, SHP_INF or SHP_CHP, size 1 or q%dim
  integer,	intent(in),	optional		::	family	SHP_LAGRANGE (default) or SHP_SERENDIPITY
  integer,	intent(in),	optional		::	nodes(:)	Node ordering: nodes(k) is the default number of node k

  Return Value type(shapefunc)

Derived Types

type, public ::  shapefunc

Shape functions and derivatives of one element type and order at the points of a cubature.

Read more…

Components

Type	Visibility	Attributes		Name		Initial
integer,	public		::	elm	=	0	Element type, CUB_*
integer,	public		::	dim	=	0	Spatial dimension
integer,	public		::	order	=	0	Polynomial order
integer,	public		::	family	=	SHP_LAGRANGE	SHP_LAGRANGE or SHP_SERENDIPITY
integer,	public		::	nnodes	=	0	Number of nodes
integer,	public		::	npoints	=	0	Number of cubature points
integer,	public		::	infin(3)	=	SHP_FIN	Infinitude per direction, SHP_*
integer,	public,	allocatable	::	lattice(:,:)			Lattice indices, shape [nfactors, nnodes], see module documentation
real(kind=rk),	public,	allocatable	::	coords(:,:)			Nodal coordinates, shape [dim, nnodes]
real(kind=rk),	public,	allocatable	::	func(:,:)			Values $N_i$, shape [nnodes, npoints]
real(kind=rk),	public,	allocatable	::	derv(:,:,:)			First derivatives $\partial N_i/\partial\xi_j$, shape [dim, nnodes, npoints]
real(kind=rk),	public,	allocatable	::	curv(:,:,:,:)			Second derivatives $\partial^2 N_i/\partial\xi_j\partial\xi_k$, shape [dim, dim, nnodes, npoints]

Constructor

Construct shape functions, e.g. shapefunc(q, 2), shapefunc(q, 1, [SHP_FIN, SHP_INF]) or shapefunc(q, 2, family=SHP_SERENDIPITY, nodes=perm)

private pure function new_shapefunc (q, order, infin, family, nodes)

Construct shape functions of order order at the points of q; see set

Type-Bound Procedures

procedure, public :: set	Build the shape functions in place
procedure, public :: eval	Evaluate at any point of the reference element
procedure, public :: is_valid	True if set and consistent
procedure, public :: summary	Write element, order, family, infinitude, nodes, points
procedure, public :: show	Write the summary, values and first derivatives
procedure, public :: numbering	Write the summary and the nodal coordinates
procedure, public :: destroy	Deallocate and reset