Skip to content

Commit

Permalink
Address review comments
Browse files Browse the repository at this point in the history
  • Loading branch information
RemDelaporteMathurin committed Oct 18, 2023
1 parent af7061a commit 1d441a6
Show file tree
Hide file tree
Showing 4 changed files with 42 additions and 81 deletions.
33 changes: 30 additions & 3 deletions festim/boundary_conditions/dirichlet_bc.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@ class DirichletBC:
species (str): the name of the species
value_fenics (fem.Function or fem.Constant): the value of the boundary condition in
fenics format
bc_expr (fem.Expression): the expression of the boundary condition that is used to
update the value_fenics
Usage:
>>> from festim import DirichletBC
Expand Down Expand Up @@ -60,21 +62,43 @@ def define_surface_subdomain_dofs(self, facet_meshtags, mesh, function_space):
condition
Args:
mesh (festim.Mesh): the domain mesh
facet_meshtags (ddolfinx.mesh.MeshTags): the mesh tags of the
surface facets
mesh (dolfinx.mesh.Mesh): the mesh
function_space (dolfinx.fem.FunctionSpace): the function space
"""
bc_facets = facet_meshtags.find(self.subdomain.id)
bc_dofs = fem.locate_dofs_topological(function_space, mesh.fdim, bc_facets)
return bc_dofs

def create_value(self, mesh, function_space, temperature, t):
def create_value(
self, mesh, function_space: fem.FunctionSpace, temperature, t: fem.Constant
):
"""Creates the value of the boundary condition as a fenics object and sets it to
self.value_fenics.
If the value is a constant, it is converted to a fenics.Constant.
If the value is a function of t, it is converted to a fenics.Constant.
Otherwise, it is converted to a fenics.Function and the
expression of the function is stored in self.bc_expr.
Args:
mesh (dolfinx.mesh.Mesh) : the mesh
function_space (dolfinx.fem.FunctionSpace): the function space
temperature (float): the temperature
t (dolfinx.fem.Constant): the time
"""
x = ufl.SpatialCoordinate(mesh)

if isinstance(self.value, (int, float)):
self.value_fenics = F.as_fenics_constant(mesh=mesh, value=self.value)

elif callable(self.value):
arguments = self.value.__code__.co_varnames

if "t" in arguments and "x" not in arguments and "T" not in arguments:
# only t is an argument
self.value_fenics = F.as_fenics_constant(
mesh=mesh, value=self.value(t=float(t))
mesh=mesh, value=self.value(t=t.value)
)
else:
self.value_fenics = fem.Function(function_space)
Expand All @@ -85,6 +109,9 @@ def create_value(self, mesh, function_space, temperature, t):
kwargs["x"] = x
if "T" in arguments:
kwargs["T"] = temperature

# store the expression of the boundary condition
# to update the value_fenics later
self.bc_expr = fem.Expression(
self.value(**kwargs),
function_space.element.interpolation_points(),
Expand Down
53 changes: 0 additions & 53 deletions festim/helpers.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,4 @@
from dolfinx import fem
import ufl
import numpy as np
from petsc4py.PETSc import ScalarType


def as_fenics_constant(value, mesh):
Expand All @@ -25,53 +22,3 @@ def as_fenics_constant(value, mesh):
raise TypeError(
f"Value must be a float, an int or a dolfinx.Constant, not {type(value)}"
)


# class SpaceTimeDependentExpression:
# def __init__(self, function, t=0):
# self.t = t
# self.function = function
# self.values = None

# def __call__(self, x):
# if self.values is None:
# self.values = np.zeros(x.shape[1], dtype=ScalarType)
# self.values = np.full(x.shape[1], self.function(x=x, t=self.t))
# return self.values


# def convert_to_appropriate_obj(object, function_space, mesh):
# """Converts a value to a dolfinx.Constant or a dolfinx.Function
# depending on the type of the value

# Args:
# object (callable or float): the value to convert
# function_space (dolfinx.fem.FunctionSpace): the function space of the domain
# mesh (dolfinx.mesh.mesh): the mesh of the domain

# Returns:
# dolfinx.Constant or dolfinx.Function: the converted value
# festim.SpaceTimeDependentExpression or None: the expression if the value is
# space and time dependent, None otherwise
# """
# x = ufl.SpatialCoordinate(mesh)
# if isinstance(object, (int, float)):
# # case 1 pressure isn't space dependent or only time dependent:
# return as_fenics_constant(mesh=mesh, value=object), None
# # case 2 pressure is space dependent
# elif callable(object):
# arguments = object.__code__.co_varnames
# if "t" in arguments and "x" in arguments:
# expr = fem.Expression(
# object(x=x, t=0), function_space.element.interpolation_points()
# )
# fenics_obj = fem.Function(function_space)
# fenics_obj.interpolate(expr)
# return fenics_obj, expr
# elif "x" in arguments:
# fenics_obj = fem.Function(function_space)
# fenics_obj.interpolate(object)
# return fenics_obj, None

# elif "t" in arguments:
# return as_fenics_constant(mesh=mesh, value=object(t=0)), None
16 changes: 8 additions & 8 deletions festim/hydrogen_transport_problem.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,8 +41,8 @@ class HydrogenTransportProblem:
ds (dolfinx.fem.ds): the surface measure of the model
function_space (dolfinx.fem.FunctionSpace): the function space 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
facet_meshtags (dolfinx.mesh.MeshTags): the facet tags of the model
volume_meshtags (dolfinx.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
Expand Down Expand Up @@ -283,23 +283,23 @@ def run(self, final_time: float):
progress = tqdm.autonotebook.tqdm(
desc="Solving H transport problem", total=final_time
)
while float(self.t) < final_time:
progress.update(float(self.dt))
self.t.value += float(self.dt)
while self.t.value < final_time:
progress.update(self.dt.value)
self.t.value += self.dt.value

# update boundary conditions
for bc in self.boundary_conditions:
bc.update(float(self.t))
bc.update(self.t.value)

self.solver.solve(self.u)

mobile_xdmf.write_function(self.u, float(self.t))
mobile_xdmf.write_function(self.u, self.t.value)

surface_flux = form(D * dot(grad(cm), n) * self.ds(2))

flux = assemble_scalar(surface_flux)
flux_values.append(flux)
times.append(float(self.t))
times.append(self.t.value)

# update previous solution
self.u_n.x.array[:] = self.u.x.array[:]
Expand Down
21 changes: 4 additions & 17 deletions test/test_dirichlet_bc.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@


def test_init():
"""Test that the attributes are set correctly"""
# create a DirichletBC object
subdomain = F.SurfaceSubdomain1D(1, x=0)
value = 1.0
Expand All @@ -29,6 +30,9 @@ def test_init():


def test_value_fenics():
"""Test that the value_fenics attribute can be set to a valid value
and that an invalid type throws an error
"""
# create a DirichletBC object
subdomain = F.SurfaceSubdomain1D(1, x=0)
value = 1.0
Expand All @@ -47,23 +51,6 @@ def test_value_fenics():
bc.value_fenics = "invalid"


# def test_define_surface_subdomain_dofs():
# # create a DirichletBC object
# subdomain = F.SurfaceSubdomain1D(1, x=0)
# value = 1.0
# species = "test"
# bc = F.DirichletBC(subdomain, value, species)

# # create facet meshtags, mesh, and function space objects
# facet_meshtags = .....
# function_space = fem.FunctionSpace(mesh, ("Lagrange", 1))

# # call the method being tested
# bc_dofs = bc.define_surface_subdomain_dofs(facet_meshtags, mesh, function_space)

# assert bc_dofs ...


def test_callable_for_value():
"""Test that the value attribute can be a callable function of x and t"""

Expand Down

0 comments on commit 1d441a6

Please sign in to comment.