test avg surface on a rhombus

Author: jhdark

test/unit/test_exports/test_derived_quantities/test_average_surface.py

@@ -1,4 +1,10 @@
-from festim import AverageSurface, AverageSurfaceCylindrical, AverageSurfaceSpherical, x
+from festim import (
+    AverageSurface,
+    AverageSurfaceCylindrical,
+    AverageSurfaceSpherical,
+    x,
+    y,
+)
 import fenics as f
 import pytest
 import numpy as np
@@ -174,3 +180,60 @@ def test_average_surface_spherical_title_no_units_temperature():
 
     my_export = AverageSurfaceSpherical("T", 9)
     assert my_export.title == "Average T surface 9 (K)"
+
+
+@pytest.mark.parametrize("j", [4, 5, 3])
+def test_average_surface_rhombus(j):
+
+    # creating a mesh with FEniCS
+    # Define the number of divisions in the x and y directions
+    nx, ny = 10, 10
+    mesh = f.RectangleMesh(f.Point(0, 0), f.Point(j, j), nx, ny)
+
+    # Access mesh coordinates
+    coordinates = mesh.coordinates()
+
+    # Rotation matrix for 45 degrees
+    theta = np.pi / 4  # 45 degrees in radians
+    rotation_matrix = np.array(
+        [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]
+    )
+
+    # Apply rotation to each coordinate
+    for i in range(len(coordinates)):
+        x_cord, y_cord = coordinates[i]
+        # Rotate the point (x, y) using the rotation matrix
+        new_coords = np.dot(rotation_matrix, np.array([x_cord, y_cord]))
+        coordinates[i] = new_coords
+
+    # Create a new mesh with the rotated coordinates
+    mesh_rotated = f.Mesh(mesh)
+
+    surface_markers = f.MeshFunction(
+        "size_t", mesh_rotated, mesh_rotated.topology().dim() - 1
+    )
+    surface_markers.set_all(0)
+
+    # find all facets along y = x and mark them
+    surf_id = 2
+    for facet in f.facets(mesh):
+        midpoint = facet.midpoint()
+        if np.isclose(midpoint[0], midpoint[1], atol=1e-10):
+            surface_markers[facet.index()] = surf_id
+
+    ds = f.Measure("ds", domain=mesh_rotated, subdomain_data=surface_markers)
+    my_export = AverageSurface("solute", surf_id)
+    V = f.FunctionSpace(mesh_rotated, "P", 1)
+
+    c_fun = lambda x, y: 2 * x + y
+    expr = f.Expression(
+        ccode(c_fun(x, y)),
+        degree=1,
+    )
+    my_export.function = f.interpolate(expr, V)
+    my_export.ds = ds
+
+    expected_value = (3 * j) / (2 * (2**0.5))
+    computed_value = my_export.compute()
+
+    assert np.isclose(computed_value, expected_value)