update mdspan (#142)

* update mdspan * Use basix.ufl (#143) * Use basix.ufl * tweak readme * Mscroggs/update ufl (#144) * Use basix.ufl * tweak readme * 1 not degree * keep mesh degree constant in Elasticity too --------- Co-authored-by: Matthew Scroggs <[email protected]>
FEniCS · Oct 7, 2023 · 770b395 · 770b395
1 parent 30a4cf9
commit 770b395
Show file tree

Hide file tree

Showing 4 changed files with 34 additions and 21 deletions.
diff --git a/README.md b/README.md
@@ -1,7 +1,7 @@
-# Performance test codes for FEniCS/DOLFIN-X
+# Performance test codes for FEniCSx/DOLFINx
 
 This repository contains solvers for testing the parallel performance of
-DOLFIN-X and the underlying linear solvers. It tests elliptic equations
+DOLFINx and the underlying linear solvers. It tests elliptic equations
 - Poisson equation and elasticity - in three dimensions.
 
 Representative performance data is available at
@@ -16,7 +16,7 @@ The source of the tests is in `src/` directory.
 
 ### Requirements
 
-- FEniCSx/DOLFINx installation (development version of dolfinx
+- FEniCSx/DOLFINx installation (development version of DOLFINx
   **required**)
 - PETSc installation
 - Boost Program Options

diff --git a/src/Elasticity.py b/src/Elasticity.py
@@ -4,8 +4,9 @@
 #
 # SPDX-License-Identifier:    MIT
 
-from ufl import (Coefficient, Identity, TestFunction, TrialFunction,
-                 VectorElement, dx, grad, inner, tetrahedron, tr)
+import basix.ufl
+from ufl import (Coefficient, Identity, FunctionSpace, Mesh, TestFunction, TrialFunction,
+                 dx, grad, inner, tetrahedron, tr)
 
 # Elasticity parameters
 E = 1.0e6
@@ -19,10 +20,12 @@
 
 forms = []
 for degree in range(1, 4):
-    element = VectorElement("Lagrange", cell, degree)
+    element = basix.ufl.element("Lagrange", "tetrahedron", degree, shape=(3, ))
+    domain = Mesh(basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3, )))
+    space = FunctionSpace(domain, element)
 
-    u, v = TrialFunction(element), TestFunction(element)
-    f = Coefficient(element)
+    u, v = TrialFunction(space), TestFunction(space)
+    f = Coefficient(space)
 
     def eps(v):
         return 0.5*(grad(v) + grad(v).T)

diff --git a/src/Poisson.py b/src/Poisson.py
@@ -4,21 +4,24 @@
 #
 # SPDX-License-Identifier:    MIT
 
-from ufl import (Coefficient, FiniteElement, TestFunction, TrialFunction, action, ds,
+import basix.ufl
+from ufl import (Coefficient, FunctionSpace, TestFunction, TrialFunction, Mesh, action, ds,
                  dx, grad, inner, tetrahedron)
 
 # Load namespace
 ns = vars()
 
 forms = []
 for degree in range(1, 4):
-    element = FiniteElement("Lagrange", tetrahedron, degree)
-
-    u = TrialFunction(element)
-    v = TestFunction(element)
-    f = Coefficient(element)
-    g = Coefficient(element)
-    un = Coefficient(element)
+    element = basix.ufl.element("Lagrange", "tetrahedron", degree)
+    domain = Mesh(basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3,)))
+    space = FunctionSpace(domain, element)
+
+    u = TrialFunction(space)
+    v = TestFunction(space)
+    f = Coefficient(space)
+    g = Coefficient(space)
+    un = Coefficient(space)
 
     aname = 'a' + str(degree)
     Lname = 'L' + str(degree)

diff --git a/src/elasticity_problem.cpp b/src/elasticity_problem.cpp
@@ -6,6 +6,7 @@
 
 #include "elasticity_problem.h"
 #include "Elasticity.h"
+#include <basix/mdspan.hpp>
 #include <dolfinx/common/Timer.h>
 #include <dolfinx/fem/DirichletBC.h>
 #include <dolfinx/fem/DofMap.h>
@@ -70,8 +71,11 @@ MatNullSpace build_near_nullspace(const fem::FunctionSpace<double>& V)
   }
 
   // Orthonormalize basis
-  la::orthonormalize(std::vector<std::reference_wrapper<la::Vector<T>>>(basis.begin(), basis.end()));
-  if (!la::is_orthonormal(std::vector<std::reference_wrapper<const la::Vector<T>>>(basis.begin(), basis.end())))
+  la::orthonormalize(std::vector<std::reference_wrapper<la::Vector<T>>>(
+      basis.begin(), basis.end()));
+  if (!la::is_orthonormal(
+          std::vector<std::reference_wrapper<const la::Vector<T>>>(
+              basis.begin(), basis.end())))
   {
     throw std::runtime_error("Space not orthonormal");
   }
@@ -145,9 +149,12 @@ elastic::problem(std::shared_ptr<mesh::Mesh<double>> mesh, int order)
       [](auto x) -> std::pair<std::vector<T>, std::vector<std::size_t>>
       {
         std::vector<T> vdata(x.extent(0) * x.extent(1));
-        namespace stdex = std::experimental;
-        stdex::mdspan<double,
-                      stdex::extents<std::size_t, 3, stdex::dynamic_extent>>
+        namespace stdex
+            = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE;
+        MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
+            T,
+            MDSPAN_IMPL_STANDARD_NAMESPACE::extents<
+                std::size_t, 3, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>
             v(vdata.data(), x.extent(0), x.extent(1));
         for (std::size_t p = 0; p < x.extent(1); ++p)
         {