Commit 36d352e8 authored by Baptiste Durand's avatar Baptiste Durand

Fix periodicity boundary conditions for cases where the bottom left corner of...

Fix periodicity boundary conditions for cases where the bottom left corner of the domain is not in (0., 0.)
parent a63fd01b
......@@ -54,7 +54,13 @@ class PeriodicDomain(fe.SubDomain):
y[i] = self.infinity[i]
@staticmethod
def pbc_dual_base(part_vectors, per_choice: str, dim=2, tol=GEO_TOLERANCE):
def pbc_dual_base(
part_vectors,
per_choice: str,
bl_corner=np.array((0.0, 0.0)),
dim=2,
tol=GEO_TOLERANCE,
):
"""Create periodic boundary only from an array
that indicate the dimensions of the part.
Appropriate for parallelepipedic domain.
......@@ -64,13 +70,17 @@ class PeriodicDomain(fe.SubDomain):
part_vectors : np.array
shape 2×2. Dimensions of the domain.
Some of them will be used as periodicity vectors.
#TODO : Ajouter précision : Vecteurs de périodicité en colonne
Each column of the array correspond to a periodicity vector.
part_vectors[:,0] -> periodicity vector for the 1st (X) direction
part_vectors[:,1] -> periodicity vector for the 2d (Y) direction
per_choice : str
Can contain X, Y (in the future : Z)
dim : int, optional
Dimension of the modeling space. (the default is 2)
tol : float, optional
geometrical tolerance for membership tests.
bl_corner : np.array
shape 2. Coordinates of the bottom-left corner of the domain.
Returns
-------
......@@ -81,15 +91,19 @@ class PeriodicDomain(fe.SubDomain):
for i in range(np.size(part_vectors, 1)):
basis.append(fe.as_vector(part_vectors[:, i]))
dualbasis.append(fe.as_vector(dual_vect[:, i]))
bl_corner = np.asarray(bl_corner)
assert len(bl_corner.shape) == 1 and bl_corner.shape[0] == dim
master_tests, slave_tests, per_vectors = list(), list(), list()
if "x" in per_choice.lower():
def left(x):
return fe.near(x.dot(dualbasis[0]), 0.0, tol)
# type x : numpy.ndarray
# shape x : (2,)
return fe.near((x - bl_corner).dot(dualbasis[0]), 0.0, tol)
# dot product return a <'ufl.constantvalue.FloatValue'>
def right(x):
return fe.near((x - basis[0]).dot(dualbasis[0]), 0.0, tol)
return fe.near((x - bl_corner - basis[0]).dot(dualbasis[0]), 0.0, tol)
master_tests.append(left)
slave_tests.append(right)
......@@ -97,10 +111,10 @@ class PeriodicDomain(fe.SubDomain):
if "y" in per_choice.lower():
def bottom(x):
return fe.near(x.dot(dualbasis[1]), 0.0, tol)
return fe.near((x - bl_corner).dot(dualbasis[1]), 0.0, tol)
def top(x):
return fe.near((x - basis[1]).dot(dualbasis[1]), 0.0, tol)
return fe.near((x - bl_corner - basis[1]).dot(dualbasis[1]), 0.0, tol)
master_tests.append(bottom)
slave_tests.append(top)
......@@ -235,7 +249,9 @@ class PeriodicExpr(fe.UserExpression):
self.base_function = base_function
self.per_vectors = np.asarray(per_vectors)
assert self.per_vectors.shape == (2, 2), "Invalid shape. Only 2D supported"
diag_test = np.allclose(self.per_vectors[[0, 1], [1, 0]], [0.0, 0.0], atol=1e-12)
diag_test = np.allclose(
self.per_vectors[[0, 1], [1, 0]], [0.0, 0.0], atol=1e-12
)
assert diag_test, "Periodicity vectors must be align with global basis."
# TODO : généralisation à des vecteurs de périodicité non alignés avec le repère.
......
......@@ -111,3 +111,54 @@ def test_pbc_from_vectors_parallelogram():
logger.debug(f"Point : {point}, map pt : {y}, expected : {map_pt}")
assert y == approx(map_pt)
# * OK, test réussi
def test_pbc_from_vectors_bottom_left_corner():
per_vect = np.array([[4.0, 2.0], [0.0, 4.0]])
bl_corner = np.array((1.0, 2.0))
pbc = periodicity.PeriodicDomain.pbc_dual_base(per_vect, "XY", bl_corner=bl_corner)
test_points = [
(0, 0),
(2, 0),
(4, 0),
(1, 2),
(5, 2),
(2, 4),
(4, 4),
(6, 4),
]
test_points = [np.array(coord, np.float64) for coord in test_points]
test_points = [p + bl_corner for p in test_points] # translation of the domain
inside_results = [
True,
True,
False,
True,
False,
False,
False,
False,
]
for point, result in zip(test_points, inside_results):
logger.debug(f"Point : {point}, master : {pbc.inside(point, True)}")
assert pbc.inside(point, on_boundary=True) == result
map_results = [
(59994, 39996),
(59994, 39996),
(0, 0),
(59994, 39996),
(1, 2),
(0, 0),
(2, 0),
(0, 0),
]
map_results = [np.array(coord, np.float64) for coord in map_results]
for i in (2, 4, 5, 6, 7):
map_results[i] += bl_corner
y = np.zeros(2)
for point, map_pt in zip(test_points, map_results):
pbc.map(point, y)
logger.debug(f"Point : {point}, map pt : {y}, expected : {map_pt}")
assert y == approx(map_pt)
# * OK, test réussi
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment