Skip to content

Models

BaseModel

Base class for defining material models.

This class provides common functionalities and utilities for material models.

Parameters:

Name Type Description Default
pars dict

Dictionary containing parameters of the material model.

 required

Attributes:

Name Type Description
la Constant

Lame coefficient lambda.
mu Constant

Lame coefficient mu.
Source code in src/fragma/models.py 
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
class BaseModel:
    """
    Base class for defining material models.

    This class provides common functionalities and utilities for material models.

    Parameters
    ----------
    pars : dict
        Dictionary containing parameters of the material model.

    Attributes
    ----------
    la : dolfinx.Constant
        Lame coefficient lambda.
    mu : dolfinx.Constant
        Lame coefficient mu.
    """

    def __init__(self, pars, domain):
        """
        Initialize the BaseModel.

        Parameters
        ----------
        pars : dict
            Dictionary containing parameters of the material model.
        domain : fragma.Domain.domain
            Domain object used to initialize heterogeneous properties.
        """
        # Get elastic parameters
        self.E = parse_parameter(pars["mechanical"]["E"], domain)
        self.nu = parse_parameter(pars["mechanical"]["nu"], domain)
        # Compute Lame coefficient
        self.la = self.E * self.nu / ((1 + self.nu) * (1 - 2 * self.nu))
        self.mu = self.E / (2 * (1 + self.nu))
        # Check the 2D assumption
        self.dim = pars["model"]["dim"]
        if self.dim == 2:
            self.assumption = pars["model"]["2D_assumption"]
            match self.assumption:
                case "plane_stress":
                    print("Plane stress assumption")
                    self.la = 2 * self.mu * self.la / (self.la + 2 * self.mu)
                case "plane_strain":
                    print("Plane strain assumption")
                case _:
                    raise ValueError(
                        f'The 2D assumption "{self.assumption}" in unknown'
                    )
        # Get the optional thermal load (to compute the thermal strain)
        self.thermal_load = pars["loading"].get("thermal_load", {})
        if self.thermal_load:
            # Get the thermal expansion coefficient
            self.a_T = parse_parameter(
                self.thermal_load["thermal_expansion_coeff"], domain
            )
            # Get the temperature field (variation)
            self.dT, self.dT_lambda = parse_parameter(
                self.thermal_load["dT"], domain, export_lambda=True
            )

    def eps_th(self):
        """
        Compute the thermal strain (for thermal loads).

        Returns
        -------
        ufl.form.Expression
            Strain tensor.
        """
        # Compute the thermal strain
        coeff = self.a_T * self.dT if self.thermal_load else 0
        return coeff * ufl.Identity(2)

    def eps(self, state):
        """
        Compute the strain tensor.

        Parameters
        ----------
        state : dict
            Dictionary containing state variables.

        Returns
        -------
        ufl.form.Expression
            Strain tensor.
        """
        return ufl.sym(ufl.grad(state["u"]))

    def eps_ela(self, state):
        """
        Compute the elastic strain tensor.

        Parameters
        ----------
        state : dict
            Dictionary containing state variables.

        Returns
        -------
        ufl.form.Expression
            Strain tensor.
        """
        return self.eps(state) - self.eps_th()

    def ela(self):
        """
        Compute the elasticity tensor.

        Returns
        -------
        ufl.form.Expression
            Elasticity tensor.
        """
        # Define index for tensorial notations
        i, j, k, l = ufl.indices(4)
        # Compute constant tensors
        Id2 = ufl.Identity(self.dim)
        Id2xId2 = ufl.outer(Id2, Id2)
        Id4 = (
            1
            / 2
            * ufl.as_tensor(Id2[i, k] * Id2[j, l] + Id2[i, l] * Id2[j, k], (i, j, k, l))
        )
        # Compute the elasticity tensor
        return 2 * self.mu * Id4 + self.la * Id2xId2

    def ela_eff(self, state):
        """
        Compute the effective elasticity tensor.

        Parameters
        ----------
        state : dict
            Dictionary containing state variables.

        Returns
        -------
        ufl.form.Expression
            Elasticity tensor.
        """
        # Compute the elasticity tensor
        return self.ela()

    def sig(self, state):
        """
        Compute the stress tensor.

        Parameters
        ----------
        state : dict
            Dictionary containing state variables.

        Returns
        -------
        ufl.form.Expression
            Stress tensor.
        """
        # Generate indices
        i, j, k, l = ufl.indices(4)
        # Get elastic parameters
        ela = self.ela()
        # Compute the strain
        eps = self.eps(state)
        # Compute the stess
        return ufl.as_tensor(ela[i, j, k, l] * eps[k, l], (i, j))

    def sig_eff(self, state):
        """
        Compute the effective stress tensor.

        Parameters
        ----------
        state : dict
            Dictionary containing state variables.

        Returns
        -------
        ufl.form.Expression
            Effective stress tensor.
        """
        # Generate indices
        i, j, k, l = ufl.indices(4)
        # Get elastic parameters
        ela_eff = self.ela_eff(state)
        # Compute the strain
        eps = self.eps(state)
        # Compute the stess
        return ufl.as_tensor(ela_eff[i, j, k, l] * eps[k, l], (i, j))

    def energy(self, state, mesh):
        """
        Compute the energy.

        This method should be implemented in the child class.

        Parameters
        ----------
        state : dict
            Dictionary containing state variables.
        mesh : dolfinx.Mesh
            The mesh representing the domain.

        Raises
        ------
        NotImplementedError
            If the method is not implemented in the child class.
        """
        raise NotImplementedError(
            "Model: The method 'energy' must be implemented in the child class."
        )

__init__(pars, domain)

Initialize the BaseModel.

Parameters:

Name Type Description Default
pars dict

Dictionary containing parameters of the material model.

 required
domain domain

Domain object used to initialize heterogeneous properties.

 required
Source code in src/fragma/models.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
def __init__(self, pars, domain):
    """
    Initialize the BaseModel.

    Parameters
    ----------
    pars : dict
        Dictionary containing parameters of the material model.
    domain : fragma.Domain.domain
        Domain object used to initialize heterogeneous properties.
    """
    # Get elastic parameters
    self.E = parse_parameter(pars["mechanical"]["E"], domain)
    self.nu = parse_parameter(pars["mechanical"]["nu"], domain)
    # Compute Lame coefficient
    self.la = self.E * self.nu / ((1 + self.nu) * (1 - 2 * self.nu))
    self.mu = self.E / (2 * (1 + self.nu))
    # Check the 2D assumption
    self.dim = pars["model"]["dim"]
    if self.dim == 2:
        self.assumption = pars["model"]["2D_assumption"]
        match self.assumption:
            case "plane_stress":
                print("Plane stress assumption")
                self.la = 2 * self.mu * self.la / (self.la + 2 * self.mu)
            case "plane_strain":
                print("Plane strain assumption")
            case _:
                raise ValueError(
                    f'The 2D assumption "{self.assumption}" in unknown'
                )
    # Get the optional thermal load (to compute the thermal strain)
    self.thermal_load = pars["loading"].get("thermal_load", {})
    if self.thermal_load:
        # Get the thermal expansion coefficient
        self.a_T = parse_parameter(
            self.thermal_load["thermal_expansion_coeff"], domain
        )
        # Get the temperature field (variation)
        self.dT, self.dT_lambda = parse_parameter(
            self.thermal_load["dT"], domain, export_lambda=True
        )

ela()

Compute the elasticity tensor.

Returns:

Type Description
Expression

Elasticity tensor.
Source code in src/fragma/models.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
def ela(self):
    """
    Compute the elasticity tensor.

    Returns
    -------
    ufl.form.Expression
        Elasticity tensor.
    """
    # Define index for tensorial notations
    i, j, k, l = ufl.indices(4)
    # Compute constant tensors
    Id2 = ufl.Identity(self.dim)
    Id2xId2 = ufl.outer(Id2, Id2)
    Id4 = (
        1
        / 2
        * ufl.as_tensor(Id2[i, k] * Id2[j, l] + Id2[i, l] * Id2[j, k], (i, j, k, l))
    )
    # Compute the elasticity tensor
    return 2 * self.mu * Id4 + self.la * Id2xId2

ela_eff(state)

Compute the effective elasticity tensor.

Parameters:

Name Type Description Default
state dict

Dictionary containing state variables.

 required

Returns:

Type Description
Expression

Elasticity tensor.
Source code in src/fragma/models.py 
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
def ela_eff(self, state):
    """
    Compute the effective elasticity tensor.

    Parameters
    ----------
    state : dict
        Dictionary containing state variables.

    Returns
    -------
    ufl.form.Expression
        Elasticity tensor.
    """
    # Compute the elasticity tensor
    return self.ela()

energy(state, mesh)

Compute the energy.

This method should be implemented in the child class.

Parameters:

Name Type Description Default
state dict

Dictionary containing state variables.

 required
mesh Mesh

The mesh representing the domain.

 required

Raises:

Type Description
NotImplementedError

If the method is not implemented in the child class.
Source code in src/fragma/models.py 
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
def energy(self, state, mesh):
    """
    Compute the energy.

    This method should be implemented in the child class.

    Parameters
    ----------
    state : dict
        Dictionary containing state variables.
    mesh : dolfinx.Mesh
        The mesh representing the domain.

    Raises
    ------
    NotImplementedError
        If the method is not implemented in the child class.
    """
    raise NotImplementedError(
        "Model: The method 'energy' must be implemented in the child class."
    )

eps(state)

Compute the strain tensor.

Parameters:

Name Type Description Default
state dict

Dictionary containing state variables.

 required

Returns:

Type Description
Expression

Strain tensor.
Source code in src/fragma/models.py 
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
def eps(self, state):
    """
    Compute the strain tensor.

    Parameters
    ----------
    state : dict
        Dictionary containing state variables.

    Returns
    -------
    ufl.form.Expression
        Strain tensor.
    """
    return ufl.sym(ufl.grad(state["u"]))

eps_ela(state)

Compute the elastic strain tensor.

Parameters:

Name Type Description Default
state dict

Dictionary containing state variables.

 required

Returns:

Type Description
Expression

Strain tensor.
Source code in src/fragma/models.py 
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def eps_ela(self, state):
    """
    Compute the elastic strain tensor.

    Parameters
    ----------
    state : dict
        Dictionary containing state variables.

    Returns
    -------
    ufl.form.Expression
        Strain tensor.
    """
    return self.eps(state) - self.eps_th()

eps_th()

Compute the thermal strain (for thermal loads).

Returns:

Type Description
Expression

Strain tensor.
Source code in src/fragma/models.py 
72
73
74
75
76
77
78
79
80
81
82
83
def eps_th(self):
    """
    Compute the thermal strain (for thermal loads).

    Returns
    -------
    ufl.form.Expression
        Strain tensor.
    """
    # Compute the thermal strain
    coeff = self.a_T * self.dT if self.thermal_load else 0
    return coeff * ufl.Identity(2)

sig(state)

Compute the stress tensor.

Parameters:

Name Type Description Default
state dict

Dictionary containing state variables.

 required

Returns:

Type Description
Expression

Stress tensor.
Source code in src/fragma/models.py 
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
def sig(self, state):
    """
    Compute the stress tensor.

    Parameters
    ----------
    state : dict
        Dictionary containing state variables.

    Returns
    -------
    ufl.form.Expression
        Stress tensor.
    """
    # Generate indices
    i, j, k, l = ufl.indices(4)
    # Get elastic parameters
    ela = self.ela()
    # Compute the strain
    eps = self.eps(state)
    # Compute the stess
    return ufl.as_tensor(ela[i, j, k, l] * eps[k, l], (i, j))

sig_eff(state)

Compute the effective stress tensor.

Parameters:

Name Type Description Default
state dict

Dictionary containing state variables.

 required

Returns:

Type Description
Expression

Effective stress tensor.
Source code in src/fragma/models.py 
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
def sig_eff(self, state):
    """
    Compute the effective stress tensor.

    Parameters
    ----------
    state : dict
        Dictionary containing state variables.

    Returns
    -------
    ufl.form.Expression
        Effective stress tensor.
    """
    # Generate indices
    i, j, k, l = ufl.indices(4)
    # Get elastic parameters
    ela_eff = self.ela_eff(state)
    # Compute the strain
    eps = self.eps(state)
    # Compute the stess
    return ufl.as_tensor(ela_eff[i, j, k, l] * eps[k, l], (i, j))

ElasticModel

Bases: BaseModel

Material model for linear elasticity.

This class implements the material model for linear elasticity.

Parameters:

Name Type Description Default
pars dict

Dictionary containing parameters of the material model.

 required
Source code in src/fragma/models.py 
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
class ElasticModel(BaseModel):
    """
    Material model for linear elasticity.

    This class implements the material model for linear elasticity.

    Parameters
    ----------
    pars : dict
        Dictionary containing parameters of the material model.
    """

    def elastic_energy(self, state, domain):
        """
        Compute the elastic energy.

        Parameters
        ----------
        state : dict
            Dictionary containing state variables.
        domain : Domain
            The domain object representing the computational domain.

        Returns
        -------
        ufl.form.Expression
            Elastic energy.
        """
        # Get the integrands
        dx = ufl.Measure("dx", domain=domain.mesh, metadata={"quadrature_degree": 12})
        # Compute the effective elasticity tensor
        ela_eff = self.ela_eff(state)
        # Compute the elastic strain
        eps_ela = self.eps_ela(state)
        # Define the total energy
        return 1 / 2 * ufl.inner(ela_eff, ufl.outer(eps_ela, eps_ela)) * dx

    def energy(self, state, domain):
        """
        Compute the total energy.

        Parameters
        ----------
        state : dict
            Dictionary containing state variables.
        domain : Domain
            The domain object representing the computational domain.

        Returns
        -------
        ufl.form.Expression
            Total energy.
        """
        # Define the energy terms
        elastic_energy = self.elastic_energy(state, domain)
        # Define the total energy
        return elastic_energy

elastic_energy(state, domain)

Compute the elastic energy.

Parameters:

Name Type Description Default
state dict

Dictionary containing state variables.

 required
domain Domain

The domain object representing the computational domain.

 required

Returns:

Type Description
Expression

Elastic energy.
Source code in src/fragma/models.py 
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
def elastic_energy(self, state, domain):
    """
    Compute the elastic energy.

    Parameters
    ----------
    state : dict
        Dictionary containing state variables.
    domain : Domain
        The domain object representing the computational domain.

    Returns
    -------
    ufl.form.Expression
        Elastic energy.
    """
    # Get the integrands
    dx = ufl.Measure("dx", domain=domain.mesh, metadata={"quadrature_degree": 12})
    # Compute the effective elasticity tensor
    ela_eff = self.ela_eff(state)
    # Compute the elastic strain
    eps_ela = self.eps_ela(state)
    # Define the total energy
    return 1 / 2 * ufl.inner(ela_eff, ufl.outer(eps_ela, eps_ela)) * dx

energy(state, domain)

Compute the total energy.

Parameters:

Name Type Description Default
state dict

Dictionary containing state variables.

 required
domain Domain

The domain object representing the computational domain.

 required

Returns:

Type Description
Expression

Total energy.
Source code in src/fragma/models.py 
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
def energy(self, state, domain):
    """
    Compute the total energy.

    Parameters
    ----------
    state : dict
        Dictionary containing state variables.
    domain : Domain
        The domain object representing the computational domain.

    Returns
    -------
    ufl.form.Expression
        Total energy.
    """
    # Define the energy terms
    elastic_energy = self.elastic_energy(state, domain)
    # Define the total energy
    return elastic_energy

FractureModel

Bases: ElasticModel

Material model for fracture mechanics.

This class implements the material model for fracture mechanics.

Parameters:

Name Type Description Default
pars dict

Dictionary containing parameters of the material model.

 required
Source code in src/fragma/models.py 
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
class FractureModel(ElasticModel):
    """
    Material model for fracture mechanics.

    This class implements the material model for fracture mechanics.

    Parameters
    ----------
    pars : dict
        Dictionary containing parameters of the material model.
    """

    def __init__(self, pars, domain):
        """
        Initialize the FractureModel.

        Parameters
        ----------
        pars : dict
            Dictionary containing parameters of the material model.
        """
        # Initialise parent class
        super().__init__(pars, domain)
        # Get the degradation model
        model_par = pars["model"]["model"]
        if model_par in ["AT1", "AT2"]:
            self.deg_model = model_par
        else:
            self.deg_model = model_par.split("-")[0]
        # Get the dissipation model
        if model_par in ["AT1", "AT2"]:
            self.dis_model = model_par
        else:
            self.dis_model = model_par.split("-")[1]
        # Get the residual crack phase
        self.alpha_res = pars["numerical"]["alpha_res"]
        # Get fracture parameters
        self.ell = parse_parameter(pars["mechanical"]["ell"], domain)
        # Check for anisotropy
        self.is_anisotropic = "theta_0" in pars["mechanical"]
        if not self.is_anisotropic:
            # Get the critical energy release rate
            self.Gc = parse_parameter(pars["mechanical"]["Gc"], domain)
        else:
            # Get the critical energy release rate (min and max)
            Gc_min = parse_parameter(pars["mechanical"]["Gc_min"], domain)
            Gc_max = parse_parameter(pars["mechanical"]["Gc_max"], domain)
            # Convert to other model parameters
            self.Gc = ufl.sqrt(1 / 2 * (Gc_min**2 + Gc_max**2))
            self.aG = 1 / 2 * (Gc_max**2 - Gc_min**2) / self.Gc**2
            # Ge the anisotropy angle
            self.theta_0 = (
                parse_parameter(pars["mechanical"]["theta_0"], domain) * np.pi / 180
                if "theta_0" in pars["mechanical"]
                else 0
            )
        # Check for model specific parameters
        if self.dis_model in ["Foc2", "Foc4", "RMBRAT1", "RMBRAT2"]:
            self.Gc = parse_parameter(pars["mechanical"]["Gc"], domain)
            self.D2 = parse_parameter(pars["mechanical"]["D2"], domain)
            self.P2 = parse_parameter(pars["mechanical"]["P2"], domain)
        if self.dis_model in ["Foc4", "RMBRAT1", "RMBRAT2"]:
            self.D4 = parse_parameter(pars["mechanical"]["D4"], domain)
            self.P4 = parse_parameter(pars["mechanical"]["P4"], domain)
        if self.dis_model in ["Foc2X"]:
            self.Gc = parse_parameter(pars["mechanical"]["Gc"], domain)
            self.th0 = parse_parameter(pars["mechanical"]["th0"], domain)
            self.dG = parse_parameter(pars["mechanical"]["dG"], domain)

    def a(self, alpha):
        """
        Degradation function.

        Parameters
        ----------
        alpha : ufl.form.Expression
            Crack phase.

        Returns
        -------
        ufl.form.Expression
            Degradation function.
        """
        # Residual crack phase
        alpha_res = self.alpha_res
        # Compute a
        match self.deg_model:
            case "AT1" | "AT2" | "Foc2" | "Foc2X" | "RMBR":
                return (1 - alpha) ** 2 + alpha_res
            case "KKL":
                return 4 * (1 - alpha) ** 3 - 4 * (1 - alpha) ** 3 + alpha_res
            case "KSM":
                return 3 * (1 - alpha) ** 2 - 3 * (1 - alpha) ** 2 + alpha_res
            case "Foc4":
                return (1 - alpha**4) ** 2 + alpha_res
            case _:
                raise ValueError(
                    f"The degradation model named '{self.deg_model}' does not exists."
                )

    def ap(self, alpha):
        """
        Derivative of the degradation function.

        Parameters
        ----------
        alpha : ufl.form.Expression
            Crack phase.

        Returns
        -------
        ufl.form.Expression
            Derivative of the degradation function.
        """
        # Define a variable
        alpha = ufl.variable(alpha)
        # Comupute the derivative
        return ufl.diff(self.a(alpha), alpha)

    def w(self, alpha):
        """
        Dissipation function.

        Parameters
        ----------
        alpha : ufl.form.Expression
            Crack phase.

        Returns
        -------
        ufl.form.Expression
            Dissipation function.
        """
        # Compute w
        match self.dis_model:
            case "AT1" | "RMBRAT1":
                return alpha
            case "AT2" | "Foc2" | "Foc2X" | "RMBRAT2":
                return alpha**2
            case "DW":
                return 16 * alpha**2 * (1 - alpha) ** 2
            case "Foc4":
                bw = 2 ** (-4 / 3)
                return 3 / bw * alpha**4
            case _:
                raise ValueError(
                    f"The dissipation model named '{self.dis_model}' does not exists."
                )

    def wp(self, alpha):
        """
        Derivative of the dissipation function.

        Parameters
        ----------
        alpha : ufl.form.Expression
            Crack phase.

        Returns
        -------
        ufl.form.Expression
            Derivative of the dissipation function.
        """
        # Define a variable
        alpha = ufl.variable(alpha)
        # Comupute the derivative
        return ufl.diff(self.w(alpha), alpha)

    def cw(self):
        """
        Normalization coefficient.

        Returns
        -------
        float
            Normalization coefficient.
        """
        match self.dis_model:
            case "AT1" | "RMBRAT1":
                return 8 / 3
            case "AT2" | "Foc2" | "Foc2X" | "RMBRAT2":
                return 2
            case "DW":
                return 4 * 2 / 3
            case "Foc4":
                return 4
            case _:
                raise ValueError(
                    f"The dissipation model named '{self.dis_model}' does not exists."
                )

    def ela_eff(self, state):
        """
        Compute the effective elasticity tensor.

        Parameters
        ----------
        state : dict
            Dictionary containing state variables.

        Returns
        -------
        ufl.form.Expression
            Elasticity tensor.
        """
        # Compute the elasticity tensor
        return self.a(state["alpha"]) * self.ela()

    def fracture_dissipation(self, state, domain):
        """
        Compute the energy dissipated by fracture.

        Parameters
        ----------
        state : dict
            Dictionary containing state variables.
        domain : Domain
            The domain object representing the computational domain.

        Returns
        -------
        ufl.form.Expression
            Energy dissipated by fracture.

        """
        # Get the integrands
        dx = ufl.Measure("dx", domain=domain.mesh, metadata={"quadrature_degree": 12})
        # Get state variables
        alpha = state["alpha"]
        # Get the fracture parameters
        Gc = self.Gc
        ell = self.ell
        cw = self.cw()
        # Check the model
        match self.dis_model:
            case "Foc2":
                # Define 2-order identity
                id2 = ufl.Identity(2)
                # Define the covariant
                h2_np = np.empty((2, 2))
                h2_np[0, 0] = -self.D2 * ufl.cos(2 * self.P2)  # Idenpendent components
                h2_np[0, 1] = -self.D2 * ufl.sin(2 * self.P2)
                h2_np[1, 1] = -h2_np[0, 0]  # Traceless condition
                h2_np[1, 0] = h2_np[0, 1]  # Symmetry conditions
                h2 = ufl.as_tensor(h2_np)
                B = id2 + h2
                # Define the anisotropy function
                grada = ufl.grad(alpha)
                grad2 = ufl.outer(grada, grada)
                phi2 = ufl.inner(B, grad2)
                # Define the dissipation terms
                dissipated_energy = Gc / cw * (self.w(alpha) / ell + ell * phi2) * dx
            case "Foc2X":
                # Define the anisotropy tensor
                A_np = np.empty((2, 2))
                A_np[0, 0] = -ufl.sin(2 * self.th0)
                A_np[0, 1] = ufl.cos(2 * self.th0)
                A_np[1, 0] = ufl.cos(2 * self.th0)
                A_np[1, 1] = ufl.sin(2 * self.th0)
                A = ufl.as_tensor(A_np)
                # Define the anisotropy function
                grada = ufl.grad(alpha)
                grad2 = ufl.outer(grada, grada)
                phi2_1 = ufl.inner(grada, grada)
                phi2_2 = self.dG / Gc * ufl.inner(A, grad2) ** 2
                phi2 = (phi2_1 + phi2_2) ** 2
                # Define the dissipation term
                dissipated_energy = Gc / cw * (self.w(alpha) / ell + ell * phi2) * dx
            case "Foc4":
                # Define constants
                id2 = ufl.Identity(2)
                id4_np = np.empty((2, 2, 2, 2))
                indices = itertools.product(range(2), range(2), range(2), range(2))
                for i, j, k, l in indices:
                    id4_np[i, j, k, l] = (
                        1 / 2 * (id2[i, k] * id2[j, l] + id2[i, l] * id2[j, k])
                    )
                id4 = ufl.as_tensor(id4_np)
                # Define the 2nd order covariant
                h2_np = np.empty((2, 2))
                h2_np[0, 0] = -self.D2 * ufl.cos(2 * self.P2)  # Idenpendent components
                h2_np[0, 1] = -self.D2 * ufl.sin(2 * self.P2)
                h2_np[1, 1] = -h2_np[0, 0]  # Traceless condition
                h2_np[1, 0] = h2_np[0, 1]  # Symmetry conditions
                h2 = ufl.as_tensor(h2_np)
                # Define the 4th order covariant
                h4_np = np.empty((2, 2, 2, 2))
                h4_np[0, 0, 0, 0] = self.D4 * ufl.cos(
                    4 * self.P4
                )  # Idenpendent components
                h4_np[0, 0, 0, 1] = self.D4 * ufl.sin(4 * self.P4)
                h4_np[0, 0, 1, 1] = -h4_np[0, 0, 0, 0]  # Traceless conditions
                h4_np[0, 1, 1, 1] = -h4_np[0, 0, 0, 1]
                h4_np[1, 1, 1, 1] = h4_np[0, 0, 0, 0]
                h4_np[0, 0, 1, 0] = h4_np[0, 0, 0, 1]  # Symmetries conditions
                h4_np[0, 1, 0, 0] = h4_np[0, 0, 0, 1]
                h4_np[1, 0, 0, 0] = h4_np[0, 0, 0, 1]
                h4_np[1, 1, 0, 0] = h4_np[0, 0, 1, 1]
                h4_np[0, 1, 0, 1] = h4_np[0, 0, 1, 1]
                h4_np[1, 0, 1, 0] = h4_np[0, 0, 1, 1]
                h4_np[1, 0, 0, 1] = h4_np[0, 0, 1, 1]
                h4_np[0, 1, 1, 0] = h4_np[0, 0, 1, 1]
                h4_np[1, 0, 1, 1] = h4_np[0, 1, 1, 1]
                h4_np[1, 1, 0, 1] = h4_np[0, 1, 1, 1]
                h4_np[1, 1, 1, 0] = h4_np[0, 1, 1, 1]
                h4 = ufl.as_tensor(h4_np)
                # Define the anistropy tensor
                B = id4 + 1.0 / 2.0 * (ufl.outer(id2, h2) + ufl.outer(h2, id2)) + h4
                # Compute the crack phase gradient
                grada = ufl.grad(alpha)
                grad2 = ufl.outer(grada, grada)
                grad4 = ufl.outer(grad2, grad2)
                # Compute the surface energy
                phi4 = ufl.inner(B, grad4)
                # Define the dissipation terms
                dissipated_energy = Gc / cw * (self.w(alpha) / ell + ell**3 * phi4) * dx
            case "RMBRAT1" | "RMBRAT2":
                alpha_0 = state["alpha0"]
                grada0 = ufl.grad(alpha_0)
                phi = ufl.conditional(
                    ufl.ge(ufl.sqrt(grada0[0] ** 2 + grada0[1] ** 2), 1e-12),
                    ufl.atan2(-grada0[1], grada0[0]),
                    np.pi / 2,
                )
                theta = phi - np.pi / 2
                Gc = self.Gc * (
                    1
                    + self.D2 * ufl.cos(2 * (theta - self.P2))
                    + self.D4 * ufl.cos(4 * (theta - self.P4))
                ) ** (1 / 4)

                wa = self.w(alpha)
                grada = ufl.grad(alpha)
                grada_grada = ufl.dot(grada, grada)
                dissipated_energy = Gc / cw * (wa / ell + ell * grada_grada) * dx
            case "AT1" | "AT2" | "DW":
                # Compute the anisotropy matrix
                A = ufl.as_tensor(np.eye(domain.mesh.geometry.dim))
                # Add the higher order terms if the model is anisotropic
                if self.is_anisotropic:
                    # Get the parameters
                    aG, theta_0 = self.aG, self.theta_0
                    #  Compute the 2nd order term of the anisotropy tensor
                    A += aG * ufl.as_tensor(
                        np.array(
                            [
                                [np.cos(2 * theta_0), np.sin(2 * theta_0)],
                                [np.sin(2 * theta_0), -np.cos(2 * theta_0)],
                            ]
                        )
                    )
                # Define the energy terms
                dissipated_energy = (
                    Gc
                    / cw
                    * (
                        self.w(alpha) / ell
                        + ell * ufl.dot(ufl.grad(alpha), A * ufl.grad(alpha))
                    )
                    * dx
                )
            case _:
                raise ValueError(
                    f"The dissipation model {self.dis_model} does not exists."
                )
        # Define the total energy
        return dissipated_energy

    def energy(self, state, domain):
        """
        Compute the energy.

        Parameters
        ----------
        state : dict
            Dictionary containing state variables.
        domain : Domain
            The domain object representing the computational domain.

        Returns
        -------
        ufl.form.Expression
            Total energy.
        """
        # Define the energy terms
        elastic_energy = self.elastic_energy(state, domain)
        dissipated_energy = self.fracture_dissipation(state, domain)
        # Define the total energy
        return elastic_energy + dissipated_energy

__init__(pars, domain)

Initialize the FractureModel.

Parameters:

Name Type Description Default
pars dict

Dictionary containing parameters of the material model.

 required
Source code in src/fragma/models.py 
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
def __init__(self, pars, domain):
    """
    Initialize the FractureModel.

    Parameters
    ----------
    pars : dict
        Dictionary containing parameters of the material model.
    """
    # Initialise parent class
    super().__init__(pars, domain)
    # Get the degradation model
    model_par = pars["model"]["model"]
    if model_par in ["AT1", "AT2"]:
        self.deg_model = model_par
    else:
        self.deg_model = model_par.split("-")[0]
    # Get the dissipation model
    if model_par in ["AT1", "AT2"]:
        self.dis_model = model_par
    else:
        self.dis_model = model_par.split("-")[1]
    # Get the residual crack phase
    self.alpha_res = pars["numerical"]["alpha_res"]
    # Get fracture parameters
    self.ell = parse_parameter(pars["mechanical"]["ell"], domain)
    # Check for anisotropy
    self.is_anisotropic = "theta_0" in pars["mechanical"]
    if not self.is_anisotropic:
        # Get the critical energy release rate
        self.Gc = parse_parameter(pars["mechanical"]["Gc"], domain)
    else:
        # Get the critical energy release rate (min and max)
        Gc_min = parse_parameter(pars["mechanical"]["Gc_min"], domain)
        Gc_max = parse_parameter(pars["mechanical"]["Gc_max"], domain)
        # Convert to other model parameters
        self.Gc = ufl.sqrt(1 / 2 * (Gc_min**2 + Gc_max**2))
        self.aG = 1 / 2 * (Gc_max**2 - Gc_min**2) / self.Gc**2
        # Ge the anisotropy angle
        self.theta_0 = (
            parse_parameter(pars["mechanical"]["theta_0"], domain) * np.pi / 180
            if "theta_0" in pars["mechanical"]
            else 0
        )
    # Check for model specific parameters
    if self.dis_model in ["Foc2", "Foc4", "RMBRAT1", "RMBRAT2"]:
        self.Gc = parse_parameter(pars["mechanical"]["Gc"], domain)
        self.D2 = parse_parameter(pars["mechanical"]["D2"], domain)
        self.P2 = parse_parameter(pars["mechanical"]["P2"], domain)
    if self.dis_model in ["Foc4", "RMBRAT1", "RMBRAT2"]:
        self.D4 = parse_parameter(pars["mechanical"]["D4"], domain)
        self.P4 = parse_parameter(pars["mechanical"]["P4"], domain)
    if self.dis_model in ["Foc2X"]:
        self.Gc = parse_parameter(pars["mechanical"]["Gc"], domain)
        self.th0 = parse_parameter(pars["mechanical"]["th0"], domain)
        self.dG = parse_parameter(pars["mechanical"]["dG"], domain)

a(alpha)

Degradation function.

Parameters:

Name Type Description Default
alpha Expression

Crack phase.

 required

Returns:

Type Description
Expression

Degradation function.
Source code in src/fragma/models.py 
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
def a(self, alpha):
    """
    Degradation function.

    Parameters
    ----------
    alpha : ufl.form.Expression
        Crack phase.

    Returns
    -------
    ufl.form.Expression
        Degradation function.
    """
    # Residual crack phase
    alpha_res = self.alpha_res
    # Compute a
    match self.deg_model:
        case "AT1" | "AT2" | "Foc2" | "Foc2X" | "RMBR":
            return (1 - alpha) ** 2 + alpha_res
        case "KKL":
            return 4 * (1 - alpha) ** 3 - 4 * (1 - alpha) ** 3 + alpha_res
        case "KSM":
            return 3 * (1 - alpha) ** 2 - 3 * (1 - alpha) ** 2 + alpha_res
        case "Foc4":
            return (1 - alpha**4) ** 2 + alpha_res
        case _:
            raise ValueError(
                f"The degradation model named '{self.deg_model}' does not exists."
            )

ap(alpha)

Derivative of the degradation function.

Parameters:

Name Type Description Default
alpha Expression

Crack phase.

 required

Returns:

Type Description
Expression

Derivative of the degradation function.
Source code in src/fragma/models.py 
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
def ap(self, alpha):
    """
    Derivative of the degradation function.

    Parameters
    ----------
    alpha : ufl.form.Expression
        Crack phase.

    Returns
    -------
    ufl.form.Expression
        Derivative of the degradation function.
    """
    # Define a variable
    alpha = ufl.variable(alpha)
    # Comupute the derivative
    return ufl.diff(self.a(alpha), alpha)

cw()

Normalization coefficient.

Returns:

Type Description
float

Normalization coefficient.
Source code in src/fragma/models.py 
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
def cw(self):
    """
    Normalization coefficient.

    Returns
    -------
    float
        Normalization coefficient.
    """
    match self.dis_model:
        case "AT1" | "RMBRAT1":
            return 8 / 3
        case "AT2" | "Foc2" | "Foc2X" | "RMBRAT2":
            return 2
        case "DW":
            return 4 * 2 / 3
        case "Foc4":
            return 4
        case _:
            raise ValueError(
                f"The dissipation model named '{self.dis_model}' does not exists."
            )

ela_eff(state)

Compute the effective elasticity tensor.

Parameters:

Name Type Description Default
state dict

Dictionary containing state variables.

 required

Returns:

Type Description
Expression

Elasticity tensor.
Source code in src/fragma/models.py 
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
def ela_eff(self, state):
    """
    Compute the effective elasticity tensor.

    Parameters
    ----------
    state : dict
        Dictionary containing state variables.

    Returns
    -------
    ufl.form.Expression
        Elasticity tensor.
    """
    # Compute the elasticity tensor
    return self.a(state["alpha"]) * self.ela()

energy(state, domain)

Compute the energy.

Parameters:

Name Type Description Default
state dict

Dictionary containing state variables.

 required
domain Domain

The domain object representing the computational domain.

 required

Returns:

Type Description
Expression

Total energy.
Source code in src/fragma/models.py 
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
def energy(self, state, domain):
    """
    Compute the energy.

    Parameters
    ----------
    state : dict
        Dictionary containing state variables.
    domain : Domain
        The domain object representing the computational domain.

    Returns
    -------
    ufl.form.Expression
        Total energy.
    """
    # Define the energy terms
    elastic_energy = self.elastic_energy(state, domain)
    dissipated_energy = self.fracture_dissipation(state, domain)
    # Define the total energy
    return elastic_energy + dissipated_energy

fracture_dissipation(state, domain)

Compute the energy dissipated by fracture.

Parameters:

Name Type Description Default
state dict

Dictionary containing state variables.

 required
domain Domain

The domain object representing the computational domain.

 required

Returns:

Type Description
Expression

Energy dissipated by fracture.
Source code in src/fragma/models.py 
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
def fracture_dissipation(self, state, domain):
    """
    Compute the energy dissipated by fracture.

    Parameters
    ----------
    state : dict
        Dictionary containing state variables.
    domain : Domain
        The domain object representing the computational domain.

    Returns
    -------
    ufl.form.Expression
        Energy dissipated by fracture.

    """
    # Get the integrands
    dx = ufl.Measure("dx", domain=domain.mesh, metadata={"quadrature_degree": 12})
    # Get state variables
    alpha = state["alpha"]
    # Get the fracture parameters
    Gc = self.Gc
    ell = self.ell
    cw = self.cw()
    # Check the model
    match self.dis_model:
        case "Foc2":
            # Define 2-order identity
            id2 = ufl.Identity(2)
            # Define the covariant
            h2_np = np.empty((2, 2))
            h2_np[0, 0] = -self.D2 * ufl.cos(2 * self.P2)  # Idenpendent components
            h2_np[0, 1] = -self.D2 * ufl.sin(2 * self.P2)
            h2_np[1, 1] = -h2_np[0, 0]  # Traceless condition
            h2_np[1, 0] = h2_np[0, 1]  # Symmetry conditions
            h2 = ufl.as_tensor(h2_np)
            B = id2 + h2
            # Define the anisotropy function
            grada = ufl.grad(alpha)
            grad2 = ufl.outer(grada, grada)
            phi2 = ufl.inner(B, grad2)
            # Define the dissipation terms
            dissipated_energy = Gc / cw * (self.w(alpha) / ell + ell * phi2) * dx
        case "Foc2X":
            # Define the anisotropy tensor
            A_np = np.empty((2, 2))
            A_np[0, 0] = -ufl.sin(2 * self.th0)
            A_np[0, 1] = ufl.cos(2 * self.th0)
            A_np[1, 0] = ufl.cos(2 * self.th0)
            A_np[1, 1] = ufl.sin(2 * self.th0)
            A = ufl.as_tensor(A_np)
            # Define the anisotropy function
            grada = ufl.grad(alpha)
            grad2 = ufl.outer(grada, grada)
            phi2_1 = ufl.inner(grada, grada)
            phi2_2 = self.dG / Gc * ufl.inner(A, grad2) ** 2
            phi2 = (phi2_1 + phi2_2) ** 2
            # Define the dissipation term
            dissipated_energy = Gc / cw * (self.w(alpha) / ell + ell * phi2) * dx
        case "Foc4":
            # Define constants
            id2 = ufl.Identity(2)
            id4_np = np.empty((2, 2, 2, 2))
            indices = itertools.product(range(2), range(2), range(2), range(2))
            for i, j, k, l in indices:
                id4_np[i, j, k, l] = (
                    1 / 2 * (id2[i, k] * id2[j, l] + id2[i, l] * id2[j, k])
                )
            id4 = ufl.as_tensor(id4_np)
            # Define the 2nd order covariant
            h2_np = np.empty((2, 2))
            h2_np[0, 0] = -self.D2 * ufl.cos(2 * self.P2)  # Idenpendent components
            h2_np[0, 1] = -self.D2 * ufl.sin(2 * self.P2)
            h2_np[1, 1] = -h2_np[0, 0]  # Traceless condition
            h2_np[1, 0] = h2_np[0, 1]  # Symmetry conditions
            h2 = ufl.as_tensor(h2_np)
            # Define the 4th order covariant
            h4_np = np.empty((2, 2, 2, 2))
            h4_np[0, 0, 0, 0] = self.D4 * ufl.cos(
                4 * self.P4
            )  # Idenpendent components
            h4_np[0, 0, 0, 1] = self.D4 * ufl.sin(4 * self.P4)
            h4_np[0, 0, 1, 1] = -h4_np[0, 0, 0, 0]  # Traceless conditions
            h4_np[0, 1, 1, 1] = -h4_np[0, 0, 0, 1]
            h4_np[1, 1, 1, 1] = h4_np[0, 0, 0, 0]
            h4_np[0, 0, 1, 0] = h4_np[0, 0, 0, 1]  # Symmetries conditions
            h4_np[0, 1, 0, 0] = h4_np[0, 0, 0, 1]
            h4_np[1, 0, 0, 0] = h4_np[0, 0, 0, 1]
            h4_np[1, 1, 0, 0] = h4_np[0, 0, 1, 1]
            h4_np[0, 1, 0, 1] = h4_np[0, 0, 1, 1]
            h4_np[1, 0, 1, 0] = h4_np[0, 0, 1, 1]
            h4_np[1, 0, 0, 1] = h4_np[0, 0, 1, 1]
            h4_np[0, 1, 1, 0] = h4_np[0, 0, 1, 1]
            h4_np[1, 0, 1, 1] = h4_np[0, 1, 1, 1]
            h4_np[1, 1, 0, 1] = h4_np[0, 1, 1, 1]
            h4_np[1, 1, 1, 0] = h4_np[0, 1, 1, 1]
            h4 = ufl.as_tensor(h4_np)
            # Define the anistropy tensor
            B = id4 + 1.0 / 2.0 * (ufl.outer(id2, h2) + ufl.outer(h2, id2)) + h4
            # Compute the crack phase gradient
            grada = ufl.grad(alpha)
            grad2 = ufl.outer(grada, grada)
            grad4 = ufl.outer(grad2, grad2)
            # Compute the surface energy
            phi4 = ufl.inner(B, grad4)
            # Define the dissipation terms
            dissipated_energy = Gc / cw * (self.w(alpha) / ell + ell**3 * phi4) * dx
        case "RMBRAT1" | "RMBRAT2":
            alpha_0 = state["alpha0"]
            grada0 = ufl.grad(alpha_0)
            phi = ufl.conditional(
                ufl.ge(ufl.sqrt(grada0[0] ** 2 + grada0[1] ** 2), 1e-12),
                ufl.atan2(-grada0[1], grada0[0]),
                np.pi / 2,
            )
            theta = phi - np.pi / 2
            Gc = self.Gc * (
                1
                + self.D2 * ufl.cos(2 * (theta - self.P2))
                + self.D4 * ufl.cos(4 * (theta - self.P4))
            ) ** (1 / 4)

            wa = self.w(alpha)
            grada = ufl.grad(alpha)
            grada_grada = ufl.dot(grada, grada)
            dissipated_energy = Gc / cw * (wa / ell + ell * grada_grada) * dx
        case "AT1" | "AT2" | "DW":
            # Compute the anisotropy matrix
            A = ufl.as_tensor(np.eye(domain.mesh.geometry.dim))
            # Add the higher order terms if the model is anisotropic
            if self.is_anisotropic:
                # Get the parameters
                aG, theta_0 = self.aG, self.theta_0
                #  Compute the 2nd order term of the anisotropy tensor
                A += aG * ufl.as_tensor(
                    np.array(
                        [
                            [np.cos(2 * theta_0), np.sin(2 * theta_0)],
                            [np.sin(2 * theta_0), -np.cos(2 * theta_0)],
                        ]
                    )
                )
            # Define the energy terms
            dissipated_energy = (
                Gc
                / cw
                * (
                    self.w(alpha) / ell
                    + ell * ufl.dot(ufl.grad(alpha), A * ufl.grad(alpha))
                )
                * dx
            )
        case _:
            raise ValueError(
                f"The dissipation model {self.dis_model} does not exists."
            )
    # Define the total energy
    return dissipated_energy

w(alpha)

Dissipation function.

Parameters:

Name Type Description Default
alpha Expression

Crack phase.

 required

Returns:

Type Description
Expression

Dissipation function.
Source code in src/fragma/models.py 
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
def w(self, alpha):
    """
    Dissipation function.

    Parameters
    ----------
    alpha : ufl.form.Expression
        Crack phase.

    Returns
    -------
    ufl.form.Expression
        Dissipation function.
    """
    # Compute w
    match self.dis_model:
        case "AT1" | "RMBRAT1":
            return alpha
        case "AT2" | "Foc2" | "Foc2X" | "RMBRAT2":
            return alpha**2
        case "DW":
            return 16 * alpha**2 * (1 - alpha) ** 2
        case "Foc4":
            bw = 2 ** (-4 / 3)
            return 3 / bw * alpha**4
        case _:
            raise ValueError(
                f"The dissipation model named '{self.dis_model}' does not exists."
            )

wp(alpha)

Derivative of the dissipation function.

Parameters:

Name Type Description Default
alpha Expression

Crack phase.

 required

Returns:

Type Description
Expression

Derivative of the dissipation function.
Source code in src/fragma/models.py 
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
def wp(self, alpha):
    """
    Derivative of the dissipation function.

    Parameters
    ----------
    alpha : ufl.form.Expression
        Crack phase.

    Returns
    -------
    ufl.form.Expression
        Derivative of the dissipation function.
    """
    # Define a variable
    alpha = ufl.variable(alpha)
    # Comupute the derivative
    return ufl.diff(self.w(alpha), alpha)