Skip to content

Commit

Permalink
Add more lazy evaluation to MappingAffine (#1106)
Browse files Browse the repository at this point in the history 
* Add more lazy evaluation to MappingAffine

* Allow prespecifying tind in MappingAffine constructor
  • Loading branch information
@kinnala
kinnala authored Mar 11, 2024
1 parent 580df1a commit 65b4a17
Show file tree
Hide file tree
Showing 2 changed files with 147 additions and 71 deletions.
202 changes: 132 additions & 70 deletions skfem/mapping/mapping_affine.py
View file
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
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)

0 comments on commit 65b4a17

Please sign in to comment.