diff --git a/festim/__init__.py b/festim/__init__.py index 10ce0105c..9a7bc7e59 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -31,7 +31,6 @@ from .meshing.mesh import Mesh from .meshing.mesh_1d import Mesh1D -from .meshing.mesh_from_refinements import MeshFromRefinements from .meshing.mesh_from_vertices import MeshFromVertices from .meshing.mesh_from_xdmf import MeshFromXDMF diff --git a/festim/meshing/mesh_from_refinements.py b/festim/meshing/mesh_from_refinements.py deleted file mode 100644 index f3bf19f3f..000000000 --- a/festim/meshing/mesh_from_refinements.py +++ /dev/null @@ -1,65 +0,0 @@ -import fenics as f -from festim import Mesh1D - - -class MeshFromRefinements(Mesh1D): - """ - 1D mesh with iterative refinements (on the left hand side of the domain) - - Args: - initial_number_of_cells (float): initial number of cells before - refinement - size (float): total size of the 1D mesh - refinements (list, optional): list of dicts - {"x": ..., "cells": ...}. For each refinement, the mesh will - have at least ["cells"] in [0, "x"]. Defaults to []. - start (float, optional): the starting point of the mesh. Defaults to - 0. - - Attributes: - initial_number_of_cells (int): initial number of cells before - refinement - size (float): total size of the 1D mesh - refinements (list): list of refinements - """ - - def __init__( - self, initial_number_of_cells, size, refinements=[], start=0.0, **kwargs - ) -> None: - super().__init__(**kwargs) - self.initial_number_of_cells = initial_number_of_cells - self.size = size - self.start = start - self.refinements = refinements - self.mesh_and_refine() - - def mesh_and_refine(self): - """Mesh and refine iteratively until meeting the refinement - conditions. - """ - - print("Meshing ...") - initial_number_of_cells = self.initial_number_of_cells - size = self.size - mesh = f.IntervalMesh(initial_number_of_cells, self.start, size) - for refinement in self.refinements: - nb_cells_ref = refinement["cells"] - refinement_point = refinement["x"] - print("Mesh size before local refinement is " + str(len(mesh.cells()))) - coarse_mesh = True - while len(mesh.cells()) < initial_number_of_cells + nb_cells_ref: - cell_markers = f.MeshFunction("bool", mesh, mesh.topology().dim()) - cell_markers.set_all(False) - for cell in f.cells(mesh): - if cell.midpoint().x() < refinement_point: - cell_markers[cell] = True - coarse_mesh = False - mesh = f.refine(mesh, cell_markers) - if coarse_mesh: - msg = ( - "Infinite loop: Initial number " + "of cells might be too small" - ) - raise ValueError(msg) - print("Mesh size after local refinement is " + str(len(mesh.cells()))) - initial_number_of_cells = len(mesh.cells()) - self.mesh = mesh diff --git a/test/heat_transfer_problem/test_solving.py b/test/heat_transfer_problem/test_solving.py index 73fe404d9..107d25702 100644 --- a/test/heat_transfer_problem/test_solving.py +++ b/test/heat_transfer_problem/test_solving.py @@ -1,6 +1,7 @@ import festim import pytest import fenics as f +import numpy as np @pytest.mark.parametrize("preconditioner", ["default", "icc"]) @@ -13,7 +14,7 @@ def test_create_functions_linear_solver_gmres(preconditioner): preconditioner (str): the preconditioning method """ - mesh = festim.MeshFromRefinements(10, size=0.1) + mesh = festim.MeshFromVertices(np.linspace(0, 0.1, num=11)) materials = festim.Materials([festim.Material(id=1, D_0=1, E_D=0, thermal_cond=1)]) mesh.define_measures(materials) @@ -60,7 +61,7 @@ def sim(self): def test_custom_solver(self): """Solves the system using the built-in solver and using the f.NewtonSolver""" - mesh = festim.MeshFromRefinements(10, size=0.1) + mesh = festim.MeshFromVertices(np.linspace(0, 0.1, num=11)) materials = festim.Materials( [festim.Material(id=1, D_0=1, E_D=0, thermal_cond=1)] ) diff --git a/test/simulation/test_postprocessing_integration.py b/test/simulation/test_postprocessing_integration.py index 51e3ba517..5897c829c 100644 --- a/test/simulation/test_postprocessing_integration.py +++ b/test/simulation/test_postprocessing_integration.py @@ -12,7 +12,7 @@ def my_sim(self): """A pytest fixture defining the festim.Simulation object""" my_sim = festim.Simulation({}) my_sim.t = 0 - my_sim.mesh = festim.MeshFromRefinements(10, 1) + my_sim.mesh = festim.MeshFromVertices(np.linspace(0, 1, num=11)) my_sim.settings = festim.Settings(None, None, final_time=10) mat1 = festim.Material(1, D_0=1, E_D=1) my_sim.materials = festim.Materials([mat1]) @@ -120,10 +120,12 @@ def test_pure_diffusion(self, my_sim): my_sim.run_post_processing() data = derived_quantities.data + print(data) + assert len(data) == i + 1 assert data[i][0] == t - assert data[i][1] == 10 - assert data[i][2] == 20 + assert data[i][1] == pytest.approx(10) + assert data[i][2] == pytest.approx(20) assert data[i][3] == pytest.approx(11) assert data[i][4] == pytest.approx(11) @@ -259,7 +261,7 @@ def test_surface_kinetics(self): """ my_sim = festim.Simulation({}) my_sim.t = 0 - my_sim.mesh = festim.MeshFromRefinements(10, 1) + my_sim.mesh = festim.MeshFromVertices(np.linspace(0, 1, num=11)) my_sim.settings = festim.Settings(None, None, final_time=10) mat1 = festim.Material(1, D_0=1, E_D=1) my_sim.materials = festim.Materials([mat1]) diff --git a/test/system/test_chemical_potential.py b/test/system/test_chemical_potential.py index 7649095e4..59ba35440 100644 --- a/test/system/test_chemical_potential.py +++ b/test/system/test_chemical_potential.py @@ -71,7 +71,7 @@ def run(h): festim.InitialCondition(field=0, value=u), ] - my_mesh = festim.MeshFromRefinements(round(size / h), size) + my_mesh = festim.MeshFromVertices(np.linspace(0, size, round(size / h) + 1)) my_bcs = [ festim.DirichletBC(surfaces=[1, 2], value=u, field=0), diff --git a/test/system/test_extrinsic_traps.py b/test/system/test_extrinsic_traps.py index 4d357ab1e..01c7365c0 100644 --- a/test/system/test_extrinsic_traps.py +++ b/test/system/test_extrinsic_traps.py @@ -1,10 +1,11 @@ import festim as F +import numpy as np def test_extrinsic_trap(): """Runs a festim sim with an extrinsic trap""" my_materials = F.Materials([F.Material(id=1, D_0=2, E_D=1, name="mat")]) - my_mesh = F.MeshFromRefinements(10, 1) + my_mesh = F.MeshFromVertices(np.linspace(0, 1, 11)) my_traps = F.ExtrinsicTrap( k_0=1, diff --git a/test/system/test_misc.py b/test/system/test_misc.py index 410d02137..a3e999806 100644 --- a/test/system/test_misc.py +++ b/test/system/test_misc.py @@ -45,7 +45,7 @@ def test_convective_flux(tmpdir): def test_error_steady_state_with_stepsize(): """Checks that an error is raised when a stepsize is given for a steady state simulation""" my_model = F.Simulation() - my_model.mesh = F.MeshFromRefinements(1000, size=1) + my_model.mesh = F.MeshFromVertices(np.linspace(0, 1, num=1001)) my_model.materials = F.Materials([F.Material(D_0=1, E_D=0, id=1)]) @@ -67,7 +67,7 @@ def test_error_steady_state_with_stepsize(): def test_error_transient_without_stepsize(): """Checks that an error is raised when a stepsize is not given for a transient simulation""" my_model = F.Simulation() - my_model.mesh = F.MeshFromRefinements(1000, size=1) + my_model.mesh = F.MeshFromVertices(np.linspace(0, 1, num=100)) my_model.materials = F.Materials([F.Material(D_0=1, E_D=0, id=1)]) diff --git a/test/system/test_system.py b/test/system/test_system.py index 98428e00c..f21be0897 100644 --- a/test/system/test_system.py +++ b/test/system/test_system.py @@ -45,7 +45,7 @@ def test_run_temperature_stationary(tmpdir): my_materials = [festim.Material(id=1, D_0=4.1e-7, E_D=0.39, thermal_cond=1)] - my_mesh = festim.MeshFromRefinements(200, size=size) + my_mesh = festim.MeshFromVertices(np.linspace(0, size, num=200)) my_boundary_conditions = [ festim.DirichletBC(value=1, field=0, surfaces=[1]), festim.DirichletBC(value=u, field="T", surfaces=[1, 2]), @@ -115,7 +115,7 @@ def test_run_temperature_transient(tmpdir): ) ] ) - my_mesh = festim.MeshFromRefinements(200, size) + my_mesh = festim.MeshFromVertices(np.linspace(0, size, num=200)) my_bcs = [ festim.DirichletBC(surfaces=[1], value=1, field=0), @@ -216,7 +216,7 @@ def run(h): festim.InitialCondition(field=1, value=v), ] - my_mesh = festim.MeshFromRefinements(round(size / h), size) + my_mesh = festim.MeshFromVertices(np.linspace(0, size, num=round(size / h) + 1)) my_bcs = [ festim.DirichletBC(surfaces=[1, 2], value=u, field=0), @@ -342,7 +342,7 @@ def run(h): festim.InitialCondition(field=1, value=v), ] - my_mesh = festim.MeshFromRefinements(round(size / h), size) + my_mesh = festim.MeshFromVertices(np.linspace(0, size, num=round(size / h) + 1)) my_bcs = [ festim.DirichletBC(surfaces=[1, 2], value=u, field=0), @@ -442,7 +442,7 @@ def test_run_chemical_pot_mass_balance(tmpdir): festim.InitialCondition(field=0, value=1), ] - my_mesh = festim.MeshFromRefinements(5, 1) + my_mesh = festim.MeshFromVertices(np.linspace(0, 1, num=6)) my_temp = festim.Temperature(700 + 210 * festim.t) @@ -512,7 +512,7 @@ def run(h): ] size = 0.1 - my_mesh = festim.MeshFromRefinements(round(size / h), size) + my_mesh = festim.MeshFromVertices(np.linspace(0, size, num=round(size / h) + 1)) my_source = festim.Source(f, 1, "solute") @@ -624,7 +624,7 @@ def run(h): ] size = 0.1 - my_mesh = festim.MeshFromRefinements(round(size / h), size) + my_mesh = festim.MeshFromVertices(np.linspace(0, size, num=round(size / h) + 1)) my_sources = [festim.Source(f, 1, "solute"), festim.Source(g, 1, "1")] @@ -718,7 +718,7 @@ def test_chemical_pot_T_solve_stationary(tmpdir): my_materials = festim.Materials( [festim.Material(id=1, D_0=1, E_D=0.1, S_0=2, E_S=0.2, thermal_cond=1)] ) - my_mesh = festim.MeshFromRefinements(10, 1) + my_mesh = festim.MeshFromVertices(np.linspace(0, 1, num=11)) my_temp = festim.HeatTransferProblem(transient=False) my_bcs = [ @@ -772,7 +772,7 @@ def test_export_particle_flux_with_chemical_pot(tmpdir): my_materials = festim.Materials( [festim.Material(id=1, D_0=2, E_D=1, S_0=2, E_S=1, thermal_cond=2)] ) - my_mesh = festim.MeshFromRefinements(10, 1) + my_mesh = festim.MeshFromVertices(np.linspace(0, 1, num=11)) my_temp = festim.Temperature(300) @@ -885,7 +885,7 @@ def test_steady_state_traps_not_everywhere(): ] ) - my_mesh = festim.MeshFromRefinements(100, 1) + my_mesh = festim.MeshFromVertices(np.linspace(0, 1, num=101)) my_trap = festim.Trap(1, 0, 1, 0, ["mat_1", "mat_3"], 1) @@ -922,7 +922,7 @@ def test_no_jacobian_update(): mat = festim.Material(id=1, D_0=1, E_D=0) my_materials = festim.Materials([mat]) - my_mesh = festim.MeshFromRefinements(10, 1) + my_mesh = festim.MeshFromVertices(np.linspace(0, 1, num=100)) my_trap = festim.Trap(1, 0, 1, 0, [mat], 1) @@ -959,7 +959,7 @@ def test_nb_iterations_bewteen_derived_quantities_compute(): def init_sim(nb_it_compute): my_materials = festim.Materials([festim.Material(id=1, D_0=1, E_D=0)]) - my_mesh = festim.MeshFromRefinements(10, 1) + my_mesh = festim.MeshFromVertices(np.linspace(0, 1, num=11)) my_temp = festim.Temperature(300) @@ -1011,7 +1011,7 @@ def test_error_steady_state_diverges(): ] ) - my_mesh = festim.MeshFromRefinements(10, 1) + my_mesh = festim.MeshFromVertices(np.linspace(0, 1, num=11)) my_temp = festim.Temperature(-1) @@ -1037,7 +1037,7 @@ def test_error_steady_state_diverges(): def test_completion_tone(): """Checks that when a completion tone is enabled, sim is not affected""" my_model = festim.Simulation(log_level=20) - my_model.mesh = festim.MeshFromRefinements(10, size=1) + my_model.mesh = festim.MeshFromVertices(np.linspace(0, 1, num=11)) my_model.materials = festim.Materials([festim.Material(id=1, D_0=1, E_D=0)]) my_model.T = festim.Temperature(100) my_model.boundary_conditions = [ diff --git a/test/unit/test_theta.py b/test/unit/test_theta.py index a00123793..db1b07a0f 100644 --- a/test/unit/test_theta.py +++ b/test/unit/test_theta.py @@ -2,6 +2,7 @@ import fenics as f from ufl.core.multiindex import Index import pytest +import numpy as np class TestInitialise: @@ -125,7 +126,7 @@ def test_get_concentration_for_a_given_material(): E_S = 0.5 my_mat = festim.Material(1, 1, 1, S_0=S_0, E_S=E_S) my_theta = festim.Theta() - my_mesh = festim.MeshFromRefinements(10, 1) + my_mesh = festim.MeshFromVertices(np.linspace(0, 1, 11)) V = f.FunctionSpace(my_mesh.mesh, "CG", 1) my_theta.solution = f.interpolate(f.Constant(100), V) my_theta.previous_solution = f.interpolate(f.Constant(200), V)