Support for internal fluxes

[RemDelaporteMathurin]

[RemDelaporteMathurin]

@jhdark you can obtain internal fluxes:

import FESTIM as F
import fenics as f
import numpy as np


class InternalFlux(F.SurfaceFlux):
    def compute(self):
        field_to_prop = {
            "0": self.D,
            "solute": self.D,
            0: self.D,
            "T": self.thermal_cond
        }
        self.prop = field_to_prop[self.field]
        flux = f.assemble(self.prop*f.dot(f.grad(self.function), self.n)('-')*dS(self.surface))
        print("Internal flux: {:.2e}".format(flux))
        return flux


mesh = F.MeshFromVertices(np.linspace(0, 1, 51))

materials = F.Materials(
    [
        F.Material(1, 1, 0, borders=[0, 0.5]),
        F.Material(2, 1, 0, borders=[0.5, 1])
    ]
)


derived_quantities = F.DerivedQuantities()
derived_quantities.derived_quantities = [
    F.SurfaceFlux(0, surface=1),
    InternalFlux(0, surface=3)
]

exports = F.Exports(
    [
        derived_quantities
    ]
)

boundary_conditions = [
    F.DirichletBC(1, 1),
    F.DirichletBC(2, 0)
]

settings = F.Settings(1e-10, 1e-10, transient=False)

temperature = F.Temperature(100)

my_sim = F.Simulation(
    mesh=mesh,
    materials=materials,
    exports=exports,
    boundary_conditions=boundary_conditions,
    temperature=temperature,
    settings=settings
)

my_sim.initialise()

# mark internal facets
id_interface = 3
x_interface = 0.5
for facet in f.facets(mesh.mesh):
    x0 = facet.midpoint()
    if f.near(x0.x(), x_interface):
        mesh.surface_markers[facet] = id_interface

# create a dS measure
dS = f.Measure('dS', domain=mesh.mesh, subdomain_data=mesh.surface_markers)


my_sim.run()