Update JOSS paper (#324)

* Update JOSS paper

* Add Makefile

* Put C++ demo in a folder

* put prints in the same order as Python

* 3

* adjust author criteria, and edit paper authors to match

* remove duplicate name

* Add second cpp demo and some improvements to test.py (#325)

* add second cpp demo and some improvements in test.py

* include cpp demo at end of python demo on website

Co-authored-by: Matthew Scroggs <[email protected]>

Co-authored-by: Igor Baratta <[email protected]>

.github/workflows/dolfin-tests.yml

@@ -41,10 +41,7 @@ jobs:
 
       - name: Run cpp demos
         run: |
-          cd demo/cpp
-          mkdir build && cd build
-          cmake .. && make
-          ./basix-demo
+          python3 -m pytest demo/cpp/test.py
 
       - name: Get DOLFINx
         uses: actions/checkout@v2

.gitignore

@@ -24,9 +24,13 @@ doc/web/web
 doc/python/source/_autosummary
 
 demo/python/*.py.rst
+demo/cpp/*/_build
 
 build*/
 
+
+joss/paper.pdf
+
 .*.swp
 
 *.egg-info

README.md

@@ -38,7 +38,7 @@ https://docs.fenicsproject.org/basix/main/.
 In Basix, the sub-entities of the reference interval are numbered as
 follows:
 
-![The numbering of a reference interval](img/interval_numbering.png)
+![The numbering of a reference interval](joss/img/interval_numbering.png)
 
 The following elements are supported on a interval:
 
@@ -53,7 +53,7 @@ The following elements are supported on a interval:
 In Basix, the sub-entities of the reference triangle are numbered as
 follows:
 
-![The numbering of a reference triangle](img/triangle_numbering.png)
+![The numbering of a reference triangle](joss/img/triangle_numbering.png)
 
 The following elements are supported on a triangle:
 
@@ -72,7 +72,7 @@ The following elements are supported on a triangle:
 In Basix, the sub-entities of the reference quadrilateral are numbered
 as follows:
 
-![The numbering of a reference quadrilateral](img/quadrilateral_numbering.png)
+![The numbering of a reference quadrilateral](joss/img/quadrilateral_numbering.png)
 
 The following elements are supported on a quadrilateral:
 
@@ -91,7 +91,7 @@ The following elements are supported on a quadrilateral:
 In Basix, the sub-entities of the reference tetrahedron are numbered as
 follows:
 
-![The numbering of a reference tetrahedron](img/tetrahedron_numbering.png)
+![The numbering of a reference tetrahedron](joss/img/tetrahedron_numbering.png)
 
 The following elements are supported on a tetrahedron:
 
@@ -110,7 +110,7 @@ The following elements are supported on a tetrahedron:
 In Basix, the sub-entities of the reference hexahedron are numbered as
 follows:
 
-![The numbering of a reference hexahedron](img/hexahedron_numbering.png)
+![The numbering of a reference hexahedron](joss/img/hexahedron_numbering.png)
 
 The following elements are supported on a hexahedron:
 
@@ -129,7 +129,7 @@ The following elements are supported on a hexahedron:
 In Basix, the sub-entities of the reference prism are numbered as
 follows:
 
-![The numbering of a reference prism](img/prism_numbering.png)
+![The numbering of a reference prism](joss/img/prism_numbering.png)
 
 The following elements are supported on a prism:
 
@@ -140,7 +140,7 @@ The following elements are supported on a prism:
 In Basix, the sub-entities of the reference pyramid are numbered as
 follows:
 
-![The numbering of a reference pyramid](img/pyramid_numbering.png)
+![The numbering of a reference pyramid](joss/img/pyramid_numbering.png)
 
 The following elements are supported on a pyramid:
 

demo/cpp/CMakeLists.txt → demo/cpp/demo_create_and_tabulate/CMakeLists.txt

@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 3.16)
-project(basix-demo)
+project(demo_create_and_tabulate)
 
 # Use C++17
 set(CMAKE_CXX_STANDARD 17)

demo/cpp/main.cpp → demo/cpp/demo_create_and_tabulate/main.cpp

@@ -23,7 +23,7 @@ int main(int argc, char* argv[])
   // the degrees of freedom (DOFs) of the element in an equally spaced lattice.
   auto variant = basix::element::lagrange_variant::equispaced;
 
-  // Create the element lagrange element
+  // Create the lagrange element
   basix::FiniteElement lagrange
       = basix::create_element(family, cell_type, k, variant);
 
@@ -38,9 +38,9 @@ int main(int argc, char* argv[])
 
   xt::xtensor<double, 4> tab = lagrange.tabulate(0, points);
 
-  std::cout << "\nTabulate data shape: " << xt::adapt(tab.shape());
   std::cout << "\nTabulate data: \n"
             << xt::view(tab, 0, xt::all(), xt::all(), 0);
+  std::cout << "\nTabulate data shape: " << xt::adapt(tab.shape());
 
   return 0;
-}
+}

demo/cpp/demo_dof_transformations/CMakeLists.txt

@@ -0,0 +1,19 @@
+cmake_minimum_required(VERSION 3.16)
+project(demo_dof_transformations)
+
+# Use C++17
+set(CMAKE_CXX_STANDARD 17)
+
+# Require C++17
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+# Do not enable compler-specific extensions
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+find_package(Basix REQUIRED)
+
+# Executable
+add_executable(${PROJECT_NAME} main.cpp)
+
+# Target libraries
+target_link_libraries(${PROJECT_NAME} Basix::basix)

demo/cpp/demo_dof_transformations/main.cpp

@@ -0,0 +1,212 @@
+// ====================================
+// DOF permutations and transformations
+// ====================================
+
+// When using high degree finite elements on general meshes, adjustments
+// may need to be made to correct for differences in the orientation of
+// mesh entities on the mesh and on the reference cell. For example, in
+// a degree 4 Lagrange element on a triangle, there are 3 degrees of
+// freedom (DOFs) associated with each edge. If two neighbouring cells in
+// a mesh disagree on the direction of the edge, they could put an
+// incorrectly combine the local basis functions to give the wrong global
+// basis function.
+
+// This issue and the use of permutations and transformations to correct
+// it is discussed in detail in `Construction of arbitrary order finite
+// element degree-of-freedom maps on polygonal and polyhedral cell
+// meshes (Scroggs, Dokken, Richardons, Wells,
+// 2021) <https://arxiv.org/abs/2102.11901>`_.
+
+// Functions to permute and transform high degree elements are
+// provided by Basix. In this demo, we show how these can be used from C++.
+
+#include <basix/finite-element.h>
+#include <basix/lattice.h>
+#include <iomanip>
+#include <xtensor/xio.hpp>
+
+static const std::map<basix::cell::type, std::string> type_to_name
+    = {{basix::cell::type::point, "point"},
+       {basix::cell::type::interval, "interval"},
+       {basix::cell::type::triangle, "triangle"},
+       {basix::cell::type::tetrahedron, "tetrahedron"},
+       {basix::cell::type::quadrilateral, "quadrilateral"},
+       {basix::cell::type::pyramid, "pyramid"},
+       {basix::cell::type::prism, "prism"},
+       {basix::cell::type::hexahedron, "hexahedron"}};
+
+int main(int argc, char* argv[])
+{
+  // Degree 5 Lagrange element
+  // =========================
+  {
+    std::cout << "Degree 5 Lagrange element \n";
+
+    // Create a degree 5 Lagrange element on a triangle
+    auto family = basix::element::family::P;
+    auto cell_type = basix::cell::type::triangle;
+    int degree = 5;
+    auto variant = basix::element::lagrange_variant::equispaced;
+
+    // Create the lagrange element
+    basix::FiniteElement lagrange
+        = basix::create_element(family, cell_type, degree, variant);
+
+    // Print bools as true/false instead of 0/1
+    std::cout << std::boolalpha;
+
+    // The value of `dof_transformations_are_identity` is `false`: this tells
+    // us that permutations or transformations are needed for this element.
+    std::cout << "Dof transformations are identity: "
+              << lagrange.dof_transformations_are_identity() << "\n";
+
+    // The value of `dof_transformations_are_permutations` is `true`: this
+    // tells us that for this element, all the corrections we need to apply
+    // permutations. This is the simpler case, and means we make the
+    // orientation corrections by applying permutations when creating the
+    // DOF map.
+    std::cout << "Dof transformations are permutations: "
+              << lagrange.dof_transformations_are_permutations() << "\n";
+
+    // We can apply permutations by using the matrices returned by the
+    // method `base_transformations`. This method will return one matrix
+    // for each edge of the cell (for 2D and 3D cells), and two matrices
+    // for each face of the cell (for 3D cells). These describe the effect
+    // of reversing the edge or reflecting and rotating the face.
+
+    // For this element, we know that the base transformations will be
+    // permutation matrices.
+    std::cout << "\nBase transformations: \n";
+    std::cout << lagrange.base_transformations() << "\n";
+
+    // The matrices returned by `base_transformations` are quite large, and
+    // are equal to the identity matrix except for a small block of the
+    // matrix. It is often easier and more efficient to use the matrices
+    // returned by the method `entity_transformations` instead.
+
+    // `entity_transformations` returns a dictionary that maps the type
+    // of entity (`"interval"`, `"triangle"`, `"quadrilateral"`) to a
+    // matrix describing the effect of permuting that entity on the DOFs
+    // on that entity.
+    std::map<basix::cell::type, xt::xtensor<double, 3>> entity__transformation
+        = lagrange.entity_transformations();
+
+    std::cout << "\nEntity transformations: \n";
+    for (auto const& [cell, transformation] : entity__transformation)
+      std::cout << " -" << type_to_name.at(cell) << ":\n"
+                << transformation << std::endl;
+
+    // For this element, we see that this method returns one matrix for
+    // an interval: this matrix reverses the order of the four DOFs
+    // associated with that edge.
+
+    // In orders to work out which DOFs are associated with each edge,
+    // we use the attribute `entity_dofs`. For example, the following can
+    // be used to see which DOF numbers are associated with edge (dim 1)
+    // number 2:
+    int dim = 1;
+    int edge_num = 2;
+    std::vector<int> entity_dofs = lagrange.entity_dofs()[dim][edge_num];
+    std::cout << "\n Entity dofs of Edge number 2: ";
+    for (const auto dof : entity_dofs)
+      std::cout << dof << " ";
+  }
+  // Degree 2 Lagrange element
+  // =========================
+  {
+    // For a degree 2 Lagrange element, no permutations or transformations
+    // are needed.
+    std::cout << "\n\nDegree 2 Lagrange element \n";
+    auto family = basix::element::family::P;
+    auto cell_type = basix::cell::type::triangle;
+    int degree = 2;
+    auto variant = basix::element::lagrange_variant::equispaced;
+
+    // Create the lagrange element
+    auto lagrange = basix::create_element(family, cell_type, degree, variant);
+
+    // We can verify this by checking that`dof_transformations_are_identity` is
+    // `True`. To confirm that the transformations are identity matrices, we
+    // also print the base transformations.
+    std::cout << "Dof transformations are identity: "
+              << lagrange.dof_transformations_are_identity() << "\n";
+
+    std::cout << "Dof transformations are permutations: "
+              << lagrange.dof_transformations_are_permutations() << "\n";
+  }
+
+  // Degree 2 Nédélec element
+  // ========================
+  {
+    std::cout << "\n\nDegree 2 Nedelec element \n";
+    // For a degree 2 Nédélec (first kind) element on a tetrahedron, the
+    // corrections are not all permutations, so both
+    // `dof_transformations_are_identity` and
+    // `dof_transformations_are_permutations` are `False`.
+
+    auto family = basix::element::family::N1E;
+    auto cell_type = basix::cell::type::tetrahedron;
+    int degree = 2;
+    auto nedelec = basix::create_element(family, cell_type, degree);
+
+    std::cout << "Dof transformations are identity: "
+              << nedelec.dof_transformations_are_identity() << "\n";
+    std::cout << "Dof transformations are permutations: "
+              << nedelec.dof_transformations_are_permutations() << "\n";
+
+    // For this element, `entity_transformations` returns a map
+    // with two entries: a matrix for an interval that describes
+    // the effect of reversing the edge; and an array of two matrices
+    // for a triangle. The first matrix for the triangle describes
+    // the effect of rotating the triangle. The second matrix describes
+    // the effect of reflecting the triangle.
+
+    // For this element, the matrix describing the effect of rotating
+    // the triangle is
+
+    // .. math::
+    //    \left(\begin{array}{cc}-1&-1\\1&0\end{array}\right).
+
+    // This is not a permutation, so this must be applied when assembling
+    // a form and cannot be applied to the DOF numbering in the DOF map.
+    std::map<basix::cell::type, xt::xtensor<double, 3>> entity__transformation
+        = nedelec.entity_transformations();
+
+    std::cout << "\nEntity transformations: \n";
+    for (auto const& [cell, transformation] : entity__transformation)
+      std::cout << " -" << type_to_name.at(cell) << ":\n"
+                << transformation << std::endl;
+
+    // To demonstrate how these transformations can be used, we create a
+    // lattice of points where we will tabulate the element.
+
+    xt::xtensor<double, 2> points = basix::lattice::create(
+        cell_type, 5, basix::lattice::type::equispaced, true);
+    int num_points = points.shape(0);
+
+    // If (for example) the direction of edge 2 in the physical cell does
+    // not match its direction on the reference, then we need to adjust the
+    // tabulated data.
+
+    xt::xtensor<double, 4> original_data = nedelec.tabulate(0, points);
+    xt::xtensor<double, 4> mod_data = nedelec.tabulate(0, points);
+    xtl::span<double> data(mod_data.data(), mod_data.size());
+
+    // If the direction of edge 2 in the physical cell is reflected, it has
+    // cell permutation info `....000010` so (from right to left):
+    //      - edge 0 is not permuted (0)
+    //      - edge 1 is reflected (1)
+    //      - edge 2 is not permuted (0)
+    //      - edge 3 is not permuted (0)
+    //      - edge 4 is not permuted (0)
+    //      - edge 5 is not permuted (0)
+
+    int cell_info = 0b000010;
+    nedelec.apply_dof_transformation(data, num_points, cell_info);
+
+    std::cout << "\n Tabulated data is equal: "
+              << xt::allclose(original_data, mod_data) << std::endl;
+  }
+
+  return 0;
+}

demo/cpp/test.py

@@ -0,0 +1,17 @@
+import os
+import pytest
+
+# Get all the demos in this folder
+path = os.path.dirname(os.path.realpath(__file__))
+demos = []
+for folder in os.listdir(path):
+    if folder.startswith("demo_"):
+        subpath = os.path.join(path, folder)
+        if os.path.isdir(subpath) and os.path.isfile(os.path.join(subpath, "main.cpp")):
+            demos.append(folder)
+
+
+@pytest.mark.parametrize("demo", demos)
+def test_demo(demo):
+    demo_build = f"{path}/{demo}/_build"
+    assert os.system(f"mkdir -p {demo_build} && cd {demo_build} && cmake .. && make && ./{demo}") == 0

demo/python/demo_create_and_tabulate.py

@@ -45,3 +45,10 @@
 # had asked for derivatives too. The second dimension (5) is the number of points.
 # The third dimension (25) is the number of DOFs. The fourth dimension (1) is the
 # value size of the element: this will be greater than 1 for vector-values elements.
+#
+# C++ demo
+# ========
+# The following C++ code runs the same demo using Basix's C++ interface:
+#
+# .. literalinclude:: ../cpp/demo_create_and_tabulate/main.cpp
+#    :language: c++

demo/python/demo_dof_transformations.py

@@ -159,3 +159,10 @@
 # More efficient functions that apply the transformations and
 # permutations directly to data can be used via Basix's C++
 # interface.
+#
+# C++ demo
+# ========
+# The following C++ code runs the same demo using Basix's C++ interface:
+#
+# .. literalinclude:: ../cpp/demo_dof_transformations/main.cpp
+#    :language: c++