diff --git a/src/bvpy/vforms/plasticity.py b/src/bvpy/vforms/plasticity.py
index a777296d711656b75726493099a14ddeb01c9f90..22d54d092ac61f5628593ce2f027da07c6bdc7f2 100644
--- a/src/bvpy/vforms/plasticity.py
+++ b/src/bvpy/vforms/plasticity.py
@@ -313,11 +313,6 @@ class AbstractGrowth(ABC):
def get_growth_increment(self):
"""
Abstract method to compute the increment of growth to add to the previous growth.
-
- Returns
- -------
- growth_increment : :class:`dolfin.cpp.function.Function`
- The computed growth increment.
"""
pass
@@ -346,7 +341,7 @@ class AbstractGrowth(ABC):
self.prev_solution = sol # previous solution
# Compute the growth increment
- self.growth_incr = self.get_growth_increment()
+ self.get_growth_increment()
# Update the previous growth with the growth increment
Id = fe.Identity(sol.geometric_dimension())
@@ -511,10 +506,10 @@ class ConstantGrowth(AbstractGrowth):
Returns
-------
- fe.Constant
- The growth increment tensor for the current time step.
+ None
+ Update `self.growth_incr` with the computed growth increment.
"""
- return fe.Constant(self.dg * self.dt)
+ self.growth_incr = fe.Constant(self.dg * self.dt)
class TimeDependentGrowth(AbstractGrowth):
@@ -609,12 +604,12 @@ class TimeDependentGrowth(AbstractGrowth):
Returns
-------
- fe.Constant
- The growth increment tensor for the current time interval `[t, t + dt]`.
+ None
+ Update `self.growth_incr` with the computed growth increment.
"""
# Integrate the growth rate function on [t, t+dt]
dg_dt, err = quad_vec(self.growth_func, self.t, self.t + self.dt, args=self.args, **self.kwargs)
- return fe.Constant(dg_dt)
+ self.growth_incr = fe.Constant(dg_dt)
class StrainDependentGrowth(AbstractGrowth):
@@ -707,8 +702,8 @@ class StrainDependentGrowth(AbstractGrowth):
Returns
-------
- fe.Function
- The growth increment tensor
+ None
+ Update `self.growth_incr` with the computed growth increment.
Notes
-----
@@ -723,7 +718,7 @@ class StrainDependentGrowth(AbstractGrowth):
# Compute the increment of growth assuming an extensibility parameter
if self.target_strain is None:
- return fe.Constant(self.extensibility) * strain * fe.Constant(self.dt)
+ self.growth_incr = fe.Constant(self.extensibility) * strain * fe.Constant(self.dt)
else:
# - Convert the current strain in numpy array
strain_obj = SolutionExplorer(strain)
@@ -743,7 +738,7 @@ class StrainDependentGrowth(AbstractGrowth):
strain_filtered.vector().set_local(strain_np_filtered.flatten())
# - Compute the increment of growth
- return fe.Constant(self.extensibility) * strain_filtered * fe.Constant(self.dt)
+ self.growth_incr = fe.Constant(self.extensibility) * strain_filtered * fe.Constant(self.dt)
class HeterogeneousGrowth(AbstractGrowth):
@@ -806,8 +801,8 @@ class HeterogeneousGrowth(AbstractGrowth):
Returns
-------
- growth_incr : fe.Function
- A FEniCS `Function` representing the final growth increment over the entire domain.
+ None
+ Update `self.growth_incr` with the computed growth increment.
"""
mesh = self.prev_solution.function_space().mesh()
@@ -830,14 +825,17 @@ class HeterogeneousGrowth(AbstractGrowth):
sub_prev_solution.interpolate(self.prev_solution)
# Compute the growth increment for the subdomain
- growth_class._prev_solution = sub_prev_solution
- sub_growth_incr = growth_class.get_growth_increment()
+ growth_class.prev_solution = sub_prev_solution
+ growth_class.get_growth_increment()
# Assert that the growth increment from the subdomain is a fe.Function
+ sub_growth_incr = growth_class.growth_incr
if not isinstance(sub_growth_incr, fe.Function):
sub_func_space = fe.FunctionSpace(sub_mesh, func_space.ufl_element())
sub_growth_incr = fe.project(sub_growth_incr, sub_func_space)
-
+ else:
+ sub_func_space = sub_growth_incr.function_space()
+ print(sub_growth_incr.compute_vertex_values())
# Assign the values computed on the subdomain to the growth increment on the whole domain
cell_idx = sub_mesh.data().array("parent_cell_indices", sub_mesh.topology().dim())
sub_dofmap = sub_func_space.dofmap()
@@ -856,7 +854,7 @@ class HeterogeneousGrowth(AbstractGrowth):
growth_incr.vector().set_local(vals)
growth_incr.vector().apply("insert")
- return growth_incr
+ self.growth_incr = growth_incr
def _propagate_property(self, prop_name, value):
"""
diff --git a/test/test_growth.py b/test/test_growth.py
index 44be5c6310af78881fd95ddc04d59c22b5e81171..ba8487e71258bd3612f14cfb8090ace3f9a81235 100644
--- a/test/test_growth.py
+++ b/test/test_growth.py
@@ -1,229 +1,251 @@
+import unittest
import numpy as np
import fenics as fe
-import pytest
-from bvpy.vforms.plasticity import MorphoElasticForm
-from bvpy.vforms.plasticity import TimeDependentGrowth, ConstantGrowth, StrainDependentGrowth, HeterogeneousGrowth
+from unittest.mock import patch
+from numpy.testing import assert_allclose
+
+from bvpy.vforms.plasticity import (
+ MorphoElasticForm,
+ TimeDependentGrowth,
+ ConstantGrowth,
+ StrainDependentGrowth,
+ HeterogeneousGrowth
+)
from bvpy.vforms.elasticity import StVenantKirchoffPotential
-@pytest.fixture
-def constant_term():
- # Example constant term for the constant growth model
- return np.array([[0.1, 0.0], [0.0, 0.1]])
-
-@pytest.fixture
-def growth_func():
- # Given a function that returns a 2x2 array
- def _growth_func(t, a, b):
- dg = np.zeros((2, 2))
- dg[0, 0] = b * np.sin(t * np.pi) + a
- return dg
-
- return _growth_func
-
-@pytest.fixture
-def fenics_setup():
- # Create a simple mesh and function space in FEniCS for testing
- mesh = fe.UnitSquareMesh(1, 1)
- V = fe.VectorFunctionSpace(mesh, 'P', 1)
- return mesh, V
-
-@pytest.fixture
-def previous_solution(fenics_setup):
- # Create a mock function representing the previous solution
- _, V = fenics_setup
- prev_sol = fe.Function(V)
- prev_sol.assign(fe.Constant([1.0, 1.0])) # Assign a constant value for simplicity
- vector = prev_sol.vector()
-
- # Assign different values to some DOFs (to create non-null strain)
- dof_indices = [0, 1, 2, 3, 4]
- new_values = [0.5, 0.8, 1.5, 1.2, 0.9]
-
- for i, value in zip(dof_indices, new_values):
- vector[i] = value
-
- return prev_sol
-
-@pytest.fixture
-def plastic_vform():
- potential_energy = StVenantKirchoffPotential()
- return MorphoElasticForm(potential_energy=potential_energy, F_g=np.eye(2), source = [0, 0])
-
-@pytest.fixture
-def growth_time_model(growth_func):
- a = 1.0
- b = 0.5
- return TimeDependentGrowth(growth_func, a, b)
-
-@pytest.fixture
-def growth_constant_model(constant_term):
- return ConstantGrowth(constant_term)
-
-@pytest.fixture
-def growth_strain_model(target_strain, extensibility, previous_solution):
- model = StrainDependentGrowth(target_strain=target_strain, extensibility=extensibility)
- model.prev_solution = previous_solution
- return model
-
-@pytest.fixture
-def target_strain():
- # Example target strain tensor for the strain-dependent growth model
- return np.array([[0.1, 0.0], [0.0, 0.1]])
-
-
-@pytest.fixture
-def extensibility():
- # Example extensibility factor for the strain-dependent growth model
- return 0.5
-
-@pytest.fixture
-def growth_heterogeneous_model(growth_constant_model, growth_strain_model, previous_solution):
- values = [1, 2, 1, 2]
- mesh = previous_solution.function_space().mesh()
- cdata = fe.MeshFunction("size_t", mesh, 2)
-
- # Assign value from cdata list to the MeshFunction's underlying array
- for cell, value in zip(fe.cells(mesh), values):
- cdata[cell.index()] = value
-
- growth_classes = {1: growth_constant_model, 2: growth_strain_model}
- model = HeterogeneousGrowth(growth_classes=growth_classes, cdata=cdata)
- return model
-
-def test_constant_growth_init(growth_constant_model, plastic_vform, constant_term):
- growth_constant_model.plastic_vform = plastic_vform
-
- assert growth_constant_model.plastic_vform == plastic_vform
- assert growth_constant_model.type == 'constant'
- assert np.array_equal(growth_constant_model.dg, constant_term)
-
-def test_constant_growth_increment(growth_constant_model, plastic_vform):
- # Ensure time and dt are set properly
- growth_constant_model.t = 0
- growth_constant_model.dt = 1
- growth_constant_model.plastic_vform == plastic_vform
-
- growth_constant_model.get_growth_increment()
- assert np.allclose(list(growth_constant_model.growth_incr.values()), [0.1, 0, 0, 0.1])
-
- # Update time
- growth_constant_model.update_time(1)
- assert growth_constant_model.t == 1
-
-def test_growth_apply(growth_constant_model, fenics_setup, previous_solution, plastic_vform):
- # Ensure time and dt are set properly
- growth_constant_model.t = 0
- growth_constant_model.dt = 1
- growth_constant_model.prev_solution = previous_solution
- growth_constant_model.prev_growth = plastic_vform.F_g
- growth_constant_model.plastic_vform = plastic_vform
-
- # Build mock problem
- _, V = fenics_setup
- u = fe.TrialFunction(V)
- v = fe.TestFunction(V)
- sol = fe.Function(V)
-
- # Call apply_growth
- growth_constant_model.apply_growth(u, v, sol)
-
- # Check the new growth values
- growth = growth_constant_model.get_growth()
- assert np.allclose(growth.compute_vertex_values(), [1.1, 1.1, 1.1, 1.1,
- 0.0, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0,
- 1.1, 1.1, 1.1, 1.1])
-
-def test_time_growth_init(growth_time_model, plastic_vform, growth_func):
- assert growth_time_model.type == 'time-dependent'
- assert growth_time_model.growth_func == growth_func
- assert growth_time_model.args == (1.0, 0.5)
-
- growth_time_model.plastic_vform = plastic_vform
- assert growth_time_model.plastic_vform == plastic_vform
-
-
-def test_time_growth_increment(monkeypatch, growth_time_model, plastic_vform):
- # Mock quad_vec to return a specific value
- dg_dt = np.zeros((2, 2))
- dg_dt[0, 0] = 1.0
- monkeypatch.setattr('scipy.integrate.quad_vec', lambda *args, **kwargs: (dg_dt, 0))
-
- # Ensure time and dt are set properly
- growth_time_model.t = 0.5
- growth_time_model.dt = 0.5
- growth_time_model.plastic_vform = plastic_vform
-
- # Call method
- growth_time_model.get_growth_increment()
-
- # Check growth increment
- assert np.allclose(growth_time_model.growth_incr.values(), [0.65915494, 0, 0, 0])
-
-def test_strain_growth_init(growth_strain_model, plastic_vform, target_strain, extensibility):
- assert growth_strain_model.type == 'strain-dependent'
- assert np.array_equal(growth_strain_model.target_strain, target_strain)
- assert growth_strain_model.extensibility == extensibility
-
- growth_strain_model.plastic_vform = plastic_vform
- assert growth_strain_model.plastic_vform == plastic_vform
-
-def test_strain_growth_increment(growth_strain_model, target_strain, plastic_vform):
- # Ensure time and dt are set properly
- growth_strain_model.t = 0
- growth_strain_model.dt = 1
- growth_strain_model.plastic_vform = plastic_vform
-
- # Call the growth increment method with non-null target strain
- growth_strain_model.get_growth_increment()
-
- # - Check if the growth increment matches the expected result
- target_strain_growth_incr = growth_strain_model.get_growth(only_increment=True)
- assert np.allclose(target_strain_growth_incr.compute_vertex_values(),
- [0.311, 0.002, 0.311, 0.311,
- -0.242, -0.006, -0.242, -0.242,
- -0.242, -0.006, -0.242, -0.242,
- 0.189, 0.020, 0.189, 0.189], atol=1e-3)
-
- # Call the growth increment method with null target strain
- growth_strain_model.target_strain = None
- growth_strain_model.get_growth_increment()
-
- # Check if the growth increment matches the expected result
- no_target_strain_growth_incr = growth_strain_model.get_growth(only_increment=True)
- assert np.allclose(no_target_strain_growth_incr.compute_vertex_values(),
- [ 0.25, -0.1775, 0.25, 0.25,
- -0.32, -0.0625, -0.32, -0.32,
- -0.32, -0.0625, -0.32, -0.32,
- 0.09, 0.0025, 0.09, 0.09])
-
-def test_heterogeneous_init(growth_heterogeneous_model, plastic_vform):
- assert (growth_heterogeneous_model.type[1] == 'constant') and (growth_heterogeneous_model.type[2] == 'strain-dependent')
- assert growth_heterogeneous_model.t == 0
-
- growth_heterogeneous_model.dt = 0
- assert growth_heterogeneous_model.dt == 0
-
- # - Propagate plastic_vform & prev_growth
- growth_heterogeneous_model.plastic_vform = plastic_vform
- assert growth_heterogeneous_model.plastic_vform is not None
-
- growth_heterogeneous_model.prev_growth = plastic_vform.F_g
- assert growth_heterogeneous_model.prev_growth is not None
-
-def test_heterogeneous_increment(growth_heterogeneous_model, previous_solution, plastic_vform):
- growth_heterogeneous_model.t = 0
- growth_heterogeneous_model.dt = 1
- growth_heterogeneous_model.prev_solution = previous_solution
- growth_heterogeneous_model.plastic_vform = plastic_vform
-
- growth_heterogeneous_model.get_growth_increment()
-
- growth_incr = growth_heterogeneous_model.get_growth(only_increment=True)
- assert np.allclose(growth_incr.compute_vertex_values(), [ 0.311, 0.100, 0.311, 0.311,
- -0.242, 0.000, -0.242, -0.242,
- -0.242, 0.000, -0.242, -0.242,
- 0.189, 0.100, 0.189, 0.189], atol=1e-3)
-
-
+class TestGrowthModels(unittest.TestCase):
+ def setUp(self):
+ # 1) Constant term for ConstantGrowth
+ self.constant_term = np.array([[0.1, 0.0],
+ [0.0, 0.1]])
+
+ # 2) Growth function for TimeDependentGrowth
+ def _growth_func(t, a, b):
+ dg = np.zeros((2, 2))
+ dg[0, 0] = b * np.sin(t * np.pi) + a
+ return dg
+ self.growth_func = _growth_func
+
+ # 3) Minimal FEniCS setup: mesh and function space
+ self.mesh = fe.UnitSquareMesh(1, 1)
+ self.V = fe.VectorFunctionSpace(self.mesh, 'P', 1)
+
+ # 4) Create a "previous solution" with some DOFs modified
+ self.prev_solution = fe.Function(self.V)
+ self.prev_solution.assign(fe.Constant([1.0, 1.0]))
+ vec = self.prev_solution.vector()
+ dof_indices = [0, 1, 2, 3, 4]
+ new_values = [0.5, 0.8, 1.5, 1.2, 0.9]
+ for idx, val in zip(dof_indices, new_values):
+ vec[idx] = val
+
+ # 5) Plastic form for all growth models
+ potential_energy = StVenantKirchoffPotential()
+ self.plastic_vform = MorphoElasticForm(
+ potential_energy=potential_energy,
+ F_g=np.eye(2),
+ source=[0.0, 0.0]
+ )
+
+ # 6) TimeDependentGrowth model
+ a = 1.0
+ b = 0.5
+ self.growth_time_model = TimeDependentGrowth(self.growth_func, a, b)
+
+ # 7) ConstantGrowth model
+ self.growth_constant_model = ConstantGrowth(self.constant_term)
+
+ # 8) StrainDependentGrowth model
+ self.target_strain = 0.1
+ self.extensibility = 0.5
+ self.growth_strain_model = StrainDependentGrowth(
+ target_strain=self.target_strain,
+ extensibility=self.extensibility
+ )
+ self.growth_strain_model.prev_solution = self.prev_solution
+
+ # 9) HeterogeneousGrowth model with subdomains
+ values = [1, 2, 1, 2]
+ cdata = fe.MeshFunction('size_t', self.mesh, self.mesh.topology().dim())
+ for cell, value in zip(fe.cells(self.mesh), values):
+ cdata[cell.index()] = value
+ growth_classes = {1: self.growth_constant_model, 2: self.growth_strain_model}
+ self.growth_heterogeneous_model = HeterogeneousGrowth(
+ growth_classes=growth_classes,
+ cdata=cdata
+ )
+
+ # -------------------------------------------------------------------------
+ # Tests for ConstantGrowth
+ # -------------------------------------------------------------------------
+ def test_constant_growth_init(self):
+ self.growth_constant_model.plastic_vform = self.plastic_vform
+
+ self.assertEqual(self.growth_constant_model.plastic_vform, self.plastic_vform)
+ self.assertEqual(self.growth_constant_model.type, 'constant')
+ self.assertTrue(np.array_equal(self.growth_constant_model.dg, self.constant_term))
+
+ def test_constant_growth_increment(self):
+ # Ensure time and dt are set properly
+ self.growth_constant_model.t = 0
+ self.growth_constant_model.dt = 1
+ self.growth_constant_model.plastic_vform = self.plastic_vform
+
+ # Get growth increment
+ self.growth_constant_model.get_growth_increment()
+ # Check the dictionary of growth_incr
+ growth_incr_values = list(self.growth_constant_model.growth_incr.values())
+ assert_allclose(growth_incr_values, [0.1, 0, 0, 0.1], atol=1e-12)
+
+ # Update time
+ self.growth_constant_model.update_time(1)
+ self.assertEqual(self.growth_constant_model.t, 1)
+
+ def test_growth_apply(self):
+ # Prepare model for apply_growth
+ self.growth_constant_model.t = 0
+ self.growth_constant_model.dt = 1
+ self.growth_constant_model.prev_solution = self.prev_solution
+ self.growth_constant_model.prev_growth = self.plastic_vform.F_g
+ self.growth_constant_model.plastic_vform = self.plastic_vform
+
+ # Build a simple problem
+ u = fe.TrialFunction(self.V)
+ v = fe.TestFunction(self.V)
+ sol = fe.Function(self.V)
+
+ # Call apply_growth
+ self.growth_constant_model.apply_growth(u, v, sol)
+
+ # Check the new growth values
+ growth = self.growth_constant_model.get_growth()
+ # Compare vertex values
+ computed_vals = growth.compute_vertex_values()
+ expected_vals = np.array([
+ 1.1, 1.1, 1.1, 1.1,
+ 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0,
+ 1.1, 1.1, 1.1, 1.1
+ ])
+ assert_allclose(computed_vals, expected_vals, atol=1e-12)
+
+ # -------------------------------------------------------------------------
+ # Tests for TimeDependentGrowth
+ # -------------------------------------------------------------------------
+ def test_time_growth_init(self):
+ self.assertEqual(self.growth_time_model.type, 'time-dependent')
+ self.assertEqual(self.growth_time_model.growth_func, self.growth_func)
+ self.assertEqual(self.growth_time_model.args, (1.0, 0.5))
+
+ self.growth_time_model.plastic_vform = self.plastic_vform
+ self.assertEqual(self.growth_time_model.plastic_vform, self.plastic_vform)
+
+ @patch('scipy.integrate.quad_vec')
+ def test_time_growth_increment(self, mock_quad_vec):
+ # Mock quad_vec to return a specific value
+ dg_dt = np.zeros((2, 2))
+ dg_dt[0, 0] = 1.0
+ mock_quad_vec.return_value = (dg_dt, 0)
+
+ # Ensure time and dt are set properly
+ self.growth_time_model.t = 0.5
+ self.growth_time_model.dt = 0.5
+ self.growth_time_model.plastic_vform = self.plastic_vform
+
+ # Call method
+ self.growth_time_model.get_growth_increment()
+
+ # Check growth increment
+ incr_vals = list(self.growth_time_model.growth_incr.values())
+ assert_allclose(incr_vals, [0.65915494, 0, 0, 0], atol=1e-7)
+
+ # -------------------------------------------------------------------------
+ # Tests for StrainDependentGrowth
+ # -------------------------------------------------------------------------
+ def test_strain_growth_init(self):
+ self.assertEqual(self.growth_strain_model.type, 'strain-dependent')
+ self.assertTrue(np.array_equal(self.growth_strain_model.target_strain, self.target_strain))
+ self.assertEqual(self.growth_strain_model.extensibility, self.extensibility)
+
+ self.growth_strain_model.plastic_vform = self.plastic_vform
+ self.assertEqual(self.growth_strain_model.plastic_vform, self.plastic_vform)
+
+ def test_strain_growth_increment(self):
+ # Ensure time and dt are set
+ self.growth_strain_model.t = 0
+ self.growth_strain_model.dt = 1
+ self.growth_strain_model.plastic_vform = self.plastic_vform
+
+ # 1) With non-null target strain
+ self.growth_strain_model.get_growth_increment()
+ target_strain_growth_incr = self.growth_strain_model.get_growth(only_increment=True)
+ assert_allclose(
+ target_strain_growth_incr.compute_vertex_values(),
+ [
+ 0.279, 0., 0.279, 0.279,
+ -0.218, 0., -0.218, -0.218,
+ -0.218, 0., -0.218, -0.218,
+ 0.170, 0., 0.170, 0.170
+ ],
+ atol=1e-3
+ )
+
+ # 2) With null target strain
+ self.growth_strain_model.target_strain = None
+ self.growth_strain_model.get_growth_increment()
+ no_target_strain_growth_incr = self.growth_strain_model.get_growth(only_increment=True)
+ assert_allclose(
+ no_target_strain_growth_incr.compute_vertex_values(),
+ [
+ 0.25, -0.1775, 0.25, 0.25,
+ -0.32, -0.0625, -0.32, -0.32,
+ -0.32, -0.0625, -0.32, -0.32,
+ 0.09, 0.0025, 0.09, 0.09
+ ],
+ atol=1e-3
+ )
+
+ # -------------------------------------------------------------------------
+ # Tests for HeterogeneousGrowth
+ # -------------------------------------------------------------------------
+ def test_heterogeneous_init(self):
+ sub_types = self.growth_heterogeneous_model.type
+ # Should be {1: 'constant', 2: 'strain-dependent'}
+ self.assertEqual(sub_types[1], 'constant')
+ self.assertEqual(sub_types[2], 'strain-dependent')
+ self.assertEqual(self.growth_heterogeneous_model.t, 0)
+
+ self.growth_heterogeneous_model.dt = 0
+ self.assertEqual(self.growth_heterogeneous_model.dt, 0)
+
+ # Plastic form
+ self.growth_heterogeneous_model.plastic_vform = self.plastic_vform
+ self.assertIsNotNone(self.growth_heterogeneous_model.plastic_vform)
+
+ # Previous growth
+ self.growth_heterogeneous_model.prev_solution = self.prev_solution
+ self.growth_heterogeneous_model.prev_growth = self.plastic_vform.F_g
+ self.assertIsNotNone(self.growth_heterogeneous_model.prev_growth)
+
+ def test_heterogeneous_increment(self):
+ self.growth_heterogeneous_model.t = 0
+ self.growth_heterogeneous_model.dt = 1
+ self.growth_heterogeneous_model.prev_solution = self.prev_solution
+ self.growth_heterogeneous_model.plastic_vform = self.plastic_vform
+
+ # Get the global increment
+ self.growth_heterogeneous_model.get_growth_increment()
+ growth_incr = self.growth_heterogeneous_model.get_growth(only_increment=True)
+
+ # Compare vertex values
+ # Expect subdomain 1 => constant model => 0.1 increment in diagonal
+ # subdomain 2 => strain model => ~ 0.179 increment (plus previous strain effect)
+ expected = [
+ 0.279, 0.100, 0.279, 0.279,
+ -0.218, 0.000, -0.218, -0.218,
+ -0.218, 0.000, -0.218, -0.218,
+ 0.170, 0.100, 0.170, 0.170
+ ]
+ computed = growth_incr.compute_vertex_values()
+ assert_allclose(computed, expected, atol=1e-3)