Add default scalar type

python/demo/demo_biharmonic.py

@@ -112,15 +112,14 @@
 
 # +
 import numpy as np
-
 import ufl
-from dolfinx import fem, io, mesh, plot
 from dolfinx.mesh import CellType, GhostMode
+from mpi4py import MPI
+from petsc4py.PETSc import ScalarType
 from ufl import (CellDiameter, FacetNormal, avg, div, dS, dx, grad, inner,
                  jump, pi, sin)
 
-from mpi4py import MPI
-from petsc4py.PETSc import ScalarType
+from dolfinx import fem, io, mesh, plot
 
 # -
 
@@ -149,15 +148,12 @@
 # function that returns `True` for points `x` on the boundary and
 # `False` otherwise.
 
-facets = mesh.locate_entities_boundary(
-    msh, dim=1,
-    marker=lambda x: np.logical_or.reduce((
-        np.isclose(x[0], 0.0),
-        np.isclose(x[0], 1.0),
-        np.isclose(x[1], 0.0),
-        np.isclose(x[1], 1.0)
-    ))
-)
+facets = mesh.locate_entities_boundary(msh, dim=1,
+                                       marker=lambda x: np.logical_or.reduce((
+                                           np.isclose(x[0], 0.0),
+                                           np.isclose(x[0], 1.0),
+                                           np.isclose(x[1], 0.0),
+                                           np.isclose(x[1], 1.0))))
 
 # We now find the degrees-of-freedom that are associated with the
 # boundary facets using {py:func}`locate_dofs_topological
@@ -217,9 +213,8 @@
 # case we use a direct (LU) solver. The {py:func}`solve
 # <dolfinx.fem.LinearProblem.solve>` will compute a solution.
 
-problem = fem.petsc.LinearProblem(a, L, bcs=[bc],
-                                  petsc_options={"ksp_type": "preonly",
-                                                 "pc_type": "lu"})
+problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly",
+                                                                 "pc_type": "lu"})
 uh = problem.solve()
 
 # The solution can be written to a  {py:class}`XDMFFile
@@ -247,7 +242,6 @@
         plotter.screenshot("uh_biharmonic.png")
     else:
         plotter.show()
-
 except ModuleNotFoundError:
     print("'pyvista' is required to visualise the solution")
     print("Install 'pyvista' with pip: 'python3 -m pip install pyvista'")

python/demo/demo_elasticity.py

@@ -24,9 +24,7 @@
 
 # +
 import numpy as np
-
 import ufl
-from dolfinx import la
 from dolfinx.fem import (Expression, Function, FunctionSpace,
                          VectorFunctionSpace, dirichletbc, form,
                          locate_dofs_topological)
@@ -35,10 +33,11 @@
 from dolfinx.io import XDMFFile
 from dolfinx.mesh import (CellType, GhostMode, create_box,
                           locate_entities_boundary)
-from ufl import dx, grad, inner
-
 from mpi4py import MPI
 from petsc4py import PETSc
+from ufl import dx, grad, inner
+
+from dolfinx import la
 
 dtype = PETSc.ScalarType
 # -

python/demo/demo_interpolation-io.py

@@ -20,13 +20,11 @@
 
 # +
 import numpy as np
-
-from dolfinx import plot
 from dolfinx.fem import Function, FunctionSpace, VectorFunctionSpace
 from dolfinx.mesh import CellType, create_rectangle, locate_entities
-
 from mpi4py import MPI
-from petsc4py.PETSc import ScalarType
+
+from dolfinx import default_scalar_type, plot
 
 # -
 
@@ -38,7 +36,7 @@
 # Create a Nédélec function space and finite element Function
 
 V = FunctionSpace(msh, ("Nedelec 1st kind H(curl)", 1))
-u = Function(V, dtype=ScalarType)
+u = Function(V, dtype=default_scalar_type)
 
 # Find cells with *all* vertices (0) $x_0 <= 0.5$ or (1) $x_0 >= 0.5$:
 
@@ -57,7 +55,7 @@
 # interpolate the $H({\rm curl})$ function `u`
 
 V0 = VectorFunctionSpace(msh, ("Discontinuous Lagrange", 1))
-u0 = Function(V0, dtype=ScalarType)
+u0 = Function(V0, dtype=default_scalar_type)
 u0.interpolate(u)
 
 # We save the interpolated function `u0` in VTX format. When visualising

python/demo/demo_lagrange_variants.py

@@ -18,19 +18,17 @@
 # We begin this demo by importing the required modules.
 
 # +
-import matplotlib.pylab as plt
-import numpy as np
-
 import basix
 import basix.ufl
+import matplotlib.pylab as plt
+import numpy as np
 import ufl
-from dolfinx import fem, mesh
+from mpi4py import MPI
 from ufl import ds, dx, grad, inner
 
-from mpi4py import MPI
-from petsc4py.PETSc import ScalarType
+from dolfinx import default_scalar_type, fem, mesh
 
-if np.issubdtype(ScalarType, np.complexfloating):
+if np.issubdtype(default_scalar_type, np.complexfloating):
     print("Demo should only be executed with DOLFINx real mode")
     exit(0)
 # -
@@ -60,7 +58,7 @@
                                basix.LagrangeVariant.equispaced)
 pts = basix.create_lattice(basix.CellType.interval, 200, basix.LatticeType.equispaced, True)
 values = element.tabulate(0, pts)[0, :, :, 0]
-if MPI.COMM_WORLD.size == 1:  # Skip this plotting in parallel
+if MPI.COMM_WORLD.size == 1:
     for i in range(values.shape[1]):
         plt.plot(pts, values[:, i])
     plt.plot(element.points, [0 for _ in element.points], "ko")
@@ -123,7 +121,7 @@
                                        marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),
                                                                       np.isclose(x[0], 2.0)))
 dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
-bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)
+bc = fem.dirichletbc(value=default_scalar_type(0), dofs=dofs, V=V)
 
 u = ufl.TrialFunction(V)
 v = ufl.TestFunction(V)

python/dolfinx/__init__.py

@@ -26,13 +26,17 @@
 
 import sys
 
-from dolfinx import common
-from dolfinx import cpp as _cpp
-from dolfinx import fem, geometry, graph, io, jit, la, log, mesh, nls, plot
+from petsc4py import PETSc as _PETSc
+
+default_scalar_type = _PETSc.ScalarType
+
 # Initialise logging
 from dolfinx.common import (TimingType, git_commit_hash, has_debug, has_kahip,
                             has_parmetis, list_timings, timing)
-# Import cpp modules
+
+from dolfinx import common
+from dolfinx import cpp as _cpp
+from dolfinx import fem, geometry, graph, io, jit, la, log, mesh, nls, plot
 from dolfinx.cpp import __version__
 
 _cpp.common.init_logging(sys.argv)
@@ -49,6 +53,7 @@ def get_include(user=False):
         return os.path.join(os.path.dirname(d), "src")
 
 
+
 __all__ = [
     "fem", "common", "geometry", "graph", "io", "jit", "la", "log", "mesh", "nls", "plot",
     "TimingType", "git_commit_hash", "has_debug", "has_kahip", "has_parmetis", "list_timings",

python/dolfinx/fem/forms.py

@@ -10,20 +10,18 @@
 import collections.abc
 import typing
 
+import numpy as np
+import ufl
 from dolfinx.fem import IntegralType
 from dolfinx.fem.function import FunctionSpace
 
+from dolfinx import cpp as _cpp
+from dolfinx import default_scalar_type, jit
+
 if typing.TYPE_CHECKING:
     from dolfinx.fem import function
     from dolfinx.mesh import Mesh
 
-import numpy as np
-import ufl
-from petsc4py import PETSc
-
-from dolfinx import cpp as _cpp
-from dolfinx import jit
-
 
 class FormMetaClass:
     def __init__(self, form, V: list[_cpp.fem.FunctionSpace], coeffs, constants,
@@ -93,7 +91,7 @@ def integral_types(self):
                           "vertex": IntegralType.vertex}
 
 
-def form(form: typing.Union[ufl.Form, typing.Iterable[ufl.Form]], dtype: np.dtype = PETSc.ScalarType,
+def form(form: typing.Union[ufl.Form, typing.Iterable[ufl.Form]], dtype: np.dtype = default_scalar_type,
          form_compiler_options: dict = {}, jit_options: dict = {}):
     """Create a Form or an array of Forms.
 

python/dolfinx/fem/function.py

@@ -22,11 +22,10 @@
 import ufl.algorithms
 import ufl.algorithms.analysis
 from dolfinx.fem import dofmap
-from petsc4py import PETSc
 from ufl.domain import extract_unique_domain
 
 from dolfinx import cpp as _cpp
-from dolfinx import jit, la
+from dolfinx import default_scalar_type, jit, la
 
 
 class Constant(ufl.Constant):
@@ -85,7 +84,7 @@ def __complex__(self):
 class Expression:
     def __init__(self, ufl_expression: ufl.core.expr.Expr, X: np.ndarray,
                  form_compiler_options: dict = {}, jit_options: dict = {},
-                 dtype=PETSc.ScalarType):
+                 dtype=default_scalar_type):
         """Create DOLFINx Expression.
 
         Represents a mathematical expression evaluated at a pre-defined
@@ -240,7 +239,7 @@ class Function(ufl.Coefficient):
     """
 
     def __init__(self, V: FunctionSpace, x: typing.Optional[la.VectorMetaClass] = None,
-                 name: typing.Optional[str] = None, dtype: np.dtype = PETSc.ScalarType):
+                 name: typing.Optional[str] = None, dtype: np.dtype = default_scalar_type):
         """Initialize a finite element Function.
 
         Args:
@@ -322,7 +321,7 @@ def eval(self, x: npt.ArrayLike, cells: npt.ArrayLike, u=None) -> np.ndarray:
         # Allocate memory for return value if not provided
         if u is None:
             value_size = ufl.product(self.ufl_element().value_shape())
-            if np.issubdtype(PETSc.ScalarType, np.complexfloating):
+            if np.issubdtype(default_scalar_type, np.complexfloating):
                 u = np.empty((num_points, value_size), dtype=np.complex128)
             else:
                 u = np.empty((num_points, value_size))

python/test/unit/fem/test_assemble_domains.py

@@ -7,20 +7,20 @@
 
 import numpy as np
 import pytest
-
 import ufl
-from dolfinx import cpp as _cpp
 from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar,
                          dirichletbc, form)
 from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                                set_bc)
 from dolfinx.mesh import (GhostMode, Mesh, create_unit_square, locate_entities,
                           locate_entities_boundary, meshtags,
                           meshtags_from_entities)
-
 from mpi4py import MPI
 from petsc4py import PETSc
 
+from dolfinx import cpp as _cpp
+from dolfinx import default_scalar_type
+
 
 @pytest.fixture
 def mesh():
@@ -192,7 +192,7 @@ def right(x):
 def test_assembly_dS_domains(mode):
     N = 10
     mesh = create_unit_square(MPI.COMM_WORLD, N, N, ghost_mode=mode)
-    one = Constant(mesh, PETSc.ScalarType(1))
+    one = Constant(mesh, default_scalar_type(1))
     val = assemble_scalar(form(one * ufl.dS))
     val = mesh.comm.allreduce(val, op=MPI.SUM)
     assert val == pytest.approx(2 * (N - 1) + N * np.sqrt(2), 1.0e-7)