Merge branch 'fenicsx' into time_as_Constant

festim/__init__.py

@@ -32,3 +32,5 @@
 
 from .subdomain.surface_subdomain import SurfaceSubdomain1D
 from .subdomain.volume_subdomain import VolumeSubdomain1D
+
+from .exports.vtx import VTXExport

festim/exports/vtx.py

@@ -0,0 +1,83 @@
+from dolfinx.io import VTXWriter
+import mpi4py
+
+import festim as F
+
+
+class VTXExport:
+    """Export functions to VTX file
+
+    Args:
+        filename (str): the name of the output file
+        field (int): the field index to export
+
+    Attributes:
+        filename (str): the name of the output file
+        writer (dolfinx.io.VTXWriter): the VTX writer
+        field (festim.Species, list of festim.Species): the field index to export
+
+    Usage:
+        >>> u = dolfinx.fem.Function(V)
+        >>> my_export = festim.VTXExport("my_export.bp")
+        >>> my_export.define_writer(mesh.comm, [u])
+        >>> for t in range(10):
+        ...    u.interpolate(lambda x: t * (x[0] ** 2 + x[1] ** 2 + x[2] ** 2))
+        ...    my_export.write(t)
+    """
+
+    def __init__(self, filename: str, field) -> None:
+        self.filename = filename
+        self.field = field
+
+    @property
+    def filename(self):
+        return self._filename
+
+    @filename.setter
+    def filename(self, value):
+        if not isinstance(value, str):
+            raise TypeError("filename must be of type str")
+        if not value.endswith(".bp"):
+            raise ValueError("filename must end with .bp")
+        self._filename = value
+
+    @property
+    def field(self):
+        return self._field
+
+    @field.setter
+    def field(self, value):
+        # check that field is festim.Species or list of festim.Species
+        if not isinstance(value, F.Species) and not isinstance(value, list):
+            raise TypeError(
+                "field must be of type festim.Species or list of festim.Species"
+            )
+        # check that all elements of list are festim.Species
+        if isinstance(value, list):
+            for element in value:
+                if not isinstance(element, F.Species):
+                    raise TypeError(
+                        "field must be of type festim.Species or list of festim.Species"
+                    )
+        # if field is festim.Species, convert to list
+        if not isinstance(value, list):
+            value = [value]
+
+        self._field = value
+
+    def define_writer(self, comm: mpi4py.MPI.Intracomm, functions: list) -> None:
+        """Define the writer
+
+        Args:
+            comm (mpi4py.MPI.Intracomm): the MPI communicator
+            functions (list): the list of functions to export
+        """
+        self.writer = VTXWriter(comm, self.filename, functions, "BP4")
+
+    def write(self, t: float):
+        """Write functions to VTX file
+
+        Args:
+            t (float): the time of export
+        """
+        self.writer.write(t)

festim/hydrogen_transport_problem.py

@@ -119,6 +119,18 @@ def initialise(self):
         self.define_boundary_conditions()
         self.create_formulation()
         self.create_solver()
+        self.defing_export_writers()
+
+    def defing_export_writers(self):
+        """Defines the export writers of the model"""
+        for export in self.exports:
+            # TODO implement when export.field is an int or str
+            # then find solution from index of species
+
+            if isinstance(export, F.VTXExport):
+                export.define_writer(
+                    MPI.COMM_WORLD, [field.solution for field in export.field]
+                )
 
     def define_function_space(self):
         elements = ufl.FiniteElement("CG", self.mesh.mesh.ufl_cell(), 1)

test/test_permeation_problem.py

@@ -38,6 +38,8 @@ def test_permeation_problem():
         ),
     ]
 
+    my_model.exports = [F.VTXExport("test.bp", field=mobile_H)]
+
     my_model.initialise()
 
     D = my_mat.get_diffusion_coefficient(my_mesh.mesh, my_model.temperature)
@@ -129,6 +131,7 @@ def test_permeation_problem_multi_volume():
             subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="H"
         ),
     ]
+    my_model.exports = [F.VTXExport("test.bp", field=mobile_H)]
 
     my_model.initialise()
 

test/test_vtx.py

@@ -0,0 +1,90 @@
+import festim as F
+import dolfinx
+from mpi4py import MPI
+import numpy as np
+
+import pytest
+
+mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
+V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
+
+
+def test_vtx_export_one_function(tmpdir):
+    """Test can add one function to a vtx export"""
+    u = dolfinx.fem.Function(V)
+
+    filename = str(tmpdir.join("my_export.bp"))
+    my_export = F.VTXExport(filename, field=F.Species("H"))
+    my_export.define_writer(mesh.comm, [u])
+
+    for t in range(10):
+        u.interpolate(lambda x: t * (x[0] ** 2 + x[1] ** 2 + x[2] ** 2))
+
+        my_export.write(t)
+
+
+def test_vtx_export_two_functions(tmpdir):
+    """Test can add two functions to a vtx export"""
+    u = dolfinx.fem.Function(V)
+    v = dolfinx.fem.Function(V)
+
+    filename = str(tmpdir.join("my_export.bp"))
+    my_export = F.VTXExport(filename, field=[F.Species("1"), F.Species("2")])
+
+    my_export.define_writer(mesh.comm, [u, v])
+
+    for t in range(10):
+        u.interpolate(lambda x: t * (x[0] ** 2 + x[1] ** 2 + x[2] ** 2))
+        v.interpolate(lambda x: t * (x[0] ** 2 + x[1] ** 2 + x[2] ** 2))
+
+        my_export.write(t)
+
+
+def test_vtx_integration_with_h_transport_problem(tmpdir):
+    my_model = F.HydrogenTransportProblem()
+    my_model.mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0]))
+    my_mat = F.Material(D_0=1, E_D=0, name="mat")
+    my_model.subdomains = [
+        F.VolumeSubdomain1D(1, borders=[0.0, 4.0], material=my_mat),
+        F.SurfaceSubdomain1D(1, x=0.0),
+        F.SurfaceSubdomain1D(2, x=4.0),
+    ]
+    my_model.species = [F.Species("H")]
+    my_model.temperature = 500
+
+    filename = str(tmpdir.join("my_export.bp"))
+    my_export = F.VTXExport(filename, field=my_model.species[0])
+    my_model.exports = [my_export]
+
+    my_model.initialise()
+
+    for t in range(10):
+        my_export.write(t)
+
+
+def test_field_attribute_is_always_list():
+    """Test that the field attribute is always a list"""
+    my_export = F.VTXExport("my_export.bp", field=F.Species("H"))
+    assert isinstance(my_export.field, list)
+
+    my_export = F.VTXExport("my_export.bp", field=[F.Species("H")])
+    assert isinstance(my_export.field, list)
+
+
+@pytest.mark.parametrize("field", ["H", 1, [F.Species("H"), 1]])
+def test_field_attribute_raises_error_when_invalid_type(field):
+    """Test that the field attribute raises an error if the type is not festim.Species or list"""
+    with pytest.raises(TypeError):
+        F.VTXExport("my_export.bp", field=field)
+
+
+def test_filename_raises_error_with_wrong_extension():
+    """Test that the filename attribute raises an error if the extension is not .bp"""
+    with pytest.raises(ValueError):
+        F.VTXExport("my_export.txt", field=[F.Species("H")])
+
+
+def test_filename_raises_error_when_wrong_type():
+    """Test that the filename attribute raises an error if the extension is not .bp"""
+    with pytest.raises(TypeError):
+        F.VTXExport(1, field=[F.Species("H")])