Mentions légales du service

Skip to content
Snippets Groups Projects
Commit e9913e74 authored by PETIT Manuel's avatar PETIT Manuel
Browse files

Fix issues related to the `get_dofs_values` reshaping. Make slight improvement...

Fix issues related to the `get_dofs_values` reshaping. Make slight improvement on Dofs gestion (_dof_per_nodes, dofs_positions).

Refactored `SolutionExplorer` to improve clarity and consistency in DoF computations, including changes to `_dof_per_nodes` logic and added support for position deduplication. Extended the `get_dofs_positions` method to handle duplicate removal and updated reshape logic for vector and tensor types in `get_dofs_values`. Added error handling for unsupported MixedElements and corresponding unit tests for better coverage.
parent 0672f229
No related branches found
No related tags found
No related merge requests found
Pipeline #1099649 passed with warnings
...@@ -65,9 +65,10 @@ class SolutionExplorer(object): ...@@ -65,9 +65,10 @@ class SolutionExplorer(object):
self.shape = self._fenics.ufl_shape self.shape = self._fenics.ufl_shape
self.geometric_dimension = self._fenics.geometric_dimension() self.geometric_dimension = self._fenics.geometric_dimension()
self._dof_per_nodes = np.prod(self.shape) if len(self.shape) == 0:
if self._dof_per_nodes == 0: self._dof_per_nodes = 1
self._dof_per_nodes = int(1) # for scalar Function else:
self._dof_per_nodes = np.prod(self.shape)
self.nb_dofs = len(self._functionSpace.dofmap().dofs()) self.nb_dofs = len(self._functionSpace.dofmap().dofs())
self.nb_nodes = int(self.nb_dofs/self._dof_per_nodes) self.nb_nodes = int(self.nb_dofs/self._dof_per_nodes)
...@@ -78,6 +79,10 @@ class SolutionExplorer(object): ...@@ -78,6 +79,10 @@ class SolutionExplorer(object):
self._type = 'vector' self._type = 'vector'
elif isinstance(self._functionSpace.ufl_element(), fe.TensorElement): elif isinstance(self._functionSpace.ufl_element(), fe.TensorElement):
self._type = 'tensor' 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): def __call__(self, *x):
return self._fenics(*x) return self._fenics(*x)
...@@ -93,23 +98,35 @@ class SolutionExplorer(object): ...@@ -93,23 +98,35 @@ class SolutionExplorer(object):
""" """
return self._mesh.coordinates() return self._mesh.coordinates()
def get_dofs_positions(self): def get_dofs_positions(self, remove_duplicates=True):
"""Generates a Fenics object containing the dofs positions. """
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 Returns
------- -------
numpy.ndarray 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 Notes
----- -----
Depending on the type of function (scalar, vector, tensor), - For scalar-valued functions, the positions correspond to unique node coordinates.
the shape of the returned array might vary and the position - For vector- or tensor-valued functions, the positions may contain duplicated
array might contain duplicated vectors. For instance, for a coordinates unless `remove_duplicates` is set to True. For example, in a 3D
3D vector-valued function, each position vector will be repeated vector-valued function, each spatial position will appear multiple times
3 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): def get_dofs_values(self):
"""Return an array of the value of the Function at the dofs. """Return an array of the value of the Function at the dofs.
...@@ -121,16 +138,14 @@ class SolutionExplorer(object): ...@@ -121,16 +138,14 @@ class SolutionExplorer(object):
""" """
values_on_cells = self._fenics.vector().get_local() values_on_cells = self._fenics.vector().get_local()
Nc = self._mesh.num_cells()
D = self.geometric_dimension
if self._type == 'scalar': if self._type == 'scalar':
return values_on_cells return values_on_cells
elif self._type == 'vector': 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': 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): def to_fenics(self):
"""Returns a fenics.Function from the considered function. """Returns a fenics.Function from the considered function.
......
import unittest import unittest
import numpy as np import numpy as np
import ufl
from bvpy.utils.post_processing import * from bvpy.utils.post_processing import *
from bvpy import BVP from bvpy import BVP
...@@ -28,7 +29,7 @@ class TestPreProcess(unittest.TestCase): ...@@ -28,7 +29,7 @@ class TestPreProcess(unittest.TestCase):
solexplo = SolutionExplorer(self.sol) solexplo = SolutionExplorer(self.sol)
solexplo.get_dofs_positions() solexplo.get_dofs_positions()
solexplo.to_vertex_values() solexplo.to_vertex_values()
# solexplo.to_dataframe() solexplo.get_dofs_values()
def test_sol_explorer_scal(self): def test_sol_explorer_scal(self):
c = fe.Constant(1) c = fe.Constant(1)
...@@ -39,11 +40,22 @@ class TestPreProcess(unittest.TestCase): ...@@ -39,11 +40,22 @@ class TestPreProcess(unittest.TestCase):
self.assertEqual(m.to_vertex_values()[mask], 0) self.assertEqual(m.to_vertex_values()[mask], 0)
def test_sol_explorer_vec(self): def test_sol_explorer_vec(self):
# - Order 1
VV = fe.VectorFunctionSpace(self.mesh, 'P', 1) VV = fe.VectorFunctionSpace(self.mesh, 'P', 1)
cc = fe.Constant([1,0]) cc = fe.Constant([1,0])
mm = SolutionExplorer(fe.project(cc, VV)) mm = SolutionExplorer(fe.project(cc, VV))
np.testing.assert_almost_equal(mm.to_vertex_values()[0], [1,0]) 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): def test_sol_explorer_tensor(self):
VVV = fe.TensorFunctionSpace(self.mesh, 'P', 1, shape=(3,3)) VVV = fe.TensorFunctionSpace(self.mesh, 'P', 1, shape=(3,3))
ccc= fe.Constant(np.diag([1,1,1])) ccc= fe.Constant(np.diag([1,1,1]))
...@@ -51,6 +63,16 @@ class TestPreProcess(unittest.TestCase): ...@@ -51,6 +63,16 @@ class TestPreProcess(unittest.TestCase):
mmm = SolutionExplorer(uuu) mmm = SolutionExplorer(uuu)
np.testing.assert_almost_equal(mmm.to_vertex_values()[0], np.diag([1,1,1])) 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): def test_sol_explorer_geom_filter(self):
c = fe.Constant(1) c = fe.Constant(1)
u = fe.project(c, self.V) u = fe.project(c, self.V)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment