From 89dd2b5cd3757fdf61416308ef1626d8716f09a5 Mon Sep 17 00:00:00 2001 From: Philipp Diercks <55233917+pdiercks@users.noreply.github.com> Date: Mon, 23 Oct 2023 12:24:18 +0200 Subject: [PATCH] [hackathon] boundary conditions (#149) --- .../boundary_conditions/bcs.py | 8 +++---- .../boundary_conditions/boundary.py | 2 +- .../test_bcs_get_boundary_dofs.py | 7 ++++-- .../test_boundary_line_at.py | 2 +- .../test_boundary_plane_at.py | 4 ++-- .../test_boundary_point_at.py | 9 ++++--- .../test_boundary_within_range.py | 4 ++-- .../boundary_conditions/test_dirichlet_bcs.py | 24 ++++++++++++------- tests/boundary_conditions/test_neumann_bcs.py | 4 +++- 9 files changed, 39 insertions(+), 25 deletions(-) diff --git a/src/fenicsxconcrete/boundary_conditions/bcs.py b/src/fenicsxconcrete/boundary_conditions/bcs.py index 4f58c4e..a0d4db2 100644 --- a/src/fenicsxconcrete/boundary_conditions/bcs.py +++ b/src/fenicsxconcrete/boundary_conditions/bcs.py @@ -9,7 +9,7 @@ from fenicsxconcrete.util import LogMixin -def get_boundary_dofs(V: dolfinx.fem.FunctionSpace, marker: Callable) -> np.ndarray: +def get_boundary_dofs(V: dolfinx.fem.FunctionSpaceBase, marker: Callable) -> np.ndarray: """Returns dofs on the boundary specified by geometrical `marker`. Args: @@ -78,7 +78,7 @@ def __init__( def add_dirichlet_bc( self, value: ( - dolfinx.fem.Function | dolfinx.fem.Constant | dolfinx.fem.DirichletBCMetaClass | np.ndarray | Callable + dolfinx.fem.Function | dolfinx.fem.Constant | dolfinx.fem.DirichletBC | np.ndarray | Callable ), boundary: int | np.ndarray | Callable | None = None, sub: int = None, @@ -103,7 +103,7 @@ def add_dirichlet_bc( topologically. Note that `entity_dim` is required if `sub` is not None and `method=geometrical`. """ - if isinstance(value, dolfinx.fem.DirichletBCMetaClass): + if isinstance(value, dolfinx.fem.DirichletBC): self._bcs.append(value) else: assert method in ("topological", "geometrical") @@ -180,7 +180,7 @@ def has_dirichlet(self) -> bool: return len(self._bcs) > 0 @property - def bcs(self) -> list[dolfinx.fem.DirichletBCMetaClass]: + def bcs(self) -> list[dolfinx.fem.DirichletBC]: """returns the list of Dirichlet BCs Returns: diff --git a/src/fenicsxconcrete/boundary_conditions/boundary.py b/src/fenicsxconcrete/boundary_conditions/boundary.py index d2e826e..b5b46f0 100644 --- a/src/fenicsxconcrete/boundary_conditions/boundary.py +++ b/src/fenicsxconcrete/boundary_conditions/boundary.py @@ -193,7 +193,7 @@ def _show_marked( if tdim in (1, 3): raise NotImplementedError(f"Not implemented for mesh of topological dimension {tdim=}.") - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1)) + V = dolfinx.fem.functionspace(domain, ("Lagrange", 1)) dofs = dolfinx.fem.locate_dofs_geometrical(V, marker) u = dolfinx.fem.Function(V) bc = dolfinx.fem.dirichletbc(u, dofs) diff --git a/tests/boundary_conditions/test_bcs_get_boundary_dofs.py b/tests/boundary_conditions/test_bcs_get_boundary_dofs.py index 92c4428..14b993d 100644 --- a/tests/boundary_conditions/test_bcs_get_boundary_dofs.py +++ b/tests/boundary_conditions/test_bcs_get_boundary_dofs.py @@ -3,6 +3,7 @@ import dolfinx import numpy as np from mpi4py import MPI +from basix.ufl import element from fenicsxconcrete.boundary_conditions.bcs import BoundaryConditions, get_boundary_dofs from fenicsxconcrete.boundary_conditions.boundary import plane_at @@ -44,7 +45,8 @@ def test_whole_boundary() -> None: dim = 2 domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n, dolfinx.mesh.CellType.quadrilateral) - V = dolfinx.fem.VectorFunctionSpace(domain, ("Lagrange", degree), dim=dim) + el = element("Lagrange", domain.basix_cell(), degree, shape=(dim,)) + V = dolfinx.fem.functionspace(domain, el) # option (a) bc_handler = BoundaryConditions(domain, V) @@ -70,7 +72,8 @@ def test_xy_plane() -> None: dim = 3 domain = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, n, n, n, dolfinx.mesh.CellType.hexahedron) - V = dolfinx.fem.VectorFunctionSpace(domain, ("Lagrange", degree), dim=dim) + el = element("Lagrange", domain.basix_cell(), degree, shape=(dim,)) + V = dolfinx.fem.functionspace(domain, el) xy_plane = plane_at(0.0, "z") # option (a) diff --git a/tests/boundary_conditions/test_boundary_line_at.py b/tests/boundary_conditions/test_boundary_line_at.py index 4235fe3..0e07e45 100644 --- a/tests/boundary_conditions/test_boundary_line_at.py +++ b/tests/boundary_conditions/test_boundary_line_at.py @@ -11,7 +11,7 @@ def test_cube() -> None: n = 4 domain = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, n, n, n, dolfinx.mesh.CellType.hexahedron) - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1)) + V = dolfinx.fem.functionspace(domain, ("Lagrange", 1)) x_axis = line_at([0, 0], ["z", "y"]) y_axis = line_at([0, 0], ["z", "x"]) diff --git a/tests/boundary_conditions/test_boundary_plane_at.py b/tests/boundary_conditions/test_boundary_plane_at.py index 9214d28..9ea926d 100644 --- a/tests/boundary_conditions/test_boundary_plane_at.py +++ b/tests/boundary_conditions/test_boundary_plane_at.py @@ -10,7 +10,7 @@ def test_square() -> None: domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, dolfinx.mesh.CellType.quadrilateral) - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1)) + V = dolfinx.fem.functionspace(domain, ("Lagrange", 1)) bottom = plane_at(0.0, "y") top = plane_at(1.0, "y") @@ -56,7 +56,7 @@ def l_shaped_boundary(x): def test_cube() -> None: nx, ny, nz = 4, 4, 4 domain = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, nx, ny, nz, dolfinx.mesh.CellType.hexahedron) - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1)) + V = dolfinx.fem.functionspace(domain, ("Lagrange", 1)) xy_plane = plane_at(0.0, "z") dofs = dolfinx.fem.locate_dofs_geometrical(V, xy_plane) diff --git a/tests/boundary_conditions/test_boundary_point_at.py b/tests/boundary_conditions/test_boundary_point_at.py index 9ed2561..27c9306 100644 --- a/tests/boundary_conditions/test_boundary_point_at.py +++ b/tests/boundary_conditions/test_boundary_point_at.py @@ -4,6 +4,7 @@ import numpy as np from mpi4py import MPI from petsc4py.PETSc import ScalarType +from basix.ufl import element from fenicsxconcrete.boundary_conditions.boundary import point_at @@ -12,7 +13,7 @@ def test_type_error() -> None: """test TypeError in conversion to float""" n = 10 domain = dolfinx.mesh.create_interval(MPI.COMM_WORLD, n, [0.0, 10.0]) - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1)) + V = dolfinx.fem.functionspace(domain, ("Lagrange", 1)) x = point_at(5) dofs = dolfinx.fem.locate_dofs_geometrical(V, x) nodal_value = 42 @@ -25,7 +26,7 @@ def test_type_error() -> None: def test_function_space() -> None: n = 101 domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n, dolfinx.mesh.CellType.quadrilateral) - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 2)) + V = dolfinx.fem.functionspace(domain, ("Lagrange", 2)) h = 1.0 / n my_point = point_at(np.array([h * 2, h * 5])) @@ -40,7 +41,9 @@ def test_function_space() -> None: def test_vector_function_space() -> None: n = 101 domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n, dolfinx.mesh.CellType.quadrilateral) - V = dolfinx.fem.VectorFunctionSpace(domain, ("Lagrange", 2)) + degree = 2 + ve = element("Lagrange", domain.basix_cell(), degree, shape=(2,)) + V = dolfinx.fem.functionspace(domain, ve) # note the inconsistency in specifying the coordinates # this is handled by `to_floats` diff --git a/tests/boundary_conditions/test_boundary_within_range.py b/tests/boundary_conditions/test_boundary_within_range.py index 176bceb..7fcd3d5 100644 --- a/tests/boundary_conditions/test_boundary_within_range.py +++ b/tests/boundary_conditions/test_boundary_within_range.py @@ -10,7 +10,7 @@ def test_1d() -> None: n = 10 domain = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, n) - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1)) + V = dolfinx.fem.functionspace(domain, ("Lagrange", 1)) # test reordering of start and end subdomain = within_range([0.75, 0.0, 0.0], [0.35, 0.0, 0.0]) dofs = dolfinx.fem.locate_dofs_geometrical(V, subdomain) @@ -22,7 +22,7 @@ def test_1d() -> None: def test_2d() -> None: n = 200 domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n, dolfinx.mesh.CellType.quadrilateral) - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1)) + V = dolfinx.fem.functionspace(domain, ("Lagrange", 1)) Δx = Δy = 1 / (n + 1) # exclude the right and top boundary, Δ must be smaller than cell size boundary = within_range([0.0, 0.0, 0.0], [1.0 - Δx, 1.0 - Δy, 0.0]) diff --git a/tests/boundary_conditions/test_dirichlet_bcs.py b/tests/boundary_conditions/test_dirichlet_bcs.py index e201c71..b711c9d 100644 --- a/tests/boundary_conditions/test_dirichlet_bcs.py +++ b/tests/boundary_conditions/test_dirichlet_bcs.py @@ -44,7 +44,8 @@ def boundary(x): def test_vector_geom() -> None: domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, dolfinx.mesh.CellType.quadrilateral) - V = dolfinx.fem.VectorFunctionSpace(domain, ("Lagrange", 2)) + ve = element("Lagrange", domain.basix_cell(), 2, shape=(2,)) + V = dolfinx.fem.functionspace(domain, ve) bc_handler = BoundaryConditions(domain, V) @@ -81,7 +82,8 @@ def left(x): def test_vector_geom_component_wise() -> None: domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, dolfinx.mesh.CellType.quadrilateral) - V = dolfinx.fem.VectorFunctionSpace(domain, ("Lagrange", 2)) + ve = element("Lagrange", domain.basix_cell(), 2, shape=(2,)) + V = dolfinx.fem.functionspace(domain, ve) bc_handler = BoundaryConditions(domain, V) @@ -105,7 +107,8 @@ def left(x): def test_scalar_geom() -> None: domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 8, 8) - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 2)) + ve = element("Lagrange", domain.basix_cell(), 2, shape=(2,)) + V = dolfinx.fem.functionspace(domain, ve) bc_handler = BoundaryConditions(domain, V) @@ -126,7 +129,7 @@ def left(x): def test_scalar_topo() -> None: n = 20 domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n) - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 2)) + V = dolfinx.fem.functionspace(domain, ("Lagrange", 2)) bc_handler = BoundaryConditions(domain, V) @@ -151,7 +154,7 @@ def test_dirichletbc() -> None: """add instance of dolfinx.fem.dirichletbc""" n = 20 domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n) - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 2)) + V = dolfinx.fem.functionspace(domain, ("Lagrange", 2)) bc_handler = BoundaryConditions(domain, V) dofs = dolfinx.fem.locate_dofs_geometrical(V, plane_at(0.0, "x")) bc = dolfinx.fem.dirichletbc(ScalarType(0), dofs, V) @@ -165,7 +168,8 @@ def test_runtimeerror_geometrical() -> None: is not None""" n = 20 domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n) - V = dolfinx.fem.VectorFunctionSpace(domain, ("Lagrange", 2)) + ve = element("Lagrange", domain.basix_cell(), 2, shape=(2,)) + V = dolfinx.fem.functionspace(domain, ve) Vsub = V.sub(0) bottom = plane_at(0.0, "y") with pytest.raises(RuntimeError): @@ -175,7 +179,9 @@ def test_runtimeerror_geometrical() -> None: def test_boundary_as_int() -> None: n = 5 domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n) - V = dolfinx.fem.VectorFunctionSpace(domain, ("Lagrange", 2)) + ve = element("Lagrange", domain.basix_cell(), 2, shape=(2,)) + V = dolfinx.fem.functionspace(domain, ve) + marker = 1011 bottom = {"bottom": (marker, plane_at(0.0, "y"))} ft, marked = create_facet_tags(domain, bottom) @@ -196,7 +202,7 @@ def test_boundary_as_int() -> None: def test_value_interpolation() -> None: n = 50 domain = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, n) - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 2)) + V = dolfinx.fem.functionspace(domain, ("Lagrange", 2)) my_value = 17.2 def expression(x): @@ -215,7 +221,7 @@ def everywhere(x): def test_clear() -> None: n = 2 domain = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, n) - V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 2)) + V = dolfinx.fem.functionspace(domain, ("Lagrange", 2)) bc_handler = BoundaryConditions(domain, V) dofs = dolfinx.fem.locate_dofs_geometrical(V, plane_at(0.0, "x")) bc = dolfinx.fem.dirichletbc(ScalarType(0), dofs, V) diff --git a/tests/boundary_conditions/test_neumann_bcs.py b/tests/boundary_conditions/test_neumann_bcs.py index 1ba16f1..42003b7 100644 --- a/tests/boundary_conditions/test_neumann_bcs.py +++ b/tests/boundary_conditions/test_neumann_bcs.py @@ -3,6 +3,7 @@ import pytest from mpi4py import MPI from petsc4py import PETSc +from basix.ufl import element from fenicsxconcrete.boundary_conditions.bcs import BoundaryConditions from fenicsxconcrete.boundary_conditions.boundary import create_facet_tags, plane_at @@ -11,7 +12,8 @@ def test_constant_traction() -> None: n = 10 domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n, dolfinx.mesh.CellType.quadrilateral) - V = dolfinx.fem.VectorFunctionSpace(domain, ("Lagrange", 1)) + ve = element("Lagrange", domain.basix_cell(), degree, shape=(2,)) + V = dolfinx.fem.functionspace(domain, ve) rmarker = 12 my_boundaries = {"right": (rmarker, plane_at(0.0, "x"))} ft, mb = create_facet_tags(domain, my_boundaries)