diff --git a/environment.yml b/environment.yml index 3654617..d24ec83 100644 --- a/environment.yml +++ b/environment.yml @@ -1,10 +1,9 @@ name: fenics-constitutive channels: - conda-forge - - defaults dependencies: - python>=3.10 - - fenics-dolfinx=0.6.0 + - fenics-dolfinx>=0.9.0 - numpy - pip - pytest diff --git a/examples/linear_elasticity/test_elasticity.py b/examples/linear_elasticity/test_elasticity.py index 2d4ef4d..1913d02 100644 --- a/examples/linear_elasticity/test_elasticity.py +++ b/examples/linear_elasticity/test_elasticity.py @@ -13,6 +13,7 @@ IncrSmallStrainProblem, PlaneStrainFrom3D, UniaxialStrainFrom3D, + norm, ) youngs_modulus = 42.0 @@ -20,8 +21,8 @@ def test_uniaxial_stress(): - mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 4) - V = df.fem.FunctionSpace(mesh, ("CG", 1)) + mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 10) + V = df.fem.functionspace(mesh, ("CG", 1)) u = df.fem.Function(V) law = LinearElasticityModel( parameters={"E": youngs_modulus, "nu": poissons_ratio}, @@ -51,25 +52,25 @@ def right_boundary(x): n, converged = solver.solve(u) # Compare the result with the analytical solution - assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.01) < 1e-10 / ( - youngs_modulus * 0.01 - ) + diff = problem.stress_1 - ufl.as_vector([youngs_modulus * 0.01,]) + assert norm(diff, problem.dxm) < 1e-10 / youngs_modulus * 0.01 problem.update() # Check that the stress is updated correctly - assert abs(problem.stress_0.x.array[0] - youngs_modulus * 0.01) < 1e-10 / ( - youngs_modulus * 0.01 - ) + diff = problem.stress_0 - ufl.as_vector([youngs_modulus * 0.01,]) + assert norm(diff, problem.dxm) < 1e-10 / youngs_modulus * 0.01 + # Check that the displacement is updated correctly - assert np.max(problem._u0.x.array) == displacement.value + max_u = MPI.COMM_WORLD.allreduce(np.max(problem._u0.x.array), MPI.MAX) + assert max_u == displacement.value displacement.value = 0.02 n, converged = solver.solve(u) # Compare the result of the updated problem with new BC with the analytical solution - assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.02) < 1e-10 / ( - youngs_modulus * 0.02 - ) + diff = problem.stress_1 - ufl.as_vector([youngs_modulus * 0.02,]) + assert norm(diff, problem.dxm) < 1e-10 / youngs_modulus * 0.02 + @pytest.mark.parametrize( @@ -83,22 +84,23 @@ def right_boundary(x): ) def test_uniaxial_stress_two_laws(factor: float): mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2) - V = df.fem.FunctionSpace(mesh, ("CG", 1)) + V = df.fem.functionspace(mesh, ("CG", 1)) u = df.fem.Function(V) + cells_local = mesh.topology.index_map(mesh.topology.dim).global_to_local(np.arange(2, dtype=np.int32)) laws = [ ( LinearElasticityModel( parameters={"E": youngs_modulus, "nu": poissons_ratio}, constraint=Constraint.UNIAXIAL_STRESS, ), - np.array([0], dtype=np.int32), + cells_local[0:1], ), ( LinearElasticityModel( parameters={"E": factor * youngs_modulus, "nu": poissons_ratio}, constraint=Constraint.UNIAXIAL_STRESS, ), - np.array([1], dtype=np.int32), + cells_local[1:2], ), ] @@ -138,7 +140,7 @@ def right_boundary(x): def test_uniaxial_strain(): mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2) - V = df.fem.FunctionSpace(mesh, ("CG", 1)) + V = df.fem.functionspace(mesh, ("CG", 1)) u = df.fem.Function(V) law = LinearElasticityModel( parameters={"E": youngs_modulus, "nu": poissons_ratio}, @@ -173,10 +175,8 @@ def right_boundary(x): * (1.0 - poissons_ratio) / ((1.0 + poissons_ratio) * (1.0 - 2.0 * poissons_ratio)) ) * displacement.value - - assert abs(problem.stress_0.x.array[0] - analytical_stress) < 1e-10 / ( - analytical_stress - ) + diff = problem.stress_0 - ufl.as_vector([analytical_stress,]) + assert norm(diff, problem.dxm) < 1e-10 / analytical_stress # test the converter from 3D model to uniaxial strain model law_3d = LinearElasticityModel( @@ -196,32 +196,28 @@ def right_boundary(x): problem_3d.update() # test that sigma_11 is the same as the analytical solution - assert abs(problem_3d.stress_0.x.array[0] - analytical_stress) < 1e-10 / ( - analytical_stress - ) + diff = problem_3d.stress_0 - ufl.as_vector([analytical_stress,]) + assert norm(diff, problem_3d.dxm) < 1e-10 / analytical_stress + # test that the stresses of the problem with uniaxial strain model # are the same as the stresses of the problem with the converted 3D model - assert ( - np.linalg.norm(problem_3d.stress_0.x.array - problem.stress_0.x.array) - / np.linalg.norm(problem.stress_0.x.array) - < 1e-14 - ) + diff = problem_3d.stress_0 - problem.stress_0 + assert norm(diff, problem_3d.dxm) < 1e-10 / norm(problem.stress_0, problem.dxm) # test that the shear stresses are zero. Since this is uniaxial strain, the # stress can have nonzero \sigma_22 and \sigma_33 components assert np.linalg.norm(wrapped_1d_law.stress_3d[3:6]) < 1e-14 # test that the displacement is the same in both versions + diff = problem_3d._u - problem._u assert ( - np.linalg.norm(problem_3d._u.x.array - problem._u.x.array) - / np.linalg.norm(problem._u.x.array) - < 1e-14 + norm(diff, problem_3d.dxm) < 1e-14 / norm(problem._u, problem.dxm) ) def test_plane_strain(): # sanity check if out of plane stress is NOT zero mesh = df.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2) - V = df.fem.VectorFunctionSpace(mesh, ("CG", 1)) + V = df.fem.functionspace(mesh, ("CG", 1,(2,))) u = df.fem.Function(V) law = LinearElasticityModel( parameters={"E": youngs_modulus, "nu": poissons_ratio}, @@ -249,14 +245,9 @@ def right_boundary(x): solver = NewtonSolver(MPI.COMM_WORLD, problem) n, converged = solver.solve(u) problem.update() - assert ( - np.linalg.norm( - problem.stress_0.x.array.reshape(-1, law.constraint.stress_strain_dim())[ - :, 2 - ] - ) - > 1e-2 - ) + # test that the stress is nonzero in 33 direction + assert norm(problem.stress_0[2], problem.dxm) > 1e-2 + # test the model conversion from 3D to plane strain law_3d = LinearElasticityModel( parameters={"E": youngs_modulus, "nu": poissons_ratio}, @@ -274,32 +265,22 @@ def right_boundary(x): n, converged = solver_3d.solve(u_3d) problem_3d.update() # test that the stress is nonzero in 33 direction - assert ( - np.linalg.norm( - problem_3d.stress_0.x.array.reshape(-1, law.constraint.stress_strain_dim())[ - :, 2 - ] - ) - > 1e-2 - ) + assert norm(problem_3d.stress_0[2], problem.dxm) > 1e-2 + # test that the displacement is the same in both versions + diff = problem_3d._u - problem._u assert ( - np.linalg.norm(problem_3d._u.x.array - problem._u.x.array) - / np.linalg.norm(problem._u.x.array) - < 1e-14 + norm(diff, problem.dxm)/ norm(problem._u, problem.dxm) < 1e-14 ) # test that the stresses are the same in both versions - assert ( - np.linalg.norm(problem_3d.stress_0.x.array - problem.stress_0.x.array) - / np.linalg.norm(problem.stress_0.x.array) - < 1e-14 - ) + diff = problem_3d.stress_0 - problem.stress_0 + assert norm(diff, problem.dxm) / norm(problem.stress_0, problem.dxm)< 1e-10 def test_plane_stress(): # just a sanity check if out of plane stress is really zero mesh = df.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2) - V = df.fem.VectorFunctionSpace(mesh, ("CG", 1)) + V = df.fem.functionspace(mesh, ("CG", 1,(2,))) u = df.fem.Function(V) law = LinearElasticityModel( parameters={"E": youngs_modulus, "nu": poissons_ratio}, @@ -328,20 +309,15 @@ def right_boundary(x): solver = NewtonSolver(MPI.COMM_WORLD, problem) n, converged = solver.solve(u) problem.update() - assert ( - np.linalg.norm( - problem.stress_0.x.array.reshape(-1, law.constraint.stress_strain_dim())[ - :, 2 - ] - ) - < 1e-10 - ) + # test that the out of plane stress is zero + assert norm(problem.stress_0[2], problem.dxm) < 1e-10 + def test_3d(): # test the 3d case against a pure fenics implementation mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 2, 2, 2) - V = df.fem.VectorFunctionSpace(mesh, ("CG", 1)) + V = df.fem.functionspace(mesh, ("CG", 1,(3,))) u = df.fem.Function(V) law = LinearElasticityModel( parameters={"E": youngs_modulus, "nu": poissons_ratio}, @@ -395,10 +371,13 @@ def sigma(v): problem_fenics.solve() # Check that the solution is the same - assert np.linalg.norm(u_fenics.x.array - u.x.array) < 1e-8 / np.linalg.norm( - u_fenics.x.array - ) + diff = u_fenics - u + assert norm(diff, problem.dxm) < 1e-8 / norm(u, problem.dxm) if __name__ == "__main__": test_uniaxial_stress() + test_uniaxial_strain() + test_plane_stress() + test_plane_strain() + test_3d() diff --git a/examples/mises_plasticity/test_plasticity.py b/examples/mises_plasticity/test_plasticity.py index 0431cbc..4a01283 100644 --- a/examples/mises_plasticity/test_plasticity.py +++ b/examples/mises_plasticity/test_plasticity.py @@ -10,8 +10,9 @@ def test_uniaxial_strain_3d(): + #no MPI, one cell only mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 1, 1, 1) - V = df.fem.VectorFunctionSpace(mesh, ("CG", 1)) + V = df.fem.functionspace(mesh, ("CG", 1, (3,))) u = df.fem.Function(V) matparam = { "p_ka": 175000, @@ -49,8 +50,8 @@ def right(x): ) dirichlet = [fix_ux_left, move_ux_right] - # - problem = IncrSmallStrainProblem(law, u, dirichlet, q_degree=2) + + problem = IncrSmallStrainProblem(law, u, dirichlet, q_degree=1) solver = NewtonSolver(MPI.COMM_WORLD, problem) @@ -62,7 +63,7 @@ def right(x): load = [0.0] for inc, time in enumerate(load_steps): - print("Load Increment:", inc) + #print("Load Increment:", inc) current_disp = time * max_disp scalar_x.value = current_disp @@ -70,7 +71,7 @@ def right(x): niter, converged = solver.solve(u) problem.update() - print(f"Converged: {converged} in {niter} iterations.") + #print(f"Converged: {converged} in {niter} iterations.") iterations = np.append(iterations, niter) stress_values = [] @@ -106,7 +107,7 @@ def right(x): def test_uniaxial_cyclic_strain_3d(): mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 1, 1, 1) - V = df.fem.VectorFunctionSpace(mesh, ("CG", 1)) + V = df.fem.functionspace(mesh, ("CG", 1, (3,))) u = df.fem.Function(V) matparam = { "p_ka": 175000, @@ -145,7 +146,7 @@ def right(x): dirichlet = [fix_ux_left, move_ux_right] # - problem = IncrSmallStrainProblem(law, u, dirichlet, q_degree=2) + problem = IncrSmallStrainProblem(law, u, dirichlet, q_degree=1) solver = NewtonSolver(MPI.COMM_WORLD, problem) @@ -157,7 +158,7 @@ def right(x): load = [0.0] for inc, time in enumerate(load_steps): - print("Load Increment:", inc) + #print("Load Increment:", inc) current_disp = np.sin(time) * max_disp scalar_x.value = current_disp @@ -165,7 +166,7 @@ def right(x): niter, converged = solver.solve(u) problem.update() - print(f"Converged: {converged} in {niter} iterations.") + #print(f"Converged: {converged} in {niter} iterations.") iterations = np.append(iterations, niter) stress_values = [] diff --git a/examples/visco_elasticity/test_viscoelasticity.py b/examples/visco_elasticity/test_viscoelasticity.py index c579026..3063c2b 100644 --- a/examples/visco_elasticity/test_viscoelasticity.py +++ b/examples/visco_elasticity/test_viscoelasticity.py @@ -24,7 +24,7 @@ def test_relaxation_uniaxial_stress(mat: IncrSmallStrainModel): """stress relaxation under uniaxial tension test for 1D, displacement controlled""" mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2) - V = df.fem.FunctionSpace(mesh, ("CG", 1)) + V = df.fem.functionspace(mesh, ("CG", 1)) u = df.fem.Function(V) law = mat( parameters={"E0": youngs_modulus, "E1": visco_modulus, "tau": relaxation_time}, @@ -89,8 +89,8 @@ def right_boundary(x): strain.append(problem._history_1[0]["strain"].x.array[-1]) viscostrain.append(problem._history_1[0]["strain_visco"].x.array[-1]) - print("0", time[0], disp[0], stress[0], strain[0], viscostrain[0]) - print("end", time[-1], disp[-1], stress[-1], strain[-1], viscostrain[-1]) + #print("0", time[0], disp[0], stress[0], strain[0], viscostrain[0]) + #print("end", time[-1], disp[-1], stress[-1], strain[-1], viscostrain[-1]) # analytic solution if isinstance(law, SpringKelvinModel): # analytic solution for 1D Kelvin model @@ -158,7 +158,7 @@ def z_boundary(x): else: raise ValueError(f"Dimension {dim} not supported") - V = df.fem.VectorFunctionSpace(mesh, ("CG", 1)) + V = df.fem.functionspace(mesh, ("CG", 1, (dim,))) u = df.fem.Function(V) def left_boundary(x): @@ -288,7 +288,7 @@ def test_kelvin_vs_maxwell(): """1D test with uniaxial stress, compare Kelvin and Maxwell model.""" mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2) - V = df.fem.FunctionSpace(mesh, ("CG", 1)) + V = df.fem.functionspace(mesh, ("CG", 1)) u = df.fem.Function(V) # Kelvin material @@ -402,7 +402,7 @@ def z_boundary(x): else: raise ValueError(f"Dimension {dim} not supported") - V = df.fem.VectorFunctionSpace(mesh, ("CG", 1)) + V = df.fem.functionspace(mesh, ("CG", 1,(dim,))) u = df.fem.Function(V) def left_boundary(x): @@ -525,7 +525,7 @@ def y_boundary(x): def create_meshtags( domain: df.mesh.Mesh, entity_dim: int, markers: dict[str, tuple[int, Callable]] -) -> tuple[df.mesh.MeshTagsMetaClass, dict[str, int]]: +) -> tuple[df.mesh.MeshTags, dict[str, int]]: """Creates meshtags for the given markers. This code is part of the FEniCSx tutorial by Jørgen S. Dokken. @@ -605,7 +605,7 @@ def z1_boundary(x): ) fixed_vector = np.array([0.0, 0.0, 0.0]) - V = df.fem.VectorFunctionSpace(mesh, ("CG", 1)) + V = df.fem.functionspace(mesh, ("CG", 1,(dim,))) u = df.fem.Function(V) # boundaries diff --git a/src/fenics_constitutive/__init__.py b/src/fenics_constitutive/__init__.py index f8c51be..e22e4c2 100644 --- a/src/fenics_constitutive/__init__.py +++ b/src/fenics_constitutive/__init__.py @@ -11,5 +11,6 @@ from .maps import * from .solver import * from .stress_strain import * +from .error_estimation import * __all__ = ["__version__"] diff --git a/src/fenics_constitutive/error_estimation.py b/src/fenics_constitutive/error_estimation.py new file mode 100644 index 0000000..c3802e0 --- /dev/null +++ b/src/fenics_constitutive/error_estimation.py @@ -0,0 +1,8 @@ +import ufl +import dolfinx as df +from mpi4py import MPI +import numpy as np + +def norm(f, dx, comm=MPI.COMM_WORLD): + norm_squared = df.fem.assemble_scalar(df.fem.form(ufl.inner(f, f) * dx)) + return np.sqrt(comm.allreduce(norm_squared, op=MPI.SUM)) diff --git a/src/fenics_constitutive/maps.py b/src/fenics_constitutive/maps.py index 4b49c50..bd88ff1 100644 --- a/src/fenics_constitutive/maps.py +++ b/src/fenics_constitutive/maps.py @@ -5,6 +5,9 @@ import dolfinx as df import numpy as np +from functools import reduce +import operator + __all__ = ["SubSpaceMap", "build_subspace_map"] @@ -41,7 +44,7 @@ def map_to_parent(self, sub: df.fem.Function, parent: df.fem.Function) -> None: ), "Parent mesh does not match" assert sub.ufl_shape == parent.ufl_shape, "Shapes do not match" - size = sub.ufl_element().value_size() + size = reduce(operator.mul, sub.ufl_shape, 1) parent_array = parent.x.array.reshape(-1, size) sub_array = sub.x.array.reshape(-1, size) @@ -63,7 +66,7 @@ def map_to_child(self, parent: df.fem.Function, sub: df.fem.Function) -> None: ), "Parent mesh does not match" assert sub.ufl_shape == parent.ufl_shape, "Shapes do not match" - size = sub.ufl_element().value_size() + size = reduce(operator.mul, sub.ufl_shape, 1) parent_array = parent.x.array.reshape(-1, size) sub_array = sub.x.array.reshape(-1, size) @@ -93,18 +96,19 @@ def build_subspace_map( mesh = V.mesh submesh, cell_map, _, _ = df.mesh.create_submesh(mesh, mesh.topology.dim, cells) fe = V.ufl_element() - V_sub = df.fem.FunctionSpace(submesh, fe) + V_sub = df.fem.functionspace(submesh, fe) submesh = V_sub.mesh view_parent = [] view_child = [] - num_sub_cells = submesh.topology.index_map(submesh.topology.dim).size_local + map_c = submesh.topology.index_map(mesh.topology.dim) + num_sub_cells = map_c.size_local + map_c.num_ghosts for cell in range(num_sub_cells): view_child.append(V_sub.dofmap.cell_dofs(cell)) view_parent.append(V.dofmap.cell_dofs(cell_map[cell])) - if view_child: + if len(view_child) > 0: map = SubSpaceMap( parent=np.hstack(view_parent), child=np.hstack(view_child), diff --git a/src/fenics_constitutive/solver.py b/src/fenics_constitutive/solver.py index 486b372..32784dc 100644 --- a/src/fenics_constitutive/solver.py +++ b/src/fenics_constitutive/solver.py @@ -1,10 +1,12 @@ from __future__ import annotations import basix +import basix.ufl import dolfinx as df import numpy as np import ufl from petsc4py import PETSc +from dolfinx.fem.petsc import NonlinearProblem from .interfaces import IncrSmallStrainModel from .maps import SubSpaceMap, build_subspace_map @@ -30,29 +32,16 @@ def build_history( history = {} for key, value in law.history_dim.items(): - match value: - case int(): - Qh = ufl.VectorElement( - "Quadrature", - mesh.ufl_cell(), - q_degree, - quad_scheme="default", - dim=value, - ) - case tuple(): - Qh = ufl.TensorElement( - "Quadrature", - mesh.ufl_cell(), - q_degree, - quad_scheme="default", - shape=value, - ) - history_space = df.fem.FunctionSpace(mesh, Qh) + value_shape = (value,) if isinstance(value, int) else value + Q = basix.ufl.quadrature_element( + mesh.topology.cell_name(), value_shape=value_shape, degree=q_degree + ) + history_space = df.fem.functionspace(mesh, Q) history[key] = df.fem.Function(history_space) return history -class IncrSmallStrainProblem(df.fem.petsc.NonlinearProblem): +class IncrSmallStrainProblem(NonlinearProblem): """ A nonlinear problem for incremental small strain models. To be used with the dolfinx NewtonSolver. @@ -79,7 +68,7 @@ def __init__( self, laws: list[tuple[IncrSmallStrainModel, np.ndarray]] | IncrSmallStrainModel, u: df.fem.Function, - bcs: list[df.fem.DirichletBCMetaClass], + bcs: list[df.fem.DirichletBC], q_degree: int, del_t: float = 1.0, form_compiler_options: dict | None = None, @@ -97,34 +86,29 @@ def __init__( law[0].constraint == constraint for law in laws ), "All laws must have the same constraint" - gdim = mesh.ufl_cell().geometric_dimension() + gdim = mesh.geometry.dim assert ( constraint.geometric_dim() == gdim ), "Geometric dimension mismatch between mesh and laws" - QVe = ufl.VectorElement( - "Quadrature", - mesh.ufl_cell(), - q_degree, - quad_scheme="default", - dim=constraint.stress_strain_dim(), + QVe = basix.ufl.quadrature_element( + mesh.topology.cell_name(), + value_shape=(constraint.stress_strain_dim(),), + degree=q_degree, ) - QTe = ufl.TensorElement( - "Quadrature", - mesh.ufl_cell(), - q_degree, - quad_scheme="default", - shape=(constraint.stress_strain_dim(), constraint.stress_strain_dim()), + QTe = basix.ufl.quadrature_element( + mesh.topology.cell_name(), + value_shape=( + constraint.stress_strain_dim(), + constraint.stress_strain_dim(), + ), + degree=q_degree, ) - Q_grad_u_e = ufl.TensorElement( - "Quadrature", - mesh.ufl_cell(), - q_degree, - quad_scheme="default", - shape=(gdim, gdim), + Q_grad_u_e = basix.ufl.quadrature_element( + mesh.topology.cell_name(), value_shape=(gdim, gdim), degree=q_degree ) - QV = df.fem.FunctionSpace(mesh, QVe) - QT = df.fem.FunctionSpace(mesh, QTe) + QV = df.fem.functionspace(mesh, QVe) + QT = df.fem.functionspace(mesh, QTe) self.laws: list[tuple[IncrSmallStrainModel, np.ndarray]] = [] self.submesh_maps: list[SubSpaceMap] = [] @@ -154,11 +138,11 @@ def __init__( self._stress.append(df.fem.Function(QV_subspace)) # subspace for grad u - Q_grad_u_subspace = df.fem.FunctionSpace(submesh, Q_grad_u_e) + Q_grad_u_subspace = df.fem.functionspace(submesh, Q_grad_u_e) self._del_grad_u.append(df.fem.Function(Q_grad_u_subspace)) # subspace for tanget - QT_subspace = df.fem.FunctionSpace(submesh, QTe) + QT_subspace = df.fem.functionspace(submesh, QTe) self._tangent.append(df.fem.Function(QT_subspace)) # subspaces for history @@ -235,8 +219,8 @@ def form(self, x: PETSc.Vec) -> None: """ super().form(x) # this copies the data from the vector x to the function _u - x.copy(self._u.vector) - self._u.vector.ghostUpdate( + x.copy(self._u.x.petsc_vec) + self._u.x.petsc_vec.ghostUpdate( addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD ) # This assertion can fail, even if everything is correct. @@ -247,11 +231,16 @@ def form(self, x: PETSc.Vec) -> None: # if len(self.laws) > 1: for k, (law, cells) in enumerate(self.laws): - # TODO: test this! - # Replace with `self._del_grad_u[k].interpolate(...)` in 0.8.0 - self.del_grad_u_expr.eval( - cells, self._del_grad_u[k].x.array.reshape(cells.size, -1) + self._del_grad_u[k].interpolate( + self.del_grad_u_expr, + cells0=cells, + cells1=np.arange(cells.size, dtype=np.int32), ) + #self.del_grad_u_expr.eval( + # self._del_grad_u[k].function_space.mesh, + # cells, + # self._del_grad_u[k].x.array.reshape(cells.size, -1), + #) self._del_grad_u[k].x.scatter_forward() if len(self.laws) > 1: self.submesh_maps[k].map_to_child(self.stress_0, self._stress[k]) diff --git a/tests/test_conversions.py b/tests/test_conversions.py index ea63355..d92ba30 100644 --- a/tests/test_conversions.py +++ b/tests/test_conversions.py @@ -71,23 +71,23 @@ def lam(x): def lam(x): return x[0] * 0.1, x[1] * 0.2, x[2] * 0.3 - P1 = df.fem.VectorFunctionSpace(mesh, ("Lagrange", 1)) + P1 = df.fem.functionspace(mesh, ("CG", 1,(constraint.geometric_dim(),))) u = df.fem.Function(P1) grad_u_ufl = ufl.grad(u) mandel_strain_ufl = ufl_mandel_strain(u, constraint) u.interpolate(lam) points, weights = basix.make_quadrature( - basix.quadrature.string_to_type("default"), - basix.cell.string_to_type(mesh.ufl_cell().cellname()), + basix.CellType[mesh.ufl_cell().cellname()], 1, + basix.quadrature.string_to_type("default"), ) expr_grad_u = df.fem.Expression(grad_u_ufl, points) expr_mandel_strain = df.fem.Expression(mandel_strain_ufl, points) n_cells = mesh.topology.index_map(mesh.topology.dim).size_local - grad_u_array = expr_grad_u.eval(np.arange(n_cells, dtype=np.int32)).flatten() - strain_array = expr_mandel_strain.eval(np.arange(n_cells, dtype=np.int32)).flatten() + grad_u_array = expr_grad_u.eval(mesh, np.arange(n_cells, dtype=np.int32)).flatten() + strain_array = expr_mandel_strain.eval(mesh, np.arange(n_cells, dtype=np.int32)).flatten() strain_array_from_grad_u = strain_from_grad_u(grad_u_array, constraint) assert np.allclose(strain_array, strain_array_from_grad_u) diff --git a/tests/test_maps.py b/tests/test_maps.py index 0ab0c4e..d19b25d 100644 --- a/tests/test_maps.py +++ b/tests/test_maps.py @@ -4,6 +4,7 @@ import numpy as np import ufl from mpi4py import MPI +import basix from fenics_constitutive import build_subspace_map @@ -12,61 +13,67 @@ def test_subspace_vector_map_vector_equals_tensor_map(): mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 5, 7, 11) map_c = mesh.topology.index_map(mesh.topology.dim) - num_cells = map_c.size_local + map_c.num_ghosts - cells = np.arange(0, num_cells, dtype=np.int32) + num_cells = map_c.size_local + size_global = map_c.size_global + cells_global = np.arange(0, size_global, dtype=np.int32) + #cells = np.arange(0, num_cells, dtype=np.int32) - Q = ufl.FiniteElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default") - QV = ufl.VectorElement( - "Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", dim=3 + Q = basix.ufl.quadrature_element( + mesh.topology.cell_name(), value_shape=(1,), degree=2 ) - QT = ufl.TensorElement( - "Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", shape=(3, 3) + QV = basix.ufl.quadrature_element( + mesh.topology.cell_name(), value_shape=(3,), degree=2 ) + QT = basix.ufl.quadrature_element( + mesh.topology.cell_name(), value_shape=(3, 3), degree=2 + ) + + Q_space = df.fem.functionspace(mesh, Q) + QV_space = df.fem.functionspace(mesh, QV) + QT_space = df.fem.functionspace(mesh, QT) - Q_space = df.fem.FunctionSpace(mesh, Q) - QV_space = df.fem.FunctionSpace(mesh, QV) - QT_space = df.fem.FunctionSpace(mesh, QT) + rng = np.random.default_rng(42) + if MPI.COMM_WORLD.rank == 0: + cell_sample_global = rng.choice(cells_global, size_global // 2, replace=False) + else: + cell_sample_global = np.empty(size_global // 2, dtype=np.int32) + MPI.COMM_WORLD.Bcast(cell_sample_global, root=0) + + cell_sample = map_c.global_to_local(cell_sample_global) + cell_sample = cell_sample[cell_sample>=0] - Q_map, _ = build_subspace_map(cells, Q_space) - QV_map, _ = build_subspace_map(cells, QV_space) - QT_map, _ = build_subspace_map(cells, QT_space) + Q_map, _ = build_subspace_map(cell_sample, Q_space) + QV_map, _ = build_subspace_map(cell_sample, QV_space) + QT_map, _ = build_subspace_map(cell_sample, QT_space) assert np.all(Q_map.parent == QV_map.parent) assert np.all(Q_map.child == QV_map.child) assert np.all(Q_map.parent == QT_map.parent) assert np.all(Q_map.child == QT_map.child) - for _ in range(10): - cell_sample = np.random.permutation(cells)[: num_cells // 2] - - Q_map, _ = build_subspace_map(cell_sample, Q_space) - QV_map, _ = build_subspace_map(cell_sample, QV_space) - QT_map, _ = build_subspace_map(cell_sample, QT_space) - - assert np.all(Q_map.parent == QV_map.parent) - assert np.all(Q_map.child == QV_map.child) - assert np.all(Q_map.parent == QT_map.parent) - assert np.all(Q_map.child == QT_map.child) - def test_map_evaluation(): - mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 5, 7, 11) + mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 5, 4, 1, cell_type=df.mesh.CellType.hexahedron) map_c = mesh.topology.index_map(mesh.topology.dim) - num_cells = map_c.size_local + map_c.num_ghosts - cells = np.arange(0, num_cells, dtype=np.int32) + num_cells_total = map_c.size_local + map_c.num_ghosts + + size_global = map_c.size_global + cells_global = np.arange(0, size_global, dtype=np.int32) - Q = ufl.FiniteElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default") - QV = ufl.VectorElement( - "Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", dim=3 + Q = basix.ufl.quadrature_element( + mesh.topology.cell_name(), value_shape=(1,), degree=1 ) - QT = ufl.TensorElement( - "Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", shape=(3, 3) + QV = basix.ufl.quadrature_element( + mesh.topology.cell_name(), value_shape=(3,), degree=1 + ) + QT = basix.ufl.quadrature_element( + mesh.topology.cell_name(), value_shape=(3, 3), degree=1 ) - Q_space = df.fem.FunctionSpace(mesh, Q) - QV_space = df.fem.FunctionSpace(mesh, QV) - QT_space = df.fem.FunctionSpace(mesh, QT) + Q_space = df.fem.functionspace(mesh, Q) + QV_space = df.fem.functionspace(mesh, QV) + QT_space = df.fem.functionspace(mesh, QT) q = df.fem.Function(Q_space) q_test = q.copy() @@ -75,22 +82,41 @@ def test_map_evaluation(): qt = df.fem.Function(QT_space) qt_test = qt.copy() - scalar_array = np.random.random(q.x.array.shape) + rng = np.random.default_rng(42) + scalar_array = rng.random(q.x.array.shape) + + q.x.array[:] = scalar_array - vector_array = np.random.random(qv.x.array.shape) + q.x.scatter_forward() + scalar_array = q.x.array + + vector_array = rng.random(qv.x.array.shape) qv.x.array[:] = vector_array - tensor_array = np.random.random(qt.x.array.shape) + qv.x.scatter_forward() + vector_array = qv.x.array + + tensor_array = rng.random(qt.x.array.shape) qt.x.array[:] = tensor_array + qt.x.scatter_forward() + tensor_array = qt.x.array for _ in range(10): - cell_sample = np.random.permutation(cells)[: num_cells // 2] + if MPI.COMM_WORLD.rank == 0: + cell_sample_global = rng.choice(cells_global, size_global // 2, replace=False) + else: + cell_sample_global = np.empty(size_global // 2, dtype=np.int32) + MPI.COMM_WORLD.Bcast(cell_sample_global, root=0) + + cell_sample = map_c.global_to_local(cell_sample_global) + cell_sample = cell_sample[cell_sample>=0] Q_sub_map, submesh = build_subspace_map( cell_sample, Q_space, return_subspace=False ) - Q_sub = df.fem.FunctionSpace(submesh, Q) - QV_sub = df.fem.FunctionSpace(submesh, QV) - QT_sub = df.fem.FunctionSpace(submesh, QT) + + Q_sub = df.fem.functionspace(submesh, Q) + QV_sub = df.fem.functionspace(submesh, QV) + QT_sub = df.fem.functionspace(submesh, QT) q_sub = df.fem.Function(Q_sub) qv_sub = df.fem.Function(QV_sub) @@ -104,14 +130,18 @@ def test_map_evaluation(): Q_sub_map.map_to_parent(qv_sub, qv_test) Q_sub_map.map_to_parent(qt_sub, qt_test) - q_view = scalar_array.reshape(num_cells, -1) - q_test_view = q_test.x.array.reshape(num_cells, -1) + q_view = scalar_array.reshape(num_cells_total, -1) + q_test_view = q_test.x.array.reshape(num_cells_total, -1) assert np.all(q_view[cell_sample] == q_test_view[cell_sample]) - qv_view = vector_array.reshape(num_cells, -1) - qv_test_view = qv_test.x.array.reshape(num_cells, -1) + qv_view = vector_array.reshape(num_cells_total, -1) + qv_test_view = qv_test.x.array.reshape(num_cells_total, -1) assert np.all(qv_view[cell_sample] == qv_test_view[cell_sample]) - qt_view = tensor_array.reshape(num_cells, -1) - qt_test_view = qt_test.x.array.reshape(num_cells, -1) + qt_view = tensor_array.reshape(num_cells_total, -1) + qt_test_view = qt_test.x.array.reshape(num_cells_total, -1) assert np.all(qt_view[cell_sample] == qt_test_view[cell_sample]) + +if __name__ == "__main__": + test_subspace_vector_map_vector_equals_tensor_map() + test_map_evaluation() \ No newline at end of file