diff --git a/discretize/curvilinear_mesh.py b/discretize/curvilinear_mesh.py index 1c6181a42..d06edb572 100644 --- a/discretize/curvilinear_mesh.py +++ b/discretize/curvilinear_mesh.py @@ -797,10 +797,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), diff --git a/discretize/unstructured_mesh.py b/discretize/unstructured_mesh.py index 408442cd0..61f707e47 100644 --- a/discretize/unstructured_mesh.py +++ b/discretize/unstructured_mesh.py @@ -519,7 +519,64 @@ def get_edge_inner_product( # NOQA D102 "The inverse of the inner product matrix with a tetrahedral mesh is not supported." ) 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 + 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 + dim = self.dim + + if dim == 2: + + if model is None: + Mu = sp.diags(face_areas) + else: + if invert_model: + model = 1.0 / model + + if (model.size == 1) | (model.size == n_faces): + Mu = sp.diags(model * face_areas) + else: + raise ValueError( + "Unrecognized size of model vector.", + "Must be scalar or have length equal to total number of faces." + ) + + return Mu + + else: + Ps = self._get_edge_surf_int_proj_mats(with_area=False) + + if model is None: + Mu = sp.diags(np.tile(face_areas, dim)) + 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)) + 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 @@ -668,10 +725,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/simplex/test_inner_products.py b/tests/simplex/test_inner_products.py index a5cae4ce6..4293b8efb 100644 --- a/tests/simplex/test_inner_products.py +++ b/tests/simplex/test_inner_products.py @@ -180,6 +180,148 @@ def test_order3_faces_invert_model(self): self.invert_model = True 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_x = lambda x, y: 2 * y + 1 # x-face properties # NOQA F841 + tau_y = lambda x, y: x + 2 # y-face properties # NOQA F841 + + 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] = eval("call(tau_{}, self.M.faces[k, :])".format(comp)) + + # 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] @@ -329,6 +471,160 @@ def test_order6_faces_invert_model(self): self.invert_model = True 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_x = lambda x, y, z: y * z + 1 # x-face properties # NOQA F841 + tau_y = lambda x, y, z: x * z + 2 # y-face properties # NOQA F841 + tau_z = lambda x, y, z: 3 + x * y # z-face properties # NOQA F841 + + 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] = eval("call(tau_{}, self.M.faces[k, :])".format(comp)) + + # 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_x = lambda x, y, z: x + 1 # x-face properties # NOQA F841 + tau_y = lambda x, y, z: y + 2 # y-face properties # NOQA F841 + tau_z = lambda x, y, z: 3 * z + 1 # z-face properties # NOQA F841 + + 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] = eval("call(tau_{}, self.M.edges[k, :])".format(comp)) + + + 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): @@ -412,6 +708,87 @@ def test_EdgeIP_3D_tensor(self): self.assertTrue(self.doTestEdge([10, 4, 5], 6)) +class TestFacePropertiesInnerProductsDerivs(unittest.TestCase): + def doTestFace(self, h, rep): + nodes, simplices = example_simplex_mesh(h) + mesh = discretize.SimplexMesh(nodes, simplices) + v = np.random.rand(mesh.n_faces) + tau = np.random.rand(1) if rep == 0 else np.random.rand(mesh.nF * rep) + + 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) + + print("Face", rep) + return discretize.tests.check_derivative(fun, tau, num=5, plotIt=False) + + def doTestEdge(self, h, rep): + nodes, simplices = example_simplex_mesh(h) + mesh = discretize.SimplexMesh(nodes, simplices) + v = np.random.rand(mesh.n_edges) + tau = np.random.rand(1) if rep == 0 else np.random.rand(mesh.nF * rep) + + 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) + + print("Edge", rep) + return discretize.tests.check_derivative(fun, tau, num=5, plotIt=False) + + def test_FaceIP_2D_float(self): + self.assertTrue(self.doTestFace([10, 4], 0)) + + def test_FaceIP_3D_float(self): + self.assertTrue(self.doTestFace([10, 4, 5], 0)) + + def test_FaceIP_2D_isotropic(self): + self.assertTrue(self.doTestFace([10, 4], 1)) + + def test_FaceIP_3D_isotropic(self): + self.assertTrue(self.doTestFace([10, 4, 5], 1)) + + def test_EdgeIP_2D_float(self): + self.assertTrue(self.doTestEdge([10, 4], 0)) + + def test_EdgeIP_3D_float(self): + self.assertTrue(self.doTestEdge([10, 4, 5], 0)) + + def test_EdgeIP_2D_isotropic(self): + self.assertTrue(self.doTestEdge([10, 4], 1)) + + def test_EdgeIP_3D_isotropic(self): + self.assertTrue(self.doTestEdge([10, 4, 5], 1)) + +class TestFacePropertiesInnerProductsDerivs(unittest.TestCase): + def doTestEdge(self, h, rep): + nodes, simplices = example_simplex_mesh(h) + mesh = discretize.SimplexMesh(nodes, simplices) + v = np.random.rand(mesh.n_edges) + tau = np.random.rand(1) if rep == 0 else np.random.rand(mesh.nE * 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) + + print("Edge", rep) + return discretize.tests.check_derivative(fun, tau, num=5, plotIt=False) + + def test_EdgeIP_2D_float(self): + self.assertTrue(self.doTestEdge([10, 4], 0)) + + def test_EdgeIP_3D_float(self): + self.assertTrue(self.doTestEdge([10, 4, 5], 0)) + + def test_EdgeIP_2D_isotropic(self): + self.assertTrue(self.doTestEdge([10, 4], 1)) + + def test_EdgeIP_3D_isotropic(self): + self.assertTrue(self.doTestEdge([10, 4, 5], 1)) + + class Test2DBoundaryIntegral(discretize.tests.OrderTest): meshSizes = [8, 16, 32] meshTypes = ["uniform simplex mesh"]