From 0e687f0bd26392c475d301689c9bd4ee6c05aa07 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 9 Dec 2024 11:31:31 +0100
Subject: [PATCH 01/43] Added sneddon testfile

---
 tests/functional/test_sneddon_2d.py | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 tests/functional/test_sneddon_2d.py

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
new file mode 100644
index 0000000000..e69de29bb2

From 96497395e75659927acf389d59d1f011b9418114 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 9 Dec 2024 11:46:10 +0100
Subject: [PATCH 02/43] TST: Added the script for sneddon 2d

---
 tests/functional/test_sneddon_2d.py | 580 ++++++++++++++++++++++++++++
 1 file changed, 580 insertions(+)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index e69de29bb2..5967a2009c 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -0,0 +1,580 @@
+
+from matplotlib import pyplot as plt
+import numpy as np
+import porepy as pp
+from porepy.applications.convergence_analysis import ConvergenceAnalysis
+
+import math
+from numba import config
+import pytest
+
+
+config.DISABLE_JIT = True
+
+
+def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray:
+    """
+    Compute the distance fracture points to the fracture centre.
+
+    Parameter::
+    pointset: Array containing coordinates on the fracture
+    center: fracture centre
+
+    Return:
+        Array of distances each point in pointset to the center
+
+    """
+    return pp.geometry.distances.point_pointset(pointset_centers, center)
+
+
+def get_bem_centers(
+    a: float, h: float, n: int, theta: float, center: np.ndarray
+) -> np.ndarray:
+    """
+    Compute coordinates of the centers of the bem segments
+
+    Parameters:
+        a: half fracture length
+        h: bem segment length
+        n: number of bem segments
+        theta: orientation of the fracture
+        center: center of the fracture
+
+    Return:
+        Array of BEM centers
+    """
+
+    # Coordinate system (4.5.1) in page 57 in in book Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics
+
+    bem_centers = np.zeros((3, n))
+    x_0 = center[0] - (a - 0.5 * h) * np.sin(theta)
+    y_0 = center[1] - (a - 0.5 * h) * np.cos(theta)
+    for i in range(0, n):
+        bem_centers[0, i] = x_0 + i * h * np.sin(theta)
+        bem_centers[1, i] = y_0 + i * h * np.cos(theta)
+
+    return bem_centers
+
+
+def analytical_displacements(
+    a: float, eta: np.ndarray, p0: float, G: float, poi: float
+) -> np.ndarray:
+    """
+    Compute Sneddon's analytical solution for the pressurized fracture.
+
+    Source: Sneddon Fourier transforms 1951 page 425 eq 92
+
+    Parameter:
+        a: half fracture length
+        eta: distance from fracture centre
+        p0: pressure
+        G: shear modulus
+        poi: poisson ratio
+
+    Return
+        List of analytical normal displacement jumps.
+    """
+    cons = (1 - poi) / G * p0 * a * 2
+    return cons * np.sqrt(1 - np.power(eta / a, 2))
+
+
+def transform(xc: np.ndarray, x: np.ndarray, alpha: float) -> np.ndarray:
+    """
+    Translation and rotation coordinate transformation of boundary face coordinates for the BEM method
+
+    Parameter
+        xc: Coordinates of BEM segment centre
+        x: Coordinates of boundary faces
+        alpha: Fracture orientation
+    Return:
+        Transformed coordinates
+    """
+    x_bar = np.zeros_like(x)
+    # Terms in (7.4.6) in Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics page 168
+    x_bar[0, :] = (x[0, :] - xc[0]) * np.cos(alpha) + (x[1, :] - xc[1]) * np.sin(alpha)
+    x_bar[1, :] = -(x[0, :] - xc[0]) * np.sin(alpha) + (x[1, :] - xc[1]) * np.cos(alpha)
+    return x_bar
+
+
+def get_bc_val(
+    g: pp.GridLike,
+    bound_faces: np.ndarray,
+    xf: np.ndarray,
+    h: float,
+    poi: float,
+    alpha: float,
+    du: float,
+) -> np.ndarray:
+    """
+    Compute semi-analytical displacement on the boundary using the BEM method for the Sneddon problem.
+
+    Parameter
+        g: Grid object
+        bound_faces: Index lists for boundary faces
+        xf: Coordinates of boundary faces
+        h: BEM segment length
+        poi: Poisson ratio
+        alpha: Fracture orientation
+        du: Sneddon's analytical relative normal displacement
+    Return:
+        Boundary values for the displacement
+    """
+
+    # Equations for f2,f3,f4.f5 can be found in book Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics pages 57, 84-92, 168
+    f2 = np.zeros(bound_faces.size)
+    f3 = np.zeros(bound_faces.size)
+    f4 = np.zeros(bound_faces.size)
+    f5 = np.zeros(bound_faces.size)
+
+    u = np.zeros((g.dim, g.num_faces))
+
+    # Constant term in (7.4.5)
+    m = 1 / (4 * np.pi * (1 - poi))
+
+    # Second term in (7.4.5)
+    f2[:] = m * (
+        np.log(np.sqrt((xf[0, :] - h) ** 2 + xf[1] ** 2))
+        - np.log(np.sqrt((xf[0, :] + h) ** 2 + xf[1] ** 2))
+    )
+
+    # First term in (7.4.5)
+    f3[:] = -m * (
+        np.arctan2(xf[1, :], (xf[0, :] - h)) - np.arctan2(xf[1, :], (xf[0, :] + h))
+    )
+
+    # The following equalities can be found on page 91 where f3,f4 as equation (5.5.3) and ux, uy in (5.5.1)
+    # Also this is in the coordinate system described in (4.5.1) at page 57, which essentially is a rotation.
+    f4[:] = m * (
+        xf[1, :] / ((xf[0, :] - h) ** 2 + xf[1, :] ** 2)
+        - xf[1, :] / ((xf[0, :] + h) ** 2 + xf[1, :] ** 2)
+    )
+
+    f5[:] = m * (
+        (xf[0, :] - h) / ((xf[0, :] - h) ** 2 + xf[1, :] ** 2)
+        - (xf[0, :] + h) / ((xf[0, :] + h) ** 2 + xf[1, :] ** 2)
+    )
+
+    u[0, bound_faces] = du * (
+        -(1 - 2 * poi) * np.cos(alpha) * f2[:]
+        - 2 * (1 - poi) * np.sin(alpha) * f3[:]
+        - xf[1, :] * (np.cos(alpha) * f4[:] + np.sin(alpha) * f5[:])
+    )
+    u[1, bound_faces] = du * (
+        -(1 - 2 * poi) * np.sin(alpha) * f2[:]
+        + 2 * (1 - poi) * np.cos(alpha) * f3[:]
+        - xf[1, :] * (np.sin(alpha) * f4[:] - np.cos(alpha) * f5[:])
+    )
+
+    return u
+
+
+def assign_bem(
+    g: pp.Grid,
+    h: float,
+    bound_faces: np.ndarray,
+    theta: float,
+    bem_centers: np.ndarray,
+    u_a: np.ndarray,
+    poi: float,
+) -> np.ndarray:
+    """
+    Compute analytical displacement using the BEM method for the Sneddon problem.
+
+    Parameter
+        g: Subdomain grid
+        h: bem segment length
+        bound_faces: boundary faces
+        theta: fracture orientation
+        bem_centers: bem segments centers
+        u_a: Sneddon's analytical relative normal displacement
+        poi: Poisson ratio
+
+    Return:
+        Semi-analytical boundary displacement values
+    """
+
+    bc_val = np.zeros((g.dim, g.num_faces))
+
+    alpha = np.pi / 2 - theta
+
+    bound_face_centers = g.face_centers[:, bound_faces]
+
+    for i in range(0, u_a.size):
+
+        new_bound_face_centers = transform(bem_centers[:, i], bound_face_centers, alpha)
+
+        u_bound = get_bc_val(
+            g, bound_faces, new_bound_face_centers, h, poi, alpha, u_a[i]
+        )
+
+        bc_val += u_bound
+
+    return bc_val
+
+
+def run_analytical_displacements(
+    gb: pp.GridLike,
+    a: float,
+    p0: float,
+    G: float,
+    poi: float,
+    height: float,
+    length: float,
+) -> tuple:
+    """
+    Compute Sneddon's analytical solution for the pressurized crack
+    problem in question.
+
+    Source: Sneddon Fourier transforms 1951 page 425 eq 92
+
+    Parameter
+        gb: Grid object
+        a: Half fracture length of fracture
+        eta: Distance from fracture centre
+        p0: Internal constant pressure inside the fracture
+        G: Shear modulus
+        poi: Poisson ratio
+        height,length: Height and length of domain
+
+    Return:
+        A tuple containing two vectors: a list of distances from the fracture center to each fracture coordinate, and the corresponding analytical apertures.
+    """
+    ambient_dim = gb.dim_max()
+    g_1 = gb.subdomains(dim=ambient_dim - 1)[0]
+    fracture_center = np.array([length / 2, height / 2, 0])
+
+    fracture_faces = g_1.cell_centers
+
+    # compute distances from fracture centre with its corresponding apertures
+    eta = compute_eta(fracture_faces, fracture_center)
+    apertures = analytical_displacements(a, eta, p0, G, poi)
+
+    return eta, apertures
+
+
+class ModifiedGeometry:
+    def set_domain(self) -> None:
+        """Defining a two-dimensional square domain with sidelength 2."""
+
+        self._domain = pp.Domain(
+            bounding_box={
+                "xmin": 0,
+                "ymin": 0,
+                "xmax": self.params["length"],
+                "ymax": self.params["height"],
+            }
+        )
+
+    def grid_type(self) -> str:
+        """Choosing the grid type for our domain."""
+        return self.params.get("grid_type", "simplex")
+
+    def set_fractures(self):
+
+        points = self.params["frac_pts"]
+        self._fractures = [pp.LineFracture(points)]
+
+
+class ModifiedBoundaryConditions:
+    def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
+        """Set boundary condition type for the problem."""
+        bounds = self.domain_boundary_sides(sd)
+
+        # Set the type of west and east boundaries to Dirichlet. North and south are
+        # Neumann by default.
+        bc = pp.BoundaryConditionVectorial(sd, bounds.all_bf, "dir")
+
+        frac_face = sd.tags["fracture_faces"]
+
+        bc.is_neu[:, frac_face] = False
+        bc.is_dir[:, frac_face] = True
+
+        return bc
+
+    def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
+        """Setting displacement boundary condition values.
+
+        This method returns an array of boundary condition values with the value 5t for
+        western boundaries and ones for the eastern boundary.
+
+        """
+
+        a, p0, G, poi = (
+            self.params["a"],
+            self.params["p0"],
+            self.params["material_constants"]["solid"].shear_modulus,
+            self.params["poi"],
+        )
+        theta, length, height = (
+            self.params["theta"],
+            self.params["length"],
+            self.params["height"],
+        )
+
+        sd = bg.parent
+        # bg = self.mdg.subdomain_to_boundary_grid(sd) # technique to extract bg via sd
+
+        if sd.dim < 2:
+            return np.zeros(self.nd * sd.num_faces)
+        box_faces = sd.get_boundary_faces()
+
+        # Set the boundary values
+        u_bc = np.zeros((sd.dim, sd.num_faces))
+
+        # apply sneddon analytical solution through BEM method
+        n = 1000
+        h = 2 * a / n
+        center = np.array([length / 2, height / 2, 0])
+        bem_centers = get_bem_centers(a, h, n, theta, center)
+        eta = compute_eta(bem_centers, center)
+
+        u_a = -analytical_displacements(a, eta, p0, G, poi)
+
+        u_bc = assign_bem(sd, h / 2, box_faces, theta, bem_centers, u_a, poi)
+
+        return bg.projection(2) @ u_bc.ravel("F")
+
+
+class PressureStressMixin:
+    def pressure(self, subdomains: pp.Grid):
+        # Return pressure in the fracture domain
+        mat = self.params["p0"] * np.ones(
+            sum(subdomain.num_cells for subdomain in subdomains)
+        )
+        return pp.ad.DenseArray(mat)
+
+    def stress_discretization(
+        self, subdomains: list[pp.Grid]
+    ) -> pp.ad.BiotAd | pp.ad.MpsaAd:
+        """Discretization of the stress tensor.
+
+        Parameters:
+            subdomains: List of subdomains where the stress is defined.
+
+        Returns:
+            Discretization operator for the stress tensor.
+
+        """
+        return pp.ad.MpsaAd(self.stress_keyword, subdomains)
+
+
+class MomentumBalanceGeometryBC(
+    PressureStressMixin,
+    ModifiedGeometry,
+    ModifiedBoundaryConditions,
+    pp.constitutive_laws.PressureStress,
+    pp.MomentumBalance,
+):
+    ...
+    """Adding geometry and modified boundary conditions to the default model."""
+
+
+def simulation(
+    frac_pts: np.ndarray,
+    theta_rad: np.ndarray,
+    h: np.ndarray,
+    a: float = 1.5,
+    height: float = 1.0,
+    length: float = 1.0,
+    p0: float = 1e-4,
+    G: float = 1.0,
+    poi: float = 0.25,
+) -> float:
+    """
+    Simulates a 2D fracture linear elasticity problem using a momentum balance model and
+    computes the relative L2 error between numerical and analytical displacement jumps.
+
+    Parameters:
+        frac_pts: Array of coordinates specifying the fracture points in the domain.
+        theta_rad: Array of angles (in radians) defining the orientation of the fractures.
+        h: Array of cell sizes used for meshing the domain.
+        a:  Half length of the fracture
+        height: Height of the domain.
+        length: Length of the domain.
+        p0: Initial constant pressure applied inside the fracture.
+        G: Shear modulus of the material.
+        poi: Poisson's ratio of the material.
+
+    Returns:
+        e: The relative L2 error between the analytical and numerical displacement jumps
+        on the fracture.
+    """
+    lam = (
+        2 * G * poi / (1 - 2 * poi)
+    )  # Convertion formula from shear modulus and poission to lame lambda parameter
+    solid = pp.SolidConstants(shear_modulus=G, lame_lambda=lam)
+
+    # Clean this up!! This is made so I can easily access the params in arbitary functions
+    params = {
+        "meshing_arguments": {"cell_size": h},
+        "prepare_simulation": True,
+        "material_constants": {"solid": solid},
+        "frac_pts": frac_pts,
+        "theta": theta_rad,
+        "p0": p0,
+        "a": a,
+        "poi": poi,
+        "length": length,  # this info can be accessed via box, so remove these as well.
+        "height": height,
+    }
+
+    model = MomentumBalanceGeometryBC(params)
+    pp.run_time_dependent_model(model, params)
+
+    frac_sd = model.mdg.subdomains(dim=model.nd - 1)
+    nd_vec_to_normal = model.normal_component(frac_sd)
+
+    # Comuting the numerical displacement jump along the fracture on the fracture cell centers.
+    u_n: pp.ad.Operator = nd_vec_to_normal @ model.displacement_jump(frac_sd)
+    u_n = u_n.value(model.equation_system)
+
+    # Checking convergence specifically on the fracture
+    u_a = run_analytical_displacements(model.mdg, a, p0, G, poi, height, length)[1]
+    e = ConvergenceAnalysis.l2_error(
+        frac_sd[0], u_a, u_n, is_scalar=False, is_cc=True, relative=True
+    )
+
+    return e
+
+
+def compute_frac_pts(
+    theta_rad: float, a: float, height: float, length: float
+) -> np.ndarray:
+    """
+    Assuming the fracture center is at the coordinate (height/2, length/2),
+    compute the endpoints of a fracture given its orientation and fracture length.
+
+    Parameters:
+    theta_rad: Angle of the fracture in radians
+    a: Half-length of the fracture.
+    height: Height of the domain.
+    length: Width of the domain.
+
+    Returns:
+
+        frac_pts : A 2x2 array where each column represents the coordinates of an end point of the fracture in 2D.
+            The first column corresponds to one end point, and the second column corresponds to the other.
+
+    """
+    # Rotate the fracture with an angle theta_rad
+    y_0 = height / 2 - a * np.cos(theta_rad)
+    x_0 = length / 2 - a * np.sin(theta_rad)
+    y_1 = height / 2 + a * np.cos(theta_rad)
+    x_1 = length / 2 + a * np.sin(theta_rad)
+
+    frac_pts = np.array([[x_0, y_0], [x_1, y_1]]).T
+    return frac_pts
+
+
+def compute_eoc(h: np.ndarray, e: np.ndarray) -> np.ndarray:
+    """
+    Compute the Error Order of Convergence (EOC) based on error and mesh size.
+
+    Parameters:
+        h: Array of mesh sizes corresponding to simulation results.
+        e: Array of error values for the corresponding mesh sizes.
+
+    Returns:
+        eoc: Array of EOC values computed for consecutive levels of refinement.
+    """
+
+    # Compute the Error Order of Convergence (EOC)
+    eoc = np.zeros(len(e) - 1)
+    for i in range(len(e) - 1):
+        eoc[i] = np.log(e[i + 1] / e[i]) / np.log(h[i + 1] / h[i])
+    return eoc
+
+
+def compute_convergence(
+    h_list: np.ndarray,
+    theta_list: np.ndarray,
+    a: float = 1,
+    height: float = 1,
+    length: float = 1,
+) -> np.ndarray:
+    """
+    Compute the convergence behavior for the Sneddon problem across mesh refinements and orientations.
+
+    Parameters:
+        h_list: Array of mesh sizes for simulations.
+        theta_list: Array of fracture orientations (in degrees).
+        a: Half-length of the fracture.
+        height: Height of the simulation domain.
+        length: Length of the simulation domain.
+
+    Returns:
+        A 2D array where rows correspond to different orientations and columns to mesh sizes.
+    """  
+    # Compute error for each orientation of the fracture 
+    err = np.zeros((len(theta_list), len(h_list)))
+    for k in range(0, len(theta_list)):
+        
+        theta_rad = math.radians(90 - theta_list[k])
+        frac_pts = compute_frac_pts(theta_rad, a, height=height, length=length)
+        for i in range(0, len(h_list)):
+            e = simulation(
+                frac_pts, theta_rad, h_list[i], a=a, height=height, length=length
+            )
+            err[k, i] = e
+
+    return err
+
+
+def test_sneddon_2d():
+    """
+    Test function to 2D Sneddon problem convergence for the MPSA method.
+    
+    This is a setup for comparing the analytical solution (also known as Sneddon's solution) with the numerical solution for linear elasticity
+    in the case of an open fracture subjected to a constant pressure p0. Fluid effects are not considered.
+    For reference about the implementation, see Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics Chapter 5.3 Pressurized crack problem.
+    Also like to note that this is a reimplementation of Sneddon solutation described in https://doi.org/10.1007/s10596-020-10002-5.
+
+    The tests performs a convergence study for a 2D fracture problem by iterating over fracture orientations and
+    mesh refinements. It computes the average error across orientations and evaluates the convergence rate (EOC)
+    through log-log linear regression.
+    
+    NOTE: We assume that the fracture half-length a = 1 and the domain is (0,1)^2 for the following reasons:
+    1. Convergence tends to be problematic if the domain is changed from (0,1) to another size configuration,
+       suggesting the presence of a bug.
+    2. While a total fracture length significantly larger than the domain shows convergence, convergence issues arise
+       when the fracture length is smaller, potentially due to numerical artifacts that occur when the fracture tip
+       approaches the boundary of the simulation domain.
+
+    Raises:
+    -------
+    AssertionError
+        If the estimated EOC is not greater than 1.0.
+    """
+
+    # Simulation settings
+    num_refs = 4  # Number of refinements
+    num_theta = 6  # Number of iteration of theta
+
+    theta0 = 0  # Inital orientation of fracture
+    delta_theta = 10  # The difference in theta between each refinement
+    theta_list = (
+        theta0 + np.arange(0, num_theta) * delta_theta
+    )  # List of all the fracture orientations we run simulation on
+
+    # Mesh refinement
+    nf0 = 6
+    nf = nf0 * 2 ** np.arange(0, num_refs)
+    h_list = 2 * 1.5 / nf  # Mesh size for simulations
+
+    # In this test we define the
+    height, length = 1, 1
+
+    # Here we define the half length of the fracture to be sufficientl big to avoid any fracture tip effects
+    a = 1
+
+    # Run convergence experiment
+    err = compute_convergence(h_list, theta_list, a=a, height=height, length=length)
+
+    # Computing average error of each theta
+    avg_err = np.mean(err, axis=0)
+
+    # Perform linear regression on the log-log scale and compute EOC
+    log_h = np.log(h_list)
+    log_avg_err = np.log(avg_err)
+    slope, _ = np.polyfit(log_h, log_avg_err, 1)
+    assert slope > 1.0, f"Estimated EOC is {slope}, which is not greater than 1"

From c483197a4b6a4769b89da8cf87b03c75884e830b Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Tue, 7 Jan 2025 11:07:57 +0100
Subject: [PATCH 03/43] Rename test to base

---
 tests/functional/{test_sneddon_2d.py => base_sneddon_2d.py} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename tests/functional/{test_sneddon_2d.py => base_sneddon_2d.py} (100%)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/base_sneddon_2d.py
similarity index 100%
rename from tests/functional/test_sneddon_2d.py
rename to tests/functional/base_sneddon_2d.py

From 806e4204e284c7cbbbfe2a3a01947499b781148a Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Tue, 7 Jan 2025 15:20:43 +0100
Subject: [PATCH 04/43] Added initial structure of the test

---
 tests/functional/test_sneddon_2d.py | 93 +++++++++++++++++++++++++++++
 1 file changed, 93 insertions(+)
 create mode 100644 tests/functional/test_sneddon_2d.py

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
new file mode 100644
index 0000000000..01cf79fd5b
--- /dev/null
+++ b/tests/functional/test_sneddon_2d.py
@@ -0,0 +1,93 @@
+from porepy.applications.convergence_analysis import ConvergenceAnalysis
+
+
+# ----> Retrieve actual order of convergence
+@pytest.fixture(scope="module")
+def actual_ooc(
+    material_constants: dict, reference_values: pp.ReferenceVariableValues
+) -> list[list[dict[str, float]]]:
+    """Retrieve actual order of convergence.
+
+    """
+    ooc = []
+    # Loop through the models
+    ooc_setup: list[dict[str, float]] = []
+    # Loop through grid type
+    # We do not perform a convergence analysis with simplices in 3d
+    grid_type = "simplex"
+    
+        
+    conv_analysis = ConvergenceAnalysis(
+        model_class=model,
+        model_params=deepcopy(params),
+        levels=4,
+        spatial_refinement_rate=2,
+        temporal_refinement_rate=4,
+    )
+    ooc_results = conv_analysis.order_of_convergence(conv_analysis.run_analysis())
+    ooc.append(ooc_results)
+
+    return ooc
+
+
+# ----> Set desired order of convergence
+@pytest.fixture(scope="module")
+def desired_ooc() -> list[list[dict[str, float]]]:
+    """Set desired order of convergence.
+
+    Returns:
+        List of lists of dictionaries, containing the desired order of convergence.
+
+    """
+    desired_ooc_2d = 1.8757689483196147        
+ 
+    return desired_ooc_2d
+
+
+@pytest.mark.skipped  # reason: slow
+@pytest.mark.parametrize("grid_type_idx", [0, 1])
+def test_order_of_convergence(
+    var: str,
+    dim_idx: int,
+    grid_type_idx: int,
+    actual_ooc: list[list[dict[str, float]]],
+    desired_ooc: list[list[dict[str, float]]],
+) -> None:
+    """Test observed order of convergence.
+
+    Note:
+        We set more flexible tolerances for simplicial grids compared to Cartesian
+        grids. This is because we would like to allow for slight changes in the
+        order of convergence if the meshes change, i.e., in newer versions of Gmsh.
+
+    Parameters:
+        var: Name of the variable to be tested.
+        grid_type_idx: Index to identify the grid type; `0` for Cartesian, and `1`
+            for simplices.
+        dim_idx: Index to identify the dimensionality of the problem; `0` for 2d, and
+            `1` for 3d.
+        actual_ooc: List of lists of dictionaries containing the actual observed
+            order of convergence.
+        desired_ooc: List of lists of dictionaries containing the desired observed
+            order of convergence.
+
+    """
+    # We require the order of convergence to always be larger than 1.0
+    if not (dim_idx == 1 and grid_type_idx == 1):  # no analysis for 3d and simplices
+        assert 1.0 < actual_ooc[dim_idx][grid_type_idx]["ooc_" + var]
+
+    if grid_type_idx == 0:  # Cartesian
+        assert np.isclose(
+            desired_ooc[dim_idx][grid_type_idx]["ooc_" + var],
+            actual_ooc[dim_idx][grid_type_idx]["ooc_" + var],
+            atol=1e-3,  # allow for an absolute difference of 0.001 in OOC
+            rtol=1e-3,  # allow for 0.1% of relative difference in OOC
+        )
+    else:  # Simplex
+        if dim_idx == 0:  # no analysis for 3d and simplices
+            assert np.isclose(
+                desired_ooc[dim_idx][grid_type_idx]["ooc_" + var],
+                actual_ooc[dim_idx][grid_type_idx]["ooc_" + var],
+                atol=1e-1,  # allow for an absolute difference of 0.1 in OOC
+                rtol=5e-1,  # allow for 5% of relative difference in OOC
+            )

From 9f650648cea3d1b65625a53e8b45b0339b757d1f Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Sun, 19 Jan 2025 21:55:02 +0100
Subject: [PATCH 05/43] Minimalistic setup

---
 tests/functional/base_sneddon_2d.py |  4 +++
 tests/functional/test_sneddon_2d.py | 54 +++++++++++------------------
 2 files changed, 25 insertions(+), 33 deletions(-)

diff --git a/tests/functional/base_sneddon_2d.py b/tests/functional/base_sneddon_2d.py
index 5967a2009c..d7898216b1 100644
--- a/tests/functional/base_sneddon_2d.py
+++ b/tests/functional/base_sneddon_2d.py
@@ -578,3 +578,7 @@ def test_sneddon_2d():
     log_avg_err = np.log(avg_err)
     slope, _ = np.polyfit(log_h, log_avg_err, 1)
     assert slope > 1.0, f"Estimated EOC is {slope}, which is not greater than 1"
+
+
+
+
diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 01cf79fd5b..dade774f73 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -1,21 +1,22 @@
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
-
+import porepy as pp
+import numpy as np
+import pytest
 
 # ----> Retrieve actual order of convergence
 @pytest.fixture(scope="module")
 def actual_ooc(
-    material_constants: dict, reference_values: pp.ReferenceVariableValues
-) -> list[list[dict[str, float]]]:
+) -> float:
     """Retrieve actual order of convergence.
 
     """
+    
     ooc = []
     # Loop through the models
     ooc_setup: list[dict[str, float]] = []
     # Loop through grid type
     # We do not perform a convergence analysis with simplices in 3d
     grid_type = "simplex"
-    
         
     conv_analysis = ConvergenceAnalysis(
         model_class=model,
@@ -24,34 +25,30 @@ def actual_ooc(
         spatial_refinement_rate=2,
         temporal_refinement_rate=4,
     )
-    ooc_results = conv_analysis.order_of_convergence(conv_analysis.run_analysis())
-    ooc.append(ooc_results)
+    order = conv_analysis.order_of_convergence(conv_analysis.run_analysis())
+    
 
-    return ooc
+    return order 
 
 
 # ----> Set desired order of convergence
 @pytest.fixture(scope="module")
-def desired_ooc() -> list[list[dict[str, float]]]:
+def desired_ooc() -> float:
     """Set desired order of convergence.
 
     Returns:
         List of lists of dictionaries, containing the desired order of convergence.
 
     """
-    desired_ooc_2d = 1.8757689483196147        
+    desired_ooc_2d = 2.0      
  
     return desired_ooc_2d
 
 
-@pytest.mark.skipped  # reason: slow
-@pytest.mark.parametrize("grid_type_idx", [0, 1])
+
 def test_order_of_convergence(
-    var: str,
-    dim_idx: int,
-    grid_type_idx: int,
-    actual_ooc: list[list[dict[str, float]]],
-    desired_ooc: list[list[dict[str, float]]],
+    actual_ooc,
+    desired_ooc,
 ) -> None:
     """Test observed order of convergence.
 
@@ -73,21 +70,12 @@ def test_order_of_convergence(
 
     """
     # We require the order of convergence to always be larger than 1.0
-    if not (dim_idx == 1 and grid_type_idx == 1):  # no analysis for 3d and simplices
-        assert 1.0 < actual_ooc[dim_idx][grid_type_idx]["ooc_" + var]
+    assert 1.0 < actual_ooc
+
 
-    if grid_type_idx == 0:  # Cartesian
-        assert np.isclose(
-            desired_ooc[dim_idx][grid_type_idx]["ooc_" + var],
-            actual_ooc[dim_idx][grid_type_idx]["ooc_" + var],
-            atol=1e-3,  # allow for an absolute difference of 0.001 in OOC
-            rtol=1e-3,  # allow for 0.1% of relative difference in OOC
-        )
-    else:  # Simplex
-        if dim_idx == 0:  # no analysis for 3d and simplices
-            assert np.isclose(
-                desired_ooc[dim_idx][grid_type_idx]["ooc_" + var],
-                actual_ooc[dim_idx][grid_type_idx]["ooc_" + var],
-                atol=1e-1,  # allow for an absolute difference of 0.1 in OOC
-                rtol=5e-1,  # allow for 5% of relative difference in OOC
-            )
+    assert np.isclose(
+        desired_ooc,
+        actual_ooc,
+        atol=1e-1,  # allow for an absolute difference of 0.1 in OOC
+        rtol=5e-1,  # allow for 5% of relative difference in OOC
+    )

From eb254550baae8229efca64cf4837c5aad93e8997 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 20 Jan 2025 10:05:28 +0100
Subject: [PATCH 06/43] Constructed setup class for manufactured solution

---
 tests/functional/base_sneddon_2d.py        | 258 +--------------------
 tests/functional/setups/manu_sneddon_2d.py | 241 +++++++++++++++++++
 2 files changed, 248 insertions(+), 251 deletions(-)
 create mode 100644 tests/functional/setups/manu_sneddon_2d.py

diff --git a/tests/functional/base_sneddon_2d.py b/tests/functional/base_sneddon_2d.py
index d7898216b1..a077e5c338 100644
--- a/tests/functional/base_sneddon_2d.py
+++ b/tests/functional/base_sneddon_2d.py
@@ -7,250 +7,10 @@
 import math
 from numba import config
 import pytest
-
-
+import  tests.functional.setups.manu_sneddon_2d as manu_sneddon_2d
 config.DISABLE_JIT = True
 
 
-def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray:
-    """
-    Compute the distance fracture points to the fracture centre.
-
-    Parameter::
-    pointset: Array containing coordinates on the fracture
-    center: fracture centre
-
-    Return:
-        Array of distances each point in pointset to the center
-
-    """
-    return pp.geometry.distances.point_pointset(pointset_centers, center)
-
-
-def get_bem_centers(
-    a: float, h: float, n: int, theta: float, center: np.ndarray
-) -> np.ndarray:
-    """
-    Compute coordinates of the centers of the bem segments
-
-    Parameters:
-        a: half fracture length
-        h: bem segment length
-        n: number of bem segments
-        theta: orientation of the fracture
-        center: center of the fracture
-
-    Return:
-        Array of BEM centers
-    """
-
-    # Coordinate system (4.5.1) in page 57 in in book Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics
-
-    bem_centers = np.zeros((3, n))
-    x_0 = center[0] - (a - 0.5 * h) * np.sin(theta)
-    y_0 = center[1] - (a - 0.5 * h) * np.cos(theta)
-    for i in range(0, n):
-        bem_centers[0, i] = x_0 + i * h * np.sin(theta)
-        bem_centers[1, i] = y_0 + i * h * np.cos(theta)
-
-    return bem_centers
-
-
-def analytical_displacements(
-    a: float, eta: np.ndarray, p0: float, G: float, poi: float
-) -> np.ndarray:
-    """
-    Compute Sneddon's analytical solution for the pressurized fracture.
-
-    Source: Sneddon Fourier transforms 1951 page 425 eq 92
-
-    Parameter:
-        a: half fracture length
-        eta: distance from fracture centre
-        p0: pressure
-        G: shear modulus
-        poi: poisson ratio
-
-    Return
-        List of analytical normal displacement jumps.
-    """
-    cons = (1 - poi) / G * p0 * a * 2
-    return cons * np.sqrt(1 - np.power(eta / a, 2))
-
-
-def transform(xc: np.ndarray, x: np.ndarray, alpha: float) -> np.ndarray:
-    """
-    Translation and rotation coordinate transformation of boundary face coordinates for the BEM method
-
-    Parameter
-        xc: Coordinates of BEM segment centre
-        x: Coordinates of boundary faces
-        alpha: Fracture orientation
-    Return:
-        Transformed coordinates
-    """
-    x_bar = np.zeros_like(x)
-    # Terms in (7.4.6) in Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics page 168
-    x_bar[0, :] = (x[0, :] - xc[0]) * np.cos(alpha) + (x[1, :] - xc[1]) * np.sin(alpha)
-    x_bar[1, :] = -(x[0, :] - xc[0]) * np.sin(alpha) + (x[1, :] - xc[1]) * np.cos(alpha)
-    return x_bar
-
-
-def get_bc_val(
-    g: pp.GridLike,
-    bound_faces: np.ndarray,
-    xf: np.ndarray,
-    h: float,
-    poi: float,
-    alpha: float,
-    du: float,
-) -> np.ndarray:
-    """
-    Compute semi-analytical displacement on the boundary using the BEM method for the Sneddon problem.
-
-    Parameter
-        g: Grid object
-        bound_faces: Index lists for boundary faces
-        xf: Coordinates of boundary faces
-        h: BEM segment length
-        poi: Poisson ratio
-        alpha: Fracture orientation
-        du: Sneddon's analytical relative normal displacement
-    Return:
-        Boundary values for the displacement
-    """
-
-    # Equations for f2,f3,f4.f5 can be found in book Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics pages 57, 84-92, 168
-    f2 = np.zeros(bound_faces.size)
-    f3 = np.zeros(bound_faces.size)
-    f4 = np.zeros(bound_faces.size)
-    f5 = np.zeros(bound_faces.size)
-
-    u = np.zeros((g.dim, g.num_faces))
-
-    # Constant term in (7.4.5)
-    m = 1 / (4 * np.pi * (1 - poi))
-
-    # Second term in (7.4.5)
-    f2[:] = m * (
-        np.log(np.sqrt((xf[0, :] - h) ** 2 + xf[1] ** 2))
-        - np.log(np.sqrt((xf[0, :] + h) ** 2 + xf[1] ** 2))
-    )
-
-    # First term in (7.4.5)
-    f3[:] = -m * (
-        np.arctan2(xf[1, :], (xf[0, :] - h)) - np.arctan2(xf[1, :], (xf[0, :] + h))
-    )
-
-    # The following equalities can be found on page 91 where f3,f4 as equation (5.5.3) and ux, uy in (5.5.1)
-    # Also this is in the coordinate system described in (4.5.1) at page 57, which essentially is a rotation.
-    f4[:] = m * (
-        xf[1, :] / ((xf[0, :] - h) ** 2 + xf[1, :] ** 2)
-        - xf[1, :] / ((xf[0, :] + h) ** 2 + xf[1, :] ** 2)
-    )
-
-    f5[:] = m * (
-        (xf[0, :] - h) / ((xf[0, :] - h) ** 2 + xf[1, :] ** 2)
-        - (xf[0, :] + h) / ((xf[0, :] + h) ** 2 + xf[1, :] ** 2)
-    )
-
-    u[0, bound_faces] = du * (
-        -(1 - 2 * poi) * np.cos(alpha) * f2[:]
-        - 2 * (1 - poi) * np.sin(alpha) * f3[:]
-        - xf[1, :] * (np.cos(alpha) * f4[:] + np.sin(alpha) * f5[:])
-    )
-    u[1, bound_faces] = du * (
-        -(1 - 2 * poi) * np.sin(alpha) * f2[:]
-        + 2 * (1 - poi) * np.cos(alpha) * f3[:]
-        - xf[1, :] * (np.sin(alpha) * f4[:] - np.cos(alpha) * f5[:])
-    )
-
-    return u
-
-
-def assign_bem(
-    g: pp.Grid,
-    h: float,
-    bound_faces: np.ndarray,
-    theta: float,
-    bem_centers: np.ndarray,
-    u_a: np.ndarray,
-    poi: float,
-) -> np.ndarray:
-    """
-    Compute analytical displacement using the BEM method for the Sneddon problem.
-
-    Parameter
-        g: Subdomain grid
-        h: bem segment length
-        bound_faces: boundary faces
-        theta: fracture orientation
-        bem_centers: bem segments centers
-        u_a: Sneddon's analytical relative normal displacement
-        poi: Poisson ratio
-
-    Return:
-        Semi-analytical boundary displacement values
-    """
-
-    bc_val = np.zeros((g.dim, g.num_faces))
-
-    alpha = np.pi / 2 - theta
-
-    bound_face_centers = g.face_centers[:, bound_faces]
-
-    for i in range(0, u_a.size):
-
-        new_bound_face_centers = transform(bem_centers[:, i], bound_face_centers, alpha)
-
-        u_bound = get_bc_val(
-            g, bound_faces, new_bound_face_centers, h, poi, alpha, u_a[i]
-        )
-
-        bc_val += u_bound
-
-    return bc_val
-
-
-def run_analytical_displacements(
-    gb: pp.GridLike,
-    a: float,
-    p0: float,
-    G: float,
-    poi: float,
-    height: float,
-    length: float,
-) -> tuple:
-    """
-    Compute Sneddon's analytical solution for the pressurized crack
-    problem in question.
-
-    Source: Sneddon Fourier transforms 1951 page 425 eq 92
-
-    Parameter
-        gb: Grid object
-        a: Half fracture length of fracture
-        eta: Distance from fracture centre
-        p0: Internal constant pressure inside the fracture
-        G: Shear modulus
-        poi: Poisson ratio
-        height,length: Height and length of domain
-
-    Return:
-        A tuple containing two vectors: a list of distances from the fracture center to each fracture coordinate, and the corresponding analytical apertures.
-    """
-    ambient_dim = gb.dim_max()
-    g_1 = gb.subdomains(dim=ambient_dim - 1)[0]
-    fracture_center = np.array([length / 2, height / 2, 0])
-
-    fracture_faces = g_1.cell_centers
-
-    # compute distances from fracture centre with its corresponding apertures
-    eta = compute_eta(fracture_faces, fracture_center)
-    apertures = analytical_displacements(a, eta, p0, G, poi)
-
-    return eta, apertures
-
 
 class ModifiedGeometry:
     def set_domain(self) -> None:
@@ -325,12 +85,10 @@ def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
         n = 1000
         h = 2 * a / n
         center = np.array([length / 2, height / 2, 0])
-        bem_centers = get_bem_centers(a, h, n, theta, center)
-        eta = compute_eta(bem_centers, center)
-
-        u_a = -analytical_displacements(a, eta, p0, G, poi)
-
-        u_bc = assign_bem(sd, h / 2, box_faces, theta, bem_centers, u_a, poi)
+        bem_centers = manu_sneddon_2d.get_bem_centers(a, h, n, theta, center)
+        eta = manu_sneddon_2d.compute_eta(bem_centers, center)
+        u_a = -manu_sneddon_2d.analytical_displacements(a, eta, p0, G, poi)
+        u_bc = manu_sneddon_2d.assign_bem(sd, h / 2, box_faces, theta, bem_centers, u_a, poi)
 
         return bg.projection(2) @ u_bc.ravel("F")
 
@@ -369,6 +127,7 @@ class MomentumBalanceGeometryBC(
     """Adding geometry and modified boundary conditions to the default model."""
 
 
+
 def simulation(
     frac_pts: np.ndarray,
     theta_rad: np.ndarray,
@@ -429,7 +188,7 @@ def simulation(
     u_n = u_n.value(model.equation_system)
 
     # Checking convergence specifically on the fracture
-    u_a = run_analytical_displacements(model.mdg, a, p0, G, poi, height, length)[1]
+    u_a = manu_sneddon_2d.run_analytical_displacements(model.mdg, a, p0, G, poi, height, length)[1]
     e = ConvergenceAnalysis.l2_error(
         frac_sd[0], u_a, u_n, is_scalar=False, is_cc=True, relative=True
     )
@@ -579,6 +338,3 @@ def test_sneddon_2d():
     slope, _ = np.polyfit(log_h, log_avg_err, 1)
     assert slope > 1.0, f"Estimated EOC is {slope}, which is not greater than 1"
 
-
-
-
diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
new file mode 100644
index 0000000000..c56715518e
--- /dev/null
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -0,0 +1,241 @@
+import numpy as np
+import porepy as pp
+
+def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray:
+    """
+    Compute the distance fracture points to the fracture centre.
+
+    Parameter::
+    pointset: Array containing coordinates on the fracture
+    center: fracture centre
+
+    Return:
+        Array of distances each point in pointset to the center
+
+    """
+    return pp.geometry.distances.point_pointset(pointset_centers, center)
+
+
+def get_bem_centers(
+    a: float, h: float, n: int, theta: float, center: np.ndarray
+) -> np.ndarray:
+    """
+    Compute coordinates of the centers of the bem segments
+
+    Parameters:
+        a: half fracture length
+        h: bem segment length
+        n: number of bem segments
+        theta: orientation of the fracture
+        center: center of the fracture
+
+    Return:
+        Array of BEM centers
+    """
+
+    # Coordinate system (4.5.1) in page 57 in in book Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics
+
+    bem_centers = np.zeros((3, n))
+    x_0 = center[0] - (a - 0.5 * h) * np.sin(theta)
+    y_0 = center[1] - (a - 0.5 * h) * np.cos(theta)
+    for i in range(0, n):
+        bem_centers[0, i] = x_0 + i * h * np.sin(theta)
+        bem_centers[1, i] = y_0 + i * h * np.cos(theta)
+
+    return bem_centers
+
+
+def analytical_displacements(
+    a: float, eta: np.ndarray, p0: float, G: float, poi: float
+) -> np.ndarray:
+    """
+    Compute Sneddon's analytical solution for the pressurized fracture.
+
+    Source: Sneddon Fourier transforms 1951 page 425 eq 92
+
+    Parameter:
+        a: half fracture length
+        eta: distance from fracture centre
+        p0: pressure
+        G: shear modulus
+        poi: poisson ratio
+
+    Return
+        List of analytical normal displacement jumps.
+    """
+    cons = (1 - poi) / G * p0 * a * 2
+    return cons * np.sqrt(1 - np.power(eta / a, 2))
+
+
+def transform(xc: np.ndarray, x: np.ndarray, alpha: float) -> np.ndarray:
+    """
+    Translation and rotation coordinate transformation of boundary face coordinates for the BEM method
+
+    Parameter
+        xc: Coordinates of BEM segment centre
+        x: Coordinates of boundary faces
+        alpha: Fracture orientation
+    Return:
+        Transformed coordinates
+    """
+    x_bar = np.zeros_like(x)
+    # Terms in (7.4.6) in Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics page 168
+    x_bar[0, :] = (x[0, :] - xc[0]) * np.cos(alpha) + (x[1, :] - xc[1]) * np.sin(alpha)
+    x_bar[1, :] = -(x[0, :] - xc[0]) * np.sin(alpha) + (x[1, :] - xc[1]) * np.cos(alpha)
+    return x_bar
+
+
+def get_bc_val(
+    g: pp.GridLike,
+    bound_faces: np.ndarray,
+    xf: np.ndarray,
+    h: float,
+    poi: float,
+    alpha: float,
+    du: float,
+) -> np.ndarray:
+    """
+    Compute semi-analytical displacement on the boundary using the BEM method for the Sneddon problem.
+
+    Parameter
+        g: Grid object
+        bound_faces: Index lists for boundary faces
+        xf: Coordinates of boundary faces
+        h: BEM segment length
+        poi: Poisson ratio
+        alpha: Fracture orientation
+        du: Sneddon's analytical relative normal displacement
+    Return:
+        Boundary values for the displacement
+    """
+
+    # Equations for f2,f3,f4.f5 can be found in book Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics pages 57, 84-92, 168
+    f2 = np.zeros(bound_faces.size)
+    f3 = np.zeros(bound_faces.size)
+    f4 = np.zeros(bound_faces.size)
+    f5 = np.zeros(bound_faces.size)
+
+    u = np.zeros((g.dim, g.num_faces))
+
+    # Constant term in (7.4.5)
+    m = 1 / (4 * np.pi * (1 - poi))
+
+    # Second term in (7.4.5)
+    f2[:] = m * (
+        np.log(np.sqrt((xf[0, :] - h) ** 2 + xf[1] ** 2))
+        - np.log(np.sqrt((xf[0, :] + h) ** 2 + xf[1] ** 2))
+    )
+
+    # First term in (7.4.5)
+    f3[:] = -m * (
+        np.arctan2(xf[1, :], (xf[0, :] - h)) - np.arctan2(xf[1, :], (xf[0, :] + h))
+    )
+
+    # The following equalities can be found on page 91 where f3,f4 as equation (5.5.3) and ux, uy in (5.5.1)
+    # Also this is in the coordinate system described in (4.5.1) at page 57, which essentially is a rotation.
+    f4[:] = m * (
+        xf[1, :] / ((xf[0, :] - h) ** 2 + xf[1, :] ** 2)
+        - xf[1, :] / ((xf[0, :] + h) ** 2 + xf[1, :] ** 2)
+    )
+
+    f5[:] = m * (
+        (xf[0, :] - h) / ((xf[0, :] - h) ** 2 + xf[1, :] ** 2)
+        - (xf[0, :] + h) / ((xf[0, :] + h) ** 2 + xf[1, :] ** 2)
+    )
+
+    u[0, bound_faces] = du * (
+        -(1 - 2 * poi) * np.cos(alpha) * f2[:]
+        - 2 * (1 - poi) * np.sin(alpha) * f3[:]
+        - xf[1, :] * (np.cos(alpha) * f4[:] + np.sin(alpha) * f5[:])
+    )
+    u[1, bound_faces] = du * (
+        -(1 - 2 * poi) * np.sin(alpha) * f2[:]
+        + 2 * (1 - poi) * np.cos(alpha) * f3[:]
+        - xf[1, :] * (np.sin(alpha) * f4[:] - np.cos(alpha) * f5[:])
+    )
+
+    return u
+
+
+def assign_bem(
+    g: pp.Grid,
+    h: float,
+    bound_faces: np.ndarray,
+    theta: float,
+    bem_centers: np.ndarray,
+    u_a: np.ndarray,
+    poi: float,
+) -> np.ndarray:
+    """
+    Compute analytical displacement using the BEM method for the Sneddon problem.
+
+    Parameter
+        g: Subdomain grid
+        h: bem segment length
+        bound_faces: boundary faces
+        theta: fracture orientation
+        bem_centers: bem segments centers
+        u_a: Sneddon's analytical relative normal displacement
+        poi: Poisson ratio
+
+    Return:
+        Semi-analytical boundary displacement values
+    """
+
+    bc_val = np.zeros((g.dim, g.num_faces))
+
+    alpha = np.pi / 2 - theta
+
+    bound_face_centers = g.face_centers[:, bound_faces]
+
+    for i in range(0, u_a.size):
+
+        new_bound_face_centers = transform(bem_centers[:, i], bound_face_centers, alpha)
+
+        u_bound = get_bc_val(
+            g, bound_faces, new_bound_face_centers, h, poi, alpha, u_a[i]
+        )
+
+        bc_val += u_bound
+
+    return bc_val
+
+
+def run_analytical_displacements(
+    gb: pp.GridLike,
+    a: float,
+    p0: float,
+    G: float,
+    poi: float,
+    height: float,
+    length: float,
+) -> tuple:
+    """
+    Compute Sneddon's analytical solution for the pressurized crack
+    problem in question.
+
+    Source: Sneddon Fourier transforms 1951 page 425 eq 92
+
+    Parameter
+        gb: Grid object
+        a: Half fracture length of fracture
+        eta: Distance from fracture centre
+        p0: Internal constant pressure inside the fracture
+        G: Shear modulus
+        poi: Poisson ratio
+        height,length: Height and length of domain
+
+    Return:
+        A tuple containing two vectors: a list of distances from the fracture center to each fracture coordinate, and the corresponding analytical apertures.
+    """
+    ambient_dim = gb.dim_max()
+    g_1 = gb.subdomains(dim=ambient_dim - 1)[0]
+    fracture_center = np.array([length / 2, height / 2, 0])
+
+    fracture_faces = g_1.cell_centers
+
+    # compute distances from fracture centre with its corresponding apertures
+    eta = compute_eta(fracture_faces, fracture_center)
+    apertures = analytical_displacements(a, eta, p0, G, poi)
+
+    return eta, apertures
\ No newline at end of file

From d9fc240525eb4ce062c430f2b8d4758c2e043544 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 20 Jan 2025 10:06:22 +0100
Subject: [PATCH 07/43] Update setup for errors

---
 tests/functional/test_sneddon_2d.py | 54 +++++++++++++++++++++++++++++
 1 file changed, 54 insertions(+)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index dade774f73..72de247579 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -3,6 +3,57 @@
 import numpy as np
 import pytest
 
+import test_sneddon_2d
+
+
+
+
+# ----> Retrieve actual L2-errors
+@pytest.fixture(scope="module")
+def actual_l2_errors(material_constants: dict) -> :
+    """Run verification setups and retrieve results for the scheduled times.
+
+    Parameters:
+        material_constants: Dictionary containing the material constant classes.
+
+    Returns:
+        List of lists of dictionaries of actual relative errors. The outer list contains
+        two items, the first contains the results for 2d and the second contains the
+        results for 3d. Both inner lists contain three items each, each of which is a
+        dictionary of results for the scheduled times, i.e., 0.5 [s] and 1.0 [s].
+
+    """
+
+    # Define model parameters
+    model_params = {
+        "grid_type": "cartesian",
+        "material_constants": material_constants,
+        "meshing_arguments": {"cell_size": 0.25},
+        "manufactured_solution": "nordbotten_2016",
+        "time_manager": pp.TimeManager([0, 0.5, 1.0], 0.5, True),
+    }
+
+    # Retrieve actual L2-relative errors.
+    errors = []
+    # Loop through models, i.e., 2d and 3d.
+    model  = SneddonSetup2d()
+    
+    # Make deep copy of params to avoid nasty bugs.
+    setup = model(deepcopy(model_params))
+    pp.run_time_dependent_model(setup)
+    
+    errors_setup: list[dict[str, float]] = []
+    # Loop through results, i.e., results for each scheduled time.
+    for result in setup.results:
+        errors_setup.append(
+            {
+                "error_pressure": getattr(result, "error_pressure"),
+            }
+        )
+    errors.append(errors_setup)
+    return errors
+
+
 # ----> Retrieve actual order of convergence
 @pytest.fixture(scope="module")
 def actual_ooc(
@@ -11,6 +62,9 @@ def actual_ooc(
 
     """
     
+    
+       
+    
     ooc = []
     # Loop through the models
     ooc_setup: list[dict[str, float]] = []

From e2edeb381b3993700551de067d760b68bd2d94cc Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 20 Jan 2025 10:51:49 +0100
Subject: [PATCH 08/43] Added class for manufactured solution

---
 tests/functional/base_sneddon_2d.py        |  8 +++-
 tests/functional/setups/manu_sneddon_2d.py | 51 +++++++++++++++++++++-
 2 files changed, 56 insertions(+), 3 deletions(-)

diff --git a/tests/functional/base_sneddon_2d.py b/tests/functional/base_sneddon_2d.py
index a077e5c338..a04e52d805 100644
--- a/tests/functional/base_sneddon_2d.py
+++ b/tests/functional/base_sneddon_2d.py
@@ -172,6 +172,7 @@ def simulation(
         "theta": theta_rad,
         "p0": p0,
         "a": a,
+        "G": G,
         "poi": poi,
         "length": length,  # this info can be accessed via box, so remove these as well.
         "height": height,
@@ -186,9 +187,12 @@ def simulation(
     # Comuting the numerical displacement jump along the fracture on the fracture cell centers.
     u_n: pp.ad.Operator = nd_vec_to_normal @ model.displacement_jump(frac_sd)
     u_n = u_n.value(model.equation_system)
-
+    
+    exact_setup = manu_sneddon_2d.ManuExactSneddon2dSetup(params)
     # Checking convergence specifically on the fracture
-    u_a = manu_sneddon_2d.run_analytical_displacements(model.mdg, a, p0, G, poi, height, length)[1]
+    #u_a = manu_sneddon_2d.run_analytical_displacements(model.mdg, a, p0, G, poi, height, length)[1]
+    u_a = exact_setup.run_analytical_displacements(model.mdg)[1]
+  
     e = ConvergenceAnalysis.l2_error(
         frac_sd[0], u_a, u_n, is_scalar=False, is_cc=True, relative=True
     )
diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index c56715518e..f0cba22b74 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -238,4 +238,53 @@ def run_analytical_displacements(
     eta = compute_eta(fracture_faces, fracture_center)
     apertures = analytical_displacements(a, eta, p0, G, poi)
 
-    return eta, apertures
\ No newline at end of file
+    return eta, apertures
+
+
+
+
+class ManuExactSneddon2dSetup:
+  
+    def __init__(self, setup):
+        # Initialize private variables from the setup dictionary
+        self.p0 = setup.get("p0")
+        self.a = setup.get("a")
+        self.G  = setup.get("G")
+        self.poi = setup.get("poi")
+        self.length = setup.get("length")
+        self.height = setup.get("height")
+         
+  
+
+
+    def run_analytical_displacements( self,
+    gb: pp.GridLike ) -> tuple:
+        """
+        Compute Sneddon's analytical solution for the pressurized crack
+        problem in question.
+
+        Source: Sneddon Fourier transforms 1951 page 425 eq 92
+
+        Parameter
+            gb: Grid object
+            a: Half fracture length of fracture
+            eta: Distance from fracture centre
+            p0: Internal constant pressure inside the fracture
+            G: Shear modulus
+            poi: Poisson ratio
+            height,length: Height and length of domain
+
+        Return:
+            A tuple containing two vectors: a list of distances from the fracture center to each fracture coordinate, and the corresponding analytical apertures.
+        """
+        ambient_dim = gb.dim_max()
+        g_1 = gb.subdomains(dim=ambient_dim - 1)[0]
+        fracture_center = np.array([self.length / 2, self.height / 2, 0])
+
+        fracture_faces = g_1.cell_centers
+
+        # compute distances from fracture centre with its corresponding apertures
+        eta = compute_eta(fracture_faces, fracture_center)
+        apertures = analytical_displacements(self.a, eta, self.p0, self.G, self.poi)
+
+        return eta, apertures

From ce68288486983008d43ce898380d13cc3ba0e6fb Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 20 Jan 2025 11:00:01 +0100
Subject: [PATCH 09/43] renaming exact sol function

---
 tests/functional/base_sneddon_2d.py        |  3 +-
 tests/functional/setups/manu_sneddon_2d.py | 60 ++++++----------------
 2 files changed, 18 insertions(+), 45 deletions(-)

diff --git a/tests/functional/base_sneddon_2d.py b/tests/functional/base_sneddon_2d.py
index a04e52d805..3180ca2abc 100644
--- a/tests/functional/base_sneddon_2d.py
+++ b/tests/functional/base_sneddon_2d.py
@@ -190,8 +190,7 @@ def simulation(
     
     exact_setup = manu_sneddon_2d.ManuExactSneddon2dSetup(params)
     # Checking convergence specifically on the fracture
-    #u_a = manu_sneddon_2d.run_analytical_displacements(model.mdg, a, p0, G, poi, height, length)[1]
-    u_a = exact_setup.run_analytical_displacements(model.mdg)[1]
+    u_a = exact_setup.exact_sol_fracture(model.mdg)[1]
   
     e = ConvergenceAnalysis.l2_error(
         frac_sd[0], u_a, u_n, is_scalar=False, is_cc=True, relative=True
diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index f0cba22b74..77c606b8f1 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -201,46 +201,6 @@ def assign_bem(
     return bc_val
 
 
-def run_analytical_displacements(
-    gb: pp.GridLike,
-    a: float,
-    p0: float,
-    G: float,
-    poi: float,
-    height: float,
-    length: float,
-) -> tuple:
-    """
-    Compute Sneddon's analytical solution for the pressurized crack
-    problem in question.
-
-    Source: Sneddon Fourier transforms 1951 page 425 eq 92
-
-    Parameter
-        gb: Grid object
-        a: Half fracture length of fracture
-        eta: Distance from fracture centre
-        p0: Internal constant pressure inside the fracture
-        G: Shear modulus
-        poi: Poisson ratio
-        height,length: Height and length of domain
-
-    Return:
-        A tuple containing two vectors: a list of distances from the fracture center to each fracture coordinate, and the corresponding analytical apertures.
-    """
-    ambient_dim = gb.dim_max()
-    g_1 = gb.subdomains(dim=ambient_dim - 1)[0]
-    fracture_center = np.array([length / 2, height / 2, 0])
-
-    fracture_faces = g_1.cell_centers
-
-    # compute distances from fracture centre with its corresponding apertures
-    eta = compute_eta(fracture_faces, fracture_center)
-    apertures = analytical_displacements(a, eta, p0, G, poi)
-
-    return eta, apertures
-
-
 
 
 class ManuExactSneddon2dSetup:
@@ -248,6 +208,8 @@ class ManuExactSneddon2dSetup:
     def __init__(self, setup):
         # Initialize private variables from the setup dictionary
         self.p0 = setup.get("p0")
+        self.theta = setup.get("theta")
+
         self.a = setup.get("a")
         self.G  = setup.get("G")
         self.poi = setup.get("poi")
@@ -255,9 +217,21 @@ def __init__(self, setup):
         self.height = setup.get("height")
          
   
-
-
-    def run_analytical_displacements( self,
+    def exact_sol_global(self, sd):
+        
+        n = 1000
+        h = 2 * self.a / n
+        box_faces = sd.get_boundary_faces()
+
+        center = np.array([self.length / 2, self.height / 2, 0])
+        bem_centers = get_bem_centers(self.a, h, n, self.theta, center)
+        eta = compute_eta(bem_centers, center)
+        u_a = -analytical_displacements(self.a, eta, self.p0, self.G, self.poi)
+        u_bc = assign_bem(sd, h / 2, box_faces, self.theta, bem_centers, u_a, poi)
+        return u_bc
+    
+
+    def exact_sol_fracture( self,
     gb: pp.GridLike ) -> tuple:
         """
         Compute Sneddon's analytical solution for the pressurized crack

From bdffff1de07ef8200f4fd6732a2c06d5059ef734 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 20 Jan 2025 11:36:21 +0100
Subject: [PATCH 10/43] Moved functionality related to sneddon to setup

---
 tests/functional/base_sneddon_2d.py        | 30 +++-------------------
 tests/functional/setups/manu_sneddon_2d.py |  3 ++-
 2 files changed, 6 insertions(+), 27 deletions(-)

diff --git a/tests/functional/base_sneddon_2d.py b/tests/functional/base_sneddon_2d.py
index 3180ca2abc..74a10b73eb 100644
--- a/tests/functional/base_sneddon_2d.py
+++ b/tests/functional/base_sneddon_2d.py
@@ -59,37 +59,15 @@ def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
 
         """
 
-        a, p0, G, poi = (
-            self.params["a"],
-            self.params["p0"],
-            self.params["material_constants"]["solid"].shear_modulus,
-            self.params["poi"],
-        )
-        theta, length, height = (
-            self.params["theta"],
-            self.params["length"],
-            self.params["height"],
-        )
-
         sd = bg.parent
-        # bg = self.mdg.subdomain_to_boundary_grid(sd) # technique to extract bg via sd
-
         if sd.dim < 2:
             return np.zeros(self.nd * sd.num_faces)
-        box_faces = sd.get_boundary_faces()
-
-        # Set the boundary values
-        u_bc = np.zeros((sd.dim, sd.num_faces))
 
+  
         # apply sneddon analytical solution through BEM method
-        n = 1000
-        h = 2 * a / n
-        center = np.array([length / 2, height / 2, 0])
-        bem_centers = manu_sneddon_2d.get_bem_centers(a, h, n, theta, center)
-        eta = manu_sneddon_2d.compute_eta(bem_centers, center)
-        u_a = -manu_sneddon_2d.analytical_displacements(a, eta, p0, G, poi)
-        u_bc = manu_sneddon_2d.assign_bem(sd, h / 2, box_faces, theta, bem_centers, u_a, poi)
-
+        exact_sol = manu_sneddon_2d.ManuExactSneddon2dSetup(self.params)
+        u_bc = exact_sol.exact_sol_global(sd)
+        
         return bg.projection(2) @ u_bc.ravel("F")
 
 
diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 77c606b8f1..7fde2493f7 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -222,12 +222,13 @@ def exact_sol_global(self, sd):
         n = 1000
         h = 2 * self.a / n
         box_faces = sd.get_boundary_faces()
+        u_bc = np.zeros((sd.dim, sd.num_faces))
 
         center = np.array([self.length / 2, self.height / 2, 0])
         bem_centers = get_bem_centers(self.a, h, n, self.theta, center)
         eta = compute_eta(bem_centers, center)
         u_a = -analytical_displacements(self.a, eta, self.p0, self.G, self.poi)
-        u_bc = assign_bem(sd, h / 2, box_faces, self.theta, bem_centers, u_a, poi)
+        u_bc = assign_bem(sd, h / 2, box_faces, self.theta, bem_centers, u_a, self.poi)
         return u_bc
     
 

From ef328f67c2da29a956fdf1a89fe60581a4f72a35 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 20 Jan 2025 12:09:44 +0100
Subject: [PATCH 11/43] Update seperation

---
 tests/functional/base_sneddon_2d.py        | 95 +--------------------
 tests/functional/setups/manu_sneddon_2d.py | 97 ++++++++++++++++++++++
 2 files changed, 98 insertions(+), 94 deletions(-)

diff --git a/tests/functional/base_sneddon_2d.py b/tests/functional/base_sneddon_2d.py
index 74a10b73eb..8b3c33b3e0 100644
--- a/tests/functional/base_sneddon_2d.py
+++ b/tests/functional/base_sneddon_2d.py
@@ -12,99 +12,6 @@
 
 
 
-class ModifiedGeometry:
-    def set_domain(self) -> None:
-        """Defining a two-dimensional square domain with sidelength 2."""
-
-        self._domain = pp.Domain(
-            bounding_box={
-                "xmin": 0,
-                "ymin": 0,
-                "xmax": self.params["length"],
-                "ymax": self.params["height"],
-            }
-        )
-
-    def grid_type(self) -> str:
-        """Choosing the grid type for our domain."""
-        return self.params.get("grid_type", "simplex")
-
-    def set_fractures(self):
-
-        points = self.params["frac_pts"]
-        self._fractures = [pp.LineFracture(points)]
-
-
-class ModifiedBoundaryConditions:
-    def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
-        """Set boundary condition type for the problem."""
-        bounds = self.domain_boundary_sides(sd)
-
-        # Set the type of west and east boundaries to Dirichlet. North and south are
-        # Neumann by default.
-        bc = pp.BoundaryConditionVectorial(sd, bounds.all_bf, "dir")
-
-        frac_face = sd.tags["fracture_faces"]
-
-        bc.is_neu[:, frac_face] = False
-        bc.is_dir[:, frac_face] = True
-
-        return bc
-
-    def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
-        """Setting displacement boundary condition values.
-
-        This method returns an array of boundary condition values with the value 5t for
-        western boundaries and ones for the eastern boundary.
-
-        """
-
-        sd = bg.parent
-        if sd.dim < 2:
-            return np.zeros(self.nd * sd.num_faces)
-
-  
-        # apply sneddon analytical solution through BEM method
-        exact_sol = manu_sneddon_2d.ManuExactSneddon2dSetup(self.params)
-        u_bc = exact_sol.exact_sol_global(sd)
-        
-        return bg.projection(2) @ u_bc.ravel("F")
-
-
-class PressureStressMixin:
-    def pressure(self, subdomains: pp.Grid):
-        # Return pressure in the fracture domain
-        mat = self.params["p0"] * np.ones(
-            sum(subdomain.num_cells for subdomain in subdomains)
-        )
-        return pp.ad.DenseArray(mat)
-
-    def stress_discretization(
-        self, subdomains: list[pp.Grid]
-    ) -> pp.ad.BiotAd | pp.ad.MpsaAd:
-        """Discretization of the stress tensor.
-
-        Parameters:
-            subdomains: List of subdomains where the stress is defined.
-
-        Returns:
-            Discretization operator for the stress tensor.
-
-        """
-        return pp.ad.MpsaAd(self.stress_keyword, subdomains)
-
-
-class MomentumBalanceGeometryBC(
-    PressureStressMixin,
-    ModifiedGeometry,
-    ModifiedBoundaryConditions,
-    pp.constitutive_laws.PressureStress,
-    pp.MomentumBalance,
-):
-    ...
-    """Adding geometry and modified boundary conditions to the default model."""
-
-
 
 def simulation(
     frac_pts: np.ndarray,
@@ -156,7 +63,7 @@ def simulation(
         "height": height,
     }
 
-    model = MomentumBalanceGeometryBC(params)
+    model = manu_sneddon_2d.MomentumBalanceGeometryBC(params)
     pp.run_time_dependent_model(model, params)
 
     frac_sd = model.mdg.subdomains(dim=model.nd - 1)
diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 7fde2493f7..b2d2b4b7af 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -263,3 +263,100 @@ def exact_sol_fracture( self,
         apertures = analytical_displacements(self.a, eta, self.p0, self.G, self.poi)
 
         return eta, apertures
+
+
+
+
+
+class ModifiedGeometry:
+    def set_domain(self) -> None:
+        """Defining a two-dimensional square domain with sidelength 2."""
+
+        self._domain = pp.Domain(
+            bounding_box={
+                "xmin": 0,
+                "ymin": 0,
+                "xmax": self.params["length"],
+                "ymax": self.params["height"],
+            }
+        )
+
+    def grid_type(self) -> str:
+        """Choosing the grid type for our domain."""
+        return self.params.get("grid_type", "simplex")
+
+    def set_fractures(self):
+
+        points = self.params["frac_pts"]
+        self._fractures = [pp.LineFracture(points)]
+
+
+class ModifiedBoundaryConditions:
+    def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
+        """Set boundary condition type for the problem."""
+        bounds = self.domain_boundary_sides(sd)
+
+        # Set the type of west and east boundaries to Dirichlet. North and south are
+        # Neumann by default.
+        bc = pp.BoundaryConditionVectorial(sd, bounds.all_bf, "dir")
+
+        frac_face = sd.tags["fracture_faces"]
+
+        bc.is_neu[:, frac_face] = False
+        bc.is_dir[:, frac_face] = True
+
+        return bc
+
+    def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
+        """Setting displacement boundary condition values.
+
+        This method returns an array of boundary condition values with the value 5t for
+        western boundaries and ones for the eastern boundary.
+
+        """
+
+        sd = bg.parent
+        if sd.dim < 2:
+            return np.zeros(self.nd * sd.num_faces)
+
+  
+        # apply sneddon analytical solution through BEM method
+        exact_sol = ManuExactSneddon2dSetup(self.params)
+        u_bc = exact_sol.exact_sol_global(sd)
+        
+        return bg.projection(2) @ u_bc.ravel("F")
+
+
+class PressureStressMixin:
+    def pressure(self, subdomains: pp.Grid):
+        # Return pressure in the fracture domain
+        mat = self.params["p0"] * np.ones(
+            sum(subdomain.num_cells for subdomain in subdomains)
+        )
+        return pp.ad.DenseArray(mat)
+
+    def stress_discretization(
+        self, subdomains: list[pp.Grid]
+    ) -> pp.ad.BiotAd | pp.ad.MpsaAd:
+        """Discretization of the stress tensor.
+
+        Parameters:
+            subdomains: List of subdomains where the stress is defined.
+
+        Returns:
+            Discretization operator for the stress tensor.
+
+        """
+        return pp.ad.MpsaAd(self.stress_keyword, subdomains)
+
+
+class MomentumBalanceGeometryBC(
+    PressureStressMixin,
+    ModifiedGeometry,
+    ModifiedBoundaryConditions,
+    pp.constitutive_laws.PressureStress,
+    pp.MomentumBalance,
+):
+    ...
+    """Adding geometry and modified boundary conditions to the default model."""
+

From b408144ef703c1f9d1716ac974d87d77c88902fe Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 20 Jan 2025 13:28:31 +0100
Subject: [PATCH 12/43] Bit closer making the Conrgence class working

---
 tests/functional/test_sneddon_2d.py | 122 +++++++++++++++-------------
 1 file changed, 64 insertions(+), 58 deletions(-)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 72de247579..603ab15ba0 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -2,56 +2,41 @@
 import porepy as pp
 import numpy as np
 import pytest
+import math
+import copy
+import  tests.functional.setups.manu_sneddon_2d as manu_sneddon_2d
 
-import test_sneddon_2d
 
-
-
-
-# ----> Retrieve actual L2-errors
-@pytest.fixture(scope="module")
-def actual_l2_errors(material_constants: dict) -> :
-    """Run verification setups and retrieve results for the scheduled times.
+def compute_frac_pts(
+    theta_rad: float, a: float, height: float, length: float
+) -> np.ndarray:
+    """
+    Assuming the fracture center is at the coordinate (height/2, length/2),
+    compute the endpoints of a fracture given its orientation and fracture length.
 
     Parameters:
-        material_constants: Dictionary containing the material constant classes.
+    theta_rad: Angle of the fracture in radians
+    a: Half-length of the fracture.
+    height: Height of the domain.
+    length: Width of the domain.
 
     Returns:
-        List of lists of dictionaries of actual relative errors. The outer list contains
-        two items, the first contains the results for 2d and the second contains the
-        results for 3d. Both inner lists contain three items each, each of which is a
-        dictionary of results for the scheduled times, i.e., 0.5 [s] and 1.0 [s].
+
+        frac_pts : A 2x2 array where each column represents the coordinates of an end point of the fracture in 2D.
+            The first column corresponds to one end point, and the second column corresponds to the other.
 
     """
+    # Rotate the fracture with an angle theta_rad
+    y_0 = height / 2 - a * np.cos(theta_rad)
+    x_0 = length / 2 - a * np.sin(theta_rad)
+    y_1 = height / 2 + a * np.cos(theta_rad)
+    x_1 = length / 2 + a * np.sin(theta_rad)
+
+    frac_pts = np.array([[x_0, y_0], [x_1, y_1]]).T
+    return frac_pts
+
 
-    # Define model parameters
-    model_params = {
-        "grid_type": "cartesian",
-        "material_constants": material_constants,
-        "meshing_arguments": {"cell_size": 0.25},
-        "manufactured_solution": "nordbotten_2016",
-        "time_manager": pp.TimeManager([0, 0.5, 1.0], 0.5, True),
-    }
 
-    # Retrieve actual L2-relative errors.
-    errors = []
-    # Loop through models, i.e., 2d and 3d.
-    model  = SneddonSetup2d()
-    
-    # Make deep copy of params to avoid nasty bugs.
-    setup = model(deepcopy(model_params))
-    pp.run_time_dependent_model(setup)
-    
-    errors_setup: list[dict[str, float]] = []
-    # Loop through results, i.e., results for each scheduled time.
-    for result in setup.results:
-        errors_setup.append(
-            {
-                "error_pressure": getattr(result, "error_pressure"),
-            }
-        )
-    errors.append(errors_setup)
-    return errors
 
 
 # ----> Retrieve actual order of convergence
@@ -59,29 +44,50 @@ def actual_l2_errors(material_constants: dict) -> :
 def actual_ooc(
 ) -> float:
     """Retrieve actual order of convergence.
-
     """
+
+    # Orientation of fracture
+    num_theta = 6  # Number of iteration of theta
+    theta0 = 0  # Inital orientation of fracture
+    delta_theta = 10  # The difference in theta between each refinement
+    theta_list = (
+        theta0 + np.arange(0, num_theta) * delta_theta
+    )  # List of all the fracture orientations we run simulation on
+    
+    poi = 0.25
+    G = 1
+    lam = (
+        2 * G * poi / (1 - 2 * poi)
+    )  # Convertion formula from shear modulus and poission to lame lambda parameter
     
+    solid = pp.SolidConstants(shear_modulus=G, lame_lambda=lam)
+    params = {
+        #"meshing_arguments": {"cell_size": h},
+        "prepare_simulation": True,
+        "material_constants": {"solid": solid},
+        #"theta": theta_rad,
+        "a" : 1.5,
+        "height": 1.0,
+        "length": 1.0,
+        "p0": 1e-4,
+        "G": G,
+        "poi": poi
+    }
     
-       
+    for k in range(0, len(theta_list)):
+        params["theta"] = math.radians(90 - theta_list[k])
+        params["frac_pts"]  = compute_frac_pts(params["theta"], params["a"], height=params["height"], length=params["length"])
+        model =  manu_sneddon_2d.MomentumBalanceGeometryBC
+   
+        conv_analysis = ConvergenceAnalysis(
+            model_class=model,
+            model_params= copy.deepcopy(params),
+            levels=4,
+            spatial_refinement_rate=2,
+            temporal_refinement_rate=4
+        )
     
-    ooc = []
-    # Loop through the models
-    ooc_setup: list[dict[str, float]] = []
-    # Loop through grid type
-    # We do not perform a convergence analysis with simplices in 3d
-    grid_type = "simplex"
-        
-    conv_analysis = ConvergenceAnalysis(
-        model_class=model,
-        model_params=deepcopy(params),
-        levels=4,
-        spatial_refinement_rate=2,
-        temporal_refinement_rate=4,
-    )
     order = conv_analysis.order_of_convergence(conv_analysis.run_analysis())
-    
-
     return order 
 
 

From 4af1de0b2694e38da321cea19510bac9b6f68a67 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Wed, 22 Jan 2025 15:31:50 +0100
Subject: [PATCH 13/43] Update result list

---
 src/porepy/viz/data_saving_model_mixin.py  |  2 +-
 tests/functional/setups/manu_sneddon_2d.py | 35 ++++++++++++++++++++--
 2 files changed, 33 insertions(+), 4 deletions(-)

diff --git a/src/porepy/viz/data_saving_model_mixin.py b/src/porepy/viz/data_saving_model_mixin.py
index b04a067ebd..b1c5e2632c 100644
--- a/src/porepy/viz/data_saving_model_mixin.py
+++ b/src/porepy/viz/data_saving_model_mixin.py
@@ -242,7 +242,7 @@ def load_data_from_pvd(
 class VerificationDataSaving(DataSavingMixin):
     """Class to store relevant data for a generic verification setup."""
 
-    results: list
+    results: list = []
     """List of objects containing the results of the verification."""
 
     def save_data_time_step(self) -> None:
diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index b2d2b4b7af..d05c7363ae 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -1,5 +1,7 @@
 import numpy as np
 import porepy as pp
+from porepy.viz.data_saving_model_mixin import VerificationDataSaving
+from dataclasses import dataclass
 
 def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray:
     """
@@ -350,13 +352,40 @@ def stress_discretization(
         return pp.ad.MpsaAd(self.stress_keyword, subdomains)
 
 
+@dataclass
+class SneddonData:
+    error_displacement: pp.number
+
+class SneddonDataSaving(VerificationDataSaving):
+    def collect_data(self):
+        print("HELLO")
+        #frac_sd = model.mdg.subdomains(dim=model.nd - 1)
+        #nd_vec_to_normal = model.normal_component(frac_sd)
+    #
+        ## Comuting the numerical displacement jump along the fracture on the fracture cell centers.
+        #u_n: pp.ad.Operator = nd_vec_to_normal @ model.displacement_jump(frac_sd)
+        #u_n = u_n.value(model.equation_system)
+        #
+        #exact_setup = manu_sneddon_2d.ManuExactSneddon2dSetup(params)
+        ## Checking convergence specifically on the fracture
+        #u_a = exact_setup.exact_sol_fracture(model.mdg)[1]
+    #
+        #e = ConvergenceAnalysis.l2_error(
+        #    frac_sd[0], u_a, u_n, is_scalar=False, is_cc=True, relative=True
+        #)
+        e = 1.0 
+        collect_data = SneddonData(error_displacement=e)
+        return collect_data
+
+
 class MomentumBalanceGeometryBC(
     PressureStressMixin,
     ModifiedGeometry,
+    SneddonDataSaving,
     ModifiedBoundaryConditions,
     pp.constitutive_laws.PressureStress,
-    pp.MomentumBalance,
+    pp.MomentumBalance
 ):
-    ...
-    """Adding geometry and modified boundary conditions to the default model."""
+    pass
+    
 

From 22cedc7597a730bba633ee474d632363aa16e120 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Thu, 23 Jan 2025 12:21:09 +0100
Subject: [PATCH 14/43] Test is passing!

---
 tests/functional/setups/manu_sneddon_2d.py | 33 +++++++++++-----------
 tests/functional/test_sneddon_2d.py        | 14 ++++-----
 2 files changed, 23 insertions(+), 24 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index d05c7363ae..707a656357 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -1,6 +1,8 @@
 import numpy as np
 import porepy as pp
 from porepy.viz.data_saving_model_mixin import VerificationDataSaving
+from porepy.applications.convergence_analysis import ConvergenceAnalysis
+
 from dataclasses import dataclass
 
 def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray:
@@ -358,22 +360,21 @@ class SneddonData:
 
 class SneddonDataSaving(VerificationDataSaving):
     def collect_data(self):
-        print("HELLO")
-        #frac_sd = model.mdg.subdomains(dim=model.nd - 1)
-        #nd_vec_to_normal = model.normal_component(frac_sd)
-    #
-        ## Comuting the numerical displacement jump along the fracture on the fracture cell centers.
-        #u_n: pp.ad.Operator = nd_vec_to_normal @ model.displacement_jump(frac_sd)
-        #u_n = u_n.value(model.equation_system)
-        #
-        #exact_setup = manu_sneddon_2d.ManuExactSneddon2dSetup(params)
-        ## Checking convergence specifically on the fracture
-        #u_a = exact_setup.exact_sol_fracture(model.mdg)[1]
-    #
-        #e = ConvergenceAnalysis.l2_error(
-        #    frac_sd[0], u_a, u_n, is_scalar=False, is_cc=True, relative=True
-        #)
-        e = 1.0 
+        frac_sd = self.mdg.subdomains(dim=self.nd - 1)
+        nd_vec_to_normal = self.normal_component(frac_sd)
+    
+        # Comuting the numerical displacement jump along the fracture on the fracture cell centers.
+        u_n: pp.ad.Operator = nd_vec_to_normal @ self.displacement_jump(frac_sd)
+        u_n = u_n.value(self.equation_system)
+        
+        exact_setup = ManuExactSneddon2dSetup(self.params)
+        # Checking convergence specifically on the fracture
+        u_a = exact_setup.exact_sol_fracture(self.mdg)[1]
+    
+        e = ConvergenceAnalysis.l2_error(
+            frac_sd[0], u_a, u_n, is_scalar=False, is_cc=True, relative=True
+        )
+        
         collect_data = SneddonData(error_displacement=e)
         return collect_data
 
diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 603ab15ba0..64897b0ca9 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -1,10 +1,10 @@
-from porepy.applications.convergence_analysis import ConvergenceAnalysis
 import porepy as pp
 import numpy as np
 import pytest
 import math
 import copy
 import  tests.functional.setups.manu_sneddon_2d as manu_sneddon_2d
+from porepy.applications.convergence_analysis import ConvergenceAnalysis
 
 
 def compute_frac_pts(
@@ -36,9 +36,6 @@ def compute_frac_pts(
     return frac_pts
 
 
-
-
-
 # ----> Retrieve actual order of convergence
 @pytest.fixture(scope="module")
 def actual_ooc(
@@ -87,8 +84,8 @@ def actual_ooc(
             temporal_refinement_rate=4
         )
     
-    order = conv_analysis.order_of_convergence(conv_analysis.run_analysis())
-    return order 
+    order_dict = conv_analysis.order_of_convergence(conv_analysis.run_analysis())
+    return order_dict
 
 
 # ----> Set desired order of convergence
@@ -129,13 +126,14 @@ def test_order_of_convergence(
             order of convergence.
 
     """
+    
     # We require the order of convergence to always be larger than 1.0
-    assert 1.0 < actual_ooc
+    assert 1.0 <   actual_ooc["ooc_displacement"]  
 
 
     assert np.isclose(
         desired_ooc,
-        actual_ooc,
+        actual_ooc["ooc_displacement"],
         atol=1e-1,  # allow for an absolute difference of 0.1 in OOC
         rtol=5e-1,  # allow for 5% of relative difference in OOC
     )

From ca62a32fd30dac631b209de4b614e3bf064bdd42 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Thu, 23 Jan 2025 15:34:52 +0100
Subject: [PATCH 15/43] Test is working

---
 tests/functional/test_sneddon_2d.py | 83 ++++++++++-------------------
 1 file changed, 27 insertions(+), 56 deletions(-)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 64897b0ca9..eb247031fe 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -6,6 +6,15 @@
 import  tests.functional.setups.manu_sneddon_2d as manu_sneddon_2d
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 
+# ----> Set up the material constants
+poi = 0.25
+G = 1
+lam = (
+    2 * G * poi / (1 - 2 * poi)
+)  # Convertion formula from shear modulus and poission to lame lambda parameter
+
+solid = pp.SolidConstants(shear_modulus=G, lame_lambda=lam)
+   
 
 def compute_frac_pts(
     theta_rad: float, a: float, height: float, length: float
@@ -43,69 +52,40 @@ def actual_ooc(
     """Retrieve actual order of convergence.
     """
 
-    # Orientation of fracture
-    num_theta = 6  # Number of iteration of theta
-    theta0 = 0  # Inital orientation of fracture
-    delta_theta = 10  # The difference in theta between each refinement
-    theta_list = (
-        theta0 + np.arange(0, num_theta) * delta_theta
-    )  # List of all the fracture orientations we run simulation on
-    
-    poi = 0.25
-    G = 1
-    lam = (
-        2 * G * poi / (1 - 2 * poi)
-    )  # Convertion formula from shear modulus and poission to lame lambda parameter
-    
-    solid = pp.SolidConstants(shear_modulus=G, lame_lambda=lam)
+    theta = 30
+
     params = {
         #"meshing_arguments": {"cell_size": h},
         "prepare_simulation": True,
         "material_constants": {"solid": solid},
-        #"theta": theta_rad,
-        "a" : 1.5,
+        "a" : 0.3,
         "height": 1.0,
         "length": 1.0,
         "p0": 1e-4,
         "G": G,
-        "poi": poi
+        "poi": poi,
+        "meshing_arguments": {"cell_size": 0.03},
     }
     
-    for k in range(0, len(theta_list)):
-        params["theta"] = math.radians(90 - theta_list[k])
-        params["frac_pts"]  = compute_frac_pts(params["theta"], params["a"], height=params["height"], length=params["length"])
-        model =  manu_sneddon_2d.MomentumBalanceGeometryBC
-   
-        conv_analysis = ConvergenceAnalysis(
-            model_class=model,
-            model_params= copy.deepcopy(params),
-            levels=4,
-            spatial_refinement_rate=2,
-            temporal_refinement_rate=4
-        )
+    params["theta"] = math.radians(90 - theta[k])
+    params["frac_pts"]  = compute_frac_pts(params["theta"], params["a"], height=params["height"], length=params["length"])
+    model =  manu_sneddon_2d.MomentumBalanceGeometryBC
+
+    conv_analysis = ConvergenceAnalysis(
+        model_class=model,
+        model_params= copy.deepcopy(params),
+        levels=2,
+        spatial_refinement_rate=2,
+        temporal_refinement_rate=1
+    )
     
-    order_dict = conv_analysis.order_of_convergence(conv_analysis.run_analysis())
+    order_dict = conv_analysis.order_of_convergence(conv_analysis.run_analysis(), data_range=slice(None, None, None))
     return order_dict
 
 
-# ----> Set desired order of convergence
-@pytest.fixture(scope="module")
-def desired_ooc() -> float:
-    """Set desired order of convergence.
-
-    Returns:
-        List of lists of dictionaries, containing the desired order of convergence.
-
-    """
-    desired_ooc_2d = 2.0      
- 
-    return desired_ooc_2d
-
-
 
 def test_order_of_convergence(
     actual_ooc,
-    desired_ooc,
 ) -> None:
     """Test observed order of convergence.
 
@@ -126,14 +106,5 @@ def test_order_of_convergence(
             order of convergence.
 
     """
-    
     # We require the order of convergence to always be larger than 1.0
-    assert 1.0 <   actual_ooc["ooc_displacement"]  
-
-
-    assert np.isclose(
-        desired_ooc,
-        actual_ooc["ooc_displacement"],
-        atol=1e-1,  # allow for an absolute difference of 0.1 in OOC
-        rtol=5e-1,  # allow for 5% of relative difference in OOC
-    )
+    assert 0.85 <   actual_ooc["ooc_displacement"]  

From c652c8826b335a659ad98876f6f1c43904f476f4 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Thu, 23 Jan 2025 15:47:22 +0100
Subject: [PATCH 16/43] Cleanup comment

---
 tests/functional/test_sneddon_2d.py | 19 +------------------
 1 file changed, 1 insertion(+), 18 deletions(-)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index eb247031fe..a055137d1f 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -88,23 +88,6 @@ def test_order_of_convergence(
     actual_ooc,
 ) -> None:
     """Test observed order of convergence.
-
-    Note:
-        We set more flexible tolerances for simplicial grids compared to Cartesian
-        grids. This is because we would like to allow for slight changes in the
-        order of convergence if the meshes change, i.e., in newer versions of Gmsh.
-
-    Parameters:
-        var: Name of the variable to be tested.
-        grid_type_idx: Index to identify the grid type; `0` for Cartesian, and `1`
-            for simplices.
-        dim_idx: Index to identify the dimensionality of the problem; `0` for 2d, and
-            `1` for 3d.
-        actual_ooc: List of lists of dictionaries containing the actual observed
-            order of convergence.
-        desired_ooc: List of lists of dictionaries containing the desired observed
-            order of convergence.
-
     """
-    # We require the order of convergence to always be larger than 1.0
+    # We require the order of convergence to always be about 1.0 
     assert 0.85 <   actual_ooc["ooc_displacement"]  

From 625a2d95803f088422fd18f22e02736011034d79 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Thu, 23 Jan 2025 15:56:28 +0100
Subject: [PATCH 17/43] Rename G to shear modulus

---
 tests/functional/setups/manu_sneddon_2d.py | 6 +++---
 tests/functional/test_sneddon_2d.py        | 9 ++++-----
 2 files changed, 7 insertions(+), 8 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 707a656357..7b74822658 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -215,7 +215,7 @@ def __init__(self, setup):
         self.theta = setup.get("theta")
 
         self.a = setup.get("a")
-        self.G  = setup.get("G")
+        self.shear_modulus  = setup.get("material_constants").get("solid").shear_modulus
         self.poi = setup.get("poi")
         self.length = setup.get("length")
         self.height = setup.get("height")
@@ -231,7 +231,7 @@ def exact_sol_global(self, sd):
         center = np.array([self.length / 2, self.height / 2, 0])
         bem_centers = get_bem_centers(self.a, h, n, self.theta, center)
         eta = compute_eta(bem_centers, center)
-        u_a = -analytical_displacements(self.a, eta, self.p0, self.G, self.poi)
+        u_a = -analytical_displacements(self.a, eta, self.p0, self.shear_modulus, self.poi)
         u_bc = assign_bem(sd, h / 2, box_faces, self.theta, bem_centers, u_a, self.poi)
         return u_bc
     
@@ -264,7 +264,7 @@ def exact_sol_fracture( self,
 
         # compute distances from fracture centre with its corresponding apertures
         eta = compute_eta(fracture_faces, fracture_center)
-        apertures = analytical_displacements(self.a, eta, self.p0, self.G, self.poi)
+        apertures = analytical_displacements(self.a, eta, self.p0, self.shear_modulus, self.poi)
 
         return eta, apertures
 
diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index a055137d1f..de0893e9db 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -8,12 +8,12 @@
 
 # ----> Set up the material constants
 poi = 0.25
-G = 1
+shear_modulus = 1
 lam = (
-    2 * G * poi / (1 - 2 * poi)
+    2 * shear_modulus * poi / (1 - 2 * poi)
 )  # Convertion formula from shear modulus and poission to lame lambda parameter
 
-solid = pp.SolidConstants(shear_modulus=G, lame_lambda=lam)
+solid = pp.SolidConstants(shear_modulus=shear_modulus, lame_lambda=lam)
    
 
 def compute_frac_pts(
@@ -62,12 +62,11 @@ def actual_ooc(
         "height": 1.0,
         "length": 1.0,
         "p0": 1e-4,
-        "G": G,
         "poi": poi,
         "meshing_arguments": {"cell_size": 0.03},
     }
     
-    params["theta"] = math.radians(90 - theta[k])
+    params["theta"] = math.radians(90 - theta)
     params["frac_pts"]  = compute_frac_pts(params["theta"], params["a"], height=params["height"], length=params["length"])
     model =  manu_sneddon_2d.MomentumBalanceGeometryBC
 

From 10f322b67f13439f3d431706a9b9f374f759fe54 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Fri, 24 Jan 2025 16:54:44 +0100
Subject: [PATCH 18/43] Update comments

---
 tests/functional/setups/manu_sneddon_2d.py | 64 +++++++++++++++++-----
 tests/functional/test_sneddon_2d.py        | 18 ++++--
 2 files changed, 63 insertions(+), 19 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 7b74822658..bb8d33b92e 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -208,6 +208,11 @@ def assign_bem(
 
 
 class ManuExactSneddon2dSetup:
+    """
+    Class for setting up the analytical solution for the pressurized fracture problem.
+    
+    
+    """
   
     def __init__(self, setup):
         # Initialize private variables from the setup dictionary
@@ -222,6 +227,13 @@ def __init__(self, setup):
          
   
     def exact_sol_global(self, sd):
+        """
+        Compute the analytical solution for the pressurized fracture problem in question.
+        
+        Parameters:
+            sd: Subdomain for which the analytical solution is to be computed.
+            
+        """
         
         n = 1000
         h = 2 * self.a / n
@@ -231,6 +243,8 @@ def exact_sol_global(self, sd):
         center = np.array([self.length / 2, self.height / 2, 0])
         bem_centers = get_bem_centers(self.a, h, n, self.theta, center)
         eta = compute_eta(bem_centers, center)
+        
+        
         u_a = -analytical_displacements(self.a, eta, self.p0, self.shear_modulus, self.poi)
         u_bc = assign_bem(sd, h / 2, box_faces, self.theta, bem_centers, u_a, self.poi)
         return u_bc
@@ -274,7 +288,7 @@ def exact_sol_fracture( self,
 
 class ModifiedGeometry:
     def set_domain(self) -> None:
-        """Defining a two-dimensional square domain with sidelength 2."""
+        """Defining a two-dimensional square domain with sidelengths length and height."""
 
         self._domain = pp.Domain(
             bounding_box={
@@ -290,6 +304,7 @@ def grid_type(self) -> str:
         return self.params.get("grid_type", "simplex")
 
     def set_fractures(self):
+        """Setting the fractures in the domain from params."""
 
         points = self.params["frac_pts"]
         self._fractures = [pp.LineFracture(points)]
@@ -297,26 +312,34 @@ def set_fractures(self):
 
 class ModifiedBoundaryConditions:
     def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
-        """Set boundary condition type for the problem."""
+        """Set boundary condition type for the problem.
+        Parameters:
+            sd: Subdomain for which the boundary conditions are to be set.
+        Returns:
+            Boundary condition type for the problem.
+        """
         bounds = self.domain_boundary_sides(sd)
 
         # Set the type of west and east boundaries to Dirichlet. North and south are
         # Neumann by default.
         bc = pp.BoundaryConditionVectorial(sd, bounds.all_bf, "dir")
-
         frac_face = sd.tags["fracture_faces"]
-
         bc.is_neu[:, frac_face] = False
         bc.is_dir[:, frac_face] = True
 
         return bc
 
     def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
-        """Setting displacement boundary condition values.
-
-        This method returns an array of boundary condition values with the value 5t for
-        western boundaries and ones for the eastern boundary.
-
+        """
+        Setting displacement boundary condition values.
+        This method sets the displacement boundary condition values for a given 
+        boundary grid using the Sneddon analytical solution through the Boundary 
+        Element Method (BEM).
+        Parameters:
+            The boundary grid for which the displacement boundary condition values 
+            are to be set.
+        Returns:
+            An array of displacement boundary condition values.
         """
 
         sd = bg.parent
@@ -324,15 +347,25 @@ def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
             return np.zeros(self.nd * sd.num_faces)
 
   
-        # apply sneddon analytical solution through BEM method
+        # Get the analytical solution for the displacement 
         exact_sol = ManuExactSneddon2dSetup(self.params)
-        u_bc = exact_sol.exact_sol_global(sd)
+        u_exact = exact_sol.exact_sol_global(sd) 
         
-        return bg.projection(2) @ u_bc.ravel("F")
+        # Project the values to the grid
+        return bg.projection(2) @ u_exact.ravel("F")
 
 
 class PressureStressMixin:
     def pressure(self, subdomains: pp.Grid):
+        """Discretization of the stress tensor.
+
+        Parameters:
+            subdomains: List of subdomains where the stress is defined.
+
+        Returns:
+            Discretization operator for the stress tensor.
+
+        """
         # Return pressure in the fracture domain
         mat = self.params["p0"] * np.ones(
             sum(subdomain.num_cells for subdomain in subdomains)
@@ -349,17 +382,22 @@ def stress_discretization(
 
         Returns:
             Discretization operator for the stress tensor.
-
         """
         return pp.ad.MpsaAd(self.stress_keyword, subdomains)
 
 
 @dataclass
 class SneddonData:
+    """Data class for storing the error in the displacement field."""
     error_displacement: pp.number
 
 class SneddonDataSaving(VerificationDataSaving):
+    """Class for saving the error in the displacement field."""
     def collect_data(self):
+        """Collecting the error in the displacement field.
+        
+        Returns: collected data dictionary
+        """
         frac_sd = self.mdg.subdomains(dim=self.nd - 1)
         nd_vec_to_normal = self.normal_component(frac_sd)
     
diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index de0893e9db..5632561633 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -24,13 +24,12 @@ def compute_frac_pts(
     compute the endpoints of a fracture given its orientation and fracture length.
 
     Parameters:
-    theta_rad: Angle of the fracture in radians
-    a: Half-length of the fracture.
-    height: Height of the domain.
-    length: Width of the domain.
+        theta_rad: Angle of the fracture in radians
+        a: Half-length of the fracture.
+        height: Height of the domain.
+        length: Width of the domain.
 
     Returns:
-
         frac_pts : A 2x2 array where each column represents the coordinates of an end point of the fracture in 2D.
             The first column corresponds to one end point, and the second column corresponds to the other.
 
@@ -50,12 +49,14 @@ def compute_frac_pts(
 def actual_ooc(
 ) -> float:
     """Retrieve actual order of convergence.
+    
+    Returns: A dictionary containing the order of convergence for the displacement.
     """
 
+    # Angle of the fracture
     theta = 30
 
     params = {
-        #"meshing_arguments": {"cell_size": h},
         "prepare_simulation": True,
         "material_constants": {"solid": solid},
         "a" : 0.3,
@@ -66,10 +67,14 @@ def actual_ooc(
         "meshing_arguments": {"cell_size": 0.03},
     }
     
+    # Construct the fracture points for the given angle and length
     params["theta"] = math.radians(90 - theta)
     params["frac_pts"]  = compute_frac_pts(params["theta"], params["a"], height=params["height"], length=params["length"])
+    
     model =  manu_sneddon_2d.MomentumBalanceGeometryBC
 
+    # Construct the convergence analysis object, which does embedd the model, parameters refinementlevels
+    # for the convergence analysis
     conv_analysis = ConvergenceAnalysis(
         model_class=model,
         model_params= copy.deepcopy(params),
@@ -78,6 +83,7 @@ def actual_ooc(
         temporal_refinement_rate=1
     )
     
+    # Dictonary containing the order of convergence for the displacement
     order_dict = conv_analysis.order_of_convergence(conv_analysis.run_analysis(), data_range=slice(None, None, None))
     return order_dict
 

From a6957d4fc87a4163d4af279727894e5cf447be7f Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Fri, 24 Jan 2025 17:10:49 +0100
Subject: [PATCH 19/43] Update whitespaces

---
 tests/functional/test_sneddon_2d.py | 40 ++++++++++++++---------------
 1 file changed, 20 insertions(+), 20 deletions(-)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 5632561633..3021f11caa 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -46,44 +46,44 @@ def compute_frac_pts(
 
 # ----> Retrieve actual order of convergence
 @pytest.fixture(scope="module")
-def actual_ooc(
-) -> float:
-    """Retrieve actual order of convergence.
+def actual_ooc() -> float:
+    """
+    Prepare parameters, fracture parameters and model for moment balance equation for the convergence analysis.
     
-    Returns: A dictionary containing the order of convergence for the displacement.
+    Returns: A dictionary containing the experimental order of convergence for the displacement.
     """
+    # Angle of the fracture in degrees
+    theta_deg = 30
 
-    # Angle of the fracture
-    theta = 30
-
+    # Simulation parameters
     params = {
         "prepare_simulation": True,
         "material_constants": {"solid": solid},
-        "a" : 0.3,
+        "a": 0.3,
         "height": 1.0,
         "length": 1.0,
         "p0": 1e-4,
         "poi": poi,
         "meshing_arguments": {"cell_size": 0.03},
     }
-    
-    # Construct the fracture points for the given angle and length
-    params["theta"] = math.radians(90 - theta)
-    params["frac_pts"]  = compute_frac_pts(params["theta"], params["a"], height=params["height"], length=params["length"])
-    
-    model =  manu_sneddon_2d.MomentumBalanceGeometryBC
 
-    # Construct the convergence analysis object, which does embedd the model, parameters refinementlevels
-    # for the convergence analysis
+    # Convert angle to radians and compute fracture points
+    params["theta"] = math.radians(90 - theta_deg)
+    params["frac_pts"] = compute_frac_pts(params["theta"], params["a"], height=params["height"], length=params["length"])
+
+    # Model for the convergence analysis
+    model = manu_sneddon_2d.MomentumBalanceGeometryBC
+
+    # Convergence analysis setup
     conv_analysis = ConvergenceAnalysis(
         model_class=model,
-        model_params= copy.deepcopy(params),
+        model_params=copy.deepcopy(params),
         levels=2,
         spatial_refinement_rate=2,
         temporal_refinement_rate=1
     )
-    
-    # Dictonary containing the order of convergence for the displacement
+
+    # Calculate and return the order of convergence for the displacement
     order_dict = conv_analysis.order_of_convergence(conv_analysis.run_analysis(), data_range=slice(None, None, None))
     return order_dict
 
@@ -95,4 +95,4 @@ def test_order_of_convergence(
     """Test observed order of convergence.
     """
     # We require the order of convergence to always be about 1.0 
-    assert 0.85 <   actual_ooc["ooc_displacement"]  
+    assert 0.85 <   actual_ooc["ooc_displacement"]

From 1ea70740c18deab906815cacb4d028ce07eb15a3 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Fri, 24 Jan 2025 17:11:35 +0100
Subject: [PATCH 20/43] Remove base sneddon, which was just a template

---
 tests/functional/base_sneddon_2d.py | 228 ----------------------------
 1 file changed, 228 deletions(-)
 delete mode 100644 tests/functional/base_sneddon_2d.py

diff --git a/tests/functional/base_sneddon_2d.py b/tests/functional/base_sneddon_2d.py
deleted file mode 100644
index 8b3c33b3e0..0000000000
--- a/tests/functional/base_sneddon_2d.py
+++ /dev/null
@@ -1,228 +0,0 @@
-
-from matplotlib import pyplot as plt
-import numpy as np
-import porepy as pp
-from porepy.applications.convergence_analysis import ConvergenceAnalysis
-
-import math
-from numba import config
-import pytest
-import  tests.functional.setups.manu_sneddon_2d as manu_sneddon_2d
-config.DISABLE_JIT = True
-
-
-
-
-def simulation(
-    frac_pts: np.ndarray,
-    theta_rad: np.ndarray,
-    h: np.ndarray,
-    a: float = 1.5,
-    height: float = 1.0,
-    length: float = 1.0,
-    p0: float = 1e-4,
-    G: float = 1.0,
-    poi: float = 0.25,
-) -> float:
-    """
-    Simulates a 2D fracture linear elasticity problem using a momentum balance model and
-    computes the relative L2 error between numerical and analytical displacement jumps.
-
-    Parameters:
-        frac_pts: Array of coordinates specifying the fracture points in the domain.
-        theta_rad: Array of angles (in radians) defining the orientation of the fractures.
-        h: Array of cell sizes used for meshing the domain.
-        a:  Half length of the fracture
-        height: Height of the domain.
-        length: Length of the domain.
-        p0: Initial constant pressure applied inside the fracture.
-        G: Shear modulus of the material.
-        poi: Poisson's ratio of the material.
-
-    Returns:
-        e: The relative L2 error between the analytical and numerical displacement jumps
-        on the fracture.
-    """
-    lam = (
-        2 * G * poi / (1 - 2 * poi)
-    )  # Convertion formula from shear modulus and poission to lame lambda parameter
-    solid = pp.SolidConstants(shear_modulus=G, lame_lambda=lam)
-
-    # Clean this up!! This is made so I can easily access the params in arbitary functions
-    params = {
-        "meshing_arguments": {"cell_size": h},
-        "prepare_simulation": True,
-        "material_constants": {"solid": solid},
-        "frac_pts": frac_pts,
-        "theta": theta_rad,
-        "p0": p0,
-        "a": a,
-        "G": G,
-        "poi": poi,
-        "length": length,  # this info can be accessed via box, so remove these as well.
-        "height": height,
-    }
-
-    model = manu_sneddon_2d.MomentumBalanceGeometryBC(params)
-    pp.run_time_dependent_model(model, params)
-
-    frac_sd = model.mdg.subdomains(dim=model.nd - 1)
-    nd_vec_to_normal = model.normal_component(frac_sd)
-
-    # Comuting the numerical displacement jump along the fracture on the fracture cell centers.
-    u_n: pp.ad.Operator = nd_vec_to_normal @ model.displacement_jump(frac_sd)
-    u_n = u_n.value(model.equation_system)
-    
-    exact_setup = manu_sneddon_2d.ManuExactSneddon2dSetup(params)
-    # Checking convergence specifically on the fracture
-    u_a = exact_setup.exact_sol_fracture(model.mdg)[1]
-  
-    e = ConvergenceAnalysis.l2_error(
-        frac_sd[0], u_a, u_n, is_scalar=False, is_cc=True, relative=True
-    )
-
-    return e
-
-
-def compute_frac_pts(
-    theta_rad: float, a: float, height: float, length: float
-) -> np.ndarray:
-    """
-    Assuming the fracture center is at the coordinate (height/2, length/2),
-    compute the endpoints of a fracture given its orientation and fracture length.
-
-    Parameters:
-    theta_rad: Angle of the fracture in radians
-    a: Half-length of the fracture.
-    height: Height of the domain.
-    length: Width of the domain.
-
-    Returns:
-
-        frac_pts : A 2x2 array where each column represents the coordinates of an end point of the fracture in 2D.
-            The first column corresponds to one end point, and the second column corresponds to the other.
-
-    """
-    # Rotate the fracture with an angle theta_rad
-    y_0 = height / 2 - a * np.cos(theta_rad)
-    x_0 = length / 2 - a * np.sin(theta_rad)
-    y_1 = height / 2 + a * np.cos(theta_rad)
-    x_1 = length / 2 + a * np.sin(theta_rad)
-
-    frac_pts = np.array([[x_0, y_0], [x_1, y_1]]).T
-    return frac_pts
-
-
-def compute_eoc(h: np.ndarray, e: np.ndarray) -> np.ndarray:
-    """
-    Compute the Error Order of Convergence (EOC) based on error and mesh size.
-
-    Parameters:
-        h: Array of mesh sizes corresponding to simulation results.
-        e: Array of error values for the corresponding mesh sizes.
-
-    Returns:
-        eoc: Array of EOC values computed for consecutive levels of refinement.
-    """
-
-    # Compute the Error Order of Convergence (EOC)
-    eoc = np.zeros(len(e) - 1)
-    for i in range(len(e) - 1):
-        eoc[i] = np.log(e[i + 1] / e[i]) / np.log(h[i + 1] / h[i])
-    return eoc
-
-
-def compute_convergence(
-    h_list: np.ndarray,
-    theta_list: np.ndarray,
-    a: float = 1,
-    height: float = 1,
-    length: float = 1,
-) -> np.ndarray:
-    """
-    Compute the convergence behavior for the Sneddon problem across mesh refinements and orientations.
-
-    Parameters:
-        h_list: Array of mesh sizes for simulations.
-        theta_list: Array of fracture orientations (in degrees).
-        a: Half-length of the fracture.
-        height: Height of the simulation domain.
-        length: Length of the simulation domain.
-
-    Returns:
-        A 2D array where rows correspond to different orientations and columns to mesh sizes.
-    """  
-    # Compute error for each orientation of the fracture 
-    err = np.zeros((len(theta_list), len(h_list)))
-    for k in range(0, len(theta_list)):
-        
-        theta_rad = math.radians(90 - theta_list[k])
-        frac_pts = compute_frac_pts(theta_rad, a, height=height, length=length)
-        for i in range(0, len(h_list)):
-            e = simulation(
-                frac_pts, theta_rad, h_list[i], a=a, height=height, length=length
-            )
-            err[k, i] = e
-
-    return err
-
-
-def test_sneddon_2d():
-    """
-    Test function to 2D Sneddon problem convergence for the MPSA method.
-    
-    This is a setup for comparing the analytical solution (also known as Sneddon's solution) with the numerical solution for linear elasticity
-    in the case of an open fracture subjected to a constant pressure p0. Fluid effects are not considered.
-    For reference about the implementation, see Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics Chapter 5.3 Pressurized crack problem.
-    Also like to note that this is a reimplementation of Sneddon solutation described in https://doi.org/10.1007/s10596-020-10002-5.
-
-    The tests performs a convergence study for a 2D fracture problem by iterating over fracture orientations and
-    mesh refinements. It computes the average error across orientations and evaluates the convergence rate (EOC)
-    through log-log linear regression.
-    
-    NOTE: We assume that the fracture half-length a = 1 and the domain is (0,1)^2 for the following reasons:
-    1. Convergence tends to be problematic if the domain is changed from (0,1) to another size configuration,
-       suggesting the presence of a bug.
-    2. While a total fracture length significantly larger than the domain shows convergence, convergence issues arise
-       when the fracture length is smaller, potentially due to numerical artifacts that occur when the fracture tip
-       approaches the boundary of the simulation domain.
-
-    Raises:
-    -------
-    AssertionError
-        If the estimated EOC is not greater than 1.0.
-    """
-
-    # Simulation settings
-    num_refs = 4  # Number of refinements
-    num_theta = 6  # Number of iteration of theta
-
-    theta0 = 0  # Inital orientation of fracture
-    delta_theta = 10  # The difference in theta between each refinement
-    theta_list = (
-        theta0 + np.arange(0, num_theta) * delta_theta
-    )  # List of all the fracture orientations we run simulation on
-
-    # Mesh refinement
-    nf0 = 6
-    nf = nf0 * 2 ** np.arange(0, num_refs)
-    h_list = 2 * 1.5 / nf  # Mesh size for simulations
-
-    # In this test we define the
-    height, length = 1, 1
-
-    # Here we define the half length of the fracture to be sufficientl big to avoid any fracture tip effects
-    a = 1
-
-    # Run convergence experiment
-    err = compute_convergence(h_list, theta_list, a=a, height=height, length=length)
-
-    # Computing average error of each theta
-    avg_err = np.mean(err, axis=0)
-
-    # Perform linear regression on the log-log scale and compute EOC
-    log_h = np.log(h_list)
-    log_avg_err = np.log(avg_err)
-    slope, _ = np.polyfit(log_h, log_avg_err, 1)
-    assert slope > 1.0, f"Estimated EOC is {slope}, which is not greater than 1"
-

From 15154f931765b6f8a812b3b19c7dc4c5102eb584 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 27 Jan 2025 10:54:04 +0100
Subject: [PATCH 21/43] Update comment

---
 tests/functional/test_sneddon_2d.py | 20 ++++++++++++--------
 1 file changed, 12 insertions(+), 8 deletions(-)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 3021f11caa..f007bdf33a 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -48,9 +48,13 @@ def compute_frac_pts(
 @pytest.fixture(scope="module")
 def actual_ooc() -> float:
     """
-    Prepare parameters, fracture parameters and model for moment balance equation for the convergence analysis.
+    Prepare parameters and model for the convergence analysis of the moment balance equation.
     
-    Returns: A dictionary containing the experimental order of convergence for the displacement.
+    This setup validates the linear elasticity model for the analytical Sneddon solution in 2D, describing the analytical displacement on the fracture. 
+    The problem consists of a 2D domain with a fracture at a given angle and internal pressure.
+    
+    Returns:
+        A dictionary containing the experimental order of convergence for the displacement.
     """
     # Angle of the fracture in degrees
     theta_deg = 30
@@ -59,11 +63,11 @@ def actual_ooc() -> float:
     params = {
         "prepare_simulation": True,
         "material_constants": {"solid": solid},
-        "a": 0.3,
-        "height": 1.0,
-        "length": 1.0,
-        "p0": 1e-4,
-        "poi": poi,
+        "a": 0.3, # Half-length of the fracture
+        "height": 1.0, # Height of the domain
+        "length": 1.0, # Length of the domain
+        "p0": 1e-4, # Internal pressure of fracture
+        "poi": poi, # Possion ratio (Not standard in solid constants)
         "meshing_arguments": {"cell_size": 0.03},
     }
 
@@ -94,5 +98,5 @@ def test_order_of_convergence(
 ) -> None:
     """Test observed order of convergence.
     """
-    # We require the order of convergence to always be about 1.0 
+    # We expect the order of L2 convergence on the fracture of dispplacement to be about 1.0 
     assert 0.85 <   actual_ooc["ooc_displacement"]

From cfdf68fc1205402a33d982f362462b06449a04bb Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 27 Jan 2025 12:14:55 +0100
Subject: [PATCH 22/43] Fixed most of mypy hints

---
 tests/functional/test_sneddon_2d.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index f007bdf33a..1d73ff5649 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -46,7 +46,7 @@ def compute_frac_pts(
 
 # ----> Retrieve actual order of convergence
 @pytest.fixture(scope="module")
-def actual_ooc() -> float:
+def actual_ooc() -> dict:
     """
     Prepare parameters and model for the convergence analysis of the moment balance equation.
     
@@ -57,7 +57,7 @@ def actual_ooc() -> float:
         A dictionary containing the experimental order of convergence for the displacement.
     """
     # Angle of the fracture in degrees
-    theta_deg = 30
+    theta_deg = 30.0
 
     # Simulation parameters
     params = {
@@ -72,8 +72,8 @@ def actual_ooc() -> float:
     }
 
     # Convert angle to radians and compute fracture points
-    params["theta"] = math.radians(90 - theta_deg)
-    params["frac_pts"] = compute_frac_pts(params["theta"], params["a"], height=params["height"], length=params["length"])
+    params["theta"] = math.radians(90.0 - theta_deg)
+    params["frac_pts"] = compute_frac_pts(theta_rad=float(params["theta"]), a=float(params["a"]), height=float(params["height"]), length=float(params["length"]))
 
     # Model for the convergence analysis
     model = manu_sneddon_2d.MomentumBalanceGeometryBC
@@ -98,5 +98,5 @@ def test_order_of_convergence(
 ) -> None:
     """Test observed order of convergence.
     """
-    # We expect the order of L2 convergence on the fracture of dispplacement to be about 1.0 
+    # We  the order of L2 convergence on the fracture of dispplacement to be about 1.0 
     assert 0.85 <   actual_ooc["ooc_displacement"]

From b14a0045eca7b89b08acc77f65bd0483f8825c26 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 27 Jan 2025 12:49:32 +0100
Subject: [PATCH 23/43] Fixed wrong name of class

---
 tests/functional/setups/manu_sneddon_2d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index bb8d33b92e..3e2cd4152c 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -251,7 +251,7 @@ def exact_sol_global(self, sd):
     
 
     def exact_sol_fracture( self,
-    gb: pp.GridLike ) -> tuple:
+    gb: pp.MixedDimensionalGrid) -> tuple:
         """
         Compute Sneddon's analytical solution for the pressurized crack
         problem in question.

From d21bc318c6285e654fc6c5ef0b7b19cc8c96ca50 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 27 Jan 2025 13:02:24 +0100
Subject: [PATCH 24/43] Fixing more mypy errors

---
 tests/functional/setups/manu_sneddon_2d.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 3e2cd4152c..f00673736b 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -2,6 +2,7 @@
 import porepy as pp
 from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
+from porepy.models.protocol import PorePyModel
 
 from dataclasses import dataclass
 
@@ -90,7 +91,7 @@ def transform(xc: np.ndarray, x: np.ndarray, alpha: float) -> np.ndarray:
 
 
 def get_bc_val(
-    g: pp.GridLike,
+    g: pp.Grid,
     bound_faces: np.ndarray,
     xf: np.ndarray,
     h: float,
@@ -112,7 +113,7 @@ def get_bc_val(
     Return:
         Boundary values for the displacement
     """
-
+    
     # Equations for f2,f3,f4.f5 can be found in book Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics pages 57, 84-92, 168
     f2 = np.zeros(bound_faces.size)
     f3 = np.zeros(bound_faces.size)
@@ -286,7 +287,7 @@ def exact_sol_fracture( self,
 
 
 
-class ModifiedGeometry:
+class ModifiedGeometry(PorePyModel):
     def set_domain(self) -> None:
         """Defining a two-dimensional square domain with sidelengths length and height."""
 
@@ -310,7 +311,7 @@ def set_fractures(self):
         self._fractures = [pp.LineFracture(points)]
 
 
-class ModifiedBoundaryConditions:
+class ModifiedBoundaryConditions(PorePyModel):
     def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
         """Set boundary condition type for the problem.
         Parameters:
@@ -355,7 +356,7 @@ def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
         return bg.projection(2) @ u_exact.ravel("F")
 
 
-class PressureStressMixin:
+class PressureStressMixin(PorePyModel):
     def pressure(self, subdomains: pp.Grid):
         """Discretization of the stress tensor.
 

From 5605ddafb996c87ec8271fa532e374dcabe773a5 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Mon, 27 Jan 2025 13:08:19 +0100
Subject: [PATCH 25/43] Improve compatibility

---
 tests/functional/setups/manu_sneddon_2d.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index f00673736b..2d6896f664 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -3,6 +3,8 @@
 from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 from porepy.models.protocol import PorePyModel
+from typing import Literal
+
 
 from dataclasses import dataclass
 
@@ -300,7 +302,7 @@ def set_domain(self) -> None:
             }
         )
 
-    def grid_type(self) -> str:
+    def grid_type(self) -> Literal["simplex", "cartesian", "tensor_grid"]:
         """Choosing the grid type for our domain."""
         return self.params.get("grid_type", "simplex")
 

From db51243d3cf0f5ca5cf8c0c73b51e9a10763ae8a Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Thu, 30 Jan 2025 14:38:07 +0100
Subject: [PATCH 26/43] Update comments

---
 tests/functional/test_sneddon_2d.py | 28 ++++++++++++----------------
 1 file changed, 12 insertions(+), 16 deletions(-)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 1d73ff5649..3021f11caa 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -46,34 +46,30 @@ def compute_frac_pts(
 
 # ----> Retrieve actual order of convergence
 @pytest.fixture(scope="module")
-def actual_ooc() -> dict:
+def actual_ooc() -> float:
     """
-    Prepare parameters and model for the convergence analysis of the moment balance equation.
+    Prepare parameters, fracture parameters and model for moment balance equation for the convergence analysis.
     
-    This setup validates the linear elasticity model for the analytical Sneddon solution in 2D, describing the analytical displacement on the fracture. 
-    The problem consists of a 2D domain with a fracture at a given angle and internal pressure.
-    
-    Returns:
-        A dictionary containing the experimental order of convergence for the displacement.
+    Returns: A dictionary containing the experimental order of convergence for the displacement.
     """
     # Angle of the fracture in degrees
-    theta_deg = 30.0
+    theta_deg = 30
 
     # Simulation parameters
     params = {
         "prepare_simulation": True,
         "material_constants": {"solid": solid},
-        "a": 0.3, # Half-length of the fracture
-        "height": 1.0, # Height of the domain
-        "length": 1.0, # Length of the domain
-        "p0": 1e-4, # Internal pressure of fracture
-        "poi": poi, # Possion ratio (Not standard in solid constants)
+        "a": 0.3,
+        "height": 1.0,
+        "length": 1.0,
+        "p0": 1e-4,
+        "poi": poi,
         "meshing_arguments": {"cell_size": 0.03},
     }
 
     # Convert angle to radians and compute fracture points
-    params["theta"] = math.radians(90.0 - theta_deg)
-    params["frac_pts"] = compute_frac_pts(theta_rad=float(params["theta"]), a=float(params["a"]), height=float(params["height"]), length=float(params["length"]))
+    params["theta"] = math.radians(90 - theta_deg)
+    params["frac_pts"] = compute_frac_pts(params["theta"], params["a"], height=params["height"], length=params["length"])
 
     # Model for the convergence analysis
     model = manu_sneddon_2d.MomentumBalanceGeometryBC
@@ -98,5 +94,5 @@ def test_order_of_convergence(
 ) -> None:
     """Test observed order of convergence.
     """
-    # We  the order of L2 convergence on the fracture of dispplacement to be about 1.0 
+    # We require the order of convergence to always be about 1.0 
     assert 0.85 <   actual_ooc["ooc_displacement"]

From c641b62ddbb373e98f03cfb5582fa9fe502482de Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Thu, 30 Jan 2025 14:39:11 +0100
Subject: [PATCH 27/43] Update comments

---
 tests/functional/test_sneddon_2d.py | 28 ++++++++++++++++------------
 1 file changed, 16 insertions(+), 12 deletions(-)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 3021f11caa..e99c71c8b4 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -46,30 +46,34 @@ def compute_frac_pts(
 
 # ----> Retrieve actual order of convergence
 @pytest.fixture(scope="module")
-def actual_ooc() -> float:
+def actual_ooc() -> dict:
     """
-    Prepare parameters, fracture parameters and model for moment balance equation for the convergence analysis.
+    Prepare parameters and model for the convergence analysis of the moment balance equation.
     
-    Returns: A dictionary containing the experimental order of convergence for the displacement.
+    This setup validates the linear elasticity model for the analytical Sneddon solution in 2D, describing the analytical displacement on the fracture. 
+    The problem consists of a 2D domain with a fracture at a given angle and internal pressure.
+    
+    Returns:
+        A dictionary containing the experimental order of convergence for the displacement.
     """
     # Angle of the fracture in degrees
-    theta_deg = 30
+    theta_deg = 30.0
 
     # Simulation parameters
     params = {
         "prepare_simulation": True,
         "material_constants": {"solid": solid},
-        "a": 0.3,
-        "height": 1.0,
-        "length": 1.0,
-        "p0": 1e-4,
-        "poi": poi,
+        "a": 0.3, # Half-length of the fracture
+        "height": 1.0, # Height of the domain
+        "length": 1.0, # Length of the domain
+        "p0": 1e-4, # Internal pressure of fracture
+        "poi": poi, # Possion ratio (Not standard in solid constants)
         "meshing_arguments": {"cell_size": 0.03},
     }
 
     # Convert angle to radians and compute fracture points
-    params["theta"] = math.radians(90 - theta_deg)
-    params["frac_pts"] = compute_frac_pts(params["theta"], params["a"], height=params["height"], length=params["length"])
+    params["theta"] = math.radians(90.0 - theta_deg)
+    params["frac_pts"] = compute_frac_pts(theta_rad=float(params["theta"]), a=float(params["a"]), height=float(params["height"]), length=float(params["length"]))
 
     # Model for the convergence analysis
     model = manu_sneddon_2d.MomentumBalanceGeometryBC
@@ -94,5 +98,5 @@ def test_order_of_convergence(
 ) -> None:
     """Test observed order of convergence.
     """
-    # We require the order of convergence to always be about 1.0 
+    # We  the order of L2 convergence on the fracture of displacement to be about 1.0 
     assert 0.85 <   actual_ooc["ooc_displacement"]

From 3b4e8f5991f37133352c9574ef1e6eed57312113 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Thu, 30 Jan 2025 15:01:04 +0100
Subject: [PATCH 28/43] test_sneddon_2d.py has successfully passed mypy

---
 tests/functional/test_sneddon_2d.py | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index e99c71c8b4..c29626880d 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -57,23 +57,28 @@ def actual_ooc() -> dict:
         A dictionary containing the experimental order of convergence for the displacement.
     """
     # Angle of the fracture in degrees
-    theta_deg = 30.0
 
     # Simulation parameters
+    theta_deg = 30.0
+    a = 0.3
+    height = 1.0
+    length = 1.0
+    theta_rad = math.radians(90 - theta_deg)
+
     params = {
         "prepare_simulation": True,
         "material_constants": {"solid": solid},
-        "a": 0.3, # Half-length of the fracture
-        "height": 1.0, # Height of the domain
-        "length": 1.0, # Length of the domain
+        "a": a, # Half-length of the fracture
+        "height": height, # Height of the domain
+        "length": length, # Length of the domain
         "p0": 1e-4, # Internal pressure of fracture
         "poi": poi, # Possion ratio (Not standard in solid constants)
         "meshing_arguments": {"cell_size": 0.03},
+        "theta": theta_rad,
     }
 
     # Convert angle to radians and compute fracture points
-    params["theta"] = math.radians(90.0 - theta_deg)
-    params["frac_pts"] = compute_frac_pts(theta_rad=float(params["theta"]), a=float(params["a"]), height=float(params["height"]), length=float(params["length"]))
+    params["frac_pts"] = compute_frac_pts(theta_rad=theta_rad, a=a, height=height, length=length)
 
     # Model for the convergence analysis
     model = manu_sneddon_2d.MomentumBalanceGeometryBC

From 164d06da87e66aad9bde9dce494d745ef4e9e557 Mon Sep 17 00:00:00 2001
From: Isak Hammer <isakhammer@gmail.com>
Date: Thu, 30 Jan 2025 15:16:34 +0100
Subject: [PATCH 29/43] Fixed obvious mypy errors

---
 tests/functional/setups/manu_sneddon_2d.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 2d6896f664..89c02689eb 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -359,7 +359,9 @@ def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
 
 
 class PressureStressMixin(PorePyModel):
-    def pressure(self, subdomains: pp.Grid):
+    stress_keyword: str
+
+    def pressure(self, subdomains: pp.SubdomainsOrBoundaries):
         """Discretization of the stress tensor.
 
         Parameters:
@@ -371,7 +373,7 @@ def pressure(self, subdomains: pp.Grid):
         """
         # Return pressure in the fracture domain
         mat = self.params["p0"] * np.ones(
-            sum(subdomain.num_cells for subdomain in subdomains)
+            sum((subdomain.num_cells for subdomain in subdomains))
         )
         return pp.ad.DenseArray(mat)
 
@@ -421,6 +423,7 @@ def collect_data(self):
 
 
 class MomentumBalanceGeometryBC(
+    PorePyModel,
     PressureStressMixin,
     ModifiedGeometry,
     SneddonDataSaving,

From 3ce3820cf93b4427d38b61c1ca1d1abfb5e6354d Mon Sep 17 00:00:00 2001
From: Veljko Lipovac <lipovac.v@outlook.com>
Date: Thu, 13 Feb 2025 12:57:18 +0100
Subject: [PATCH 30/43] TST: Removing obsolete Verification data saving.

---
 tests/functional/setups/manu_sneddon_2d.py | 13 +++++--------
 1 file changed, 5 insertions(+), 8 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 89c02689eb..62e92cedc8 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -1,6 +1,5 @@
 import numpy as np
 import porepy as pp
-from porepy.viz.data_saving_model_mixin import VerificationDataSaving
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 from porepy.models.protocol import PorePyModel
 from typing import Literal
@@ -8,6 +7,7 @@
 
 from dataclasses import dataclass
 
+
 def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray:
     """
     Compute the distance fracture points to the fracture centre.
@@ -208,8 +208,6 @@ def assign_bem(
     return bc_val
 
 
-
-
 class ManuExactSneddon2dSetup:
     """
     Class for setting up the analytical solution for the pressurized fracture problem.
@@ -286,9 +284,6 @@ def exact_sol_fracture( self,
         return eta, apertures
 
 
-
-
-
 class ModifiedGeometry(PorePyModel):
     def set_domain(self) -> None:
         """Defining a two-dimensional square domain with sidelengths length and height."""
@@ -396,9 +391,11 @@ class SneddonData:
     """Data class for storing the error in the displacement field."""
     error_displacement: pp.number
 
-class SneddonDataSaving(VerificationDataSaving):
+
+class SneddonDataSaving(pp.PorePyModel):
     """Class for saving the error in the displacement field."""
-    def collect_data(self):
+
+    def collect_data(self) -> SneddonData:
         """Collecting the error in the displacement field.
         
         Returns: collected data dictionary

From 57716a2261c57ebedc427311b05a440ff96dd573 Mon Sep 17 00:00:00 2001
From: Veljko Lipovac <lipovac.v@outlook.com>
Date: Thu, 13 Feb 2025 13:00:27 +0100
Subject: [PATCH 31/43] STY: isort, ruff

---
 tests/functional/setups/manu_sneddon_2d.py | 71 +++++++++++-----------
 tests/functional/test_sneddon_2d.py        | 46 +++++++-------
 2 files changed, 59 insertions(+), 58 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 62e92cedc8..69c2dc380a 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -1,11 +1,11 @@
+from dataclasses import dataclass
+from typing import Literal
+
 import numpy as np
+
 import porepy as pp
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 from porepy.models.protocol import PorePyModel
-from typing import Literal
-
-
-from dataclasses import dataclass
 
 
 def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray:
@@ -115,7 +115,7 @@ def get_bc_val(
     Return:
         Boundary values for the displacement
     """
-    
+
     # Equations for f2,f3,f4.f5 can be found in book Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics pages 57, 84-92, 168
     f2 = np.zeros(bound_faces.size)
     f3 = np.zeros(bound_faces.size)
@@ -196,7 +196,6 @@ def assign_bem(
     bound_face_centers = g.face_centers[:, bound_faces]
 
     for i in range(0, u_a.size):
-
         new_bound_face_centers = transform(bem_centers[:, i], bound_face_centers, alpha)
 
         u_bound = get_bc_val(
@@ -211,31 +210,30 @@ def assign_bem(
 class ManuExactSneddon2dSetup:
     """
     Class for setting up the analytical solution for the pressurized fracture problem.
-    
-    
+
+
     """
-  
+
     def __init__(self, setup):
         # Initialize private variables from the setup dictionary
         self.p0 = setup.get("p0")
         self.theta = setup.get("theta")
 
         self.a = setup.get("a")
-        self.shear_modulus  = setup.get("material_constants").get("solid").shear_modulus
+        self.shear_modulus = setup.get("material_constants").get("solid").shear_modulus
         self.poi = setup.get("poi")
         self.length = setup.get("length")
         self.height = setup.get("height")
-         
-  
+
     def exact_sol_global(self, sd):
         """
         Compute the analytical solution for the pressurized fracture problem in question.
-        
+
         Parameters:
             sd: Subdomain for which the analytical solution is to be computed.
-            
+
         """
-        
+
         n = 1000
         h = 2 * self.a / n
         box_faces = sd.get_boundary_faces()
@@ -244,15 +242,14 @@ def exact_sol_global(self, sd):
         center = np.array([self.length / 2, self.height / 2, 0])
         bem_centers = get_bem_centers(self.a, h, n, self.theta, center)
         eta = compute_eta(bem_centers, center)
-        
-        
-        u_a = -analytical_displacements(self.a, eta, self.p0, self.shear_modulus, self.poi)
+
+        u_a = -analytical_displacements(
+            self.a, eta, self.p0, self.shear_modulus, self.poi
+        )
         u_bc = assign_bem(sd, h / 2, box_faces, self.theta, bem_centers, u_a, self.poi)
         return u_bc
-    
 
-    def exact_sol_fracture( self,
-    gb: pp.MixedDimensionalGrid) -> tuple:
+    def exact_sol_fracture(self, gb: pp.MixedDimensionalGrid) -> tuple:
         """
         Compute Sneddon's analytical solution for the pressurized crack
         problem in question.
@@ -279,7 +276,9 @@ def exact_sol_fracture( self,
 
         # compute distances from fracture centre with its corresponding apertures
         eta = compute_eta(fracture_faces, fracture_center)
-        apertures = analytical_displacements(self.a, eta, self.p0, self.shear_modulus, self.poi)
+        apertures = analytical_displacements(
+            self.a, eta, self.p0, self.shear_modulus, self.poi
+        )
 
         return eta, apertures
 
@@ -330,11 +329,11 @@ def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
     def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
         """
         Setting displacement boundary condition values.
-        This method sets the displacement boundary condition values for a given 
-        boundary grid using the Sneddon analytical solution through the Boundary 
+        This method sets the displacement boundary condition values for a given
+        boundary grid using the Sneddon analytical solution through the Boundary
         Element Method (BEM).
         Parameters:
-            The boundary grid for which the displacement boundary condition values 
+            The boundary grid for which the displacement boundary condition values
             are to be set.
         Returns:
             An array of displacement boundary condition values.
@@ -344,11 +343,10 @@ def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
         if sd.dim < 2:
             return np.zeros(self.nd * sd.num_faces)
 
-  
-        # Get the analytical solution for the displacement 
+        # Get the analytical solution for the displacement
         exact_sol = ManuExactSneddon2dSetup(self.params)
-        u_exact = exact_sol.exact_sol_global(sd) 
-        
+        u_exact = exact_sol.exact_sol_global(sd)
+
         # Project the values to the grid
         return bg.projection(2) @ u_exact.ravel("F")
 
@@ -389,6 +387,7 @@ def stress_discretization(
 @dataclass
 class SneddonData:
     """Data class for storing the error in the displacement field."""
+
     error_displacement: pp.number
 
 
@@ -397,24 +396,24 @@ class SneddonDataSaving(pp.PorePyModel):
 
     def collect_data(self) -> SneddonData:
         """Collecting the error in the displacement field.
-        
+
         Returns: collected data dictionary
         """
         frac_sd = self.mdg.subdomains(dim=self.nd - 1)
         nd_vec_to_normal = self.normal_component(frac_sd)
-    
+
         # Comuting the numerical displacement jump along the fracture on the fracture cell centers.
         u_n: pp.ad.Operator = nd_vec_to_normal @ self.displacement_jump(frac_sd)
         u_n = u_n.value(self.equation_system)
-        
+
         exact_setup = ManuExactSneddon2dSetup(self.params)
         # Checking convergence specifically on the fracture
         u_a = exact_setup.exact_sol_fracture(self.mdg)[1]
-    
+
         e = ConvergenceAnalysis.l2_error(
             frac_sd[0], u_a, u_n, is_scalar=False, is_cc=True, relative=True
         )
-        
+
         collect_data = SneddonData(error_displacement=e)
         return collect_data
 
@@ -426,8 +425,6 @@ class MomentumBalanceGeometryBC(
     SneddonDataSaving,
     ModifiedBoundaryConditions,
     pp.constitutive_laws.PressureStress,
-    pp.MomentumBalance
+    pp.MomentumBalance,
 ):
     pass
-    
-
diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index c29626880d..a6090dbf2a 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -1,9 +1,11 @@
-import porepy as pp
+import copy
+import math
+
 import numpy as np
 import pytest
-import math
-import copy
-import  tests.functional.setups.manu_sneddon_2d as manu_sneddon_2d
+
+import porepy as pp
+import tests.functional.setups.manu_sneddon_2d as manu_sneddon_2d
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 
 # ----> Set up the material constants
@@ -14,7 +16,7 @@
 )  # Convertion formula from shear modulus and poission to lame lambda parameter
 
 solid = pp.SolidConstants(shear_modulus=shear_modulus, lame_lambda=lam)
-   
+
 
 def compute_frac_pts(
     theta_rad: float, a: float, height: float, length: float
@@ -49,10 +51,10 @@ def compute_frac_pts(
 def actual_ooc() -> dict:
     """
     Prepare parameters and model for the convergence analysis of the moment balance equation.
-    
-    This setup validates the linear elasticity model for the analytical Sneddon solution in 2D, describing the analytical displacement on the fracture. 
+
+    This setup validates the linear elasticity model for the analytical Sneddon solution in 2D, describing the analytical displacement on the fracture.
     The problem consists of a 2D domain with a fracture at a given angle and internal pressure.
-    
+
     Returns:
         A dictionary containing the experimental order of convergence for the displacement.
     """
@@ -68,17 +70,19 @@ def actual_ooc() -> dict:
     params = {
         "prepare_simulation": True,
         "material_constants": {"solid": solid},
-        "a": a, # Half-length of the fracture
-        "height": height, # Height of the domain
-        "length": length, # Length of the domain
-        "p0": 1e-4, # Internal pressure of fracture
-        "poi": poi, # Possion ratio (Not standard in solid constants)
+        "a": a,  # Half-length of the fracture
+        "height": height,  # Height of the domain
+        "length": length,  # Length of the domain
+        "p0": 1e-4,  # Internal pressure of fracture
+        "poi": poi,  # Possion ratio (Not standard in solid constants)
         "meshing_arguments": {"cell_size": 0.03},
         "theta": theta_rad,
     }
 
     # Convert angle to radians and compute fracture points
-    params["frac_pts"] = compute_frac_pts(theta_rad=theta_rad, a=a, height=height, length=length)
+    params["frac_pts"] = compute_frac_pts(
+        theta_rad=theta_rad, a=a, height=height, length=length
+    )
 
     # Model for the convergence analysis
     model = manu_sneddon_2d.MomentumBalanceGeometryBC
@@ -89,19 +93,19 @@ def actual_ooc() -> dict:
         model_params=copy.deepcopy(params),
         levels=2,
         spatial_refinement_rate=2,
-        temporal_refinement_rate=1
+        temporal_refinement_rate=1,
     )
 
     # Calculate and return the order of convergence for the displacement
-    order_dict = conv_analysis.order_of_convergence(conv_analysis.run_analysis(), data_range=slice(None, None, None))
+    order_dict = conv_analysis.order_of_convergence(
+        conv_analysis.run_analysis(), data_range=slice(None, None, None)
+    )
     return order_dict
 
 
-
 def test_order_of_convergence(
     actual_ooc,
 ) -> None:
-    """Test observed order of convergence.
-    """
-    # We  the order of L2 convergence on the fracture of displacement to be about 1.0 
-    assert 0.85 <   actual_ooc["ooc_displacement"]
+    """Test observed order of convergence."""
+    # We  the order of L2 convergence on the fracture of displacement to be about 1.0
+    assert 0.85 < actual_ooc["ooc_displacement"]

From 36a760678b03c6bbf8ae0016be3acde1c92aed92 Mon Sep 17 00:00:00 2001
From: Veljko Lipovac <lipovac.v@outlook.com>
Date: Thu, 13 Feb 2025 13:30:55 +0100
Subject: [PATCH 32/43] MOD: Refactoring Sneddon setup.

---
 tests/functional/setups/manu_sneddon_2d.py | 199 +++++++++++----------
 tests/functional/test_sneddon_2d.py        |   2 +-
 2 files changed, 104 insertions(+), 97 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 69c2dc380a..3b4bc2f0f2 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -1,11 +1,10 @@
 from dataclasses import dataclass
-from typing import Literal
+from typing import Callable, Literal
 
 import numpy as np
 
 import porepy as pp
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
-from porepy.models.protocol import PorePyModel
 
 
 def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray:
@@ -207,15 +206,10 @@ def assign_bem(
     return bc_val
 
 
-class ManuExactSneddon2dSetup:
-    """
-    Class for setting up the analytical solution for the pressurized fracture problem.
-
-
-    """
+class ManuSneddonExactSolution2d:
+    """Class representing the analytical solution for the pressurized fracture problem."""
 
-    def __init__(self, setup):
-        # Initialize private variables from the setup dictionary
+    def __init__(self, setup: dict):
         self.p0 = setup.get("p0")
         self.theta = setup.get("theta")
 
@@ -225,9 +219,9 @@ def __init__(self, setup):
         self.length = setup.get("length")
         self.height = setup.get("height")
 
-    def exact_sol_global(self, sd):
-        """
-        Compute the analytical solution for the pressurized fracture problem in question.
+    def exact_sol_global(self, sd: pp.Grid) -> np.ndarray:
+        """Compute the analytical solution for the pressurized fracture problem in
+        question.
 
         Parameters:
             sd: Subdomain for which the analytical solution is to be computed.
@@ -249,27 +243,26 @@ def exact_sol_global(self, sd):
         u_bc = assign_bem(sd, h / 2, box_faces, self.theta, bem_centers, u_a, self.poi)
         return u_bc
 
-    def exact_sol_fracture(self, gb: pp.MixedDimensionalGrid) -> tuple:
-        """
-        Compute Sneddon's analytical solution for the pressurized crack
+    def exact_sol_fracture(
+        self, mdg: pp.MixedDimensionalGrid
+    ) -> tuple[np.ndarray, np.ndarray]:
+        """Compute Sneddon's analytical solution for the pressurized crack
         problem in question.
 
-        Source: Sneddon Fourier transforms 1951 page 425 eq 92
+        References:
+            Sneddon Fourier transforms 1951 page 425 eq 92
 
-        Parameter
-            gb: Grid object
-            a: Half fracture length of fracture
-            eta: Distance from fracture centre
-            p0: Internal constant pressure inside the fracture
-            G: Shear modulus
-            poi: Poisson ratio
-            height,length: Height and length of domain
+        Parameters:
+            mdg: Mixed-dimensional domain of the setup.
 
         Return:
-            A tuple containing two vectors: a list of distances from the fracture center to each fracture coordinate, and the corresponding analytical apertures.
+            A 2-tuple of numpy arrays, where the first vector contains the distances
+            between fracture center to each fracture coordinate, and the second vector
+            the respective analytical apertures.
+
         """
-        ambient_dim = gb.dim_max()
-        g_1 = gb.subdomains(dim=ambient_dim - 1)[0]
+        ambient_dim = mdg.dim_max()
+        g_1 = mdg.subdomains(dim=ambient_dim - 1)[0]
         fracture_center = np.array([self.length / 2, self.height / 2, 0])
 
         fracture_faces = g_1.cell_centers
@@ -283,9 +276,10 @@ def exact_sol_fracture(self, gb: pp.MixedDimensionalGrid) -> tuple:
         return eta, apertures
 
 
-class ModifiedGeometry(PorePyModel):
+class ManuSneddonGeometry2d(pp.PorePyModel):
     def set_domain(self) -> None:
-        """Defining a two-dimensional square domain with sidelengths length and height."""
+        """Defining a two-dimensional rectangle with parametrizable ``'length'`` and
+        ``'height'`` in the model parameters."""
 
         self._domain = pp.Domain(
             bounding_box={
@@ -301,19 +295,24 @@ def grid_type(self) -> Literal["simplex", "cartesian", "tensor_grid"]:
         return self.params.get("grid_type", "simplex")
 
     def set_fractures(self):
-        """Setting the fractures in the domain from params."""
+        """Setting a single line fracture in the domain using ``params['frac_pts']``."""
 
         points = self.params["frac_pts"]
         self._fractures = [pp.LineFracture(points)]
 
 
-class ModifiedBoundaryConditions(PorePyModel):
+class ManuSneddonBoundaryConditions(pp.PorePyModel):
+    exact_sol: ManuSneddonExactSolution2d
+
     def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
         """Set boundary condition type for the problem.
+
         Parameters:
             sd: Subdomain for which the boundary conditions are to be set.
+
         Returns:
             Boundary condition type for the problem.
+
         """
         bounds = self.domain_boundary_sides(sd)
 
@@ -327,104 +326,112 @@ def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
         return bc
 
     def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
-        """
-        Setting displacement boundary condition values.
-        This method sets the displacement boundary condition values for a given
+        """This method sets the displacement boundary condition values for a given
         boundary grid using the Sneddon analytical solution through the Boundary
         Element Method (BEM).
-        Parameters:
-            The boundary grid for which the displacement boundary condition values
-            are to be set.
-        Returns:
-            An array of displacement boundary condition values.
-        """
-
-        sd = bg.parent
-        if sd.dim < 2:
-            return np.zeros(self.nd * sd.num_faces)
-
-        # Get the analytical solution for the displacement
-        exact_sol = ManuExactSneddon2dSetup(self.params)
-        u_exact = exact_sol.exact_sol_global(sd)
-
-        # Project the values to the grid
-        return bg.projection(2) @ u_exact.ravel("F")
-
-
-class PressureStressMixin(PorePyModel):
-    stress_keyword: str
-
-    def pressure(self, subdomains: pp.SubdomainsOrBoundaries):
-        """Discretization of the stress tensor.
 
         Parameters:
-            subdomains: List of subdomains where the stress is defined.
+            bg: The boundary grid for which the displacement boundary condition values
+                are to be set.
 
         Returns:
-            Discretization operator for the stress tensor.
+            An array of displacement values on the boundary.
 
         """
-        # Return pressure in the fracture domain
-        mat = self.params["p0"] * np.ones(
-            sum((subdomain.num_cells for subdomain in subdomains))
-        )
-        return pp.ad.DenseArray(mat)
-
-    def stress_discretization(
-        self, subdomains: list[pp.Grid]
-    ) -> pp.ad.BiotAd | pp.ad.MpsaAd:
-        """Discretization of the stress tensor.
 
-        Parameters:
-            subdomains: List of subdomains where the stress is defined.
+        sd = bg.parent
+        if sd.dim < 2:
+            return np.zeros(self.nd * sd.num_faces)
+        else:
+            # Get the analytical solution for the displacement
+            u_exact = self.exact_sol.exact_sol_global(sd)
 
-        Returns:
-            Discretization operator for the stress tensor.
-        """
-        return pp.ad.MpsaAd(self.stress_keyword, subdomains)
+            # Project the values to the grid
+            return bg.projection(2) @ u_exact.ravel("F")
 
 
 @dataclass
-class SneddonData:
+class ManuSneddonSaveData:
     """Data class for storing the error in the displacement field."""
 
     error_displacement: pp.number
 
 
-class SneddonDataSaving(pp.PorePyModel):
+class ManuSneddonDataSaving(pp.PorePyModel):
     """Class for saving the error in the displacement field."""
 
-    def collect_data(self) -> SneddonData:
-        """Collecting the error in the displacement field.
+    exact_sol: ManuSneddonExactSolution2d
 
-        Returns: collected data dictionary
-        """
+    displacement_jump: Callable[[list[pp.Grid]], pp.ad.Operator]
+
+    def collect_data(self) -> ManuSneddonSaveData:
+        """Collecting the error in the displacement field."""
         frac_sd = self.mdg.subdomains(dim=self.nd - 1)
         nd_vec_to_normal = self.normal_component(frac_sd)
 
-        # Comuting the numerical displacement jump along the fracture on the fracture cell centers.
+        # Computing the numerical displacement jump along the fracture on the fracture
+        # cell centers.
         u_n: pp.ad.Operator = nd_vec_to_normal @ self.displacement_jump(frac_sd)
         u_n = u_n.value(self.equation_system)
 
-        exact_setup = ManuExactSneddon2dSetup(self.params)
         # Checking convergence specifically on the fracture
-        u_a = exact_setup.exact_sol_fracture(self.mdg)[1]
+        u_a = self.exact_sol.exact_sol_fracture(self.mdg)[1]
 
-        e = ConvergenceAnalysis.l2_error(
+        e = ConvergenceAnalysis.lp_error(
             frac_sd[0], u_a, u_n, is_scalar=False, is_cc=True, relative=True
         )
 
-        collect_data = SneddonData(error_displacement=e)
+        collect_data = ManuSneddonSaveData(error_displacement=e)
         return collect_data
 
 
-class MomentumBalanceGeometryBC(
-    PorePyModel,
-    PressureStressMixin,
-    ModifiedGeometry,
-    SneddonDataSaving,
-    ModifiedBoundaryConditions,
-    pp.constitutive_laws.PressureStress,
+class ManuSneddonConstitutiveLaws(pp.constitutive_laws.PressureStress):
+    """Inherits the relevant pressure-stress relations, defines a constant pressure
+    and sets the stress discretization to MPSA.
+
+    Note that the latter is relevant because the inherited constitutive law
+    uses the Biot discretization.
+
+    """
+
+    def pressure(self, domains: pp.SubdomainsOrBoundaries):
+        """Defines a constant pressure given by ``params['p0']``.
+
+        Parameters:
+            subdomains: List of grids or boundary grids where the pressure is defined.
+
+        Returns:
+            A dense array broadcasting the constant pressure into a suitable array.
+
+        """
+        p = self.params["p0"] * np.ones(sum((grid.num_cells for grid in domains)))
+        return pp.ad.DenseArray(p)
+
+    def stress_discretization(self, subdomains: list[pp.Grid]) -> pp.ad.MpsaAd:
+        """Discretization class for the stress tensor.
+
+        Parameters:
+            subdomains: List of subdomains where the stress is defined.
+
+        Returns:
+            The MPSA discretization in Ad operator form.
+
+        """
+        return pp.ad.MpsaAd(self.stress_keyword, subdomains)
+
+
+class ManuSneddonSetup2d(
+    ManuSneddonGeometry2d,
+    ManuSneddonDataSaving,
+    ManuSneddonBoundaryConditions,
+    ManuSneddonConstitutiveLaws,
     pp.MomentumBalance,
 ):
-    pass
+    """Complete Sneddon setup, including data collection and the analytical solution."""
+
+    exact_sol: ManuSneddonExactSolution2d
+
+    def set_materials(self):
+        """Initiating additionally the exact solution."""
+        super().set_materials()
+        self.exact_sol = ManuSneddonExactSolution2d(self.params)
diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index a6090dbf2a..069755dd74 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -85,7 +85,7 @@ def actual_ooc() -> dict:
     )
 
     # Model for the convergence analysis
-    model = manu_sneddon_2d.MomentumBalanceGeometryBC
+    model = manu_sneddon_2d.ManuSneddonSetup2d
 
     # Convergence analysis setup
     conv_analysis = ConvergenceAnalysis(

From 19a3a27e0ce0d624b75cf3c6275a46224e0b5670 Mon Sep 17 00:00:00 2001
From: Veljko Lipovac <lipovac.v@outlook.com>
Date: Thu, 13 Feb 2025 13:42:15 +0100
Subject: [PATCH 33/43] STY: Reworking documentation of BEM methods for tests.

---
 tests/functional/setups/manu_sneddon_2d.py | 134 +++++++++++----------
 1 file changed, 71 insertions(+), 63 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 3b4bc2f0f2..6c99de2529 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -5,41 +5,42 @@
 
 import porepy as pp
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
+from porepy.geometry.distances import point_pointset
 
 
 def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray:
-    """
-    Compute the distance fracture points to the fracture centre.
+    """Compute the distance from fracture points to the fracture centre.
 
-    Parameter::
-    pointset: Array containing coordinates on the fracture
-    center: fracture centre
+    Parameter:
+        pointset_centers: Coordinates of fracture points.
+        center: Coordinates of the fracture center.
 
     Return:
-        Array of distances each point in pointset to the center
+        Array of distances each point in pointset to the center.
 
     """
-    return pp.geometry.distances.point_pointset(pointset_centers, center)
+    return point_pointset(pointset_centers, center)
 
 
 def get_bem_centers(
     a: float, h: float, n: int, theta: float, center: np.ndarray
 ) -> np.ndarray:
-    """
-    Compute coordinates of the centers of the bem segments
+    """Compute coordinates of the centers of the BEM segments.
 
     Parameters:
-        a: half fracture length
-        h: bem segment length
-        n: number of bem segments
-        theta: orientation of the fracture
-        center: center of the fracture
+        a: Half of the fracture length.
+        h: BEM segment length.
+        n: Number of BEM segments.
+        theta: Orientation of the fracture.
+        center: Coordinates of the fracture center.
 
     Return:
-        Array of BEM centers
+        Array containing the BEM segment centers.
+
     """
 
-    # Coordinate system (4.5.1) in page 57 in in book Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics
+    # Coordinate system (4.5.1) in page 57 in in book
+    # Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics
 
     bem_centers = np.zeros((3, n))
     x_0 = center[0] - (a - 0.5 * h) * np.sin(theta)
@@ -54,45 +55,48 @@ def get_bem_centers(
 def analytical_displacements(
     a: float, eta: np.ndarray, p0: float, G: float, poi: float
 ) -> np.ndarray:
-    """
-    Compute Sneddon's analytical solution for the pressurized fracture.
+    """Compute Sneddon's analytical solution for the pressurized fracture displacement.
 
-    Source: Sneddon Fourier transforms 1951 page 425 eq 92
+    References:
+        Sneddon Fourier transforms 1951 page 425 eq 92
 
-    Parameter:
-        a: half fracture length
-        eta: distance from fracture centre
-        p0: pressure
-        G: shear modulus
-        poi: poisson ratio
+    Parameters:
+        a: Half of the fracture length.
+        eta: Distances of fracture points to the fracture center.
+        p0: Constant, fixed pressure.
+        G: Shear modulus.
+        poi: Poisson ratio.
 
     Return
-        List of analytical normal displacement jumps.
+        Array containing the analytical normal displacement jumps.
+
     """
     cons = (1 - poi) / G * p0 * a * 2
     return cons * np.sqrt(1 - np.power(eta / a, 2))
 
 
 def transform(xc: np.ndarray, x: np.ndarray, alpha: float) -> np.ndarray:
-    """
-    Translation and rotation coordinate transformation of boundary face coordinates for the BEM method
+    """Translation and rotation transform of boundary face coordinates for the BEM.
+
+    Parameters:
+        xc: Coordinates of BEM segment centre.
+        x: Coordinates of boundary faces.
+        alpha: Fracture orientation.
 
-    Parameter
-        xc: Coordinates of BEM segment centre
-        x: Coordinates of boundary faces
-        alpha: Fracture orientation
     Return:
-        Transformed coordinates
+        Transformed coordinates.
+
     """
     x_bar = np.zeros_like(x)
-    # Terms in (7.4.6) in Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics page 168
+    # Terms in (7.4.6) in
+    # Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics page 168
     x_bar[0, :] = (x[0, :] - xc[0]) * np.cos(alpha) + (x[1, :] - xc[1]) * np.sin(alpha)
     x_bar[1, :] = -(x[0, :] - xc[0]) * np.sin(alpha) + (x[1, :] - xc[1]) * np.cos(alpha)
     return x_bar
 
 
 def get_bc_val(
-    g: pp.Grid,
+    grid: pp.Grid,
     bound_faces: np.ndarray,
     xf: np.ndarray,
     h: float,
@@ -100,28 +104,31 @@ def get_bc_val(
     alpha: float,
     du: float,
 ) -> np.ndarray:
-    """
-    Compute semi-analytical displacement on the boundary using the BEM method for the Sneddon problem.
+    """Computes semi-analytical displacement values on the boundary using the BEM for
+    the Sneddon problem.
 
     Parameter
-        g: Grid object
-        bound_faces: Index lists for boundary faces
-        xf: Coordinates of boundary faces
-        h: BEM segment length
-        poi: Poisson ratio
-        alpha: Fracture orientation
-        du: Sneddon's analytical relative normal displacement
+        grid: The matrix grid.
+        bound_faces: Array of indices of boundary faces of ``grid``.
+        xf: Coordinates of boundary faces.
+        h: BEM segment length.
+        poi: Poisson ratio.
+        alpha: Fracture orientation.
+        du: Sneddon's analytical relative normal displacement.
+
     Return:
-        Boundary values for the displacement
+        Boundary values for the displacement.
+
     """
 
-    # Equations for f2,f3,f4.f5 can be found in book Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics pages 57, 84-92, 168
+    # Equations for f2,f3,f4.f5 can be found in book
+    # Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics pages 57, 84-92, 168
     f2 = np.zeros(bound_faces.size)
     f3 = np.zeros(bound_faces.size)
     f4 = np.zeros(bound_faces.size)
     f5 = np.zeros(bound_faces.size)
 
-    u = np.zeros((g.dim, g.num_faces))
+    u = np.zeros((grid.dim, grid.num_faces))
 
     # Constant term in (7.4.5)
     m = 1 / (4 * np.pi * (1 - poi))
@@ -164,41 +171,42 @@ def get_bc_val(
 
 
 def assign_bem(
-    g: pp.Grid,
+    grid: pp.Grid,
     h: float,
     bound_faces: np.ndarray,
-    theta: float,
+    alpha: float,
     bem_centers: np.ndarray,
     u_a: np.ndarray,
     poi: float,
 ) -> np.ndarray:
-    """
-    Compute analytical displacement using the BEM method for the Sneddon problem.
+    """Computes semi-analytical displacement values using the BEM for the Sneddon
+    problem.
 
     Parameter
-        g: Subdomain grid
-        h: bem segment length
-        bound_faces: boundary faces
-        theta: fracture orientation
-        bem_centers: bem segments centers
-        u_a: Sneddon's analytical relative normal displacement
-        poi: Poisson ratio
+        grid: The matrix grid.
+        h: BEM segment length.
+        bound_faces: Array of indices of boundary faces of ``grid``.
+        alpha: Fracture orientation.
+        bem_centers: BEM segments centers.
+        u_a: Sneddon's analytical relative normal displacement.
+        poi: Poisson ratio.
 
     Return:
-        Semi-analytical boundary displacement values
+        Semi-analytical boundary displacement values.
+
     """
 
-    bc_val = np.zeros((g.dim, g.num_faces))
+    bc_val = np.zeros((grid.dim, grid.num_faces))
 
-    alpha = np.pi / 2 - theta
+    alpha = np.pi / 2 - alpha
 
-    bound_face_centers = g.face_centers[:, bound_faces]
+    bound_face_centers = grid.face_centers[:, bound_faces]
 
     for i in range(0, u_a.size):
         new_bound_face_centers = transform(bem_centers[:, i], bound_face_centers, alpha)
 
         u_bound = get_bc_val(
-            g, bound_faces, new_bound_face_centers, h, poi, alpha, u_a[i]
+            grid, bound_faces, new_bound_face_centers, h, poi, alpha, u_a[i]
         )
 
         bc_val += u_bound

From e3085a48154ea03e4647fa684ebe8356f9c3fc57 Mon Sep 17 00:00:00 2001
From: Veljko Lipovac <lipovac.v@outlook.com>
Date: Thu, 13 Feb 2025 13:44:12 +0100
Subject: [PATCH 34/43] MIN: Renaming init arg of Sneddon exact solution.

---
 tests/functional/setups/manu_sneddon_2d.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 6c99de2529..a106e08832 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -217,15 +217,15 @@ def assign_bem(
 class ManuSneddonExactSolution2d:
     """Class representing the analytical solution for the pressurized fracture problem."""
 
-    def __init__(self, setup: dict):
-        self.p0 = setup.get("p0")
-        self.theta = setup.get("theta")
-
-        self.a = setup.get("a")
-        self.shear_modulus = setup.get("material_constants").get("solid").shear_modulus
-        self.poi = setup.get("poi")
-        self.length = setup.get("length")
-        self.height = setup.get("height")
+    def __init__(self, params: dict):
+        self.p0 = params.get("p0")
+        self.theta = params.get("theta")
+
+        self.a = params.get("a")
+        self.shear_modulus = params.get("material_constants").get("solid").shear_modulus
+        self.poi = params.get("poi")
+        self.length = params.get("length")
+        self.height = params.get("height")
 
     def exact_sol_global(self, sd: pp.Grid) -> np.ndarray:
         """Compute the analytical solution for the pressurized fracture problem in

From 5552527bd0ef57d4941924aa1629cd6067b4a53a Mon Sep 17 00:00:00 2001
From: Veljko Lipovac <lipovac.v@outlook.com>
Date: Thu, 13 Feb 2025 13:52:46 +0100
Subject: [PATCH 35/43] MIN: Consistent naming of fracture orientation.

---
 tests/functional/setups/manu_sneddon_2d.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index a106e08832..294b6b52a0 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -23,7 +23,7 @@ def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray:
 
 
 def get_bem_centers(
-    a: float, h: float, n: int, theta: float, center: np.ndarray
+    a: float, h: float, n: int, alpha: float, center: np.ndarray
 ) -> np.ndarray:
     """Compute coordinates of the centers of the BEM segments.
 
@@ -31,7 +31,7 @@ def get_bem_centers(
         a: Half of the fracture length.
         h: BEM segment length.
         n: Number of BEM segments.
-        theta: Orientation of the fracture.
+        alpha: Orientation of the fracture.
         center: Coordinates of the fracture center.
 
     Return:
@@ -43,11 +43,11 @@ def get_bem_centers(
     # Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics
 
     bem_centers = np.zeros((3, n))
-    x_0 = center[0] - (a - 0.5 * h) * np.sin(theta)
-    y_0 = center[1] - (a - 0.5 * h) * np.cos(theta)
+    x_0 = center[0] - (a - 0.5 * h) * np.sin(alpha)
+    y_0 = center[1] - (a - 0.5 * h) * np.cos(alpha)
     for i in range(0, n):
-        bem_centers[0, i] = x_0 + i * h * np.sin(theta)
-        bem_centers[1, i] = y_0 + i * h * np.cos(theta)
+        bem_centers[0, i] = x_0 + i * h * np.sin(alpha)
+        bem_centers[1, i] = y_0 + i * h * np.cos(alpha)
 
     return bem_centers
 

From 7992f9df8cbd1d9aa32b4b8f00cb91ad827f5237 Mon Sep 17 00:00:00 2001
From: Veljko Lipovac <lipovac.v@outlook.com>
Date: Thu, 13 Feb 2025 14:00:37 +0100
Subject: [PATCH 36/43] MIN: Formatting Sneddon test module.

---
 tests/functional/test_sneddon_2d.py | 27 +++++++++++++--------------
 1 file changed, 13 insertions(+), 14 deletions(-)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 069755dd74..01c0459bc4 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -8,7 +8,7 @@
 import tests.functional.setups.manu_sneddon_2d as manu_sneddon_2d
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 
-# ----> Set up the material constants
+# Set up the material constants
 poi = 0.25
 shear_modulus = 1
 lam = (
@@ -21,8 +21,7 @@
 def compute_frac_pts(
     theta_rad: float, a: float, height: float, length: float
 ) -> np.ndarray:
-    """
-    Assuming the fracture center is at the coordinate (height/2, length/2),
+    """Assuming the fracture center is at the coordinate (height/2, length/2),
     compute the endpoints of a fracture given its orientation and fracture length.
 
     Parameters:
@@ -32,8 +31,9 @@ def compute_frac_pts(
         length: Width of the domain.
 
     Returns:
-        frac_pts : A 2x2 array where each column represents the coordinates of an end point of the fracture in 2D.
-            The first column corresponds to one end point, and the second column corresponds to the other.
+        A 2x2 array where each column represents the coordinates of an end point of the
+        fracture in 2D. The first column corresponds to one end point, and the second
+        column corresponds to the other.
 
     """
     # Rotate the fracture with an angle theta_rad
@@ -46,17 +46,18 @@ def compute_frac_pts(
     return frac_pts
 
 
-# ----> Retrieve actual order of convergence
 @pytest.fixture(scope="module")
 def actual_ooc() -> dict:
-    """
-    Prepare parameters and model for the convergence analysis of the moment balance equation.
+    """Performs convergence analysis of the Sneddon setup.
 
-    This setup validates the linear elasticity model for the analytical Sneddon solution in 2D, describing the analytical displacement on the fracture.
-    The problem consists of a 2D domain with a fracture at a given angle and internal pressure.
+    This setup validates the linear elasticity model for the analytical Sneddon solution
+    in 2D, describing the analytical displacement on the fracture. The problem consists
+    of a 2D domain with a fracture at a given angle and internal pressure.
 
     Returns:
-        A dictionary containing the experimental order of convergence for the displacement.
+        A dictionary containing the experimental order of convergence for the
+        displacement.
+
     """
     # Angle of the fracture in degrees
 
@@ -103,9 +104,7 @@ def actual_ooc() -> dict:
     return order_dict
 
 
-def test_order_of_convergence(
-    actual_ooc,
-) -> None:
+def test_order_of_convergence(actual_ooc: dict) -> None:
     """Test observed order of convergence."""
     # We  the order of L2 convergence on the fracture of displacement to be about 1.0
     assert 0.85 < actual_ooc["ooc_displacement"]

From 54de9293fc2aad3c355aa76dbd72a4d06ccb79ef Mon Sep 17 00:00:00 2001
From: Veljko Lipovac <lipovac.v@outlook.com>
Date: Thu, 13 Feb 2025 15:34:46 +0100
Subject: [PATCH 37/43] MOD: Applying suggestions from review.

---
 tests/functional/setups/manu_sneddon_2d.py | 116 +++++++++------------
 tests/functional/test_sneddon_2d.py        |  11 +-
 2 files changed, 52 insertions(+), 75 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 294b6b52a0..f7aa7d1d77 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -1,10 +1,24 @@
+"""Model setup for the fracturized pressure problem using Sneddon's analytical solution.
+
+References:
+
+    - [1] Starfield, C.: Boundary Element Methods in Solid Mechanics (1983)
+    - [2] Sneddon, I.N.: Fourier transforms (1951)
+
+"""
+
+from __future__ import annotations
+
 from dataclasses import dataclass
-from typing import Callable, Literal
+from typing import Callable
 
 import numpy as np
 
 import porepy as pp
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
+from porepy.applications.md_grids.model_geometries import (
+    SquareDomainOrthogonalFractures,
+)
 from porepy.geometry.distances import point_pointset
 
 
@@ -39,9 +53,7 @@ def get_bem_centers(
 
     """
 
-    # Coordinate system (4.5.1) in page 57 in in book
-    # Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics
-
+    # See [1], p. 57 coordinate system (4.5.1)
     bem_centers = np.zeros((3, n))
     x_0 = center[0] - (a - 0.5 * h) * np.sin(alpha)
     y_0 = center[1] - (a - 0.5 * h) * np.cos(alpha)
@@ -57,9 +69,6 @@ def analytical_displacements(
 ) -> np.ndarray:
     """Compute Sneddon's analytical solution for the pressurized fracture displacement.
 
-    References:
-        Sneddon Fourier transforms 1951 page 425 eq 92
-
     Parameters:
         a: Half of the fracture length.
         eta: Distances of fracture points to the fracture center.
@@ -71,6 +80,7 @@ def analytical_displacements(
         Array containing the analytical normal displacement jumps.
 
     """
+    # See [2], p. 425 equ. 92
     cons = (1 - poi) / G * p0 * a * 2
     return cons * np.sqrt(1 - np.power(eta / a, 2))
 
@@ -88,15 +98,14 @@ def transform(xc: np.ndarray, x: np.ndarray, alpha: float) -> np.ndarray:
 
     """
     x_bar = np.zeros_like(x)
-    # Terms in (7.4.6) in
-    # Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics page 168
+    # See [1], p. 168 equ. 7.4.6
     x_bar[0, :] = (x[0, :] - xc[0]) * np.cos(alpha) + (x[1, :] - xc[1]) * np.sin(alpha)
     x_bar[1, :] = -(x[0, :] - xc[0]) * np.sin(alpha) + (x[1, :] - xc[1]) * np.cos(alpha)
     return x_bar
 
 
 def get_bc_val(
-    grid: pp.Grid,
+    sd: pp.Grid,
     bound_faces: np.ndarray,
     xf: np.ndarray,
     h: float,
@@ -108,8 +117,8 @@ def get_bc_val(
     the Sneddon problem.
 
     Parameter
-        grid: The matrix grid.
-        bound_faces: Array of indices of boundary faces of ``grid``.
+        sd: The matrix grid.
+        bound_faces: Array of indices of boundary faces of ``sd``.
         xf: Coordinates of boundary faces.
         h: BEM segment length.
         poi: Poisson ratio.
@@ -120,32 +129,25 @@ def get_bc_val(
         Boundary values for the displacement.
 
     """
-
-    # Equations for f2,f3,f4.f5 can be found in book
-    # Crouch Starfield 1983 Boundary Element Methods in Solid Mechanics pages 57, 84-92, 168
+    # See [1], pages 57, 84-92, 168
     f2 = np.zeros(bound_faces.size)
     f3 = np.zeros(bound_faces.size)
     f4 = np.zeros(bound_faces.size)
     f5 = np.zeros(bound_faces.size)
 
-    u = np.zeros((grid.dim, grid.num_faces))
+    u = np.zeros((sd.dim, sd.num_faces))
 
-    # Constant term in (7.4.5)
+    # See [1], equ. (7.4.5)
     m = 1 / (4 * np.pi * (1 - poi))
-
-    # Second term in (7.4.5)
     f2[:] = m * (
         np.log(np.sqrt((xf[0, :] - h) ** 2 + xf[1] ** 2))
         - np.log(np.sqrt((xf[0, :] + h) ** 2 + xf[1] ** 2))
     )
-
-    # First term in (7.4.5)
     f3[:] = -m * (
         np.arctan2(xf[1, :], (xf[0, :] - h)) - np.arctan2(xf[1, :], (xf[0, :] + h))
     )
 
-    # The following equalities can be found on page 91 where f3,f4 as equation (5.5.3) and ux, uy in (5.5.1)
-    # Also this is in the coordinate system described in (4.5.1) at page 57, which essentially is a rotation.
+    # See [1], equations 5.5.3, 5.5.1, and p. 57, equ. 4.5.1
     f4[:] = m * (
         xf[1, :] / ((xf[0, :] - h) ** 2 + xf[1, :] ** 2)
         - xf[1, :] / ((xf[0, :] + h) ** 2 + xf[1, :] ** 2)
@@ -214,18 +216,19 @@ def assign_bem(
     return bc_val
 
 
-class ManuSneddonExactSolution2d:
+class SneddonExactSolution2d:
     """Class representing the analytical solution for the pressurized fracture problem."""
 
-    def __init__(self, params: dict):
-        self.p0 = params.get("p0")
-        self.theta = params.get("theta")
+    def __init__(self, model: "ManuSneddonSetup2d"):
+        self.p0 = model.params.get("p0")
+        self.theta = model.params.get("theta")
 
-        self.a = params.get("a")
-        self.shear_modulus = params.get("material_constants").get("solid").shear_modulus
-        self.poi = params.get("poi")
-        self.length = params.get("length")
-        self.height = params.get("height")
+        self.a = model.params.get("a")
+        self.shear_modulus = (
+            model.params.get("material_constants").get("solid").shear_modulus
+        )
+        self.poi = model.params.get("poi")
+        self.length = model.domain_size
 
     def exact_sol_global(self, sd: pp.Grid) -> np.ndarray:
         """Compute the analytical solution for the pressurized fracture problem in
@@ -241,7 +244,7 @@ def exact_sol_global(self, sd: pp.Grid) -> np.ndarray:
         box_faces = sd.get_boundary_faces()
         u_bc = np.zeros((sd.dim, sd.num_faces))
 
-        center = np.array([self.length / 2, self.height / 2, 0])
+        center = np.array([self.length / 2, self.length / 2, 0])
         bem_centers = get_bem_centers(self.a, h, n, self.theta, center)
         eta = compute_eta(bem_centers, center)
 
@@ -257,9 +260,6 @@ def exact_sol_fracture(
         """Compute Sneddon's analytical solution for the pressurized crack
         problem in question.
 
-        References:
-            Sneddon Fourier transforms 1951 page 425 eq 92
-
         Parameters:
             mdg: Mixed-dimensional domain of the setup.
 
@@ -269,13 +269,13 @@ def exact_sol_fracture(
             the respective analytical apertures.
 
         """
+
         ambient_dim = mdg.dim_max()
         g_1 = mdg.subdomains(dim=ambient_dim - 1)[0]
-        fracture_center = np.array([self.length / 2, self.height / 2, 0])
-
+        fracture_center = np.array([self.length / 2, self.length / 2, 0])
         fracture_faces = g_1.cell_centers
 
-        # compute distances from fracture centre with its corresponding apertures
+        # See [2], p. 425, equ. 92
         eta = compute_eta(fracture_faces, fracture_center)
         apertures = analytical_displacements(
             self.a, eta, self.p0, self.shear_modulus, self.poi
@@ -284,33 +284,17 @@ def exact_sol_fracture(
         return eta, apertures
 
 
-class ManuSneddonGeometry2d(pp.PorePyModel):
-    def set_domain(self) -> None:
-        """Defining a two-dimensional rectangle with parametrizable ``'length'`` and
-        ``'height'`` in the model parameters."""
-
-        self._domain = pp.Domain(
-            bounding_box={
-                "xmin": 0,
-                "ymin": 0,
-                "xmax": self.params["length"],
-                "ymax": self.params["height"],
-            }
-        )
-
-    def grid_type(self) -> Literal["simplex", "cartesian", "tensor_grid"]:
-        """Choosing the grid type for our domain."""
-        return self.params.get("grid_type", "simplex")
+class ManuSneddonGeometry2d(SquareDomainOrthogonalFractures):
+    """Square domain but with single line fracture."""
 
     def set_fractures(self):
         """Setting a single line fracture in the domain using ``params['frac_pts']``."""
 
-        points = self.params["frac_pts"]
-        self._fractures = [pp.LineFracture(points)]
+        self._fractures = [pp.LineFracture(self.params["frac_pts"])]
 
 
 class ManuSneddonBoundaryConditions(pp.PorePyModel):
-    exact_sol: ManuSneddonExactSolution2d
+    exact_sol: SneddonExactSolution2d
 
     def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
         """Set boundary condition type for the problem.
@@ -335,8 +319,7 @@ def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
 
     def bc_values_displacement(self, bg: pp.BoundaryGrid) -> np.ndarray:
         """This method sets the displacement boundary condition values for a given
-        boundary grid using the Sneddon analytical solution through the Boundary
-        Element Method (BEM).
+        boundary grid using the Sneddon analytical solution through the BEM.
 
         Parameters:
             bg: The boundary grid for which the displacement boundary condition values
@@ -368,7 +351,7 @@ class ManuSneddonSaveData:
 class ManuSneddonDataSaving(pp.PorePyModel):
     """Class for saving the error in the displacement field."""
 
-    exact_sol: ManuSneddonExactSolution2d
+    exact_sol: SneddonExactSolution2d
 
     displacement_jump: Callable[[list[pp.Grid]], pp.ad.Operator]
 
@@ -437,9 +420,6 @@ class ManuSneddonSetup2d(
 ):
     """Complete Sneddon setup, including data collection and the analytical solution."""
 
-    exact_sol: ManuSneddonExactSolution2d
-
-    def set_materials(self):
-        """Initiating additionally the exact solution."""
-        super().set_materials()
-        self.exact_sol = ManuSneddonExactSolution2d(self.params)
+    def __init__(self, params=None):
+        super().__init__(params)
+        self.exact_sol = SneddonExactSolution2d(self)
diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 01c0459bc4..03ad81f15e 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -5,8 +5,8 @@
 import pytest
 
 import porepy as pp
-import tests.functional.setups.manu_sneddon_2d as manu_sneddon_2d
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
+from tests.functional.setups.manu_sneddon_2d import ManuSneddonSetup2d
 
 # Set up the material constants
 poi = 0.25
@@ -72,11 +72,11 @@ def actual_ooc() -> dict:
         "prepare_simulation": True,
         "material_constants": {"solid": solid},
         "a": a,  # Half-length of the fracture
-        "height": height,  # Height of the domain
-        "length": length,  # Length of the domain
+        "domain_size": height,  # Length of square domain
         "p0": 1e-4,  # Internal pressure of fracture
         "poi": poi,  # Possion ratio (Not standard in solid constants)
         "meshing_arguments": {"cell_size": 0.03},
+        "grid_type": "simplex",
         "theta": theta_rad,
     }
 
@@ -85,12 +85,9 @@ def actual_ooc() -> dict:
         theta_rad=theta_rad, a=a, height=height, length=length
     )
 
-    # Model for the convergence analysis
-    model = manu_sneddon_2d.ManuSneddonSetup2d
-
     # Convergence analysis setup
     conv_analysis = ConvergenceAnalysis(
-        model_class=model,
+        model_class=ManuSneddonSetup2d,
         model_params=copy.deepcopy(params),
         levels=2,
         spatial_refinement_rate=2,

From 011e03966cbfccc023ad81ea21459cb20429dc22 Mon Sep 17 00:00:00 2001
From: Yury Zabegaev <yuri.zabegaev@uib.no>
Date: Fri, 14 Feb 2025 10:11:59 +0000
Subject: [PATCH 38/43] MAINT: Suggestiongs from the review

---
 tests/functional/setups/manu_sneddon_2d.py | 70 ++++++++++++++++------
 tests/functional/test_sneddon_2d.py        | 55 ++++-------------
 2 files changed, 61 insertions(+), 64 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index f7aa7d1d77..ec0cc81602 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -57,9 +57,9 @@ def get_bem_centers(
     bem_centers = np.zeros((3, n))
     x_0 = center[0] - (a - 0.5 * h) * np.sin(alpha)
     y_0 = center[1] - (a - 0.5 * h) * np.cos(alpha)
-    for i in range(0, n):
-        bem_centers[0, i] = x_0 + i * h * np.sin(alpha)
-        bem_centers[1, i] = y_0 + i * h * np.cos(alpha)
+    i = np.arange(n)
+    bem_centers[0] = x_0 + i * h * np.sin(alpha)
+    bem_centers[1] = y_0 + i * h * np.cos(alpha)
 
     return bem_centers
 
@@ -220,14 +220,13 @@ class SneddonExactSolution2d:
     """Class representing the analytical solution for the pressurized fracture problem."""
 
     def __init__(self, model: "ManuSneddonSetup2d"):
-        self.p0 = model.params.get("p0")
-        self.theta = model.params.get("theta")
+        self.p0 = model.params["p0"]
+        self.theta_rad = model.params["theta_rad"]
 
-        self.a = model.params.get("a")
-        self.shear_modulus = (
-            model.params.get("material_constants").get("solid").shear_modulus
-        )
-        self.poi = model.params.get("poi")
+        self.a = model.params["a"]
+        self.shear_modulus = model.params["material_constants"]["solid"].shear_modulus
+        self.poi = model.params["poi"]
+        self.n = model.params["num_bem_segments"]
         self.length = model.domain_size
 
     def exact_sol_global(self, sd: pp.Grid) -> np.ndarray:
@@ -239,19 +238,20 @@ def exact_sol_global(self, sd: pp.Grid) -> np.ndarray:
 
         """
 
-        n = 1000
-        h = 2 * self.a / n
+        h = 2 * self.a / self.n
         box_faces = sd.get_boundary_faces()
         u_bc = np.zeros((sd.dim, sd.num_faces))
 
         center = np.array([self.length / 2, self.length / 2, 0])
-        bem_centers = get_bem_centers(self.a, h, n, self.theta, center)
+        bem_centers = get_bem_centers(self.a, h, self.n, self.theta_rad, center)
         eta = compute_eta(bem_centers, center)
 
         u_a = -analytical_displacements(
             self.a, eta, self.p0, self.shear_modulus, self.poi
         )
-        u_bc = assign_bem(sd, h / 2, box_faces, self.theta, bem_centers, u_a, self.poi)
+        u_bc = assign_bem(
+            sd, h / 2, box_faces, self.theta_rad, bem_centers, u_a, self.poi
+        )
         return u_bc
 
     def exact_sol_fracture(
@@ -288,9 +288,43 @@ class ManuSneddonGeometry2d(SquareDomainOrthogonalFractures):
     """Square domain but with single line fracture."""
 
     def set_fractures(self):
-        """Setting a single line fracture in the domain using ``params['frac_pts']``."""
+        """Setting a single line fracture which runs through the domain center."""
+
+        theta_rad = self.params["theta_rad"]
+        a = self.params["a"]  # Half-length of the fracture
+        height = length = self.params["domain_size"]
+        frac_pts = self.compute_frac_pts(
+            theta_rad=theta_rad, a=a, height=height, length=length
+        )
+        self._fractures = [pp.LineFracture(frac_pts)]
+
+    @staticmethod
+    def compute_frac_pts(
+        theta_rad: float, a: float, height: float, length: float
+    ) -> np.ndarray:
+        """Assuming the fracture center is at the coordinate (height/2, length/2),
+        compute the endpoints of a fracture given its orientation and fracture length.
+
+        Parameters:
+            theta_rad: Angle of the fracture in radians
+            a: Half-length of the fracture.
+            height: Height of the domain.
+            length: Width of the domain.
+
+        Returns:
+            A 2x2 array where each column represents the coordinates of an end point of
+            the fracture in 2D. The first column corresponds to one end point, and the
+            second column corresponds to the other.
+
+        """
+        # Rotate the fracture with an angle theta_rad
+        y_0 = height / 2 - a * np.cos(theta_rad)
+        x_0 = length / 2 - a * np.sin(theta_rad)
+        y_1 = height / 2 + a * np.cos(theta_rad)
+        x_1 = length / 2 + a * np.sin(theta_rad)
 
-        self._fractures = [pp.LineFracture(self.params["frac_pts"])]
+        frac_pts = np.array([[x_0, y_0], [x_1, y_1]]).T
+        return frac_pts
 
 
 class ManuSneddonBoundaryConditions(pp.PorePyModel):
@@ -311,9 +345,7 @@ def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
         # Set the type of west and east boundaries to Dirichlet. North and south are
         # Neumann by default.
         bc = pp.BoundaryConditionVectorial(sd, bounds.all_bf, "dir")
-        frac_face = sd.tags["fracture_faces"]
-        bc.is_neu[:, frac_face] = False
-        bc.is_dir[:, frac_face] = True
+        bc.internal_to_dirichlet()
 
         return bc
 
diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 03ad81f15e..a53ab28108 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -1,50 +1,12 @@
 import copy
 import math
 
-import numpy as np
 import pytest
 
 import porepy as pp
 from porepy.applications.convergence_analysis import ConvergenceAnalysis
 from tests.functional.setups.manu_sneddon_2d import ManuSneddonSetup2d
 
-# Set up the material constants
-poi = 0.25
-shear_modulus = 1
-lam = (
-    2 * shear_modulus * poi / (1 - 2 * poi)
-)  # Convertion formula from shear modulus and poission to lame lambda parameter
-
-solid = pp.SolidConstants(shear_modulus=shear_modulus, lame_lambda=lam)
-
-
-def compute_frac_pts(
-    theta_rad: float, a: float, height: float, length: float
-) -> np.ndarray:
-    """Assuming the fracture center is at the coordinate (height/2, length/2),
-    compute the endpoints of a fracture given its orientation and fracture length.
-
-    Parameters:
-        theta_rad: Angle of the fracture in radians
-        a: Half-length of the fracture.
-        height: Height of the domain.
-        length: Width of the domain.
-
-    Returns:
-        A 2x2 array where each column represents the coordinates of an end point of the
-        fracture in 2D. The first column corresponds to one end point, and the second
-        column corresponds to the other.
-
-    """
-    # Rotate the fracture with an angle theta_rad
-    y_0 = height / 2 - a * np.cos(theta_rad)
-    x_0 = length / 2 - a * np.sin(theta_rad)
-    y_1 = height / 2 + a * np.cos(theta_rad)
-    x_1 = length / 2 + a * np.sin(theta_rad)
-
-    frac_pts = np.array([[x_0, y_0], [x_1, y_1]]).T
-    return frac_pts
-
 
 @pytest.fixture(scope="module")
 def actual_ooc() -> dict:
@@ -65,8 +27,16 @@ def actual_ooc() -> dict:
     theta_deg = 30.0
     a = 0.3
     height = 1.0
-    length = 1.0
     theta_rad = math.radians(90 - theta_deg)
+    
+    # Set up the material constants
+    poi = 0.25
+    shear_modulus = 1
+    lam = (
+        2 * shear_modulus * poi / (1 - 2 * poi)
+    )  # Convertion formula from shear modulus and poission to lame lambda parameter
+
+    solid = pp.SolidConstants(shear_modulus=shear_modulus, lame_lambda=lam)
 
     params = {
         "prepare_simulation": True,
@@ -77,14 +47,9 @@ def actual_ooc() -> dict:
         "poi": poi,  # Possion ratio (Not standard in solid constants)
         "meshing_arguments": {"cell_size": 0.03},
         "grid_type": "simplex",
-        "theta": theta_rad,
+        "theta_rad": theta_rad,
     }
 
-    # Convert angle to radians and compute fracture points
-    params["frac_pts"] = compute_frac_pts(
-        theta_rad=theta_rad, a=a, height=height, length=length
-    )
-
     # Convergence analysis setup
     conv_analysis = ConvergenceAnalysis(
         model_class=ManuSneddonSetup2d,

From e3d94cb2c1d20a9204c591d34e868aa1497e416e Mon Sep 17 00:00:00 2001
From: Yury Zabegaev <yuri.zabegaev@uib.no>
Date: Fri, 14 Feb 2025 10:13:20 +0000
Subject: [PATCH 39/43] MAINT: Suggestiongs from the review

---
 tests/functional/test_sneddon_2d.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index a53ab28108..3b9c2ac329 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -48,6 +48,7 @@ def actual_ooc() -> dict:
         "meshing_arguments": {"cell_size": 0.03},
         "grid_type": "simplex",
         "theta_rad": theta_rad,
+        "num_bem_segments": 1000,
     }
 
     # Convergence analysis setup

From 325b6afa811ee2f359534e733333a3dcd1e405f0 Mon Sep 17 00:00:00 2001
From: Veljko Lipovac <lipovac.v@outlook.com>
Date: Thu, 20 Feb 2025 14:09:55 +0100
Subject: [PATCH 40/43] MIN: Some suggestions from review and fixing test
 model.

---
 tests/functional/setups/manu_sneddon_2d.py | 2 +-
 tests/functional/test_sneddon_2d.py        | 5 +----
 2 files changed, 2 insertions(+), 5 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index ec0cc81602..d64205005d 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -345,7 +345,7 @@ def bc_type_mechanics(self, sd: pp.Grid) -> pp.BoundaryConditionVectorial:
         # Set the type of west and east boundaries to Dirichlet. North and south are
         # Neumann by default.
         bc = pp.BoundaryConditionVectorial(sd, bounds.all_bf, "dir")
-        bc.internal_to_dirichlet()
+        bc.internal_to_dirichlet(sd)
 
         return bc
 
diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 3b9c2ac329..7c89b0f399 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -57,13 +57,10 @@ def actual_ooc() -> dict:
         model_params=copy.deepcopy(params),
         levels=2,
         spatial_refinement_rate=2,
-        temporal_refinement_rate=1,
     )
 
     # Calculate and return the order of convergence for the displacement
-    order_dict = conv_analysis.order_of_convergence(
-        conv_analysis.run_analysis(), data_range=slice(None, None, None)
-    )
+    order_dict = conv_analysis.order_of_convergence(conv_analysis.run_analysis())
     return order_dict
 
 

From 9ea17b17b8745179ff971e9df6b91b2693b44df2 Mon Sep 17 00:00:00 2001
From: Veljko Lipovac <lipovac.v@outlook.com>
Date: Thu, 20 Feb 2025 15:13:43 +0100
Subject: [PATCH 41/43] TST: Refactoring Sneddon test setup methods for BEM.

---
 tests/functional/setups/manu_sneddon_2d.py | 377 ++++++++++-----------
 tests/functional/test_sneddon_2d.py        |   2 +-
 2 files changed, 173 insertions(+), 206 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index d64205005d..342f6dda8c 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -22,200 +22,6 @@
 from porepy.geometry.distances import point_pointset
 
 
-def compute_eta(pointset_centers: np.ndarray, center: np.ndarray) -> np.ndarray:
-    """Compute the distance from fracture points to the fracture centre.
-
-    Parameter:
-        pointset_centers: Coordinates of fracture points.
-        center: Coordinates of the fracture center.
-
-    Return:
-        Array of distances each point in pointset to the center.
-
-    """
-    return point_pointset(pointset_centers, center)
-
-
-def get_bem_centers(
-    a: float, h: float, n: int, alpha: float, center: np.ndarray
-) -> np.ndarray:
-    """Compute coordinates of the centers of the BEM segments.
-
-    Parameters:
-        a: Half of the fracture length.
-        h: BEM segment length.
-        n: Number of BEM segments.
-        alpha: Orientation of the fracture.
-        center: Coordinates of the fracture center.
-
-    Return:
-        Array containing the BEM segment centers.
-
-    """
-
-    # See [1], p. 57 coordinate system (4.5.1)
-    bem_centers = np.zeros((3, n))
-    x_0 = center[0] - (a - 0.5 * h) * np.sin(alpha)
-    y_0 = center[1] - (a - 0.5 * h) * np.cos(alpha)
-    i = np.arange(n)
-    bem_centers[0] = x_0 + i * h * np.sin(alpha)
-    bem_centers[1] = y_0 + i * h * np.cos(alpha)
-
-    return bem_centers
-
-
-def analytical_displacements(
-    a: float, eta: np.ndarray, p0: float, G: float, poi: float
-) -> np.ndarray:
-    """Compute Sneddon's analytical solution for the pressurized fracture displacement.
-
-    Parameters:
-        a: Half of the fracture length.
-        eta: Distances of fracture points to the fracture center.
-        p0: Constant, fixed pressure.
-        G: Shear modulus.
-        poi: Poisson ratio.
-
-    Return
-        Array containing the analytical normal displacement jumps.
-
-    """
-    # See [2], p. 425 equ. 92
-    cons = (1 - poi) / G * p0 * a * 2
-    return cons * np.sqrt(1 - np.power(eta / a, 2))
-
-
-def transform(xc: np.ndarray, x: np.ndarray, alpha: float) -> np.ndarray:
-    """Translation and rotation transform of boundary face coordinates for the BEM.
-
-    Parameters:
-        xc: Coordinates of BEM segment centre.
-        x: Coordinates of boundary faces.
-        alpha: Fracture orientation.
-
-    Return:
-        Transformed coordinates.
-
-    """
-    x_bar = np.zeros_like(x)
-    # See [1], p. 168 equ. 7.4.6
-    x_bar[0, :] = (x[0, :] - xc[0]) * np.cos(alpha) + (x[1, :] - xc[1]) * np.sin(alpha)
-    x_bar[1, :] = -(x[0, :] - xc[0]) * np.sin(alpha) + (x[1, :] - xc[1]) * np.cos(alpha)
-    return x_bar
-
-
-def get_bc_val(
-    sd: pp.Grid,
-    bound_faces: np.ndarray,
-    xf: np.ndarray,
-    h: float,
-    poi: float,
-    alpha: float,
-    du: float,
-) -> np.ndarray:
-    """Computes semi-analytical displacement values on the boundary using the BEM for
-    the Sneddon problem.
-
-    Parameter
-        sd: The matrix grid.
-        bound_faces: Array of indices of boundary faces of ``sd``.
-        xf: Coordinates of boundary faces.
-        h: BEM segment length.
-        poi: Poisson ratio.
-        alpha: Fracture orientation.
-        du: Sneddon's analytical relative normal displacement.
-
-    Return:
-        Boundary values for the displacement.
-
-    """
-    # See [1], pages 57, 84-92, 168
-    f2 = np.zeros(bound_faces.size)
-    f3 = np.zeros(bound_faces.size)
-    f4 = np.zeros(bound_faces.size)
-    f5 = np.zeros(bound_faces.size)
-
-    u = np.zeros((sd.dim, sd.num_faces))
-
-    # See [1], equ. (7.4.5)
-    m = 1 / (4 * np.pi * (1 - poi))
-    f2[:] = m * (
-        np.log(np.sqrt((xf[0, :] - h) ** 2 + xf[1] ** 2))
-        - np.log(np.sqrt((xf[0, :] + h) ** 2 + xf[1] ** 2))
-    )
-    f3[:] = -m * (
-        np.arctan2(xf[1, :], (xf[0, :] - h)) - np.arctan2(xf[1, :], (xf[0, :] + h))
-    )
-
-    # See [1], equations 5.5.3, 5.5.1, and p. 57, equ. 4.5.1
-    f4[:] = m * (
-        xf[1, :] / ((xf[0, :] - h) ** 2 + xf[1, :] ** 2)
-        - xf[1, :] / ((xf[0, :] + h) ** 2 + xf[1, :] ** 2)
-    )
-
-    f5[:] = m * (
-        (xf[0, :] - h) / ((xf[0, :] - h) ** 2 + xf[1, :] ** 2)
-        - (xf[0, :] + h) / ((xf[0, :] + h) ** 2 + xf[1, :] ** 2)
-    )
-
-    u[0, bound_faces] = du * (
-        -(1 - 2 * poi) * np.cos(alpha) * f2[:]
-        - 2 * (1 - poi) * np.sin(alpha) * f3[:]
-        - xf[1, :] * (np.cos(alpha) * f4[:] + np.sin(alpha) * f5[:])
-    )
-    u[1, bound_faces] = du * (
-        -(1 - 2 * poi) * np.sin(alpha) * f2[:]
-        + 2 * (1 - poi) * np.cos(alpha) * f3[:]
-        - xf[1, :] * (np.sin(alpha) * f4[:] - np.cos(alpha) * f5[:])
-    )
-
-    return u
-
-
-def assign_bem(
-    grid: pp.Grid,
-    h: float,
-    bound_faces: np.ndarray,
-    alpha: float,
-    bem_centers: np.ndarray,
-    u_a: np.ndarray,
-    poi: float,
-) -> np.ndarray:
-    """Computes semi-analytical displacement values using the BEM for the Sneddon
-    problem.
-
-    Parameter
-        grid: The matrix grid.
-        h: BEM segment length.
-        bound_faces: Array of indices of boundary faces of ``grid``.
-        alpha: Fracture orientation.
-        bem_centers: BEM segments centers.
-        u_a: Sneddon's analytical relative normal displacement.
-        poi: Poisson ratio.
-
-    Return:
-        Semi-analytical boundary displacement values.
-
-    """
-
-    bc_val = np.zeros((grid.dim, grid.num_faces))
-
-    alpha = np.pi / 2 - alpha
-
-    bound_face_centers = grid.face_centers[:, bound_faces]
-
-    for i in range(0, u_a.size):
-        new_bound_face_centers = transform(bem_centers[:, i], bound_face_centers, alpha)
-
-        u_bound = get_bc_val(
-            grid, bound_faces, new_bound_face_centers, h, poi, alpha, u_a[i]
-        )
-
-        bc_val += u_bound
-
-    return bc_val
-
-
 class SneddonExactSolution2d:
     """Class representing the analytical solution for the pressurized fracture problem."""
 
@@ -243,14 +49,12 @@ def exact_sol_global(self, sd: pp.Grid) -> np.ndarray:
         u_bc = np.zeros((sd.dim, sd.num_faces))
 
         center = np.array([self.length / 2, self.length / 2, 0])
-        bem_centers = get_bem_centers(self.a, h, self.n, self.theta_rad, center)
-        eta = compute_eta(bem_centers, center)
+        bem_centers = self.get_bem_centers(h, center)
+        eta = point_pointset(bem_centers, center)
 
-        u_a = -analytical_displacements(
-            self.a, eta, self.p0, self.shear_modulus, self.poi
-        )
-        u_bc = assign_bem(
-            sd, h / 2, box_faces, self.theta_rad, bem_centers, u_a, self.poi
+        u_a = -self.analytical_displacements(eta)
+        u_bc = self.displacement_bc_values_from_bem(
+            sd, h / 2, box_faces, bem_centers, u_a
         )
         return u_bc
 
@@ -276,13 +80,176 @@ def exact_sol_fracture(
         fracture_faces = g_1.cell_centers
 
         # See [2], p. 425, equ. 92
-        eta = compute_eta(fracture_faces, fracture_center)
-        apertures = analytical_displacements(
-            self.a, eta, self.p0, self.shear_modulus, self.poi
-        )
+        eta = point_pointset(fracture_faces, fracture_center)
+        apertures = self.analytical_displacements(eta)
 
         return eta, apertures
 
+    def get_bem_centers(self, h: float, center: np.ndarray) -> np.ndarray:
+        """Compute coordinates of the centers of the BEM segments.
+
+        Parameters:
+            h: BEM segment length.
+            center: Coordinates of the fracture center.
+
+        Return:
+            Array containing the BEM segment centers.
+
+        """
+
+        # See [1], p. 57 coordinate system (4.5.1)
+        bem_centers = np.zeros((3, self.n))
+        x_0 = center[0] - (self.a - 0.5 * h) * np.sin(self.theta_rad)
+        y_0 = center[1] - (self.a - 0.5 * h) * np.cos(self.theta_rad)
+        i = np.arange(self.n)
+        bem_centers[0] = x_0 + i * h * np.sin(self.theta_rad)
+        bem_centers[1] = y_0 + i * h * np.cos(self.theta_rad)
+
+        return bem_centers
+
+    def analytical_displacements(self, eta: np.ndarray) -> np.ndarray:
+        """Compute Sneddon's analytical solution for the pressurized fracture displacement.
+
+        Parameters:
+            eta: Distances of fracture points to the fracture center.
+
+        Return
+            Array containing the analytical normal displacement jumps.
+
+        """
+        # See [2], p. 425 equ. 92
+        cons = (1 - self.poi) / self.shear_modulus * self.p0 * self.a * 2
+        return cons * np.sqrt(1 - np.power(eta / self.a, 2))
+
+    def tranform_boundary_face_coordinates_bem(
+        self, xc: np.ndarray, x: np.ndarray
+    ) -> np.ndarray:
+        """Translation and rotation transform of boundary face coordinates for the BEM.
+
+        Parameters:
+            xc: Coordinates of BEM segment centre.
+            x: Coordinates of boundary faces.
+
+        Return:
+            Transformed coordinates.
+
+        """
+        alpha = np.pi / 2 - self.theta_rad
+        x_bar = np.zeros_like(x)
+        # See [1], p. 168 equ. 7.4.6
+        x_bar[0, :] = (x[0, :] - xc[0]) * np.cos(alpha) + (x[1, :] - xc[1]) * np.sin(
+            alpha
+        )
+        x_bar[1, :] = -(x[0, :] - xc[0]) * np.sin(alpha) + (x[1, :] - xc[1]) * np.cos(
+            alpha
+        )
+        return x_bar
+
+    def _get_bc_val_from_normal_displacement(
+        self,
+        sd: pp.Grid,
+        bound_faces: np.ndarray,
+        xf: np.ndarray,
+        h: float,
+        du: float,
+    ) -> np.ndarray:
+        """Computes semi-analytical displacement values on the boundary using the BEM for
+        the Sneddon problem.
+
+        Parameter
+            sd: The matrix grid.
+            bound_faces: Array of indices of boundary faces of ``sd``.
+            xf: Coordinates of boundary faces.
+            h: BEM segment length.
+            du: Sneddon's analytical relative normal displacement.
+
+        Return:
+            Boundary values for the displacement.
+
+        """
+        # See [1], pages 57, 84-92, 168
+        f2 = np.zeros(bound_faces.size)
+        f3 = np.zeros(bound_faces.size)
+        f4 = np.zeros(bound_faces.size)
+        f5 = np.zeros(bound_faces.size)
+
+        u = np.zeros((sd.dim, sd.num_faces))
+
+        alpha = np.pi / 2 - self.theta_rad
+
+        # See [1], equ. (7.4.5)
+        m = 1 / (4 * np.pi * (1 - self.poi))
+        f2[:] = m * (
+            np.log(np.sqrt((xf[0, :] - h) ** 2 + xf[1] ** 2))
+            - np.log(np.sqrt((xf[0, :] + h) ** 2 + xf[1] ** 2))
+        )
+        f3[:] = -m * (
+            np.arctan2(xf[1, :], (xf[0, :] - h)) - np.arctan2(xf[1, :], (xf[0, :] + h))
+        )
+
+        # See [1], equations 5.5.3, 5.5.1, and p. 57, equ. 4.5.1
+        f4[:] = m * (
+            xf[1, :] / ((xf[0, :] - h) ** 2 + xf[1, :] ** 2)
+            - xf[1, :] / ((xf[0, :] + h) ** 2 + xf[1, :] ** 2)
+        )
+
+        f5[:] = m * (
+            (xf[0, :] - h) / ((xf[0, :] - h) ** 2 + xf[1, :] ** 2)
+            - (xf[0, :] + h) / ((xf[0, :] + h) ** 2 + xf[1, :] ** 2)
+        )
+
+        u[0, bound_faces] = du * (
+            -(1 - 2 * self.poi) * np.cos(alpha) * f2[:]
+            - 2 * (1 - self.poi) * np.sin(alpha) * f3[:]
+            - xf[1, :] * (np.cos(alpha) * f4[:] + np.sin(alpha) * f5[:])
+        )
+        u[1, bound_faces] = du * (
+            -(1 - 2 * self.poi) * np.sin(alpha) * f2[:]
+            + 2 * (1 - self.poi) * np.cos(alpha) * f3[:]
+            - xf[1, :] * (np.sin(alpha) * f4[:] - np.cos(alpha) * f5[:])
+        )
+
+        return u
+
+    def displacement_bc_values_from_bem(
+        self,
+        grid: pp.Grid,
+        h: float,
+        bound_faces: np.ndarray,
+        bem_centers: np.ndarray,
+        u_a: np.ndarray,
+    ) -> np.ndarray:
+        """Computes semi-analytical displacement values using the BEM for the Sneddon
+        problem.
+
+        Parameter
+            grid: The matrix grid.
+            h: BEM segment length.
+            bound_faces: Array of indices of boundary faces of ``grid``.
+            bem_centers: BEM segments centers.
+            u_a: Sneddon's analytical relative normal displacement.
+
+        Return:
+            Semi-analytical boundary displacement values.
+
+        """
+
+        bc_val = np.zeros((grid.dim, grid.num_faces))
+        bound_face_centers = grid.face_centers[:, bound_faces]
+
+        for i in range(0, u_a.size):
+            new_bound_face_centers = self.tranform_boundary_face_coordinates_bem(
+                bem_centers[:, i], bound_face_centers
+            )
+
+            u_bound = self._get_bc_val_from_normal_displacement(
+                grid, bound_faces, new_bound_face_centers, h, u_a[i]
+            )
+
+            bc_val += u_bound
+
+        return bc_val
+
 
 class ManuSneddonGeometry2d(SquareDomainOrthogonalFractures):
     """Square domain but with single line fracture."""
diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index 7c89b0f399..d989f0e198 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -28,7 +28,7 @@ def actual_ooc() -> dict:
     a = 0.3
     height = 1.0
     theta_rad = math.radians(90 - theta_deg)
-    
+
     # Set up the material constants
     poi = 0.25
     shear_modulus = 1

From af71f2a5e7c6c63590a8181aabcf7c5546b0bfa2 Mon Sep 17 00:00:00 2001
From: Yury Zabegaev <yuri.zabegaev@uib.no>
Date: Thu, 20 Feb 2025 14:52:23 +0000
Subject: [PATCH 42/43] TEST: Truncation of fracture tips

---
 tests/functional/setups/manu_sneddon_2d.py | 58 +++++++++++++++++-----
 tests/functional/test_sneddon_2d.py        | 11 +++-
 2 files changed, 55 insertions(+), 14 deletions(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index 342f6dda8c..fdda0780d6 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -58,9 +58,7 @@ def exact_sol_global(self, sd: pp.Grid) -> np.ndarray:
         )
         return u_bc
 
-    def exact_sol_fracture(
-        self, mdg: pp.MixedDimensionalGrid
-    ) -> tuple[np.ndarray, np.ndarray]:
+    def exact_sol_fracture(self, mdg: pp.MixedDimensionalGrid) -> np.ndarray:
         """Compute Sneddon's analytical solution for the pressurized crack
         problem in question.
 
@@ -68,9 +66,7 @@ def exact_sol_fracture(
             mdg: Mixed-dimensional domain of the setup.
 
         Return:
-            A 2-tuple of numpy arrays, where the first vector contains the distances
-            between fracture center to each fracture coordinate, and the second vector
-            the respective analytical apertures.
+            A numpy array with the analytical apertures.
 
         """
 
@@ -83,7 +79,7 @@ def exact_sol_fracture(
         eta = point_pointset(fracture_faces, fracture_center)
         apertures = self.analytical_displacements(eta)
 
-        return eta, apertures
+        return apertures
 
     def get_bem_centers(self, h: float, center: np.ndarray) -> np.ndarray:
         """Compute coordinates of the centers of the BEM segments.
@@ -359,13 +355,21 @@ def collect_data(self) -> ManuSneddonSaveData:
         frac_sd = self.mdg.subdomains(dim=self.nd - 1)
         nd_vec_to_normal = self.normal_component(frac_sd)
 
-        # Computing the numerical displacement jump along the fracture on the fracture
+        # Computing the numerical displacement jump along the fracture on the fracture.
         # cell centers.
-        u_n: pp.ad.Operator = nd_vec_to_normal @ self.displacement_jump(frac_sd)
-        u_n = u_n.value(self.equation_system)
+        u_n = (nd_vec_to_normal @ self.displacement_jump(frac_sd)).value(
+            self.equation_system
+        )
+
+        # Checking convergence specifically on the fracture.
+        u_a = self.exact_sol.exact_sol_fracture(self.mdg)
 
-        # Checking convergence specifically on the fracture
-        u_a = self.exact_sol.exact_sol_fracture(self.mdg)[1]
+        # Truncate the solution at the fracture cells close to the boundary.
+        # How close is determined by the parameter eps (percent).
+        eps = self.params["error_exclusion_zone_fracture_tips"]
+        close_to_fracture_tips = self._index_fracture_cells_close_to_fracture_tips(eps)
+        u_a[close_to_fracture_tips] = 0
+        u_n[close_to_fracture_tips] = 0
 
         e = ConvergenceAnalysis.lp_error(
             frac_sd[0], u_a, u_n, is_scalar=False, is_cc=True, relative=True
@@ -374,6 +378,36 @@ def collect_data(self) -> ManuSneddonSaveData:
         collect_data = ManuSneddonSaveData(error_displacement=e)
         return collect_data
 
+    def _index_fracture_cells_close_to_fracture_tips(self, eps: float) -> np.ndarray:
+        """Compute the boolean mask array with `True` values for the fracture cells
+        close to the fracture tips.
+
+        Parameters:
+            eps: The relative threshold to determine how close the fracture cells
+                are to the fracture tips to be truncated, must be in `[0, 1]`.
+
+        Returns:
+            A boolean array of shape `(num_fracture_cells,)`.
+
+        """
+        frac_sd = self.mdg.subdomains(dim=self.nd - 1)
+        assert len(frac_sd) == 1, "Must be only one fracture."
+
+        # Find the domain center.
+        xmax, xmin = self.domain.bounding_box["xmax"], self.domain.bounding_box["xmin"]
+        ymax, ymin = self.domain.bounding_box["ymax"], self.domain.bounding_box["ymin"]
+        xmean = (xmax - xmin) * 0.5
+        ymean = (ymax - ymin) * 0.5
+
+        # Compute distance from the domain center.
+        x, y, _ = frac_sd[0].cell_centers
+        distance_from_center = np.sqrt((x - xmean) ** 2 + (y - ymean) ** 2)
+
+        # Compute the percentage of how close the cells are to the fracture tips.
+        a = self.exact_sol.a  # Fracture half-length.
+        relative_distance_from_center = (a - distance_from_center) / a
+        return relative_distance_from_center < eps
+
 
 class ManuSneddonConstitutiveLaws(pp.constitutive_laws.PressureStress):
     """Inherits the relevant pressure-stress relations, defines a constant pressure
diff --git a/tests/functional/test_sneddon_2d.py b/tests/functional/test_sneddon_2d.py
index d989f0e198..b85661a618 100644
--- a/tests/functional/test_sneddon_2d.py
+++ b/tests/functional/test_sneddon_2d.py
@@ -1,6 +1,7 @@
 import copy
 import math
 
+import numpy as np
 import pytest
 
 import porepy as pp
@@ -49,6 +50,8 @@ def actual_ooc() -> dict:
         "grid_type": "simplex",
         "theta_rad": theta_rad,
         "num_bem_segments": 1000,
+        # Truncate results from the cells closer than 10% to the fracture tips.
+        "error_exclusion_zone_fracture_tips": 0.1,
     }
 
     # Convergence analysis setup
@@ -66,5 +69,9 @@ def actual_ooc() -> dict:
 
 def test_order_of_convergence(actual_ooc: dict) -> None:
     """Test observed order of convergence."""
-    # We  the order of L2 convergence on the fracture of displacement to be about 1.0
-    assert 0.85 < actual_ooc["ooc_displacement"]
+    # The `error_exclusion_zone_fracture_tips`` is set to 10% to balance the observation
+    # of boundary effects and good convergence. This test evaluates changes in fracture
+    # tip treatment and the inner domain. Raising the threshold to 15% yielded a
+    # convergence order of ~2. Decreasing the threshold to 0% yielded a convergence
+    # order of ~0.85.
+    assert np.isclose(1.66752, actual_ooc["ooc_displacement"])

From 2438952cf3877d0620f9a557b8bc2de611d5eacf Mon Sep 17 00:00:00 2001
From: Veljko Lipovac <68282540+vlipovac@users.noreply.github.com>
Date: Fri, 21 Feb 2025 10:28:48 +0100
Subject: [PATCH 43/43] Update tests/functional/setups/manu_sneddon_2d.py

Co-authored-by: Ivar Stefansson <ivar.stefansson@uib.no>
---
 tests/functional/setups/manu_sneddon_2d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/functional/setups/manu_sneddon_2d.py b/tests/functional/setups/manu_sneddon_2d.py
index fdda0780d6..eac7eacf00 100644
--- a/tests/functional/setups/manu_sneddon_2d.py
+++ b/tests/functional/setups/manu_sneddon_2d.py
@@ -248,7 +248,7 @@ def displacement_bc_values_from_bem(
 
 
 class ManuSneddonGeometry2d(SquareDomainOrthogonalFractures):
-    """Square domain but with single line fracture."""
+    """Square domain  with a single, arbitrarily oriented fracture at the center."""
 
     def set_fractures(self):
         """Setting a single line fracture which runs through the domain center."""