diff --git a/festim/__init__.py b/festim/__init__.py index ca22327fe..de4a2652c 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -21,4 +21,7 @@ from .species import Species, Trap +from .subdomain.surface_subdomain import SurfaceSubdomain1D +from .subdomain.volume_subdomain import VolumeSubdomain1D + from .hydrogen_transport_problem import HydrogenTransportProblem diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 20a037fec..e5dfd56b2 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -3,12 +3,9 @@ import ufl from mpi4py import MPI from dolfinx.fem import Function -from ufl import ( - TestFunction, - dot, - grad, - exp, -) +from dolfinx.mesh import meshtags, locate_entities +from ufl import TestFunction, dot, grad, exp, Measure +import numpy as np import festim as F @@ -38,8 +35,8 @@ class HydrogenTransportProblem: dx (dolfinx.fem.dx): the volume measure of the model ds (dolfinx.fem.ds): the surface measure of the model function_space (dolfinx.fem.FunctionSpace): the function space of the model - facet_tags (dolfinx.cpp.mesh.MeshTags): the facet tags of the model - volume_tags (dolfinx.cpp.mesh.MeshTags): the volume tags of the model + facet_meshtags (dolfinx.cpp.mesh.MeshTags): the facet tags of the model + volume_meshtags (dolfinx.cpp.mesh.MeshTags): the volume tags of the model formulation (ufl.form.Form): the formulation of the model solver (dolfinx.nls.newton.NewtonSolver): the solver of the model @@ -88,8 +85,8 @@ def __init__( self.dx = None self.ds = None self.function_space = None - self.facet_tags = None - self.volume_tags = None + self.facet_meshtags = None + self.volume_meshtags = None self.formulation = None def initialise(self): @@ -98,12 +95,7 @@ def initialise(self): """ self.define_function_space() - ( - self.facet_tags, - self.volume_tags, - self.dx, - self.ds, - ) = self.mesh.create_measures_and_tags(self.function_space) + self.define_markers_and_measures() self.assign_functions_to_species() self.create_formulation() @@ -111,6 +103,46 @@ def define_function_space(self): elements = ufl.FiniteElement("CG", self.mesh.mesh.ufl_cell(), 1) self.function_space = fem.FunctionSpace(self.mesh.mesh, elements) + def define_markers_and_measures(self): + dofs_facets, tags_facets = [], [] + + # TODO this should be a property of mesh + fdim = self.mesh.mesh.topology.dim - 1 + vdim = self.mesh.mesh.topology.dim + + # find all cells in domain and mark them as 0 + num_cells = self.mesh.mesh.topology.index_map(vdim).size_local + mesh_cell_indices = np.arange(num_cells, dtype=np.int32) + tags_volumes = np.full(num_cells, 0, dtype=np.int32) + + for sub_dom in self.subdomains: + if isinstance(sub_dom, F.SurfaceSubdomain1D): + dof = sub_dom.locate_dof(self.function_space) + dofs_facets.append(dof) + tags_facets.append(sub_dom.id) + if isinstance(sub_dom, F.VolumeSubdomain1D): + # find all cells in subdomain and mark them as sub_dom.id + entities = sub_dom.locate_subdomain_entities(self.mesh.mesh, vdim) + tags_volumes[entities] = sub_dom.id + + # dofs and tags need to be in np.in32 format for meshtags + dofs_facets = np.array(dofs_facets, dtype=np.int32) + tags_facets = np.array(tags_facets, dtype=np.int32) + + # define mesh tags + self.facet_meshtags = meshtags(self.mesh.mesh, fdim, dofs_facets, tags_facets) + self.volume_meshtags = meshtags( + self.mesh.mesh, vdim, mesh_cell_indices, tags_volumes + ) + + # define measures + self.ds = Measure( + "ds", domain=self.mesh.mesh, subdomain_data=self.facet_meshtags + ) + self.dx = Measure( + "dx", domain=self.mesh.mesh, subdomain_data=self.volume_meshtags + ) + def assign_functions_to_species(self): """Creates for each species the solution, prev solution and test function""" if len(self.species) > 1: @@ -124,8 +156,6 @@ def create_formulation(self): """Creates the formulation of the model""" if len(self.sources) > 1: raise NotImplementedError("Sources not implemented yet") - if len(self.subdomains) > 1: - raise NotImplementedError("Multiple subdomains not implemented yet") if len(self.species) > 1: raise NotImplementedError("Multiple species not implemented yet") diff --git a/festim/mesh/mesh.py b/festim/mesh/mesh.py index 8d6a47d43..90f7746bc 100644 --- a/festim/mesh/mesh.py +++ b/festim/mesh/mesh.py @@ -26,17 +26,3 @@ def __init__(self, mesh=None): self.mesh.topology.create_connectivity( self.mesh.topology.dim - 1, self.mesh.topology.dim ) - - def create_measures_and_tags(self, function_space): - """Creates the ufl.measure.Measure objects for self.ds and - self.dx, also passes the facet and volume tags - """ - facet_tags, volume_tags = self.create_meshtags(function_space) - dx = ufl.Measure("dx", domain=self.mesh, subdomain_data=volume_tags) - ds = ufl.Measure("ds", domain=self.mesh, subdomain_data=facet_tags) - return ( - facet_tags, - volume_tags, - dx, - ds, - ) diff --git a/festim/mesh/mesh_1d.py b/festim/mesh/mesh_1d.py index 84f9a9faa..6a29e4c54 100644 --- a/festim/mesh/mesh_1d.py +++ b/festim/mesh/mesh_1d.py @@ -36,34 +36,3 @@ def generate_mesh(self): cells = np.stack((indexes[:-1], indexes[1:]), axis=-1) return mesh.create_mesh(MPI.COMM_WORLD, cells, mesh_points, domain) - - def create_meshtags(self, function_space): - """Creates the meshtags for a given function space - Args: - function_space (dolfinx.fem.function.FunctionSpace): the function - space of the model - - Returns: - dolfinx.mesh.MeshTagsMetaClass: the tags containing the facet - and volume tags - """ - - dofs_L = fem.locate_dofs_geometrical( - function_space, lambda x: np.isclose(x[0], min(self.vertices)) - ) - dofs_R = fem.locate_dofs_geometrical( - function_space, lambda x: np.isclose(x[0], max(self.vertices)) - ) - - dofs_facets = np.array([dofs_L[0], dofs_R[0]], dtype=np.int32) - tags_facets = np.array([1, 2], dtype=np.int32) - - facet_dimension = self.mesh.topology.dim - 1 - mesh_tags_facets = mesh.meshtags( - self.mesh, facet_dimension, dofs_facets, tags_facets - ) - - # TODO implement this - mesh_tags_volumes = None - - return mesh_tags_facets, mesh_tags_volumes diff --git a/festim/subdomain/__init__.py b/festim/subdomain/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/festim/subdomain/surface_subdomain.py b/festim/subdomain/surface_subdomain.py new file mode 100644 index 000000000..409ca6671 --- /dev/null +++ b/festim/subdomain/surface_subdomain.py @@ -0,0 +1,39 @@ +from dolfinx import fem +import numpy as np + + +class SurfaceSubdomain1D: + """ + Surface subdomain class for 1D cases + + Args: + id (int): the id of the surface subdomain + x (float): the x coordinate of the surface subdomain + + Attributes: + id (int): the id of the surface subdomain + x (float): the x coordinate of the surface subdomain + + Usage: + >>> surf_subdomain = F.SurfaceSubdomain1D(id=1, x=1) + """ + + def __init__(self, id, x) -> None: + self.id = id + self.x = x + + def locate_dof(self, function_space): + """Locates the dof of the surface subdomain within the function space + + Args: + function_space (dolfinx.fem.FunctionSpace): the function space of + the model + + Returns: + dof (np.array): the first value in the list of dofs of the surface + subdomain + """ + dofs = fem.locate_dofs_geometrical( + function_space, lambda x: np.isclose(x[0], self.x) + ) + return dofs[0] diff --git a/festim/subdomain/volume_subdomain.py b/festim/subdomain/volume_subdomain.py new file mode 100644 index 000000000..a3398993c --- /dev/null +++ b/festim/subdomain/volume_subdomain.py @@ -0,0 +1,45 @@ +from dolfinx.mesh import locate_entities +import numpy as np + + +class VolumeSubdomain1D: + """ + Volume subdomain class for 1D cases + + Args: + id (int): the id of the volume subdomain + borders (list of float): the borders of the volume subdomain + material (festim.Material): the material of the volume subdomain + + Attributes: + id (int): the id of the volume subdomain + borders (list of float): the borders of the volume subdomain + material (festim.Material): the material of the volume subdomain + + Usage: + >>> vol_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 1], + ... material=F.Material(...)) + """ + + def __init__(self, id, borders, material) -> None: + self.borders = borders + self.material = material + self.id = id + + def locate_subdomain_entities(self, mesh, vdim): + """Locates all cells in subdomain borders within domain + + Args: + mesh (dolfinx.cpp.mesh.Mesh): the mesh of the model + vdim (int): the dimension of the volumes of the mesh, + for 1D this is always 1 + + Returns: + entities (np.array): the entities of the subdomain + """ + entities = locate_entities( + mesh, + vdim, + lambda x: np.logical_and(x[0] >= self.borders[0], x[0] <= self.borders[1]), + ) + return entities diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 3179684ed..87492dec4 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -26,6 +26,11 @@ def test_permeation_problem(): my_model = F.HydrogenTransportProblem() my_model.mesh = my_mesh + my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=None) + left_surface = F.SurfaceSubdomain1D(id=1, x=0) + right_surface = F.SurfaceSubdomain1D(id=2, x=L) + my_model.subdomains = [my_subdomain, left_surface, right_surface] + mobile_H = F.Species("H") my_model.species = [mobile_H] @@ -45,9 +50,9 @@ def siverts_law(T, S_0, E_S, pressure): return S * pressure**0.5 fdim = my_mesh.mesh.topology.dim - 1 - left_facets = my_model.facet_tags.find(1) + left_facets = my_model.facet_meshtags.find(1) left_dofs = locate_dofs_topological(V, fdim, left_facets) - right_facets = my_model.facet_tags.find(2) + right_facets = my_model.facet_meshtags.find(2) right_dofs = locate_dofs_topological(V, fdim, right_facets) S_0 = 4.02e21 @@ -116,10 +121,9 @@ def siverts_law(T, S_0, E_S, pressure): analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1) analytical_flux = np.abs(analytical_flux) - flux_values = np.array(np.abs(flux_values)) - relative_error = (flux_values - analytical_flux) / analytical_flux + relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) relative_error = relative_error[ np.where(analytical_flux > 0.01 * np.max(analytical_flux)) diff --git a/test/test_subdomains.py b/test/test_subdomains.py new file mode 100644 index 000000000..c446804fc --- /dev/null +++ b/test/test_subdomains.py @@ -0,0 +1,48 @@ +import numpy as np +import festim as F + + +def test_different_surface_ids(): + my_model = F.HydrogenTransportProblem() + + L = 3e-02 + my_model.mesh = F.Mesh1D(np.linspace(0, L, num=3)) + + surfacec_subdomains_ids = [3, 8] + surface_subdomain_1 = F.SurfaceSubdomain1D(id=surfacec_subdomains_ids[0], x=0) + surface_subdomain_2 = F.SurfaceSubdomain1D(id=surfacec_subdomains_ids[1], x=L) + my_model.subdomains = [surface_subdomain_1, surface_subdomain_2] + + my_model.define_function_space() + my_model.define_markers_and_measures() + + for surf_id in surfacec_subdomains_ids: + assert surf_id in np.array(my_model.facet_meshtags.values) + + +def test_different_volume_ids(): + my_model = F.HydrogenTransportProblem() + + sub_dom_1 = np.linspace(0, 1e-04, num=3) + sub_dom_2 = np.linspace(1e-04, 2e-04, num=4) + sub_dom_3 = np.linspace(2e-04, 3e-04, num=5) + my_model.mesh = F.Mesh1D( + np.unique(np.concatenate([sub_dom_1, sub_dom_2, sub_dom_3])) + ) + + vol_subdomains_ids = [2, 16, 7] + vol_subdomain_1 = F.VolumeSubdomain1D( + id=vol_subdomains_ids[0], borders=[sub_dom_1[0], sub_dom_1[-1]], material=None + ) + vol_subdomain_2 = F.VolumeSubdomain1D( + id=vol_subdomains_ids[1], borders=[sub_dom_2[0], sub_dom_2[-1]], material=None + ) + vol_subdomain_3 = F.VolumeSubdomain1D( + id=vol_subdomains_ids[2], borders=[sub_dom_3[0], sub_dom_3[-1]], material=None + ) + my_model.subdomains = [vol_subdomain_1, vol_subdomain_2, vol_subdomain_3] + + my_model.define_markers_and_measures() + + for vol_id in vol_subdomains_ids: + assert vol_id in np.array(my_model.volume_meshtags.values)