hyperion/cylindrical_flux.py at main · festim-dev/hyperion · GitHub

1
2
3
4
5
6
7
8
9
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
import ufl
import numpy as np
import dolfinx
import festim as F


# here we define our own SurfaceFlux class that takes into account the cylindrical coordinate system


class CylindricalFlux(F.SurfaceFlux):
    azimuth_range: tuple[float, float] = (0.0, 2 * np.pi)

    def compute(
        self,
        u: dolfinx.fem.Function | ufl.indexed.Indexed,
        ds: ufl.Measure,
        entity_maps=None,
        restriction: str | None = None,
        subdomain_id: int | None = None,
    ):
        """Computes the value of the flux at the surface

        J =int(- D * grad(c) . n * r dS)

        Args:
            u: field for which the flux is computed
            ds: surface measure of the model
            entity_maps: entity maps relating parent mesh and submesh
            restriction: side of an interior facet to evaluate the flux on
            subdomain_id: id to index ds with, when it is not the surface's own
        """
        from scifem import assemble_scalar

        # obtain mesh normal from field
        # if case multispecies, solution is an index, use sub_function_space
        if isinstance(u, ufl.indexed.Indexed):
            mesh = self.field.sub_function_space.mesh
        else:
            mesh = u.function_space.mesh
        n = ufl.FacetNormal(mesh)
        x = ufl.SpatialCoordinate(mesh)
        r = x[0]

        if subdomain_id is None:
            subdomain_id = self.surface.id
        integrand = -self.D * r * ufl.dot(ufl.grad(u), n)
        if restriction is not None:
            integrand = ufl.as_ufl(integrand)(restriction)

        flux = assemble_scalar(
            dolfinx.fem.form(
                integrand * ds(subdomain_id),
                entity_maps=entity_maps,
            )
        )

        flux *= self.azimuth_range[1] - self.azimuth_range[0]

        self.value = flux
        self.data.append(self.value)