1r"""Implementations of material model formulations and stress tensors
  2
  3# Material Model Formulations
  4
  5Valid material model formulations are listed in `FORMULATIONS`, below
  6are the strain-energy densities they correspond to. It is supplied to
  7`Mechanics` in the `formulation` argument of its constructor.
  8
  9## "beta"
 10
 11This is the official best formulation.
 12`mod_repr` is $\beta(\mathbf{x}_0)$.
 13
 14$$\Psi=c_1e^{\beta(\mathbf{x}_0)}\left(I_1-3-2\ln{J}\right)+D_1(\ln{J})^2$$
 15
 16## "alpha"
 17
 18This is the formulation described in our previous approach.
 19`mod_repr` is $\alpha(\mathbf{x}_0)$
 20
 21$$\Psi=\alpha(\mathbf{x}_0)c_1\left(I_1-3-2\ln{J}\right)+D_1(\ln{J})^2$$
 22
 23## "alpha_on_all"
 24
 25This formulation models both material parameters as varying in the same
 26way.
 27`mod_repr` is $\alpha(\mathbf{x}_0)$
 28
 29$$\Psi=\alpha(\mathbf{x}_0)\left[c_1\left(I_1-3-2\ln{J}\right)+D_1(\ln{J})^2\right]$$
 30
 31## "exclude_all_penalty"
 32
 33This formulation has only been used for testing purposes; testing
 34revealed that this formulation predicts nonzero stress in some
 35undeformed configurations, and is therefore aphysical.
 36`mod_repr` is $\alpha(\mathbf{x}_0)$
 37
 38$$\Psi=\alpha(\mathbf{x}_0)c_1\left(I_1-3\right)-2c_1\ln{J}+D_1(\ln{J})^2$$
 39
 40## "beta_tilde"
 41
 42This formulation was designed to utilize the benefits of the exponential
 43formulation while restricting modulus prediction to lie within a
 44prescribed range. The range must be supplied to `kwargs` in `Mechanics`
 45through float entries "beta_min" $\beta_{min}$ and "beta_max"
 46$\beta_{max}$.
 47`mod_repr` is $\beta(\mathbf{x}_0)$. A helper field is defined
 48
 49$$\tilde{\beta}(\mathbf{x}_0)=a\tanh\left(m\beta(\mathbf{x}_0)+b\right)+c$$
 50
 51where
 52
 53$$a=\frac{\beta_{max}-\beta_{min}}{2}$$
 54$$c=\frac{\beta_{max}+\beta_{min}}{2}$$
 55$$b=-\text{arctanh}\left(\frac{c}{a}\right)$$
 56$$m=\frac{1}{a\left(1-\left(\frac{c}{a}\right)^2\right)}$$
 57
 58and then strain energy density is
 59
 60$$\Psi=c_1e^{\tilde{\beta}(\mathbf{x}_0)}\left(I_1-3-2\ln{J}\right)+D_1(\ln{J})^2$$
 61
 62# Modulus Representation "mod_repr"
 63
 64The `mod_repr` argument to `Mechanics` is an appropriate FEniCS
 65function. In contrast, the `create_mod_repr` function takes a str desc
 66and generates an appropriate FEniCS 1st-order Lagrange scalar field on
 67the mesh.
 68
 69If `desc` matches an option in the dict `MOD_REPR_FIELDS`, then the
 70field will be assembled correspondingly. Otherwise, it will be
 71interpreted as the path to a full-shape/write-checkpoint .xdmf file with
 72the field in the attribute "mod_repr"
 73
 74# API
 75"""
 76from .header import *
 77from .geometry import *
 78from .kinematics import *
 79
 80from ufl import tanh
 81
 82
 83FORMULATIONS = [
 84    "beta",
 85    "alpha",
 86    "alpha_on_all",
 87    "exclude_all_penalty", # Note: stress not always 0 at no deformation
 88    "beta_tilde"
 89]
 90"""List of str valid formulation names"""
 91
 92# 
 93ZERO_FIX_F = [
 94    "beta",
 95    "beta_tilde"
 96]
 97"""List of formulations for which mod_repr==0 is unmodified gel
 98
 99Rest have mod_repr==1 is unmodified gel.
100"""
101
102MU_FF_DEFAULT = 108 # Pa
103r"""Rheometric measurement of shear modulus for $c_1=\frac{\mu}{2}$ in Pa"""
104
105
106# Input validation helpers
107def validate_formulation(f):
108    """Raises `ValueError` if `f` str is not a valid formulation."""
109    if f not in FORMULATIONS:
110        raise ValueError(f"{f} not a valid formulation")
111
112
113def _get_neohookean_psi(
114        formulation,
115        kinematics,
116        mu_ff,
117        nu,
118        mod_repr,
119        **kwargs
120    ):
121    geo = kinematics.geo
122
123    # Far field values
124    lmbda = 2*mu_ff*nu/(1-2*nu) # 1st Lame' parameter
125    c1 = mu_ff/2 * 1e-6
126    d1 = lmbda/2 * 1e-6
127
128    I1 = kinematics.Ic
129
130    lnJ = ln(kinematics.Ju)
131
132    # Stored strain energy density (compressible neo-Hookean model)
133    if formulation == "beta":
134        psi = exp(mod_repr)*c1*(I1 - 3 - 2*lnJ) + d1*(lnJ)**2
135    elif formulation == "alpha":
136        psi = mod_repr*c1*(I1 - 3 - 2*lnJ) + d1*(lnJ)**2
137    elif formulation == "alpha_on_all":
138        psi = mod_repr*(c1*(I1 - 3 - 2*lnJ) + d1*(lnJ)**2)
139    elif formulation == "exclude_all_penalty":
140        psi = mod_repr*c1*(I1 - 3) - 2*c1*lnJ + d1*(lnJ)**2
141    elif formulation == "beta_tilde":
142        beta_max = kwargs["beta_max"]
143        beta_min = kwargs["beta_min"]
144
145        a = 0.5*(beta_max - beta_min)
146        c = 0.5*(beta_max + beta_min)
147        b = -np.arctanh(c/a)
148        m = 1/(a*(1-((c/a)**2)))
149
150        beta = a*tanh(m*mod_repr + b) + c
151
152        psi = exp(beta)*c1*(I1 - 3 - 2*lnJ) + d1*(lnJ)**2
153    else:
154        # Assume no modification
155        psi = c1*(I1 - 3) - 2*c1*lnJ + d1*(lnJ)**2
156
157    return psi
158
159
160def get_homogeneous_field(geo, value):
161    """Returns FEniCS function with a homogeneous scalar value
162
163    * `geo`: type `gel.geometry.Geometry` with mesh whereupon the scalar
164    field will be defined
165    * `value`: float value the field will take
166    """
167    return interpolate(Constant(value), geo.V0)
168
169
170MOD_REPR_FIELDS = {
171    "zero" : (lambda geo : get_homogeneous_field(geo, 0.0)),
172    "one" : (lambda geo : get_homogeneous_field(geo, 1.0))
173}
174"""dict of names of pre-defined modulus representation fields to FEniCS
175function factories that take in a `gel.geometry.Geometry` as a single
176argument
177"""
178
179
180def validate_mod_repr_field(mod_repr_field):
181    """Raises `ValueError` if `mod_repr_field` cannot be interpreted as
182    a valid scalar field (see top of `gel.mechanics`)
183    """
184    if mod_repr_field not in MOD_REPR_FIELDS:
185        if not os.path.exists(mod_repr_field):
186            raise ValueError(
187                f"{mod_repr_field} not a known valid mod_repr field"
188            )
189
190
191def create_mod_repr(geo, desc):
192    """Returns FEniCS FE scalar function according to the description.
193
194    * `geo`: `gel.geometry.Geometry` with the mesh and function space
195    * `desc`: str description of the field, see the intro to
196    `gel.mechanics`.
197
198    Returns: FEniCS field named "ctl" due to primary use as control
199    variable
200    """
201    if isinstance(desc, str):
202        if desc in MOD_REPR_FIELDS:
203            # Case lookup
204            mod_repr = MOD_REPR_FIELDS[desc](geo)
205        else:
206            # Case read from file
207            mod_repr = load_shape(geo.V0, desc, "mod_repr")
208            mod_repr.set_allow_extrapolation(True)
209    else:
210        # Case given the function
211        mod_repr = desc
212    mod_repr.rename("ctl", "ctl")
213
214    return mod_repr
215
216
217def _get_nu_for_target_d1c1(d1c1):
218    """Returns Poisson ratio float from D1/C1 `d1c1` float."""
219    nu = 0.5*d1c1/(1+d1c1)
220    return nu
221
222
223def _get_d1c1_for_target_nu(nu):
224    """Returns D1/C1 float from Poisson ratio `nu` float."""
225    twice_nu = 2*nu
226    return twice_nu / (1 - twice_nu)
227
228
229class Mechanics:
230    """Information on and functions to retrieve material properties"""
231
232    def __init__(
233            self,
234            kinematics,
235            formulation="beta",
236            mu_ff=MU_FF_DEFAULT,
237            d1c1=1.0,
238            mod_repr=None,
239            bx_T=Geometry.CELL,
240            **kwargs
241        ):
242        r"""Object for convenient handling of mechanical quantities
243
244        * `kinematics`: `gel.kinematics.Kinematics` object with
245        displacements to use by default
246        * `formulation`: str the name of the formulation (see
247        `gel.mechanics`)
248        * `mu_ff`: float unmodified gel shear modulus from rheometry,
249        used to compute $c_1=\frac{\mu}{2}$
250        * `d1c1`: float ratio $\frac{D_1}{c_1}$ encoding compressibility
251        * `mod_repr`: FEniCS scalar field
252        * `bx_T`: surface tag used by `kinematics.geo` (see
253        `gel.geometry.Geometry`) on which to apply traction force when
254        provided to subsequent function calls.
255        * `kwargs`: dict additional material model formulation
256        parameters, see `gel.mechanics` intro.
257        """
258        validate_formulation(formulation)
259
260        self.kinematics = kinematics
261        """Type `gel.kinematics.Kinematics` with deformation data with
262        which to compute mechanical quantities by default
263        """
264
265        self.formulation = formulation
266        """str name of the formulation from `FORMULATIONS`"""
267        self.formulation_kwargs = kwargs
268        """dict additional arguments to formulation, see intro to
269        `gel.mechanics`
270        """
271
272        self.mu_ff = mu_ff
273        """Far-field (unmodified) infinitesimal strain shear modulus"""
274
275        self.d1c1 = d1c1
276        r"""Measure of incompressibility $\frac{D_1}{c_1}$"""
277
278        self.nu = _get_nu_for_target_d1c1(self.d1c1)
279        """Infinitesimal strain Poisson ratio for unmodified gel"""
280
281        self.bx_T = bx_T
282        """Boundary on which external traction would be applied"""
283
284        if mod_repr is None:
285            mod_repr = MOD_REPR_FIELDS["zero"](kinematics.geo)
286
287        self.mod_repr = mod_repr
288        """FEniCS scalar function for modulus control/representation"""
289
290    def get_psi(self, kinematics=None):
291        """Returns FEniCS expression for strain-energy density
292
293        Units are MPa
294
295        * `kinematics`: `gel.kinematics.Kinematics` optionally use a
296        different instance than the one originally provided to this
297        object for purposes such as computing stress tensors
298        """
299        if kinematics is None:
300            kinematics = self.kinematics
301
302        psi = _get_neohookean_psi(
303            self.formulation,
304            kinematics,
305            self.mu_ff,
306            self.nu,
307            self.mod_repr,
308            **self.formulation_kwargs
309        )
310
311        return psi
312
313    def get_energy(self, kinematics=None, B=None, T=None):
314        """
315        Returns FEniCS expression for total energy
316
317        Units are pJ
318
319        * `kinematics`: `gel.kinematics.Kinematics` optionally use a
320        different instance than the one originally provided to this
321        object for purposes such as computing stress tensors
322        * `B`: FEniCS expression for body force in hydrogel domain, by
323        default 0
324        * `T`: FEniCS expression for traction force applied on `bx_T`,
325        by default 0
326        """
327        if kinematics is None:
328            kinematics = self.kinematics
329
330        if B is None:
331            B = Constant((0.0, 0.0, 0.0))
332
333        if T is None:
334            T = Constant((0.0, 0.0, 0.0))
335
336        psi = self.get_psi(kinematics=kinematics)
337
338        Pi = (
339            psi*kinematics.geo.dx
340            - dot(B, kinematics.u)*kinematics.geo.dx
341            - dot(T, kinematics.u)*kinematics.geo.ds(self.bx_T)
342        )
343
344        return Pi
345
346    def get_pk1_stress(self):
347        """Returns 1st Piola-Kirchoff stress tensor FEniCS expression.
348
349        Units are MPa
350
351        Not ideal for computing surface tractions due to inaccuracy
352        in numerical integration, see `get_nodal_traction` for that
353        instead.
354        """
355        kinematics_probe = Kinematics(
356            self.kinematics.geo,
357            u=self.kinematics.u,
358            differentiable_F=True
359        )
360
361        psi = self.get_psi(kinematics_probe)
362
363        P = diff(psi, kinematics_probe.F) # 1st PK stress tensor
364
365        return P
366    
367    def get_cauchy_stress(self):
368        """Returns Cauchy stress tensor FEniCS expression.
369
370        Units are MPa
371        """
372        P = self.get_pk1_stress()
373        cauchy_stress = P*self.kinematics.F.T / self.kinematics.Ju
374        return cauchy_stress
375
376    def get_total_strain_energy(self):
377        """Returns float total strain energy in units of pJ"""
378        psi = self.get_psi()
379        Pi = psi*self.kinematics.geo.dx
380        return assemble(Pi) # In MPa*um^3 = J * 10^{-12}
381
382    def get_nodal_traction(self, surf_marker=Geometry.CELL):
383        """Returns FEniCS function with nodal surface traction forces.
384
385        Units are MPa. Only the case for `formulation` "beta" allowed.
386
387        * `surf_marker`: surface tag from `gel.geometry.Geometry` on
388        which to compute traction force.
389
390        Returns: a function in the vector 1st order Lagrange basis on
391        the full hydrogel mesh. DoFs on the requested surface are the
392        projected traction forces.
393
394        Raises: `ValueError` is formulation not implemented
395        """
396        if self.formulation != "beta":
397            raise ValueError(
398                "Nodal traction computation only implemented for beta "
399                "formulation due to requirement to manually define form"
400                " of PK tensor for full accuracy"
401            )
402
403        geo = self.kinematics.geo
404        kin = self.kinematics
405
406        alpha = exp(self.mod_repr)
407
408        F = kin.F
409        J = kin.Ju
410        kin.u.set_allow_extrapolation(True)
411
412        Finv = inv(F)
413
414        P = self.mu_ff*1e-6*(
415            alpha*(F.T - Finv.T)
416            + self.d1c1*ln(J)*Finv.T
417        )
418
419        Vproj = geo.V
420        t = Function(Vproj)
421
422        n = FacetNormal(geo.mesh)
423        T = dot(P, n)
424
425        dt = TrialFunction(Vproj)
426        v = TestFunction(Vproj)
427
428        L = assemble(dot(T, v)*geo.ds(surf_marker))
429        A = assemble(dot(dt, v)*geo.ds(surf_marker), keep_diagonal=True)
430        A.ident_zeros()
431
432        solve(A, t.vector(), L)
433
434        return t