Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add more lazy evaluation to MappingAffine #1106

Merged
merged 6 commits into from
Mar 11, 2024
Merged

Add more lazy evaluation to MappingAffine #1106

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
b86c1c8
Add more lazy evaluation to MappingAffine
kinnala Mar 11, 2024
1a87458
Fix flake8
kinnala Mar 11, 2024
8e7cdc8
Allow prespecifying tind in MappingAffine constructor
kinnala Mar 11, 2024
2cdc47f
Add test for memory optimization
kinnala Mar 11, 2024
29c2d9d
Fix import
kinnala Mar 11, 2024
2006946
Fix assertion
kinnala Mar 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
202 changes: 132 additions & 70 deletions skfem/mapping/mapping_affine.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,70 +6,150 @@
class MappingAffine(Mapping):
"""An affine mapping for simplical elements."""

def __init__(self, mesh):
dim = mesh.p.shape[0]
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

if mesh.t.shape[0] > 0:
nt = mesh.t.shape[1]
# initialize the affine mapping
self.A = np.empty((dim, dim, nt))
self.b = np.empty((dim, nt))
@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

for i in range(dim):
self.b[i] = mesh.p[i, 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]])
@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] 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))

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] if self.tind is None else len(self.tind)
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."""
Expand Down Expand Up @@ -100,26 +180,8 @@ 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:
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]
Expand All @@ -131,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]
Expand All @@ -141,15 +203,15 @@ 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]

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]
Expand All @@ -161,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]
Expand Down
16 changes: 15 additions & 1 deletion tests/test_mapping.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,9 +1,10 @@
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
from skfem.mapping import MappingAffine


class TestIsoparamNormals(unittest.TestCase):
Expand Down Expand Up @@ -87,3 +88,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 len(orig.detA) > len(opt.detA)
Loading