fenics/basix

FEniCSx finite element basis evaluation library

basix's Issues

Set Debug mode in CI

Implement base transformations for Regge element

Expose cell::str_to_type and cell::type_to_str to python layer

To make it as easy as possible to use basix along-side dolfin-x, I suggest that cell:str_to_type and cell::type_to_str is expose, such that one can easily convert between dolfinx::mesh::CellType and basix::CellType.
This is for instance relevant when creating custom quadrature rules, as make_quadrature require a basix cell type.
In dolfin-x, one can convert the dolfinx::mesh::CellType to a string using dolfinx.cpp.mesh.to_string(mesh.topology.cell_type)

Remove "auto" with Eigen wherever possible

e.g. in polyset.cpp

Implement combining of base transformations

It would be useful to have a function that takes a uint32_t that decodes the cell orientation, then returns the permutation. Dolfinx could then use this instead of the version in the ffc generated code (see FEniCS/ffcx#300)

Create Spack package

We need to create a Spack package rather urgently as the DOLFINX Spack package is broken and will remain so until basix is added to Spack.

Add permutation information to cell

Improve documentation

Some specific suggestions:
FEniCS/dolfinx#1324 (comment)
FEniCS/dolfinx#1324 (comment)

We should also write some demos

ffcx interfacing

Work in progress at https://github.com/FEniCS/ffcx/tree/chris/libtab

Improve quadrature interface

Currently we have make_quadrature which generates Gauss-Jacobi rules on all the cell types.

• We need a way to choose different schemes.
  • User defined
  • Pre-defined with a name, e.g. "Gauss-Jacobi", "Gauss-Radau", "Strang-Fix" or just "default"

Create INSTALL instructions

for python and cmake imstalls.

Update install instructions

Reflect changes after restructuring of build systems.

Tidy up entity_dofs and topology

Stop using int& references to create entity_dofs

pip install fails with new eigen versions on AMD machines

$ pip3 install git+https://github.com/FEniCS/basix.git --upgrade
...
...
-- Build files have been written to: /tmp/pip-req-build-xvapwdaa/_skbuild/linux-x86_64-3.8/cmake-build
  [1/18] Building CXX object cpp/CMakeFiles/basix.dir/dof-permutations.cpp.o
  FAILED: cpp/CMakeFiles/basix.dir/dof-permutations.cpp.o
  /usr/bin/c++ -DBASIX_VERSION=0.0.1 -DEIGEN_MAX_ALIGN_BYTES=32 -Dbasix_EXPORTS -I/usr/local/include/eigen3 -Wno-comment -Wall -Werror -std=c++17  -O3 -DNDEBUG -fPIC -std=c++17 -MD -MT cpp/CMakeFiles/basix.dir/dof-permutations.cpp.o -MF cpp/CMakeFiles/basix.dir/dof-permutations.cpp.o.d -o cpp/CMakeFiles/basix.dir/dof-permutations.cpp.o -c ../../../cpp/dof-permutations.cpp
  In file included from /usr/local/include/eigen3/Eigen/Core:159,
                   from /usr/local/include/eigen3/Eigen/Dense:1,
                   from ../../../cpp/dof-permutations.h:7,
                   from ../../../cpp/dof-permutations.cpp:5:
  /usr/local/include/eigen3/Eigen/src/Core/util/Memory.h: In function ‘void Eigen::internal::queryCacheSizes_amd(int&, int&, int&)’:
  /usr/local/include/eigen3/Eigen/src/Core/util/Memory.h:1064:14: error: comparison of integer expressions of different signedness: ‘int’ and ‘unsigned int’ [-Werror=sign-compare]
   1064 |   if(abcd[0] >= 0x80000006)
        |      ~~~~~~~~^~~~~~~~~~~~~

The issue can be solved by removing the -Werror flag from release mode.

It happens with eigen 3.9 and the development versions.

Simplify integral moment and interpolation matrix generation

Currently, there is a lot of very similar code in these functions, eg see https://github.com/FEniCS/libtab/blob/723b70a72a52d4867af93321303aafc44a589aee/cpp/moments.cpp#L62-L100 and https://github.com/FEniCS/libtab/blob/723b70a72a52d4867af93321303aafc44a589aee/cpp/moments.cpp#L134-L170.

Use GLL points for tensor product Lagrange elements

Default to GLL points rather than equispaced points.

Use Eigen::Matrix when an object is logically a matrix

Probably need to go through the code quite carefully. The different behaviour of Array and Matrix multiplication is used quite extensively.

Header and library installation and finder

Probably should have something like get_prefix() to find the prefix dir for headers and library.

Remove tp.cpp

Put "TensorProduct" into Lagrange elements

Rename base permutations to base transformations

Documentation with doxygen and sphinx

• Go through all documentation checking it works OK with doxygen
• Add docstrings to python using separate file e.g. wrapper_docs.h and a lookup map.

Check array order consistency - RowMajor vs ColumnMajor

Example:
basix::tabulate receives a pointer to the memory of a set of points [npoints, tdim], which is expected to be in row-major order, maps it into a row-major eigen array, but in the subsequent function calls it's treated as a column-major array.

Add interface for performing interpolation

It would be helpful to support interpolation.

• For an element, return the points on the reference cell where it requires function evaluations
• Given function evaluations on the reference element, return dofs.

Fix mapping name

Should be "identity" as per UFL, not "affine".

value_shape

Some elements may need a value shape as well as a value size

Transpose the data returned by tabulate

Currently, tabulate returns (row-major) data with shape [(nd+tdim)!/nd!tdim!, value_size, dim, npoints]. However data in dolfinx is stored (row-major) as [npoints, dim, value_size].

This means that the same implementation of (eg) map_push_forward cannot be applied to both the result of FiniteElement::tabulate and data passed in from dolfinx.

Make pip installable

See some discussion at pybind/pybind11#2326 on pip, CMake, and building dependencies.

Branch https://github.com/FEniCS/libtab/tree/skbuild now installs with pip and scikit-build.
It will also build and install with cmake, which is nice. It could probably be improved, and we still need to export the public C++ API headers somehow.

Restructure the function mappings

• basix::map_push_forward and basix::map_pull_back should work on multiple points at ones
• FiniteElement::map_push_forward and FiniteElement::map_pull_back should use in/out arguments, rather than returning arrays
• A pointer to the correct mapping function should be stored in FiniteElement instead of querying the mapping type for every point and every value

Enable editable pip install (-e)

Support pip install -e ..

Interface to quadrature rules

FFCx wants to use quadrature rules from libtab (https://github.com/FEniCS/ffcx/blob/3c7727c75afcb4882214f1678d005eb92e2eb737/ffcx/libtab_interface.py#L15), but currently I think they are only available externally for simplices

Use shorter family names

Use names which are shorter, e.g. from FEEC.

Add "create_element" in C++

Interface to get an element from string ("family", "celltype", degree)

Support "variants" of elements

Do we have "variants" of elements, e.g. Lagrange at GLL points? Discontinuous Lagrange? Point eval vs. integral moment?

Or should each variant just have a new name?

Implement Piola mappings

These can then be used directly by dolfinx (see FEniCS/ffcx#299)

Inform whether a call to compute_reference_geometry is needed for interpolation

Simple regression example:

Number of degrees of freedom: 1030301
element = FiniteElement("Lagrange", tetrahedron, 2)
BoxMesh {50, 50, 50}

Build BoxMesh               |     1  1.140848  1.140848
Build dofmap data           |     2  0.985655  1.971311
Assemble Vector             |     1  0.907553  0.907553


interpolate                 |     1  5.931514  5.931514

Without compute_reference_geometry:

interpolate                 |     1  2.030478 2.030478

Store entity-dof information

There needs to be a way of associating dofs with topological entities.

At first, we can try just storing the dofs in increasing topological order i.e. vertex, edge, face, cell.
In that case, we just need to know how many dofs there are on each dimension.
We can make the assumption that all similar entities of a given dimension must have the same number of dofs (otherwise continuity is difficult). This is a simplification, and may be revised later if necessary.

For face and cell dofs, the layout is not simple (unlike vertex and edge), especially when a cell may have both triangular and quadrilateral faces. In this case, we could specify the number of "edge_dofs", e.g. for a triangle edge_dofs = 3 implies 6 dofs. For a quadrialteral, edge_dofs = 3 implies 9 dofs.

Add all elements used in FEniCSx tests

Used in dolfinx tests:

• Brezzi-Douglas-Marini
• RTCE
• RTCF
• NCE
• NCF
• DQ

Used in ffc demos:

• #379
• #261
• Bubble (#71)

Missing spaces that probably need to be dealt with in FFCx, not in Basix:

• Hdiv trace
• EnrichedElement
• NodalEnrichedElement
• Quadrature
• Real

Custom element in dolfin

In #62, we make finite element class store a family name string. This will be used by dolfin to recontruct the element for tabulation.

For custom element (eg the elements in test_create.py), this will not work as the element cannot be constructed.

Make quadrature points C-contiguous

One of the UFL file demos in dolfin was throwing an error ("ndarray is not C-contiguous") when trying to hash an array of quadrature points obtained from libtab. I added this to get around it: https://github.com/FEniCS/ffcx/blob/90055afa735312cde8544f42d411c4c540660a8a/ffcx/ir/representationutils.py#L21

We may be able to avoid reordering there by changing the return type of create_quadrature?

Check topology/geometry

Especially for quad, hex, prism, pyramid.

Implement all elements from the periodic table

These are:

• Simplices
  • First family
    • Lagrange
    • Nedlec first kind
    • Raviart-Thomas
  • Second family
    • Lagrange
    • Nedelec second kind
    • BDM
• TP
  • First family
    • Q
    • RTCe / NCe
    • RTCf / NCf
  • Second family
    • Serendipity
    • BDMCe / AAe
    • BDMCf / AAf

Add unique integer identifier for each element type

It would be useful to have an integer identifier for each element type, and embed in this ID the basix version number (or release short git hash). Generated FFCX code could then just store the unique element ID, and the element can be then be easily and reliably fetched from basix.

Remove EIgen from public interface and use plain C API

Exposing Eigen in the public interface can cause all sort of problems when interfacing with other libraries that may use a different alignment. Simple fix is to not expose Eigen.

Add Bubble element

Should be easy as it is just based on Lagrange. Needed for "mini" element demos.

Test the cpp interface

Currently all the tests use the Python interface (via wrapper.cpp).

It would be good to add tests that directly interface through cpp so we can test this without having to use the dolfin tests

Integration with dolfinx

• get libtab from Python
• Find header files of libtab

Folder structure

As more elements are added, we're going to end up with a lot of files in the cpp/ directory. Would it be neater to put them in subdirectories, eg:

• cpp contains basix and wrapper
• cpp/core containing cell, dof-permutations, finite-element, indexing, moment, polyset, quadrature
• cpp/elements containing lagrange, crouzeix-raviart, raviart-thomas, nedelec, etc

Nedelec 2nd kind H(curl) fails on 32-bit systems

Judging by the build errors at https://buildd.debian.org/status/package.php?p=basix, looks like Nedelec 2nd kind H(curl) is not compatible with 32 bit systems. e.g.
i386 https://buildd.debian.org/status/fetch.php?pkg=basix&arch=i386&ver=0.0.1%7Egit20210122.4f10ef2-1&stamp=1612391676&raw=0
armel https://buildd.debian.org/status/fetch.php?pkg=basix&arch=armel&ver=0.0.1%7Egit20210122.4f10ef2-1&stamp=1612398463&raw=0
armhf https://buildd.debian.org/status/fetch.php?pkg=basix&arch=armhf&ver=0.0.1%7Egit20210122.4f10ef2-1&stamp=1612392128&raw=0

Is this problem expected? Is it something that can be fixed, or should Nedelec 2nd kind H(curl) just be considered not supported on 32 bit systems?

The error message is

__ test_permutation_of_tabulated_data_tetrahedron[5-Nedelec 2nd kind H(curl)] __

element_name = 'Nedelec 2nd kind H(curl)', order = 5

    @pytest.mark.parametrize("element_name", tetrahedron_elements)
    @pytest.mark.parametrize("order", range(1, 6))
    def test_permutation_of_tabulated_data_tetrahedron(element_name, order):
        if element_name == "Crouzeix-Raviart" and order != 1:
            pytest.xfail()
    
        e = basix.create_element(element_name, "tetrahedron", order)
    
        N = 4
        points = np.array([[i / N, j / N, k / N]
                           for i in range(N + 1) for j in range(N + 1 - i) for k in range(N + 1 - i - j)])
        values = e.tabulate(0, points)[0]
    
        start = sum(e.entity_dofs[0])
        ndofs = e.entity_dofs[1][0]
        if ndofs != 0:
            # Check that the 0th permutation undoes the effect of reflecting edge 0
            reflected_points = np.array([[p[0], p[2], p[1]] for p in points])
            reflected_values = e.tabulate(0, reflected_points)[0]
    
            if e.mapping_name == "affine":
                pass
            elif e.mapping_name == "covariant piola":
                for i, r in enumerate(reflected_values):
                    for j in range(e.dim):
                        reflected_values[i][j + e.dim::e.dim] = r[j + e.dim::e.dim][::-1]
            elif e.mapping_name == "contravariant piola":
                for i, r in enumerate(reflected_values):
                    for j in range(e.dim):
                        reflected_values[i][j] = -r[j]
                        reflected_values[i][j + e.dim::e.dim] = -r[j + e.dim::e.dim][::-1]
            elif e.mapping_name == "double covariant piola":
                pytest.skip()  # TODO: implement double covariant piola
            else:
                raise ValueError(f"Unknown mapping: {e.mapping_name}")
    
            for i, j in zip(values, reflected_values):
                for d in range(e.value_size):
                    i_slice = i[d * e.dim:(d + 1) * e.dim]
                    j_slice = j[d * e.dim:(d + 1) * e.dim]
>                   assert np.allclose((e.base_permutations[0].dot(i_slice))[start: start + ndofs],
                                       j_slice[start: start + ndofs])
E                   assert False
E                    +  where False = <function allclose at 0xf6d2a658>(array([ 0., -0., -0., -0., -0.,  0.]), array([-1.34, -1.12, -0.28,  0.68,  0.8 , -1.5 ]))
E                    +    where <function allclose at 0xf6d2a658> = np.allclose

test/test_permutations.py:227: AssertionError

In the meantime since only this one test fails, I'll just skip this test in the 32-bit Debian builds. So the error won't show in future builds (but will still be there).

Create a public API

• The library should provide single header API that should not expose any other libtab headers
• Provide a plain C interface
• Wrap with pybind11 just the public API

Would it be convenient to provide an API that using Eigen? Would this be portable?

Check that ndofs is consistent

Should check ndofs (i.e. coeffs.rows()) against sum(entity_dofs) in constructor.