diff --git a/src/bvpy/utils/post_processing.py b/src/bvpy/utils/post_processing.py index c8bd88f646432adeb3e172540654191ba7cffa9f..34a36314125f4755fdaab69380b2730433521800 100644 --- a/src/bvpy/utils/post_processing.py +++ b/src/bvpy/utils/post_processing.py @@ -65,9 +65,10 @@ class SolutionExplorer(object): self.shape = self._fenics.ufl_shape self.geometric_dimension = self._fenics.geometric_dimension() - self._dof_per_nodes = np.prod(self.shape) - if self._dof_per_nodes == 0: - self._dof_per_nodes = int(1) # for scalar Function + if len(self.shape) == 0: + self._dof_per_nodes = 1 + else: + self._dof_per_nodes = np.prod(self.shape) self.nb_dofs = len(self._functionSpace.dofmap().dofs()) self.nb_nodes = int(self.nb_dofs/self._dof_per_nodes) @@ -78,6 +79,10 @@ class SolutionExplorer(object): self._type = 'vector' elif isinstance(self._functionSpace.ufl_element(), fe.TensorElement): self._type = 'tensor' + elif isinstance(self._functionSpace.ufl_element(), fe.MixedElement): + raise ValueError("MixedElement is not supported yet.") + else: + raise ValueError("Unknown element type.") def __call__(self, *x): return self._fenics(*x) @@ -93,23 +98,35 @@ class SolutionExplorer(object): """ return self._mesh.coordinates() - def get_dofs_positions(self): - """Generates a Fenics object containing the dofs positions. + def get_dofs_positions(self, remove_duplicates=True): + """ + Retrieves the degree of freedom (DoF) positions from the function space. + + Parameters + ---------- + remove_duplicates : bool, optional + If True, returns unique DoF positions by removing duplicates. This is + useful for vector- or tensor-valued functions where the same position + vector appears multiple times (default is True). Returns ------- numpy.ndarray - The positions vectors of the dofs. - + A 2D array containing the positions of the degrees of freedom (DoFs). + Each row in the array represents a position vector of a DoF. + Notes ----- - Depending on the type of function (scalar, vector, tensor), - the shape of the returned array might vary and the position - array might contain duplicated vectors. For instance, for a - 3D vector-valued function, each position vector will be repeated - 3 times. + - For scalar-valued functions, the positions correspond to unique node coordinates. + - For vector- or tensor-valued functions, the positions may contain duplicated + coordinates unless `remove_duplicates` is set to True. For example, in a 3D + vector-valued function, each spatial position will appear multiple times + (once for each vector component). """ - return self._functionSpace.tabulate_dof_coordinates() + if remove_duplicates: + return self._functionSpace.tabulate_dof_coordinates()[::self._dof_per_nodes, :] + else: + return self._functionSpace.tabulate_dof_coordinates() def get_dofs_values(self): """Return an array of the value of the Function at the dofs. @@ -121,16 +138,14 @@ class SolutionExplorer(object): """ values_on_cells = self._fenics.vector().get_local() - - Nc = self._mesh.num_cells() - D = self.geometric_dimension if self._type == 'scalar': return values_on_cells elif self._type == 'vector': - return values_on_cells.reshape((Nc, D)) + return values_on_cells.reshape((self.nb_nodes, *self.shape)) elif self._type == 'tensor': - return values_on_cells.reshape((Nc, D, D)) + return values_on_cells.reshape((self.nb_nodes, *self.shape)) + def to_fenics(self): """Returns a fenics.Function from the considered function. diff --git a/test/test_post_process.py b/test/test_post_process.py index 905de60ebbba932ea410c634eb8ff80a28de5b43..018f6b3aaa65d48a0748976ae2feba30a57603db 100644 --- a/test/test_post_process.py +++ b/test/test_post_process.py @@ -1,6 +1,7 @@ import unittest import numpy as np +import ufl from bvpy.utils.post_processing import * from bvpy import BVP @@ -28,7 +29,7 @@ class TestPreProcess(unittest.TestCase): solexplo = SolutionExplorer(self.sol) solexplo.get_dofs_positions() solexplo.to_vertex_values() - # solexplo.to_dataframe() + solexplo.get_dofs_values() def test_sol_explorer_scal(self): c = fe.Constant(1) @@ -39,11 +40,22 @@ class TestPreProcess(unittest.TestCase): self.assertEqual(m.to_vertex_values()[mask], 0) def test_sol_explorer_vec(self): + # - Order 1 VV = fe.VectorFunctionSpace(self.mesh, 'P', 1) cc = fe.Constant([1,0]) mm = SolutionExplorer(fe.project(cc, VV)) np.testing.assert_almost_equal(mm.to_vertex_values()[0], [1,0]) + # - Order 2: assert that the number of dofs changed + VV = fe.VectorFunctionSpace(self.mesh, 'P', 2) + cc = fe.Constant([1,0]) + mm = SolutionExplorer(fe.project(cc, VV)) + self.assertEqual(mm.get_dofs_values().shape, mm.get_dofs_positions().shape) + self.assertEqual( + mm.get_dofs_positions(remove_duplicates=False).shape[0], + 2 * mm.get_dofs_positions(remove_duplicates=True).shape[0] + ) + def test_sol_explorer_tensor(self): VVV = fe.TensorFunctionSpace(self.mesh, 'P', 1, shape=(3,3)) ccc= fe.Constant(np.diag([1,1,1])) @@ -51,6 +63,16 @@ class TestPreProcess(unittest.TestCase): mmm = SolutionExplorer(uuu) np.testing.assert_almost_equal(mmm.to_vertex_values()[0], np.diag([1,1,1])) + def test_sol_explorer_mixed_error(self): + P = ufl.FiniteElement("Lagrange", self.mesh.ufl_cell(), 2) + U = ufl.VectorElement("Lagrange", self.mesh.ufl_cell(), 2) + element = ufl.MixedElement([P, U]) + PU = FunctionSpace(self.mesh, element) + cc = fe.Constant(np.array([1, 0, 2])) + + with self.assertRaises(ValueError, msg="MixedElement is not supported yet."): + SolutionExplorer(fe.project(cc, PU)) + def test_sol_explorer_geom_filter(self): c = fe.Constant(1) u = fe.project(c, self.V)