From b86c1c813703bda84e77482ecba2d5178c5ccabc Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 11 Mar 2024 12:41:10 +0200 Subject: [PATCH 1/6] Add more lazy evaluation to MappingAffine --- skfem/mapping/mapping_affine.py | 161 ++++++++++++++++++++------------ 1 file changed, 99 insertions(+), 62 deletions(-) diff --git a/skfem/mapping/mapping_affine.py b/skfem/mapping/mapping_affine.py index 26379ef11..2508b0e0c 100644 --- a/skfem/mapping/mapping_affine.py +++ b/skfem/mapping/mapping_affine.py @@ -7,69 +7,124 @@ class MappingAffine(Mapping): """An affine mapping for simplical elements.""" def __init__(self, mesh): - dim = mesh.p.shape[0] + self.dim = mesh.p.shape[0] + self.mesh = mesh - if mesh.t.shape[0] > 0: - nt = mesh.t.shape[1] + @property + def A(self): + """Affine mapping matrix.""" + if not hasattr(self, '_A'): + self._init_Ab() + return self._A + + @property + def b(self): + """Affine mapping constant.""" + if not hasattr(self, '_b'): + self._init_Ab() + return self._b + + @property + def detA(self): + """Affine mapping determinant.""" + if not hasattr(self, '_detA'): + self._init_invA() + return self._detA + + @property + def invA(self): + """Affine mapping matrix inverse.""" + if not hasattr(self, '_invA'): + self._init_invA() + return self._invA + + def _init_Ab(self): + """For lazy evaluation of interior mapping.""" + if self.mesh.t.shape[0] > 0: # todo: why this guard exists? + nt = self.mesh.t.shape[1] + dim = self.dim # initialize the affine mapping - self.A = np.empty((dim, dim, nt)) - self.b = np.empty((dim, nt)) + self._A = np.empty((dim, dim, nt)) + self._b = np.empty((dim, nt)) for i in range(dim): - self.b[i] = mesh.p[i, mesh.t[0]] + self._b[i] = self.mesh.p[i, self.mesh.t[0]] for j in range(dim): - self.A[i, j] = (mesh.p[i, mesh.t[j + 1]] - - mesh.p[i, mesh.t[0]]) - + self._A[i, j] = (self.mesh.p[i, self.mesh.t[j + 1]] - + self.mesh.p[i, self.mesh.t[0]]) + + def _init_invA(self): + """For lazy evaluation of interior mapping.""" + if self.mesh.t.shape[0] > 0: # todo: why this guard exists? + nt = self.mesh.t.shape[1] + dim = self.dim # determinants if dim == 1: - self.detA = self.A[0, 0] + self._detA = self.A[0, 0] elif dim == 2: - self.detA = (self.A[0, 0] * self.A[1, 1] - - self.A[0, 1] * self.A[1, 0]) + self._detA = (self.A[0, 0] * self.A[1, 1] - + self.A[0, 1] * self.A[1, 0]) elif dim == 3: - self.detA = self.A[0, 0] * (self.A[1, 1] * self.A[2, 2] - - self.A[1, 2] * self.A[2, 1]) -\ - self.A[0, 1] * (self.A[1, 0] * self.A[2, 2] - - self.A[1, 2] * self.A[2, 0]) +\ - self.A[0, 2] * (self.A[1, 0] * self.A[2, 1] - - self.A[1, 1] * self.A[2, 0]) + self._detA = (self.A[0, 0] * (self.A[1, 1] * self.A[2, 2] - + self.A[1, 2] * self.A[2, 1]) - + self.A[0, 1] * (self.A[1, 0] * self.A[2, 2] - + self.A[1, 2] * self.A[2, 0]) + + self.A[0, 2] * (self.A[1, 0] * self.A[2, 1] - + self.A[1, 1] * self.A[2, 0])) else: raise Exception("Not implemented for the given dimension.") # affine mapping inverses - self.invA = np.empty((dim, dim, nt)) + self._invA = np.empty((dim, dim, nt)) if dim == 1: - self.invA[0, 0] = 1. / self.A[0, 0] + self._invA[0, 0] = 1. / self.A[0, 0] elif dim == 2: - self.invA[0, 0] = self.A[1, 1] / self.detA # noqa - self.invA[0, 1] = -self.A[0, 1] / self.detA - self.invA[1, 0] = -self.A[1, 0] / self.detA - self.invA[1, 1] = self.A[0, 0] / self.detA # noqa + self._invA[0, 0] = self.A[1, 1] / self.detA # noqa + self._invA[0, 1] = -self.A[0, 1] / self.detA + self._invA[1, 0] = -self.A[1, 0] / self.detA + self._invA[1, 1] = self.A[0, 0] / self.detA # noqa elif dim == 3: - self.invA[0, 0] = (-self.A[1, 2] * self.A[2, 1] + - self.A[1, 1] * self.A[2, 2]) / self.detA - self.invA[1, 0] = (self.A[1, 2] * self.A[2, 0] - - self.A[1, 0] * self.A[2, 2]) / self.detA - self.invA[2, 0] = (-self.A[1, 1] * self.A[2, 0] + - self.A[1, 0] * self.A[2, 1]) / self.detA - self.invA[0, 1] = (self.A[0, 2] * self.A[2, 1] - - self.A[0, 1] * self.A[2, 2]) / self.detA - self.invA[1, 1] = (-self.A[0, 2] * self.A[2, 0] + - self.A[0, 0] * self.A[2, 2]) / self.detA - self.invA[2, 1] = (self.A[0, 1] * self.A[2, 0] - - self.A[0, 0] * self.A[2, 1]) / self.detA - self.invA[0, 2] = (-self.A[0, 2] * self.A[1, 1] + - self.A[0, 1] * self.A[1, 2]) / self.detA - self.invA[1, 2] = (self.A[0, 2] * self.A[1, 0] - - self.A[0, 0] * self.A[1, 2]) / self.detA - self.invA[2, 2] = (-self.A[0, 1] * self.A[1, 0] + - self.A[0, 0] * self.A[1, 1]) / self.detA + self._invA[0, 0] = (-self.A[1, 2] * self.A[2, 1] + + self.A[1, 1] * self.A[2, 2]) / self.detA + self._invA[1, 0] = (self.A[1, 2] * self.A[2, 0] - + self.A[1, 0] * self.A[2, 2]) / self.detA + self._invA[2, 0] = (-self.A[1, 1] * self.A[2, 0] + + self.A[1, 0] * self.A[2, 1]) / self.detA + self._invA[0, 1] = (self.A[0, 2] * self.A[2, 1] - + self.A[0, 1] * self.A[2, 2]) / self.detA + self._invA[1, 1] = (-self.A[0, 2] * self.A[2, 0] + + self.A[0, 0] * self.A[2, 2]) / self.detA + self._invA[2, 1] = (self.A[0, 1] * self.A[2, 0] - + self.A[0, 0] * self.A[2, 1]) / self.detA + self._invA[0, 2] = (-self.A[0, 2] * self.A[1, 1] + + self.A[0, 1] * self.A[1, 2]) / self.detA + self._invA[1, 2] = (self.A[0, 2] * self.A[1, 0] - + self.A[0, 0] * self.A[1, 2]) / self.detA + self._invA[2, 2] = (-self.A[0, 1] * self.A[1, 0] + + self.A[0, 0] * self.A[1, 1]) / self.detA else: raise Exception("Not implemented for the given dimension.") - self.dim = dim - self.mesh = mesh # this is required in ElementH2 + @property + def B(self): + """Boundary affine mapping matrix.""" + if not hasattr(self, '_B'): + self._init_boundary_mapping() + return self._B + + @property + def c(self): + """Boundary affine mapping constant.""" + if not hasattr(self, '_c'): + self._init_boundary_mapping() + return self._c + + @property + def detB(self): + """Boundary affine mapping determinant.""" + if not hasattr(self, '_detB'): + self._init_boundary_mapping() + return self._detB def _init_boundary_mapping(self): """For lazy evaluation of boundary mapping.""" @@ -100,24 +155,6 @@ def _init_boundary_mapping(self): else: raise Exception("Not implemented for the given dimension.") - @property - def B(self): - if not hasattr(self, '_B'): - self._init_boundary_mapping() - return self._B - - @property - def c(self): - if not hasattr(self, '_c'): - self._init_boundary_mapping() - return self._c - - @property - def detB(self): - if not hasattr(self, '_detB'): - self._init_boundary_mapping() - return self._detB - def F(self, X, tind=None): if tind is None: A, b = self.A, self.b From 1a8745844e77f5e419ee80eb06ab375923f63c54 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 11 Mar 2024 12:43:29 +0200 Subject: [PATCH 2/6] Fix flake8 --- skfem/mapping/mapping_affine.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skfem/mapping/mapping_affine.py b/skfem/mapping/mapping_affine.py index 2508b0e0c..c60388f36 100644 --- a/skfem/mapping/mapping_affine.py +++ b/skfem/mapping/mapping_affine.py @@ -66,7 +66,7 @@ def _init_invA(self): self.A[0, 1] * self.A[1, 0]) elif dim == 3: self._detA = (self.A[0, 0] * (self.A[1, 1] * self.A[2, 2] - - self.A[1, 2] * self.A[2, 1]) - + self.A[1, 2] * self.A[2, 1]) - self.A[0, 1] * (self.A[1, 0] * self.A[2, 2] - self.A[1, 2] * self.A[2, 0]) + self.A[0, 2] * (self.A[1, 0] * self.A[2, 1] - From 8e7cdc8b83289d6bc6416e0e63d8a2950ceb1077 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 11 Mar 2024 13:09:08 +0200 Subject: [PATCH 3/6] Allow prespecifying tind in MappingAffine constructor --- skfem/mapping/mapping_affine.py | 51 ++++++++++++++++++++++++--------- 1 file changed, 38 insertions(+), 13 deletions(-) diff --git a/skfem/mapping/mapping_affine.py b/skfem/mapping/mapping_affine.py index c60388f36..41255983a 100644 --- a/skfem/mapping/mapping_affine.py +++ b/skfem/mapping/mapping_affine.py @@ -6,9 +6,23 @@ class MappingAffine(Mapping): """An affine mapping for simplical elements.""" - def __init__(self, mesh): + def __init__(self, mesh, tind=None): + """Local-to-global mappings for simplex elements. + + Parameters + ---------- + mesh + An object of type :class:`~skfem.mesh.Mesh`. + tind + An optional array of element indices for memory + optimizations. If provided, only the included element + indices are used and all methods will ignore the + respective tind parameter. + + """ self.dim = mesh.p.shape[0] self.mesh = mesh + self.tind = tind @property def A(self): @@ -41,22 +55,33 @@ def invA(self): def _init_Ab(self): """For lazy evaluation of interior mapping.""" if self.mesh.t.shape[0] > 0: # todo: why this guard exists? - nt = self.mesh.t.shape[1] + nt = self.mesh.t.shape[1] if self.tind is None else len(self.tind) dim = self.dim # initialize the affine mapping self._A = np.empty((dim, dim, nt)) self._b = np.empty((dim, nt)) - for i in range(dim): - self._b[i] = self.mesh.p[i, self.mesh.t[0]] - for j in range(dim): - self._A[i, j] = (self.mesh.p[i, self.mesh.t[j + 1]] - - self.mesh.p[i, self.mesh.t[0]]) + if self.tind is None: + for i in range(dim): + self._b[i] = self.mesh.p[i, self.mesh.t[0]] + for j in range(dim): + self._A[i, j] = ( + self.mesh.p[i, self.mesh.t[j + 1]] - + self.mesh.p[i, self.mesh.t[0]] + ) + else: + for i in range(dim): + self._b[i] = self.mesh.p[i, self.mesh.t[0, self.tind]] + for j in range(dim): + self._A[i, j] = ( + self.mesh.p[i, self.mesh.t[j + 1, self.tind]] - + self.mesh.p[i, self.mesh.t[0, self.tind]] + ) def _init_invA(self): """For lazy evaluation of interior mapping.""" if self.mesh.t.shape[0] > 0: # todo: why this guard exists? - nt = self.mesh.t.shape[1] + nt = self.mesh.t.shape[1] if self.tind is None else len(self.tind) dim = self.dim # determinants if dim == 1: @@ -156,7 +181,7 @@ def _init_boundary_mapping(self): raise Exception("Not implemented for the given dimension.") def F(self, X, tind=None): - if tind is None: + if tind is None or self.tind is not None: A, b = self.A, self.b else: A, b = self.A[:, :, tind], self.b[:, tind] @@ -168,7 +193,7 @@ def F(self, X, tind=None): raise NotImplementedError def invF(self, x, tind=None): - if tind is None: + if tind is None or self.tind is not None: invA, b = self.invA, self.b else: invA, b = self.invA[:, :, tind], self.b[:, tind] @@ -178,7 +203,7 @@ def invF(self, x, tind=None): return np.einsum('ijk,jkl->ikl', invA, y) def detDF(self, X, tind=None): - if tind is None: + if tind is None or self.tind is not None: detDF = self.detA else: detDF = self.detA[tind] @@ -186,7 +211,7 @@ def detDF(self, X, tind=None): return np.tile(detDF, (X.shape[-1], 1)).T def DF(self, X, tind=None): - if tind is None: + if tind is None or self.tind is not None: DF = self.A else: DF = self.A[:, :, tind] @@ -198,7 +223,7 @@ def DF(self, X, tind=None): raise NotImplementedError def invDF(self, X, tind=None): - if tind is None: + if tind is None or self.tind is not None: invDF = self.invA else: invDF = self.invA[:, :, tind] From 2cdc47fce1685b488ea34801b4033255602bcc74 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 11 Mar 2024 13:38:44 +0200 Subject: [PATCH 4/6] Add test for memory optimization --- tests/test_mapping.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/tests/test_mapping.py b/tests/test_mapping.py index a22de28c7..6dc1d3570 100644 --- a/tests/test_mapping.py +++ b/tests/test_mapping.py @@ -1,7 +1,7 @@ import unittest import numpy as np -from skfem.mesh import MeshHex, MeshQuad, MeshTri +from skfem.mesh import MeshHex, MeshQuad, MeshTri, MeshTet from skfem.element import ElementHex1, ElementQuad1, ElementHex2 from skfem.assembly import FacetBasis @@ -87,3 +87,16 @@ class TestInverseMappingHex2(TestInverseMappingHex): """This should be equivalent to TestInverseMappingHex.""" element = ElementHex2 + + +def test_mapping_memory_optimization(): + + m = MeshTet.init_tensor(np.linspace(0, 1, 100), + np.linspace(0, 1, 100), + np.linspace(0, 1, 100)) + m = m.with_subdomains({ + 'omega0': [0, 1, 2, 3, 4, 5], + }) + orig = MappingAffine(m) + opt = MappingAffine(m, tind=m.subdomains['omega0']) + assert orig.detA > opt.detA From 29c2d9d5959a742fa4ac0075854de6686b0b27de Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 11 Mar 2024 13:49:25 +0200 Subject: [PATCH 5/6] Fix import --- tests/test_mapping.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_mapping.py b/tests/test_mapping.py index 6dc1d3570..6b98d3f08 100644 --- a/tests/test_mapping.py +++ b/tests/test_mapping.py @@ -4,6 +4,7 @@ from skfem.mesh import MeshHex, MeshQuad, MeshTri, MeshTet from skfem.element import ElementHex1, ElementQuad1, ElementHex2 from skfem.assembly import FacetBasis +from skfem.mapping import MappingAffine class TestIsoparamNormals(unittest.TestCase): From 2006946ab2b19959b49285d8b5ac01b8c228d2db Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 11 Mar 2024 14:03:24 +0200 Subject: [PATCH 6/6] Fix assertion --- tests/test_mapping.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_mapping.py b/tests/test_mapping.py index 6b98d3f08..7c1375f52 100644 --- a/tests/test_mapping.py +++ b/tests/test_mapping.py @@ -100,4 +100,4 @@ def test_mapping_memory_optimization(): }) orig = MappingAffine(m) opt = MappingAffine(m, tind=m.subdomains['omega0']) - assert orig.detA > opt.detA + assert len(orig.detA) > len(opt.detA)