Skip to content

Commit

Permalink
Merge pull request #283 from kinnala/optimize-mapping-isoparam-by-cac…
Browse files Browse the repository at this point in the history 
…hing

Optimize mapping isoparam by caching
  • Loading branch information
@gdmcbain
gdmcbain committed Dec 18, 2019
2 parents 48c836c + ca3d86d commit e77096e
Show file tree
Hide file tree
Showing 2 changed files with 30 additions and 58 deletions.
86 changes: 29 additions & 57 deletions skfem/mapping/mapping_isoparametric.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,20 +145,17 @@ def invF(self, x, tind=None, newton_max_iters = 50, newton_tol = 1e-8):
def F(self, X, tind=None):
return np.array([self.map(i, X, tind) for i in range(X.shape[0])])

def detDF(self, X, tind=None):
def detDF(self, X, tind=None, J=None):
if J is None:
J = [[self.J(i, j, X, tind=tind) for j in range(self.dim)]
for i in range(self.dim)]

if self.dim == 2:
detDF = (self.J(0, 0, X, tind=tind) * self.J(1, 1, X, tind=tind) -
self.J(0, 1, X, tind=tind) * self.J(1, 0, X, tind=tind))
detDF = J[0][0] * J[1][1] - J[0][1] * J[1][0]
elif self.dim == 3:
detDF = (self.J(0, 0, X, tind=tind) * \
(self.J(1, 1, X, tind=tind) * self.J(2, 2, X, tind=tind) -
self.J(1, 2, X, tind=tind) * self.J(2, 1, X, tind=tind))
- self.J(0, 1, X, tind=tind) * \
(self.J(1, 0, X, tind=tind) * self.J(2, 2, X, tind=tind) -
self.J(1, 2, X, tind=tind) * self.J(2, 0, X, tind=tind))
+ self.J(0, 2, X, tind=tind) * \
(self.J(1, 0, X, tind=tind) * self.J(2, 1, X, tind=tind) -
self.J(1, 1, X, tind=tind) * self.J(2, 0, X, tind=tind)))
detDF = (J[0][0] * (J[1][1] * J[2][2] - J[1][2] * J[2][1]) -
J[0][1] * (J[1][0] * J[2][2] - J[1][2] * J[2][0]) +
J[0][2] * (J[1][0] * J[2][1] - J[1][1] * J[2][0]))
else:
raise Exception("Not implemented for the given dimension.")

Expand All @@ -168,57 +165,32 @@ def detDF(self, X, tind=None):
return detDF

def invDF(self, X, tind=None):
detDF = self.detDF(X, tind)
J = [[self.J(i, j, X, tind=tind) for j in range(self.dim)]
for i in range(self.dim)]
detDF = self.detDF(X, tind, J=J)
invDF = np.empty((self.dim, self.dim) + J[0][0].shape)

if self.dim == 2:
invDF = np.empty((2, 2) + self.J(0, 0, X, tind=tind).shape)
invDF[0, 0] = self.J(1, 1, X, tind=tind) / detDF
invDF[0, 1] = -self.J(0, 1, X, tind=tind) / detDF
invDF[1, 0] = -self.J(1, 0, X, tind=tind) / detDF
invDF[1, 1] = self.J(0, 0, X, tind=tind) / detDF
detDF = self.detDF(X, tind)
invDF[0, 0] = J[1][1]
invDF[0, 1] = -J[0][1]
invDF[1, 0] = -J[1][0]
invDF[1, 1] = J[0][0]
elif self.dim == 3:
invDF = np.empty((3, 3) + self.J(0, 0, X, tind=tind).shape)
invDF[0, 0] = (-self.J(1, 2, X, tind=tind) * \
self.J(2, 1, X, tind=tind) +
self.J(1, 1, X, tind=tind) * \
self.J(2, 2, X, tind=tind)) / detDF
invDF[1, 0] = (self.J(1, 2, X, tind=tind) * \
self.J(2, 0, X, tind=tind) -
self.J(1, 0, X, tind=tind) * \
self.J(2, 2, X, tind=tind)) / detDF
invDF[2, 0] = (-self.J(1, 1, X, tind=tind) * \
self.J(2, 0, X, tind=tind) +
self.J(1, 0, X, tind=tind) * \
self.J(2, 1, X, tind=tind)) / detDF
invDF[0, 1] = (self.J(0, 2, X, tind=tind) * \
self.J(2, 1, X, tind=tind) +
-self.J(0, 1, X, tind=tind) * \
self.J(2, 2, X, tind=tind)) / detDF
invDF[1, 1] = (-self.J(0, 2, X, tind=tind) * \
self.J(2, 0, X, tind=tind) +
self.J(0, 0, X, tind=tind) * \
self.J(2, 2, X, tind=tind)) / detDF
invDF[2, 1] = (self.J(0, 1, X, tind=tind) * \
self.J(2, 0, X, tind=tind) -
self.J(0, 0, X, tind=tind) * \
self.J(2, 1, X, tind=tind)) / detDF
invDF[0, 2] = (-self.J(0, 2, X, tind=tind) * \
self.J(1, 1, X, tind=tind) +
self.J(0, 1, X, tind=tind) * \
self.J(1, 2, X, tind=tind)) / detDF
invDF[1, 2] = (self.J(0, 2, X, tind=tind) * \
self.J(1, 0, X, tind=tind) -
self.J(0, 0, X, tind=tind) * \
self.J(1, 2, X, tind=tind)) / detDF
invDF[2, 2] = (-self.J(0, 1, X, tind=tind) * \
self.J(1, 0, X, tind=tind) +
self.J(0, 0, X, tind=tind) * \
self.J(1, 1, X, tind=tind)) / detDF
invDF[0, 0] = -J[1][2] * J[2][1] + J[1][1] * J[2][2]
invDF[1, 0] = J[1][2] * J[2][0] - J[1][0] * J[2][2]
invDF[2, 0] = -J[1][1] * J[2][0] + J[1][0] * J[2][1]
invDF[0, 1] = J[0][2] * J[2][1] - J[0][1] * J[2][2]
invDF[1, 1] = -J[0][2] * J[2][0] + J[0][0] * J[2][2]
invDF[2, 1] = J[0][1] * J[2][0] - J[0][0] * J[2][1]
invDF[0, 2] = -J[0][2] * J[1][1] + J[0][1] * J[1][2]
invDF[1, 2] = J[0][2] * J[1][0] - J[0][0] * J[1][2]
invDF[2, 2] = -J[0][1] * J[1][0] + J[0][0] * J[1][1]
else:
raise Exception("Not implemented for the given dimension.")

return invDF
return invDF / detDF

def normals(self, X, tind, find, t2f):
if self.dim == 1:
Nref = np.array([[-1.0],
Expand Down
2 changes: 1 addition & 1 deletion tests/test_mapping.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ def runTest(self):
if itr == jtr:
self.assertTrue((fb.normals[jtr][case] == -1).all())
else:
self.assertTrue((fb.normals[jtr][case] == 0).all())
self.assertTrue((np.abs(fb.normals[jtr][case]) < 1e-14).all())


class TestIsoparamNormalsQuad(TestIsoparamNormals):
Expand Down

0 comments on commit e77096e

Please sign in to comment.