diff --git a/discretize/curvilinear_mesh.py b/discretize/curvilinear_mesh.py index 1c6181a42..d2336fc5b 100644 --- a/discretize/curvilinear_mesh.py +++ b/discretize/curvilinear_mesh.py @@ -13,6 +13,7 @@ from discretize.base import BaseRectangularMesh from discretize.operators import DiffOperators, InnerProducts from discretize.mixins import InterfaceMixins +from discretize.utils import spzeros # Some helper functions. @@ -33,7 +34,10 @@ def _normalize3D(x): class CurvilinearMesh( - DiffOperators, InnerProducts, BaseRectangularMesh, InterfaceMixins + DiffOperators, + InnerProducts, + BaseRectangularMesh, + InterfaceMixins, ): """Curvilinear mesh class. @@ -797,10 +801,13 @@ def _get_edge_surf_int_proj_mats(self, only_boundary=False, with_area=True): edge_dirs = self.edge_tangents[node_edges] t_for = np.concatenate((edge_dirs, face_normals[:, None, :]), axis=1) t_inv = np.linalg.inv(t_for) - t_inv = t_inv[:, :, :-1] / 4 # n_edges_per_thing + t_inv = t_inv[:, :, :-1] if with_area: t_inv *= face_areas[:, None, None] + t_inv /= 4 # n_edges_per_thing + else: + t_inv /= 2 # sqrt n_edges_per_thing T = C2F @ sp.csr_matrix( (t_inv.reshape(-1), T_col_inds, T_ind_ptr), @@ -809,6 +816,132 @@ def _get_edge_surf_int_proj_mats(self, only_boundary=False, with_area=True): Ps.append((T @ P)) return Ps + def get_edge_inner_product_surface( # NOQA D102 + self, model=None, invert_model=False, invert_matrix=False + ): + # Documentation inherited from discretize.base.BaseMesh + dim = self.dim + if dim == 2: + # in 2D faces_x -> edges_y and edges_x -> faces_y, so need to permute the x and y faces to x and y edges. + P_e2f = sp.diags( + [1, 1], + (-self.n_edges_y, self.n_edges_x), + shape=(self.n_edges, self.n_edges), + format="csr", + ) + return ( + P_e2f.T + @ super().get_face_inner_product_surface( + model=model, invert_model=invert_model, invert_matrix=invert_matrix + ) + @ P_e2f + ) + + if invert_matrix: + raise NotImplementedError( + "The inverse of the inner product matrix with a tetrahedral mesh is not supported." + ) + + # Edge inner product surface projection matrices + n_faces = self.n_faces + face_areas = self.face_areas + + Ps = self._get_edge_surf_int_proj_mats(with_area=False) + + if model is None: + Mu = sp.diags(np.tile(face_areas, dim)) # Number of edges per face + else: + if invert_model: + model = 1.0 / model + + if (model.size == 1) | (model.size == n_faces): + Mu = sp.diags( + np.tile(model * face_areas, dim) + ) # Number of edges per face + else: + raise ValueError( + "Unrecognized size of model vector.", + "Must be scalar or have length equal to total number of faces.", + ) + + A = np.sum([P.T @ Mu @ P for P in Ps]) + + return A + + def get_edge_inner_product_surface_deriv( # NOQA D102 + self, + model, + invert_model=False, + invert_matrix=False, + ): + # Documentation inherited from discretize.base.BaseMesh + dim = self.dim + if dim == 2: + # in 2D faces_x -> edges_y and edges_x -> faces_y, so need to permute the x and y faces to x and y edges. + P_e2f = sp.diags( + [1, 1], + (-self.n_edges_y, self.n_edges_x), + shape=(self.n_edges, self.n_edges), + format="csr", + ) + face_func = self.get_face_inner_product_surface_deriv( + model=model, invert_model=invert_model, invert_matrix=invert_matrix + ) + + def func(v): + return P_e2f.T @ (face_func(P_e2f @ v)) + + return func + + if invert_model: + raise NotImplementedError( + "Inverted model derivatives are not supported here" + ) + if invert_matrix: + raise NotImplementedError( + "The inverse of the inner product matrix with a tetrahedral mesh is not supported." + ) + model = np.asarray(model) + # Edge inner product surface projection matrices + n_faces = self.n_faces + n_edges = self.n_edges + face_areas = self.face_areas + + Ps = self._get_edge_surf_int_proj_mats(with_area=False) + area = sp.diags(np.tile(np.sqrt(face_areas), dim)) + Ps = list([area @ P for P in Ps]) + + if model.size == 1: + + def func(v): + dMdm = spzeros(n_edges, 1) + for P in Ps: + dMdm = dMdm + sp.csr_matrix( + (P.T @ (P @ v), (range(n_edges), np.zeros(n_edges))), + shape=(n_edges, 1), + ) + return dMdm + + elif model.size == n_faces: + col_inds = np.tile(np.arange(n_faces), dim) + ind_ptr = np.arange(n_faces * dim + 1) + + def func(v): + dMdm = spzeros(n_edges, n_faces) + for P in Ps: + ys = P @ v + dMdm = dMdm + P.T @ sp.csr_matrix( + (ys, col_inds, ind_ptr), shape=(n_faces * dim, n_faces) + ) + return dMdm + + else: + raise ValueError( + "Unrecognized size of model vector.", + "Must be scalar or have length equal to total number of faces.", + ) + return func + @property def boundary_edge_vector_integral(self): # NOQA D102 # Documentation inherited from discretize.base.BaseMesh diff --git a/discretize/operators/__init__.py b/discretize/operators/__init__.py index 70c324760..b371f10c1 100644 --- a/discretize/operators/__init__.py +++ b/discretize/operators/__init__.py @@ -4,7 +4,7 @@ ================================================ .. currentmodule:: discretize.operators -The ``operators`` package contains the classes discretize meshes with regular structure +The ``operators`` package contains the classes discretize meshes use to construct discrete versions of the differential operators. Operator Classes diff --git a/discretize/unstructured_mesh.py b/discretize/unstructured_mesh.py index 408442cd0..faec51c10 100644 --- a/discretize/unstructured_mesh.py +++ b/discretize/unstructured_mesh.py @@ -520,6 +520,47 @@ def get_edge_inner_product( # NOQA D102 ) return self.__get_inner_product("E", model, invert_model) + def get_edge_inner_product_surface( # NOQA D102 + self, model=None, invert_model=False, invert_matrix=False + ): + # Documentation inherited from discretize.base.BaseMesh + dim = self.dim + if dim == 2: + return self.get_face_inner_product_surface( + model=model, invert_model=invert_model, invert_matrix=invert_matrix + ) + + if invert_matrix: + raise NotImplementedError( + "The inverse of the inner product matrix with a tetrahedral mesh is not supported." + ) + + # Edge inner product surface projection matrices + n_faces = self.n_faces + face_areas = self.face_areas + + Ps = self._get_edge_surf_int_proj_mats(with_area=False) + + if model is None: + Mu = sp.diags(np.tile(face_areas, dim)) # Number of edges per face + else: + if invert_model: + model = 1.0 / model + + if (model.size == 1) | (model.size == n_faces): + Mu = sp.diags( + np.tile(model * face_areas, dim) + ) # Number of edges per face + else: + raise ValueError( + "Unrecognized size of model vector.", + "Must be scalar or have length equal to total number of faces.", + ) + + A = np.sum([P.T @ Mu @ P for P in Ps]) + + return A + def __get_inner_product_deriv_func(self, i_type, model): Ps, _ = self.__get_inner_product_projection_matrices(i_type) dim = self.dim @@ -611,6 +652,68 @@ def get_edge_inner_product_deriv( # NOQA D102 raise NotImplementedError("Inverted matrix derivatives are not supported") return self.__get_inner_product_deriv_func("E", model) + def get_edge_inner_product_surface_deriv( # NOQA D102 + self, + model, + invert_model=False, + invert_matrix=False, + ): + # Documentation inherited from discretize.base.BaseMesh + dim = self.dim + if dim == 2: + return super().get_face_inner_product_surface_deriv( + model=model, invert_model=invert_model, invert_matrix=invert_matrix + ) + + if invert_model: + raise NotImplementedError( + "Inverted model derivatives are not supported here" + ) + if invert_matrix: + raise NotImplementedError( + "The inverse of the inner product matrix with a tetrahedral mesh is not supported." + ) + model = np.asarray(model) + # Edge inner product surface projection matrices + n_faces = self.n_faces + n_edges = self.n_edges + face_areas = self.face_areas + + Ps = self._get_edge_surf_int_proj_mats(with_area=False) + area = sp.diags(np.tile(np.sqrt(face_areas), dim)) + Ps = list([area @ P for P in Ps]) + + if model.size == 1: + + def func(v): + dMdm = spzeros(n_edges, 1) + for P in Ps: + dMdm = dMdm + sp.csr_matrix( + (P.T @ (P @ v), (range(n_edges), np.zeros(n_edges))), + shape=(n_edges, 1), + ) + return dMdm + + elif model.size == n_faces: + col_inds = np.tile(np.arange(n_faces), dim) + ind_ptr = np.arange(n_faces * dim + 1) + + def func(v): + dMdm = spzeros(n_edges, n_faces) + for P in Ps: + ys = P @ v + dMdm = dMdm + P.T @ sp.csr_matrix( + (ys, col_inds, ind_ptr), shape=(n_faces * dim, n_faces) + ) + return dMdm + + else: + raise ValueError( + "Unrecognized size of model vector.", + "Must be scalar or have length equal to total number of faces.", + ) + return func + def _get_edge_surf_int_proj_mats(self, only_boundary=False, with_area=True): """Return the projection operators for integrating edges on each face. @@ -668,10 +771,13 @@ def _get_edge_surf_int_proj_mats(self, only_boundary=False, with_area=True): edge_dirs = self.edge_tangents[node_edges] t_for = np.concatenate((edge_dirs, face_normals[:, None, :]), axis=1) t_inv = invert_blocks(t_for) - t_inv = t_inv[:, :, :-1] / 3 # n_edges_per_thing + t_inv = t_inv[:, :, :-1] if with_area: t_inv *= face_areas[:, None, None] + t_inv /= 3 # n_edges_per_thing + else: + t_inv /= np.sqrt(3) # sqrt n_edges_per_thing T = C2F @ sp.csr_matrix( (t_inv.reshape(-1), T_col_inds, T_ind_ptr), diff --git a/tests/base/test_curvilinear.py b/tests/base/test_curvilinear.py index 65cf5c9d2..c71c401cb 100644 --- a/tests/base/test_curvilinear.py +++ b/tests/base/test_curvilinear.py @@ -1,7 +1,10 @@ import numpy as np import unittest +import pytest + from discretize import TensorMesh, CurvilinearMesh from discretize.utils import ndgrid +from discretize.tests import setup_mesh, check_derivative, assert_expected_order class BasicCurvTests(unittest.TestCase): @@ -287,5 +290,222 @@ def test_grid(self): self.assertTrue(np.all(self.Curv3.gridEz == self.TM3.gridEz)) +@pytest.mark.parametrize("u_type", ["edge", "face"]) +@pytest.mark.parametrize("dim", [2, 3], ids=["2D", "3D"]) +@pytest.mark.parametrize("rep", [0, 1], ids=["uniform", "isotropic"]) +def test_surface_inner_product_prop_deriv(u_type, dim, rep): + rng = np.random.default_rng(6732) + mesh, _ = setup_mesh("rotateCurv", 20, dim) + tau = rng.uniform(1, 2, 1) if rep == 0 else rng.uniform(1, 2, mesh.n_faces * rep) + + match u_type: + case "edge": + v = rng.uniform(1, 2, mesh.n_edges) + + def fun(tau): + M = mesh.get_edge_inner_product_surface(tau) + Md = mesh.get_edge_inner_product_surface_deriv(tau) + return M * v, Md(v) + + case "face": + v = rng.uniform(1, 2, mesh.n_faces) + + def fun(tau): + M = mesh.get_face_inner_product_surface(tau) + Md = mesh.get_face_inner_product_surface_deriv(tau) + return M * v, Md(v) + + case _: + raise Exception("Invalid test parameter.") + + check_derivative(fun, tau, num=5, random_seed=rng) + + +@pytest.mark.parametrize("dim", [2, 3], ids=["2D", "3D"]) +@pytest.mark.parametrize("rep", [0, 1], ids=["uniform", "isotropic"]) +def test_line_inner_product_prop_deriv(dim, rep): + rng = np.random.default_rng(6732) + mesh, _ = setup_mesh("rotateCurv", 20, dim) + v = rng.uniform(1, 2, mesh.n_edges) + tau = rng.uniform(1, 2, 1) if rep == 0 else rng.uniform(1, 2, mesh.n_edges * rep) + + def fun(tau): + M = mesh.get_edge_inner_product_line(tau) + Md = mesh.get_edge_inner_product_line_deriv(tau) + return M * v, Md(v) + + check_derivative(fun, tau, num=5, random_seed=rng) + + +def test_edge_surface_integral_2d(): + """ + testing the line integral here of: + vector field w: [y**2, x**2] + physical property u: (1 - x) * (1 + x) + over the path x=t, y=(t - 1) * (t + 1) + from t=-1 to 1 + """ + + def error_eval(nx): + ny = 4 # not really important to this test, as only one "surface" is non-zeros + xlocs = np.linspace(-1, 1, nx + 1) + nodes_x = xlocs[:, None] * np.ones(ny + 1)[None, :] + + ylocs = (xlocs - 1) * (xlocs + 1) + nodes_y = ylocs[:, None] + np.linspace(-1, 1, ny + 1) + + mesh = CurvilinearMesh((nodes_x, nodes_y)) + + faces = (np.arange(mesh.n_faces_y) + mesh.n_faces_x).reshape(ny + 1, nx) + face_inds = faces[2, :] + + xcs = (xlocs[1:] + xlocs[:-1]) * 0.5 + prop = (1 - xcs) * (xcs + 1) + 1 + face_props = np.zeros(mesh.n_faces) + face_props[face_inds] = prop + + edge_vectors = np.c_[mesh.edges[:, 1] ** 2, mesh.edges[:, 0] ** 2] + edge_vals = mesh.project_edge_vector(edge_vectors) + + M = mesh.get_edge_inner_product_surface(model=face_props) + + discrete_val = np.sum(M @ edge_vals) + reference_value = 208 / 105 + return np.abs(discrete_val - reference_value), xlocs[1] - xlocs[0] + + assert_expected_order(error_eval, [20.0, 30.0, 40.0, 50.0]) + + +def test_edge_surface_integral_3d(): + """ + For this test, we will use a single surface within the curvilinear mesh + parameterized as: + x = u + y = v + z = u**2 + v**2 + (a simple parabola) + + we want to test the integral of: + mu * w.dot(w) ds + over this surface, for some a vector field that is always tangent to the surface + the two vectors that are perpendicular and tangent to this surface are: + r_u = [1, 0, 2*u] + r_v = [0, 1, 2*v] + + so any linear combination of these vectors will be tangent to the surface, let's use: + w = x * y * z * (r_u + r_v) + + the ds integrand is sqrt(|| r_u.cross(r_v) ||) or: + sqrt(4 * x**2 + 4 * y**2 + 1) + + and a property field: mu = x**2 + y**2 + z**2 + + the analytic integral over the range u=[-1, 1], v=[-1,1] found by numeric integration: + + >>> def int_f(x, y): + ... xsq = x**2 + ... ysq = y**2 + ... z = x**2 + y**2 + ... zsq = z**2 + ... v1 = 2 * xsq * ysq * zsq + (2 * xsq * y * z + 2 * x * ysq * z)**2 + ... v2 = (xsq + ysq + zsq) * np.sqrt(4 * xsq + 4 * ysq + 1) + ... return v1 * v2 + >>> scint.dblquad(int_f, -1, 1, -1, 1, epsrel=1e-20, epsabs=1E-20) + (51.8667132595898, 2.241102757162104e-12) + """ + + def error_eval(nx): + ny = nx + nz = 4 # not really important to this test, as only one "surface" is non-zero + xlocs = np.linspace(-1, 1, nx + 1) + ylocs = np.linspace(-1, 1, ny + 1) + nodes_x = xlocs[:, None, None] * np.ones((1, ny + 1, nz + 1)) + nodes_y = ylocs[None, :, None] * np.ones((nx + 1, 1, nz + 1)) + + nodes_z = nodes_x**2 + nodes_y**2 + np.linspace(-1, 1, nz + 1)[None, None, :] + mesh = CurvilinearMesh((nodes_x, nodes_y, nodes_z)) + + faces = (np.arange(mesh.n_faces_z) + mesh.n_faces_x + mesh.n_faces_y).reshape( + nz + 1, nx * ny + ) + face_inds = faces[2, :] + + prop = mesh.faces[:, 0] ** 2 + mesh.faces[:, 1] ** 2 + mesh.faces[:, 2] ** 2 + face_props = np.zeros(mesh.n_faces) + face_props[face_inds] = prop[face_inds] + + ex = mesh.edges[:, 0] + ey = mesh.edges[:, 1] + w_vec = (ex * ey * (ex**2 + ey**2))[:, None] * np.c_[ + np.ones_like(ex), + np.ones_like(ex), + 2 * ex + 2 * ey, + ] + w = mesh.project_edge_vector(w_vec) + + M = mesh.get_edge_inner_product_surface(model=face_props) + + discrete_val = w @ M @ w + reference_value = 51.8667132595898 + return np.abs(discrete_val - reference_value), xlocs[1] - xlocs[0] + + assert_expected_order(error_eval, [20.0, 30.0, 40.0, 50.0]) + + +def test_edge_line_integral_3d(): + """ + testing the line integral here of: + scalar field: v = x**2 + y**2 + times the property: mu = (1 - x) * (x + 1) + 1 + over the path x=t, y = 2 * t, z=(t - 1) * (t + 1) + from t=-1 to 1 + + >>> t = sy.Symbol('t') + >>> x, y, z = t, 2 * t, (t-1)*(t+1) + >>> mu = (1 - t) * (t+1) + 1 + >>> v = x**2 + y**2 + >>> dx = sy.diff(x, t) + >>> dy = sy.diff(y, t) + >>> dz = sy.diff(z, t) + >>> dl = sy.sqrt(dx **2 + dy**2 + dz**2) + >>> integrand = mu * v * dl + >>> val = sy.integrate(integrand, (t, -1, 1)) + >>> float(val) + 12.490674765352512 + + """ + + def error_eval(nx): + nz = ny = 4 # not really important to this test, as only one "line" is non-zero + xlocs = np.linspace(-1, 1, nx + 1) + nodes_x = xlocs[:, None, None] * np.ones((1, ny + 1, nz + 1)) + nodes_y = 2 * nodes_x + np.linspace(-1, 1, ny + 1)[None, :, None] + nodes_z = (nodes_x - 1) * (nodes_x + 1) + np.linspace(-1, 1, nz + 1)[ + None, None, : + ] + + mesh = CurvilinearMesh((nodes_x, nodes_y, nodes_z)) + + # edge indices along the path... are x-edges: + edge_inds = np.arange(mesh.n_edges_x).reshape((nz + 1, ny + 1, nx))[2, 2, :] + + prop = (1 - mesh.edges[:, 0]) * (mesh.edges[:, 0] + 1) + 1 + edge_props = np.zeros(mesh.n_edges) + edge_props[edge_inds] = prop[edge_inds] + + ex = mesh.edges[:, 0] + ey = mesh.edges[:, 1] + w = ex**2 + ey**2 + + M = mesh.get_edge_inner_product_line(model=edge_props) + + discrete_val = np.sum(M @ w) + reference_value = 12.490674765352512 + print(discrete_val) + return np.abs(discrete_val - reference_value), xlocs[1] - xlocs[0] + + assert_expected_order(error_eval, [20.0, 40.0, 80.0]) + + if __name__ == "__main__": unittest.main() diff --git a/tests/simplex/test_inner_products.py b/tests/simplex/test_inner_products.py index a5cae4ce6..80550afbf 100644 --- a/tests/simplex/test_inner_products.py +++ b/tests/simplex/test_inner_products.py @@ -1,5 +1,6 @@ import numpy as np import unittest +import pytest import discretize import scipy.sparse as sp from discretize.utils import example_simplex_mesh @@ -181,6 +182,149 @@ def test_order3_faces_invert_model(self): self.orderTest() +class TestInnerProductsFaceProperties2D(discretize.tests.OrderTest): + meshSizes = [8, 16, 32] + meshTypes = ["uniform simplex mesh"] + + def setupMesh(self, n): + points, simplices = example_simplex_mesh((n, n)) + self.M = discretize.SimplexMesh(points, simplices) + return 1.0 / n + + def getError(self): + + call = lambda fun, xy: fun(xy[:, 0], xy[:, 1]) # NOQA F841 + + ex = lambda x, y: x**2 + y + ey = lambda x, y: (y**2) * x + + tau_funcs = { + "x": lambda x, y: 2 * y + 1, # x-face properties + "y": lambda x, y: x + 2, # y-face properties + } + + mesh = self.M + + tau = 1e-8 * np.ones(mesh.n_faces) + for ii, comp in enumerate(["x", "y"]): + k = np.isclose(self.M.faces[:, ii], 0.5) # x, or y location for each plane + tau[k] = call(tau_funcs[comp], self.M.faces[k, :]) + + # integrate components parallel to the plane of integration + if self.location == "edges": + analytic = 2.24166666666667 # Found using sympy. + + p = mesh.edges + Ec = np.c_[ex(*p.T), ey(*p.T)] + E = mesh.project_edge_vector(Ec) + + if self.invert_model: + A = self.M.get_edge_inner_product_surface(1 / tau, invert_model=True) + else: + A = self.M.get_edge_inner_product_surface(tau) + + numeric = E.T.dot(A.dot(E)) + + # integrate component normal to the plane of integration + elif self.location == "faces": + analytic = 1.59895833333333 # Found using sympy. + + p = mesh.faces + Fc = np.c_[ex(*p.T), ey(*p.T)] + F = mesh.project_face_vector(Fc) + + if self.invert_model: + A = self.M.get_face_inner_product_surface(1 / tau, invert_model=True) + else: + A = self.M.get_face_inner_product_surface(tau) + + numeric = F.T.dot(A.dot(F)) + + err = np.abs(numeric - analytic) + + return err + + def test_order1_edges(self): + self.name = "Edge Inner Product - Isotropic" + self.location = "edges" + self.invert_model = False + self.orderTest() + + def test_order1_edges_invert_model(self): + self.name = "Edge Inner Product - Isotropic - invert_model" + self.location = "edges" + self.invert_model = True + self.orderTest() + + def test_order1_faces(self): + self.name = "Face Inner Product - Isotropic" + self.location = "faces" + self.invert_model = False + self.orderTest() + + def test_order1_faces_invert_model(self): + self.name = "Face Inner Product - Isotropic - invert_model" + self.location = "faces" + self.invert_model = True + self.orderTest() + + +class TestInnerProductsEdgeProperties2D(discretize.tests.OrderTest): + """Integrate a function over a line within a unit cube domain + using edgeInnerProducts.""" + + meshTypes = ["uniformTree", "notatreeTree"] + meshDimension = 2 + meshSizes = [16, 32] + + def getError(self): + call = lambda fun, xy: fun(xy[:, 0], xy[:, 1]) # NOQA F841 + + ex = lambda x, y: x**2 + y + ey = lambda x, y: (x**2) * y + + tau_x = lambda x, y: x + 1 # x-face properties # NOQA F841 + tau_y = lambda x, y: y + 2 # y-face properties # NOQA F841 + + mesh = self.M + + tau = 1e-8 * np.ones(mesh.n_edges) + for ii, comp in enumerate(["x", "y"]): + k = np.isclose(self.M.edges[:, ii - 1], 0.5) & np.isclose( + self.M.edges[:, ii - 2], 0.5 + ) # x, y or z location for each line + tau[k] = eval("call(tau_{}, self.M.edges[k, :])".format(comp)) + + analytic = 1.38229166666667 # Found using sympy. + + p = mesh.edges + Ec = np.c_[ex(*p.T), ey(*p.T)] + E = mesh.project_edge_vector(Ec) + + if self.invert_model: + A = self.M.get_edge_inner_product_line(1 / tau, invert_model=True) + else: + A = self.M.get_edge_inner_product_line(tau) + + numeric = E.T.dot(A.dot(E)) + + err = np.abs(numeric - analytic) + + return err + + def test_order1_edges(self): + self.name = "Edge Inner Product - Isotropic" + self.location = "edges" + self.invert_model = False + self.orderTest() + + def test_order1_edges_invert_model(self): + self.name = "Edge Inner Product - Isotropic - invert_model" + self.location = "edges" + self.invert_model = True + self.orderTest() + + class TestInnerProducts3D(discretize.tests.OrderTest): meshSizes = [8, 16, 32] meshTypes = ["uniform simplex mesh"] @@ -330,6 +474,163 @@ def test_order6_faces_invert_model(self): self.orderTest() +class TestInnerProductsFaceProperties3D(discretize.tests.OrderTest): + meshSizes = [8, 16, 32] + meshTypes = ["uniform simplex mesh"] + + def setupMesh(self, n): + points, simplices = example_simplex_mesh((n, n, n)) + self.M = discretize.SimplexMesh(points, simplices) + return 1.0 / n + + def getError(self): + + call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) + + ex = lambda x, y, z: x**2 + y * z + ey = lambda x, y, z: (z**2) * x + y * z + ez = lambda x, y, z: y**2 + x * z + + tau_funcs = { + "x": lambda x, y, z: y * z + 1, # x-face properties + "y": lambda x, y, z: x * z + 2, # y-face properties + "z": lambda x, y, z: 3 + x * y, # z-face properties + } + + mesh = self.M + + tau = 1e-8 * np.ones(mesh.n_faces) + for ii, comp in enumerate(["x", "y", "z"]): + k = np.isclose(self.M.faces[:, ii], 0.5) # x, or y location for each plane + tau[k] = call(tau_funcs[comp], self.M.faces[k, :]) + + # integrate components parallel to the plane of integration + if self.location == "edges": + analytic = 5.02760416666667 # Found using sympy. + + p = mesh.edges + Ec = np.c_[ex(*p.T), ey(*p.T), ez(*p.T)] + E = mesh.project_edge_vector(Ec) + + if self.invert_model: + A = self.M.get_edge_inner_product_surface(1 / tau, invert_model=True) + else: + A = self.M.get_edge_inner_product_surface(tau) + + numeric = E.T.dot(A.dot(E)) + + # integrate component normal to the plane of integration + elif self.location == "faces": + analytic = 2.66979166666667 # Found using sympy. + + p = mesh.faces + Fc = np.c_[ex(*p.T), ey(*p.T), ez(*p.T)] + F = mesh.project_face_vector(Fc) + + if self.invert_model: + A = self.M.get_face_inner_product_surface(1 / tau, invert_model=True) + else: + A = self.M.get_face_inner_product_surface(tau) + + numeric = F.T.dot(A.dot(F)) + + print(analytic) + print(numeric) + print(analytic / numeric) + + err = np.abs(numeric - analytic) + + return err + + def test_order1_edges(self): + self.name = "Edge Inner Product - Isotropic" + self.location = "edges" + self.invert_model = False + self.orderTest() + + def test_order1_edges_invert_model(self): + self.name = "Edge Inner Product - Isotropic - invert_model" + self.location = "edges" + self.invert_model = True + self.orderTest() + + def test_order1_faces(self): + self.name = "Face Inner Product - Isotropic" + self.location = "faces" + self.invert_model = False + self.orderTest() + + def test_order1_faces_invert_model(self): + self.name = "Face Inner Product - Isotropic - invert_model" + self.location = "faces" + self.invert_model = True + self.orderTest() + + +class TestInnerProductsEdgeProperties3D(discretize.tests.OrderTest): + """Integrate a function over a line within a unit cube domain + using edgeInnerProducts.""" + + meshSizes = [8, 16, 32] + meshTypes = ["uniform simplex mesh"] + + def setupMesh(self, n): + points, simplices = example_simplex_mesh((n, n, n)) + self.M = discretize.SimplexMesh(points, simplices) + return 1.0 / n + + def getError(self): + call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) + + ex = lambda x, y, z: x**2 + y * z + ey = lambda x, y, z: (z**2) * x + y * z + ez = lambda x, y, z: y**2 + x * z + + tau_funcs = { + "x": lambda x, y, z: x + 1, # x-face properties + "y": lambda x, y, z: y + 2, # y-face properties + "z": lambda x, y, z: 3 * z + 1, # z-face properties + } + + mesh = self.M + + tau = 1e-8 * np.ones(mesh.n_edges) + for ii, comp in enumerate(["x", "y", "z"]): + k = np.isclose(self.M.edges[:, ii - 1], 0.5) & np.isclose( + self.M.edges[:, ii - 2], 0.5 + ) # x, y or z location for each line + tau[k] = call(tau_funcs[comp], self.M.edges[k, :]) + + analytic = 1.98906250000000 # Found using sympy. + + p = mesh.edges + Ec = np.c_[ex(*p.T), ey(*p.T), ez(*p.T)] + E = mesh.project_edge_vector(Ec) + + if self.invert_model: + A = self.M.get_edge_inner_product_line(1 / tau, invert_model=True) + else: + A = self.M.get_edge_inner_product_line(tau) + + numeric = E.T.dot(A.dot(E)) + + err = np.abs(numeric - analytic) + + return err + + def test_order1_edges(self): + self.name = "Edge Inner Product - Isotropic" + self.location = "edges" + self.invert_model = False + self.orderTest() + + def test_order1_edges_invert_model(self): + self.name = "Edge Inner Product - Isotropic - invert_model" + self.location = "edges" + self.invert_model = True + self.orderTest() + + class TestInnerProductsDerivs(unittest.TestCase): def doTestFace(self, h, rep): nodes, simplices = example_simplex_mesh(h) @@ -412,6 +713,55 @@ def test_EdgeIP_3D_tensor(self): self.assertTrue(self.doTestEdge([10, 4, 5], 6)) +@pytest.mark.parametrize("u_type", ["edge", "face"]) +@pytest.mark.parametrize("h", [(10, 4), (10, 4, 5)], ids=["2D", "3D"]) +@pytest.mark.parametrize("rep", [0, 1], ids=["uniform", "isotropic"]) +def test_surface_inner_product_prop_deriv(u_type, h, rep): + rng = np.random.default_rng(6732) + nodes, simplices = example_simplex_mesh(h) + mesh = discretize.SimplexMesh(nodes, simplices) + tau = rng.uniform(1, 2, 1) if rep == 0 else rng.uniform(1, 2, mesh.n_faces * rep) + + match u_type: + case "edge": + v = rng.uniform(1, 2, mesh.n_edges) + + def fun(tau): + M = mesh.get_edge_inner_product_surface(tau) + Md = mesh.get_edge_inner_product_surface_deriv(tau) + return M * v, Md(v) + + case "face": + v = rng.uniform(1, 2, mesh.n_faces) + + def fun(tau): + M = mesh.get_face_inner_product_surface(tau) + Md = mesh.get_face_inner_product_surface_deriv(tau) + return M * v, Md(v) + + case _: + raise Exception("Invalid test parameter.") + + discretize.tests.check_derivative(fun, tau, num=5, random_seed=rng) + + +@pytest.mark.parametrize("h", [(10, 4), (10, 4, 5)], ids=["2D", "3D"]) +@pytest.mark.parametrize("rep", [0, 1], ids=["uniform", "isotropic"]) +def test_line_inner_product_prop_deriv(h, rep): + rng = np.random.default_rng(6732) + nodes, simplices = example_simplex_mesh(h) + mesh = discretize.SimplexMesh(nodes, simplices) + v = rng.uniform(1, 2, mesh.n_edges) + tau = rng.uniform(1, 2, 1) if rep == 0 else rng.uniform(1, 2, mesh.n_edges * rep) + + def fun(tau): + M = mesh.get_edge_inner_product_line(tau) + Md = mesh.get_edge_inner_product_line_deriv(tau) + return M * v, Md(v) + + discretize.tests.check_derivative(fun, tau, num=5, random_seed=rng) + + class Test2DBoundaryIntegral(discretize.tests.OrderTest): meshSizes = [8, 16, 32] meshTypes = ["uniform simplex mesh"]