==============================================
Defining conforming Crouzeix--Raviart elements
==============================================

In this demo, we show how Basix's custom element functionality can be used to
create a conforming Crouzeix--Raviart element.

First, we import Basix and Numpy.

::

  import basix
  import numpy as np
  from basix import CellType, MapType, PolynomialType, LatticeType

  from mpl_toolkits import mplot3d  # noqa: F401
  import matplotlib.pyplot as plt

Conforming CR element on a triangle
===================================

We begin by implementing this element on a triangle. The following function
implements this element for an arbitrary degree. Details of the definition of
this element can be found at
https://defelement.com/elements/conforming-crouzeix-raviart.html.

For degrees 1 and 2, this element is equal to a degree 1 Lagrange element.
For higher degrees, we define the element. Note that in Basix, an element's degree
is necessarily equal to the highest degree polynomial in the space is spans, so
the degree we use here is one higher than that shown on DefElement.


::

  def create_ccr_triangle(degree):

      if degree in [1, 2]:
          wcoeffs = np.eye(3)

          x = [[], [], [], []]
          x[0].append(np.array([[0.0, 0.0]]))
          x[0].append(np.array([[1.0, 0.0]]))
          x[0].append(np.array([[0.0, 1.0]]))

          for _ in range(3):
              x[1].append(np.zeros((0, 2)))
          x[2].append(np.zeros((0, 2)))

          M = [[], [], [], []]
          for _ in range(3):
              M[0].append(np.array([[[1.]]]))

          for _ in range(3):
              M[1].append(np.zeros((0, 1, 0)))
          M[2].append(np.zeros((0, 1, 0)))

          return basix.create_custom_element(
              CellType.triangle, 1, [], wcoeffs, x, M, MapType.identity, False, 1)

      npoly = (degree + 1) * (degree + 2) // 2
      ndofs = (degree - 1) * (degree + 4) // 2
      wcoeffs = np.zeros((ndofs, npoly))

      dof_n = 0
      for i in range(degree * (degree + 1) // 2):
          wcoeffs[dof_n, dof_n] = 1
          dof_n += 1

      pts, wts = basix.make_quadrature(CellType.triangle, 2 * degree)
      poly = basix.tabulate_polynomials(PolynomialType.legendre, CellType.triangle, degree, pts)
      for i in range(1, degree - 1):
          x = pts[:, 0]
          y = pts[:, 1]
          f = x ** i * y ** (degree - 1 - i) * (x + y)

          for j in range(npoly):
              wcoeffs[dof_n, j] = sum(f * poly[:, j] * wts)
          dof_n += 1

      geometry = basix.geometry(CellType.triangle)
      topology = basix.topology(CellType.triangle)
      x = [[], [], [], []]
      M = [[], [], [], []]
      for v in topology[0]:
          x[0].append(np.array(geometry[v]))
          M[0].append(np.array([[[1.]]]))
      pts = basix.create_lattice(CellType.interval, degree - 1, LatticeType.equispaced, False)
      mat = np.zeros((len(pts), 1, len(pts)))
      mat[:, 0, :] = np.eye(len(pts))
      for e in topology[1]:
          edge_pts = []
          v0 = geometry[e[0]]
          v1 = geometry[e[1]]
          for p in pts:
              edge_pts.append(v0 + p * (v1 - v0))
          x[1].append(np.array(edge_pts))
          M[1].append(mat)
      pts = basix.create_lattice(CellType.triangle, degree, LatticeType.equispaced, False)
      x[2].append(pts)
      mat = np.zeros((len(pts), 1, len(pts)))
      mat[:, 0, :] = np.eye(len(pts))
      M[2].append(mat)

      return basix.create_custom_element(
          CellType.triangle, degree, [], wcoeffs, x, M, MapType.identity, False, degree - 1)


We can then create a degree 3 conforming CR element.

::

  e = create_ccr_triangle(3)

We now visualise the basis functions of the element we have created.

::

  pts = basix.create_lattice(CellType.triangle, 30, LatticeType.equispaced, True)

  x = pts[:, 0]
  y = pts[:, 1]
  z = e.tabulate(0, pts)[0]

  fig = plt.figure(figsize=(8, 8))

  for n in range(7):
      if n == 6:
          ax = plt.subplot(3, 3, n + 2, projection='3d')
      else:
          ax = plt.subplot(3, 3, n + 1, projection='3d')
      ax.plot([0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 0], "k-")
      ax.scatter(x, y, z[:, n, 0])

  plt.savefig("ccr_triangle_3.png")

.. image:: ccr_triangle_3.png
 :width: 100%
 :alt: The basis functions of a degree 3 conforming CR element

We also visualise the basis functions of a degree 4 conforming CR element.

::

  e = create_ccr_triangle(4)

  pts = basix.create_lattice(CellType.triangle, 30, LatticeType.equispaced, True)

  x = pts[:, 0]
  y = pts[:, 1]
  z = e.tabulate(0, pts)[0]

  fig = plt.figure(figsize=(11, 8))

  for n in range(12):
      ax = plt.subplot(3, 4, n + 1, projection='3d')
      ax.plot([0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 0], "k-")
      ax.scatter(x, y, z[:, n, 0])

  plt.savefig("ccr_triangle_4.png")

.. image:: ccr_triangle_4.png
 :width: 100%
 :alt: The basis functions of a degree 3 conforming CR element