Coverage for physioblocks / library / model_components / rheology.py: 100%

1# SPDX-FileCopyrightText: Copyright INRIA 

2# 

3# SPDX-License-Identifier: LGPL-3.0-only 

4# 

5# Copyright INRIA 

6# 

7# This file is part of PhysioBlocks, a library mostly developed by the 

8# [Ananke project-team](https://team.inria.fr/ananke) at INRIA. 

9# 

10# Authors: 

11# - Colin Drieu 

12# - Dominique Chapelle 

13# - François Kimmig 

14# - Philippe Moireau 

15# 

16# PhysioBlocks is free software: you can redistribute it and/or modify it under the 

17# terms of the GNU Lesser General Public License as published by the Free Software 

18# Foundation, version 3 of the License. 

19# 

20# PhysioBlocks is distributed in the hope that it will be useful, but WITHOUT ANY 

21# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 

22# PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. 

23# 

24# You should have received a copy of the GNU Lesser General Public License along with 

25# PhysioBlocks. If not, see <https://www.gnu.org/licenses/>. 

26 

27""" 

28Describes rheology models 

29""" 

30 

31from dataclasses import dataclass 

32from typing import Any 

33 

34import numpy as np 

35 

36from physioblocks.computing import Expression, ModelComponent, Quantity, diff, mid_point 

37from physioblocks.registers import register_type 

38from physioblocks.simulation import Time 

39 

40# Constant for the rheology model type id 

41RHEOLOGY_FIBER_ADDITIVE_TYPE_ID = "rheology_fiber_additive" 

42 

43 

44@register_type(RHEOLOGY_FIBER_ADDITIVE_TYPE_ID) 

45@dataclass 

46class RheologyFiberAdditiveModelComponent(ModelComponent): 

47 r""" 

48 Describes a Fiber Additive Rheology model implementation. 

49 

50 **Internal Equation:** 

51 

52 .. math:: 

53 

54 \dot{e}_c + T_c - k_s \Bigl( \frac{y}{R_0}-e_c \Bigr) = 0 

55 

56 

57 **Discretised:** 

58 

59 .. math:: 

60 

61 \mu \frac{e_c^{n+1}-e_c^n}{\Delta t^n} 

62 - k_s \Bigl( \frac{y^{n + \frac{1}{2}}}{R_0}-e_c^{n + \frac{1}{2}} \Bigr) 

63 + T_c^{{n + \frac{1}{2}}\sharp} = 0 

64 

65 """ 

66 

67 fib_deform: Quantity[np.float64] 

68 """:math:`e_c` the fiber deformation""" 

69 

70 disp: Quantity[np.float64] 

71 """:math:`y` the displacement""" 

72 

73 active_tension_discr: Quantity[np.float64] 

74 """:math:`T_c^{{n + \frac{1}{2}}\\sharp}` the active tension discretisation""" 

75 

76 radius: Quantity[np.float64] 

77 """:math:`R_0` sphere initial radius""" 

78 

79 damping_parallel: Quantity[np.float64] 

80 """:math:`\\mu` the damping parallel""" 

81 

82 series_stiffness: Quantity[np.float64] 

83 """:math:`k_s` the series stiffness""" 

84 

85 time: Time 

86 """:math:`t` the simulation time""" 

87 

88 def initialize(self) -> None: 

89 """ 

90 Initialize model's radius inverse 

91 """ 

92 self._inv_radius = 1 / self.radius.current 

93 

94 def fib_deform_equation(self) -> Any: 

95 """ 

96 Compute the equation representing the fiber deformation. 

97 

98 :return: the relation value 

99 :rtype: np.float64 

100 """ 

101 

102 disp_mid_pt = mid_point(self.disp) 

103 fib_deform_mid_pt = mid_point(self.fib_deform) 

104 

105 return ( 

106 self.active_tension_discr.new 

107 + (self.damping_parallel.current * self.time.inv_dt * diff(self.fib_deform)) 

108 - ( 

109 self.series_stiffness.current 

110 * (disp_mid_pt * self._inv_radius - fib_deform_mid_pt) 

111 ) 

112 ) 

113 

114 def dfib_deform_equation_dfib_deform(self) -> Any: 

115 """ 

116 Compute the equation partial derivative for ``fib_deform`` 

117 

118 :return: the partial derivative value 

119 :rtype: np.float64 

120 """ 

121 

122 return ( 

123 self.damping_parallel.current * self.time.inv_dt 

124 + self.series_stiffness.current * 0.5 

125 ) 

126 

127 def dfib_deform_equation_ddisp(self) -> np.float64: 

128 """ 

129 Compute the equation partial derivative for ``disp`` 

130 

131 :return: the partial derivative value 

132 :rtype: np.float64 

133 """ 

134 return -self.series_stiffness.current * self._inv_radius * 0.5 

135 

136 def dfib_deform_equation_dactive_tension_discr(self) -> Any: 

137 """ 

138 Compute the equation partial derivative for the ``active_tension_discr`` 

139 

140 :return: the partial derivative value 

141 :rtype: np.float64 

142 """ 

143 return 1.0 

144 

145 

146# Define the expression of the equation of the model. 

147_fib_deform_equation_expr = Expression( 

148 1, 

149 RheologyFiberAdditiveModelComponent.fib_deform_equation, 

150 { 

151 "fib_deform": RheologyFiberAdditiveModelComponent.dfib_deform_equation_dfib_deform, # noqa: E501 

152 "disp": RheologyFiberAdditiveModelComponent.dfib_deform_equation_ddisp, 

153 "active_tension_discr": RheologyFiberAdditiveModelComponent.dfib_deform_equation_dactive_tension_discr, # noqa: E501 

154 }, 

155) 

156 

157RheologyFiberAdditiveModelComponent.declares_internal_expression( 

158 "fib_deform", 

159 _fib_deform_equation_expr, 

160)