From da7fbbb461833bf4eb00f0ab2d1c0b7c4de2480f Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 9 Nov 2015 19:51:46 -0800 Subject: [PATCH] edgeCurl is O(h) not O(h^2) ???? --- SimPEG/Mesh/PointerTree.py | 302 ++++++++++++++++++++++--------- SimPEG/Tests.py | 27 +++ tests/mesh/test_TreeOperators.py | 103 +++++++++++ tests/mesh/test_pointerMesh.py | 61 ++++++- 4 files changed, 397 insertions(+), 96 deletions(-) create mode 100644 tests/mesh/test_TreeOperators.py diff --git a/SimPEG/Mesh/PointerTree.py b/SimPEG/Mesh/PointerTree.py index d3262fee..2c71bd75 100644 --- a/SimPEG/Mesh/PointerTree.py +++ b/SimPEG/Mesh/PointerTree.py @@ -1,8 +1,13 @@ from SimPEG import np, sp, Utils, Solver, Mesh import matplotlib.pyplot as plt import matplotlib +from mpl_toolkits.mplot3d import Axes3D +import matplotlib.colors as colors +import matplotlib.cm as cmx + import TreeUtils from SimPEG.Mesh.InnerProducts import InnerProducts +from SimPEG.Mesh.BaseMesh import BaseMesh import time MAX_BITS = 20 @@ -49,8 +54,8 @@ def SortGrid(grid, offset=0): class NotBalancedException(Exception): pass -class Tree(InnerProducts): - def __init__(self, h_in, levels=3): +class Tree(BaseMesh, InnerProducts): + def __init__(self, h_in, x0_in=None, levels=3): assert type(h_in) is list, 'h_in must be a list' assert len(h_in) > 1, "len(h_in) must be greater than 1" @@ -67,6 +72,24 @@ class Tree(InnerProducts): h[i] = h_i[:] # make a copy. self.h = h + x0 = np.zeros(len(h)) + if x0_in is not None: + assert len(h) == len(x0_in), "Dimension mismatch. x0 != len(h)" + for i in range(len(h)): + x_i, h_i = x0_in[i], h[i] + if Utils.isScalar(x_i): + x0[i] = x_i + elif x_i == '0': + x0[i] = 0.0 + elif x_i == 'C': + x0[i] = -h_i.sum()*0.5 + elif x_i == 'N': + x0[i] = -h_i.sum() + else: + raise Exception("x0[%i] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative." % i) + + BaseMesh.__init__(self, [len(_) for _ in h], x0) + self._levels = levels self._levelBits = int(np.ceil(np.sqrt(levels)))+1 @@ -295,11 +318,11 @@ class Tree(InnerProducts): return v in self._cells return self._index(v) in self._cells - def refine(self, function=None, recursive=True, cells=None, balance=True, _inRecursion=False): + def refine(self, function=None, recursive=True, cells=None, balance=True, verbose=False, _inRecursion=False): if not _inRecursion: self.__dirty__ = True - print 'Refining Mesh' + if verbose: print 'Refining Mesh' cells = cells if cells is not None else sorted(self._cells) recurse = [] @@ -310,7 +333,7 @@ class Tree(InnerProducts): if do: recurse += self._refineCell(cell) - print ' ', time.time() - tic + if verbose: print ' ', time.time() - tic if recursive and len(recurse) > 0: recurse += self.refine(function=function, recursive=True, cells=recurse, balance=balance, _inRecursion=True) @@ -457,12 +480,12 @@ class Tree(InnerProducts): return self._getNextCell(self._parentPointer(pointer), direction=direction, positive=positive) - def balance(self, recursive=True, cells=None, _inRecursion=False): + def balance(self, recursive=True, cells=None, verbose=False, _inRecursion=False): tic = time.time() if not _inRecursion: self.__dirty__ = True - print 'Balancing Mesh:' + if verbose: print 'Balancing Mesh:' cells = cells if cells is not None else sorted(self._cells) @@ -497,7 +520,7 @@ class Tree(InnerProducts): recurse.update([_ for _ in cs if type(_) in [int, long]]) # only add the bigger ones! recurse.update(newCells) - print ' ', len(cells), time.time() - tic + if verbose: print ' ', len(cells), time.time() - tic if recursive and len(recurse) > 0: self.balance(cells=sorted(recurse), _inRecursion=True) @@ -792,15 +815,28 @@ class Tree(InnerProducts): w = self._levelWidth(sl) if self.dim == 2: - self._hangingFx[self._fx2i[test ]] = ([self._fx2i[fx], 0.5], ) - self._hangingFx[self._fx2i[self._index([p[0] , p[1] + w, sl])]] = ([self._fx2i[fx], 0.5], ) + chy0 = self._cellH([p[0] , p[1] , sl])[1] + chy1 = self._cellH([p[0] , p[1] + w, sl])[1] + A = (chy0 + chy1) + + self._hangingFx[self._fx2i[test ]] = ([self._fx2i[fx], chy0 / A], ) + self._hangingFx[self._fx2i[self._index([p[0] , p[1] + w, sl])]] = ([self._fx2i[fx], chy1 / A], ) n0, n1 = fx, self._index([p[0], p[1] + 2*w, p[-1]]) self._hangingN[self._n2i[test ]] = ([self._n2i[n0], 1.0], ) - self._hangingN[self._n2i[self._index([p[0] , p[1] + w, sl])]] = ([self._n2i[n0], 0.5], [self._n2i[n1], 0.5]) + self._hangingN[self._n2i[self._index([p[0] , p[1] + w, sl])]] = ([self._n2i[n0], 1.0 - chy0 / A], [self._n2i[n1], 1.0 - chy1 / A]) self._hangingN[self._n2i[self._index([p[0] , p[1] + 2*w, sl])]] = ([self._n2i[n1], 1.0], ) elif self.dim == 3: + + chy0 = self._cellH([p[0] , p[1] , p[2] , sl])[1] + chy1 = self._cellH([p[0] , p[1] + w, p[2] , sl])[1] + chz0 = self._cellH([p[0] , p[1] , p[2] , sl])[2] + chz1 = self._cellH([p[0] , p[1] , p[2] + w, sl])[2] + lenY = chy0 + chy1 + lenZ = chz0 + chz1 + A = lenY * lenZ + ey0 = fx ey1 = self._index([p[0], p[1] , p[2] + 2*w, p[-1]]) ez0 = fx @@ -811,24 +847,38 @@ class Tree(InnerProducts): n2 = self._index([p[0], p[1] , p[2] + 2*w, p[-1]]) n3 = self._index([p[0], p[1] + 2*w, p[2] + 2*w, p[-1]]) - self._hangingFx[self._fx2i[test ]] = ([self._fx2i[fx], 0.25], ) - self._hangingFx[self._fx2i[self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._fx2i[fx], 0.25], ) - self._hangingFx[self._fx2i[self._index([p[0], p[1] , p[2] + w, sl])]] = ([self._fx2i[fx], 0.25], ) - self._hangingFx[self._fx2i[self._index([p[0], p[1] + w, p[2] + w, sl])]] = ([self._fx2i[fx], 0.25], ) + self._hangingFx[self._fx2i[test ]] = ([self._fx2i[fx], chy0*chz0 / A ], ) + self._hangingFx[self._fx2i[self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._fx2i[fx], chy1*chz0 / A ], ) + self._hangingFx[self._fx2i[self._index([p[0], p[1] , p[2] + w, sl])]] = ([self._fx2i[fx], chy0*chz1 / A ], ) + self._hangingFx[self._fx2i[self._index([p[0], p[1] + w, p[2] + w, sl])]] = ([self._fx2i[fx], chy1*chz1 / A ], ) - self._hangingEy[self._ey2i[test ]] = ([self._ey2i[ey0], 0.5], ) - self._hangingEy[self._ey2i[self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._ey2i[ey0], 0.5], ) - self._hangingEy[self._ey2i[self._index([p[0], p[1] , p[2] + w, sl])]] = ([self._ey2i[ey0], 0.25], [self._ey2i[ey1], 0.25]) - self._hangingEy[self._ey2i[self._index([p[0], p[1] + w, p[2] + w, sl])]] = ([self._ey2i[ey0], 0.25], [self._ey2i[ey1], 0.25]) - self._hangingEy[self._ey2i[self._index([p[0], p[1] , p[2] + 2*w, sl])]] = ([self._ey2i[ey1], 0.5], ) - self._hangingEy[self._ey2i[self._index([p[0], p[1] + w, p[2] + 2*w, sl])]] = ([self._ey2i[ey1], 0.5], ) + self._hangingEy[self._ey2i[test ]] = ([self._ey2i[ey0], 1.0], ) + self._hangingEy[self._ey2i[self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._ey2i[ey0], 1.0], ) + self._hangingEy[self._ey2i[self._index([p[0], p[1] , p[2] + w, sl])]] = ([self._ey2i[ey0], 0.5], [self._ey2i[ey1], 0.5]) + self._hangingEy[self._ey2i[self._index([p[0], p[1] + w, p[2] + w, sl])]] = ([self._ey2i[ey0], 0.5], [self._ey2i[ey1], 0.5]) + self._hangingEy[self._ey2i[self._index([p[0], p[1] , p[2] + 2*w, sl])]] = ([self._ey2i[ey1], 1.0], ) + self._hangingEy[self._ey2i[self._index([p[0], p[1] + w, p[2] + 2*w, sl])]] = ([self._ey2i[ey1], 1.0], ) - self._hangingEz[self._ez2i[test ]] = ([self._ez2i[ez0], 0.5], ) - self._hangingEz[self._ez2i[self._index([p[0], p[1] , p[2] + w, sl])]] = ([self._ez2i[ez0], 0.5], ) - self._hangingEz[self._ez2i[self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._ez2i[ez0], 0.25], [self._ez2i[ez1], 0.25]) - self._hangingEz[self._ez2i[self._index([p[0], p[1] + w, p[2] + w, sl])]] = ([self._ez2i[ez0], 0.25], [self._ez2i[ez1], 0.25]) - self._hangingEz[self._ez2i[self._index([p[0], p[1] + 2*w, p[2] , sl])]] = ([self._ez2i[ez1], 0.5], ) - self._hangingEz[self._ez2i[self._index([p[0], p[1] + 2*w, p[2] + w, sl])]] = ([self._ez2i[ez1], 0.5], ) + self._hangingEz[self._ez2i[test ]] = ([self._ez2i[ez0], 1.0], ) + self._hangingEz[self._ez2i[self._index([p[0], p[1] , p[2] + w, sl])]] = ([self._ez2i[ez0], 1.0], ) + self._hangingEz[self._ez2i[self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._ez2i[ez0], 0.5], [self._ez2i[ez1], 0.5]) + self._hangingEz[self._ez2i[self._index([p[0], p[1] + w, p[2] + w, sl])]] = ([self._ez2i[ez0], 0.5], [self._ez2i[ez1], 0.5]) + self._hangingEz[self._ez2i[self._index([p[0], p[1] + 2*w, p[2] , sl])]] = ([self._ez2i[ez1], 1.0], ) + self._hangingEz[self._ez2i[self._index([p[0], p[1] + 2*w, p[2] + w, sl])]] = ([self._ez2i[ez1], 1.0], ) + + # self._hangingEy[self._ey2i[test ]] = ([self._ey2i[ey0], chy0 / lenY], ) + # self._hangingEy[self._ey2i[self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._ey2i[ey0], chy1 / lenY], ) + # self._hangingEy[self._ey2i[self._index([p[0], p[1] , p[2] + w, sl])]] = ([self._ey2i[ey0], chy0 / lenY / 2.0], [self._ey2i[ey1], chy0 / lenY / 2.0]) + # self._hangingEy[self._ey2i[self._index([p[0], p[1] + w, p[2] + w, sl])]] = ([self._ey2i[ey0], chy1 / lenY / 2.0], [self._ey2i[ey1], chy1 / lenY / 2.0]) + # self._hangingEy[self._ey2i[self._index([p[0], p[1] , p[2] + 2*w, sl])]] = ([self._ey2i[ey1], chy0 / lenY], ) + # self._hangingEy[self._ey2i[self._index([p[0], p[1] + w, p[2] + 2*w, sl])]] = ([self._ey2i[ey1], chy1 / lenY], ) + + # self._hangingEz[self._ez2i[test ]] = ([self._ez2i[ez0], chz0 / lenZ], ) + # self._hangingEz[self._ez2i[self._index([p[0], p[1] , p[2] + w, sl])]] = ([self._ez2i[ez0], chz1 / lenZ], ) + # self._hangingEz[self._ez2i[self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._ez2i[ez0], chz0 / lenZ / 2.0], [self._ez2i[ez1], chz0 / lenZ / 2.0]) + # self._hangingEz[self._ez2i[self._index([p[0], p[1] + w, p[2] + w, sl])]] = ([self._ez2i[ez0], chz1 / lenZ / 2.0], [self._ez2i[ez1], chz1 / lenZ / 2.0]) + # self._hangingEz[self._ez2i[self._index([p[0], p[1] + 2*w, p[2] , sl])]] = ([self._ez2i[ez1], chz0 / lenZ], ) + # self._hangingEz[self._ez2i[self._index([p[0], p[1] + 2*w, p[2] + w, sl])]] = ([self._ez2i[ez1], chz1 / lenZ], ) self._hangingN[ self._n2i[ test ]] = ([self._n2i[n0], 1.0], ) self._hangingN[ self._n2i[ self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._n2i[n0], 0.5], [self._n2i[n1], 0.5]) @@ -852,8 +902,11 @@ class Tree(InnerProducts): w = self._levelWidth(sl) if self.dim == 2: - self._hangingFy[self._fy2i[test ]] = ([self._fy2i[fy], 0.5], ) - self._hangingFy[self._fy2i[self._index([p[0] + w, p[1] , sl])]] = ([self._fy2i[fy], 0.5], ) + chx0 = self._cellH([p[0] , p[1] , sl])[0] + chx1 = self._cellH([p[0] + w, p[1] , sl])[0] + + self._hangingFy[self._fy2i[test ]] = ([self._fy2i[fy], chx0 / (chx0 + chx1)], ) + self._hangingFy[self._fy2i[self._index([p[0] + w, p[1] , sl])]] = ([self._fy2i[fy], chx1 / (chx0 + chx1)], ) n0, n1 = fy, self._index([p[0] + 2*w, p[1], p[-1]]) self._hangingN[self._n2i[test ]] = ([self._n2i[n0], 1.0], ) @@ -861,6 +914,15 @@ class Tree(InnerProducts): self._hangingN[self._n2i[self._index([p[0] + 2*w, p[1] , sl])]] = ([self._n2i[n1], 1.0], ) elif self.dim == 3: + + chx0 = self._cellH([p[0] , p[1] , p[2] , sl])[0] + chx1 = self._cellH([p[0] + w, p[1] , p[2] , sl])[0] + chz0 = self._cellH([p[0] , p[1] , p[2] , sl])[2] + chz1 = self._cellH([p[0] , p[1] , p[2] + w, sl])[2] + lenX = chx0 + chx1 + lenZ = chz0 + chz1 + A = lenX * lenZ + ex0 = fy ex1 = self._index([p[0] , p[1], p[2] + 2*w, p[-1]]) ez0 = fy @@ -871,24 +933,38 @@ class Tree(InnerProducts): n2 = self._index([p[0] , p[1], p[2] + 2*w, p[-1]]) n3 = self._index([p[0] + 2*w, p[1], p[2] + 2*w, p[-1]]) - self._hangingFy[self._fy2i[test ]] = ([self._fy2i[fy], 0.25], ) - self._hangingFy[self._fy2i[self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._fy2i[fy], 0.25], ) - self._hangingFy[self._fy2i[self._index([p[0] , p[1], p[2] + w, sl])]] = ([self._fy2i[fy], 0.25], ) - self._hangingFy[self._fy2i[self._index([p[0] + w, p[1], p[2] + w, sl])]] = ([self._fy2i[fy], 0.25], ) + self._hangingFy[self._fy2i[test ]] = ([self._fy2i[fy], chx0*chz0 / A ], ) + self._hangingFy[self._fy2i[self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._fy2i[fy], chx1*chz0 / A ], ) + self._hangingFy[self._fy2i[self._index([p[0] , p[1], p[2] + w, sl])]] = ([self._fy2i[fy], chx0*chz1 / A ], ) + self._hangingFy[self._fy2i[self._index([p[0] + w, p[1], p[2] + w, sl])]] = ([self._fy2i[fy], chx1*chz1 / A ], ) - self._hangingEx[self._ex2i[test ]] = ([self._ex2i[ex0], 0.5], ) - self._hangingEx[self._ex2i[self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._ex2i[ex0], 0.5], ) - self._hangingEx[self._ex2i[self._index([p[0] , p[1], p[2] + w, sl])]] = ([self._ex2i[ex0], 0.25], [self._ex2i[ex1], 0.25]) - self._hangingEx[self._ex2i[self._index([p[0] + w, p[1], p[2] + w, sl])]] = ([self._ex2i[ex0], 0.25], [self._ex2i[ex1], 0.25]) - self._hangingEx[self._ex2i[self._index([p[0] , p[1], p[2] + 2*w, sl])]] = ([self._ex2i[ex1], 0.5], ) - self._hangingEx[self._ex2i[self._index([p[0] + w, p[1], p[2] + 2*w, sl])]] = ([self._ex2i[ex1], 0.5], ) + self._hangingEx[self._ex2i[test ]] = ([self._ex2i[ex0], 1.0], ) + self._hangingEx[self._ex2i[self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._ex2i[ex0], 1.0], ) + self._hangingEx[self._ex2i[self._index([p[0] , p[1], p[2] + w, sl])]] = ([self._ex2i[ex0], 0.5], [self._ex2i[ex1], 0.5]) + self._hangingEx[self._ex2i[self._index([p[0] + w, p[1], p[2] + w, sl])]] = ([self._ex2i[ex0], 0.5], [self._ex2i[ex1], 0.5]) + self._hangingEx[self._ex2i[self._index([p[0] , p[1], p[2] + 2*w, sl])]] = ([self._ex2i[ex1], 1.0], ) + self._hangingEx[self._ex2i[self._index([p[0] + w, p[1], p[2] + 2*w, sl])]] = ([self._ex2i[ex1], 1.0], ) - self._hangingEz[self._ez2i[test ]] = ([self._ez2i[ez0], 0.5], ) - self._hangingEz[self._ez2i[self._index([p[0] , p[1], p[2] + w, sl])]] = ([self._ez2i[ez0], 0.5], ) - self._hangingEz[self._ez2i[self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._ez2i[ez0], 0.25], [self._ez2i[ez1], 0.25]) - self._hangingEz[self._ez2i[self._index([p[0] + w, p[1], p[2] + w, sl])]] = ([self._ez2i[ez0], 0.25], [self._ez2i[ez1], 0.25]) - self._hangingEz[self._ez2i[self._index([p[0] + 2*w, p[1], p[2] , sl])]] = ([self._ez2i[ez1], 0.5], ) - self._hangingEz[self._ez2i[self._index([p[0] + 2*w, p[1], p[2] + w, sl])]] = ([self._ez2i[ez1], 0.5], ) + self._hangingEz[self._ez2i[test ]] = ([self._ez2i[ez0], 1.0], ) + self._hangingEz[self._ez2i[self._index([p[0] , p[1], p[2] + w, sl])]] = ([self._ez2i[ez0], 1.0], ) + self._hangingEz[self._ez2i[self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._ez2i[ez0], 0.5], [self._ez2i[ez1], 0.5]) + self._hangingEz[self._ez2i[self._index([p[0] + w, p[1], p[2] + w, sl])]] = ([self._ez2i[ez0], 0.5], [self._ez2i[ez1], 0.5]) + self._hangingEz[self._ez2i[self._index([p[0] + 2*w, p[1], p[2] , sl])]] = ([self._ez2i[ez1], 1.0], ) + self._hangingEz[self._ez2i[self._index([p[0] + 2*w, p[1], p[2] + w, sl])]] = ([self._ez2i[ez1], 1.0], ) + + # self._hangingEx[self._ex2i[test ]] = ([self._ex2i[ex0], chx0 / lenX], ) + # self._hangingEx[self._ex2i[self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._ex2i[ex0], chx1 / lenX], ) + # self._hangingEx[self._ex2i[self._index([p[0] , p[1], p[2] + w, sl])]] = ([self._ex2i[ex0], chx0 / lenX / 2.0], [self._ex2i[ex1], chx0 / lenX / 2.0]) + # self._hangingEx[self._ex2i[self._index([p[0] + w, p[1], p[2] + w, sl])]] = ([self._ex2i[ex0], chx1 / lenX / 2.0], [self._ex2i[ex1], chx1 / lenX / 2.0]) + # self._hangingEx[self._ex2i[self._index([p[0] , p[1], p[2] + 2*w, sl])]] = ([self._ex2i[ex1], chx0 / lenX], ) + # self._hangingEx[self._ex2i[self._index([p[0] + w, p[1], p[2] + 2*w, sl])]] = ([self._ex2i[ex1], chx1 / lenX], ) + + # self._hangingEz[self._ez2i[test ]] = ([self._ez2i[ez0], chz0 / lenZ], ) + # self._hangingEz[self._ez2i[self._index([p[0] , p[1], p[2] + w, sl])]] = ([self._ez2i[ez0], chz1 / lenZ], ) + # self._hangingEz[self._ez2i[self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._ez2i[ez0], chz0 / lenZ / 2.0], [self._ez2i[ez1], chz0 / lenZ / 2.0]) + # self._hangingEz[self._ez2i[self._index([p[0] + w, p[1], p[2] + w, sl])]] = ([self._ez2i[ez0], chz1 / lenZ / 2.0], [self._ez2i[ez1], chz1 / lenZ / 2.0]) + # self._hangingEz[self._ez2i[self._index([p[0] + 2*w, p[1], p[2] , sl])]] = ([self._ez2i[ez1], chz0 / lenZ], ) + # self._hangingEz[self._ez2i[self._index([p[0] + 2*w, p[1], p[2] + w, sl])]] = ([self._ez2i[ez1], chz1 / lenZ], ) self._hangingN[ self._n2i[ test ]] = ([self._n2i[n0], 1.0], ) self._hangingN[ self._n2i[ self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._n2i[n0], 0.5], [self._n2i[n1], 0.5]) @@ -915,6 +991,14 @@ class Tree(InnerProducts): continue w = self._levelWidth(sl) + chx0 = self._cellH([p[0] , p[1] , p[2] , sl])[0] + chx1 = self._cellH([p[0] + w, p[1] , p[2] , sl])[0] + chy0 = self._cellH([p[0] , p[1] , p[2] , sl])[1] + chy1 = self._cellH([p[0] , p[1] + w, p[2] , sl])[1] + lenX = chx0 + chx1 + lenY = chy0 + chy1 + A = lenX * lenY + ex0 = fz ex1 = self._index([p[0] , p[1] + 2*w, p[2], p[-1]]) ey0 = fz @@ -925,24 +1009,38 @@ class Tree(InnerProducts): n2 = self._index([p[0] , p[1] + 2*w, p[2], p[-1]]) n3 = self._index([p[0] + 2*w, p[1] + 2*w, p[2], p[-1]]) - self._hangingFz[self._fz2i[test ]] = ([self._fz2i[fz], 0.25], ) - self._hangingFz[self._fz2i[self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._fz2i[fz], 0.25], ) - self._hangingFz[self._fz2i[self._index([p[0] , p[1] + w, p[2], sl])]] = ([self._fz2i[fz], 0.25], ) - self._hangingFz[self._fz2i[self._index([p[0] + w, p[1] + w, p[2], sl])]] = ([self._fz2i[fz], 0.25], ) + self._hangingFz[self._fz2i[test ]] = ([self._fz2i[fz], chx0*chy0 / A ], ) + self._hangingFz[self._fz2i[self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._fz2i[fz], chx1*chy0 / A ], ) + self._hangingFz[self._fz2i[self._index([p[0] , p[1] + w, p[2], sl])]] = ([self._fz2i[fz], chx0*chy1 / A ], ) + self._hangingFz[self._fz2i[self._index([p[0] + w, p[1] + w, p[2], sl])]] = ([self._fz2i[fz], chx1*chy1 / A ], ) - self._hangingEx[self._ex2i[test ]] = ([self._ex2i[ex0], 0.5], ) - self._hangingEx[self._ex2i[self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._ex2i[ex0], 0.5], ) - self._hangingEx[self._ex2i[self._index([p[0] , p[1] + w, p[2], sl])]] = ([self._ex2i[ex0], 0.25], [self._ex2i[ex1], 0.25]) - self._hangingEx[self._ex2i[self._index([p[0] + w, p[1] + w, p[2], sl])]] = ([self._ex2i[ex0], 0.25], [self._ex2i[ex1], 0.25]) - self._hangingEx[self._ex2i[self._index([p[0] , p[1] + 2*w, p[2], sl])]] = ([self._ex2i[ex1], 0.5], ) - self._hangingEx[self._ex2i[self._index([p[0] + w, p[1] + 2*w, p[2], sl])]] = ([self._ex2i[ex1], 0.5], ) + self._hangingEx[self._ex2i[test ]] = ([self._ex2i[ex0], 1.0], ) + self._hangingEx[self._ex2i[self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._ex2i[ex0], 1.0], ) + self._hangingEx[self._ex2i[self._index([p[0] , p[1] + w, p[2], sl])]] = ([self._ex2i[ex0], 0.5], [self._ex2i[ex1], 0.5]) + self._hangingEx[self._ex2i[self._index([p[0] + w, p[1] + w, p[2], sl])]] = ([self._ex2i[ex0], 0.5], [self._ex2i[ex1], 0.5]) + self._hangingEx[self._ex2i[self._index([p[0] , p[1] + 2*w, p[2], sl])]] = ([self._ex2i[ex1], 1.0], ) + self._hangingEx[self._ex2i[self._index([p[0] + w, p[1] + 2*w, p[2], sl])]] = ([self._ex2i[ex1], 1.0], ) - self._hangingEy[self._ey2i[test ]] = ([self._ey2i[ey0], 0.5], ) - self._hangingEy[self._ey2i[self._index([p[0] , p[1] + w, p[2], sl])]] = ([self._ey2i[ey0], 0.5], ) - self._hangingEy[self._ey2i[self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._ey2i[ey0], 0.25], [self._ey2i[ey1], 0.25]) - self._hangingEy[self._ey2i[self._index([p[0] + w, p[1] + w, p[2], sl])]] = ([self._ey2i[ey0], 0.25], [self._ey2i[ey1], 0.25]) - self._hangingEy[self._ey2i[self._index([p[0] + 2*w, p[1] , p[2], sl])]] = ([self._ey2i[ey1], 0.5], ) - self._hangingEy[self._ey2i[self._index([p[0] + 2*w, p[1] + w, p[2], sl])]] = ([self._ey2i[ey1], 0.5], ) + self._hangingEy[self._ey2i[test ]] = ([self._ey2i[ey0], 1.0], ) + self._hangingEy[self._ey2i[self._index([p[0] , p[1] + w, p[2], sl])]] = ([self._ey2i[ey0], 1.0], ) + self._hangingEy[self._ey2i[self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._ey2i[ey0], 0.5], [self._ey2i[ey1], 0.5]) + self._hangingEy[self._ey2i[self._index([p[0] + w, p[1] + w, p[2], sl])]] = ([self._ey2i[ey0], 0.5], [self._ey2i[ey1], 0.5]) + self._hangingEy[self._ey2i[self._index([p[0] + 2*w, p[1] , p[2], sl])]] = ([self._ey2i[ey1], 1.0], ) + self._hangingEy[self._ey2i[self._index([p[0] + 2*w, p[1] + w, p[2], sl])]] = ([self._ey2i[ey1], 1.0], ) + + # self._hangingEx[self._ex2i[test ]] = ([self._ex2i[ex0], chx0 / lenX], ) + # self._hangingEx[self._ex2i[self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._ex2i[ex0], chx1 / lenX], ) + # self._hangingEx[self._ex2i[self._index([p[0] , p[1] + w, p[2], sl])]] = ([self._ex2i[ex0], chx0 / lenX / 2.0], [self._ex2i[ex1], chx0 / lenX / 2.0]) + # self._hangingEx[self._ex2i[self._index([p[0] + w, p[1] + w, p[2], sl])]] = ([self._ex2i[ex0], chx1 / lenX / 2.0], [self._ex2i[ex1], chx1 / lenX / 2.0]) + # self._hangingEx[self._ex2i[self._index([p[0] , p[1] + 2*w, p[2], sl])]] = ([self._ex2i[ex1], chx0 / lenX], ) + # self._hangingEx[self._ex2i[self._index([p[0] + w, p[1] + 2*w, p[2], sl])]] = ([self._ex2i[ex1], chx1 / lenX], ) + + # self._hangingEy[self._ey2i[test ]] = ([self._ey2i[ey0], chy0 / lenY], ) + # self._hangingEy[self._ey2i[self._index([p[0] , p[1] + w, p[2], sl])]] = ([self._ey2i[ey0], chy1 / lenY], ) + # self._hangingEy[self._ey2i[self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._ey2i[ey0], chy0 / lenY / 2.0], [self._ey2i[ey1], chy0 / lenY / 2.0]) + # self._hangingEy[self._ey2i[self._index([p[0] + w, p[1] + w, p[2], sl])]] = ([self._ey2i[ey0], chy1 / lenY / 2.0], [self._ey2i[ey1], chy1 / lenY / 2.0]) + # self._hangingEy[self._ey2i[self._index([p[0] + 2*w, p[1] , p[2], sl])]] = ([self._ey2i[ey1], chy0 / lenY], ) + # self._hangingEy[self._ey2i[self._index([p[0] + 2*w, p[1] + w, p[2], sl])]] = ([self._ey2i[ey1], chy1 / lenY], ) self._hangingN[ self._n2i[ test ]] = ([self._n2i[n0], 1.0], ) self._hangingN[ self._n2i[ self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._n2i[n0], 0.5], [self._n2i[n1], 0.5]) @@ -956,12 +1054,12 @@ class Tree(InnerProducts): self.__dirtyHanging__ = False - def number(self, force=False): + def number(self, balance=True, force=False): if not self.__dirty__ and not force: return - self.balance() + if balance: self.balance() self._hanging(force=force) - def _deflationMatrix(self, location, withHanging=True): + def _deflationMatrix(self, location, withHanging=True, asOnes=False): assert location in ['N','F','Fx','Fy'] + (['Fz','E','Ex','Ey','Ez'] if self.dim == 3 else []) args = dict() @@ -974,11 +1072,11 @@ class Tree(InnerProducts): args['Ey'] = (self._edgesY, self._hangingEy, self._ey2i) args['Ez'] = (self._edgesZ, self._hangingEz, self._ez2i) if location in ['F', 'E']: - Rlist = [self._deflationMatrix(location + subLoc, withHanging=withHanging) for subLoc in ['x','y','z'][:self.dim]] + Rlist = [self._deflationMatrix(location + subLoc, withHanging=withHanging, asOnes=asOnes) for subLoc in ['x','y','z'][:self.dim]] return sp.block_diag(Rlist) - return self.__deflationMatrix(*args[location], withHanging=withHanging) + return self.__deflationMatrix(*args[location], withHanging=withHanging, asOnes=asOnes) - def __deflationMatrix(self, theSet, theHang, theIndex, withHanging=True): + def __deflationMatrix(self, theSet, theHang, theIndex, withHanging=True, asOnes=False): reducedInd = dict() # final reduced index ii = 0 I,J,V = [],[],[] @@ -994,7 +1092,10 @@ class Tree(InnerProducts): hf = theHang[hfkey] I += [hfkey]*len(hf) J += [reducedInd[_[0]] for _ in hf] - V += [_[1] for _ in hf] + if asOnes: + V += [1.0]*len(hf) + else: + V += [_[1] for _ in hf] return sp.csr_matrix((V,(I,J)), shape=(len(theSet), len(reducedInd))) @property @@ -1037,10 +1138,13 @@ class Tree(InnerProducts): V += [pm] D = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntF)) - R = self._deflationMatrix('F') + R = self._deflationMatrix('F',asOnes=True) VOL = self.vol - S = self.area - self._faceDiv = Utils.sdiag(1.0/VOL)*D*R*Utils.sdiag(S) + if self.dim == 2: + S = np.r_[self._areaFxFull, self._areaFyFull] + elif self.dim == 3: + S = np.r_[self._areaFxFull, self._areaFyFull, self._areaFzFull] + self._faceDiv = Utils.sdiag(1.0/VOL)*D*Utils.sdiag(S)*R return self._faceDiv @property @@ -1072,7 +1176,7 @@ class Tree(InnerProducts): J += [edge + off] V += [pm] - faceOffset = self.nFx + faceOffset = self.ntFx offset = [0]*2 + [self.ntEx+self.ntEy]*2 PM = [-1, 1, 1, -1] for ii, fy in enumerate(sorted(self._facesY)): @@ -1092,7 +1196,7 @@ class Tree(InnerProducts): J += [edge + off] V += [pm] - faceOffset = self.nFx + self.nFy + faceOffset = self.ntFx + self.ntFy offset = [0]*2 + [self.ntEx]*2 PM = [1, -1, -1, 1] for ii, fz in enumerate(sorted(self._facesZ)): @@ -1112,15 +1216,19 @@ class Tree(InnerProducts): J += [edge + off] V += [pm] - Rf = self._deflationMatrix('F') + Rf = self._deflationMatrix('F', withHanging=False, asOnes=False) Re = self._deflationMatrix('E') - C = sp.csr_matrix((V,(I,J)), shape=(self.ntF, self.ntE)) - S = self.area - L = self.edge - self._edgeCurl = Utils.sdiag(1/S)*Rf.T*C*Re*Utils.sdiag(L) - return self._edgeCurl + Rf_ave = Utils.sdiag(1./Rf.sum(axis=0)) * Rf.T + # print Rf_ave + + C = sp.csr_matrix((V,(I,J)), shape=(self.ntF, self.ntE)) + S = np.r_[self._areaFxFull, self._areaFyFull, self._areaFzFull] + L = np.r_[self._edgeExFull, self._edgeEyFull, self._edgeEzFull] + # self._edgeCurl = Rf.T*Utils.sdiag(1.0/S)*C*Utils.sdiag(L)*Re + self._edgeCurl = Rf_ave*Utils.sdiag(1.0/S)*C*Utils.sdiag(L)*Re + return self._edgeCurl def _getFaceP(self, xFace, yFace, zFace): ind1, ind2, ind3 = [], [], [] @@ -1145,7 +1253,7 @@ class Tree(InnerProducts): PXXX = sp.coo_matrix((np.ones(self.dim*self.nC), (range(self.dim*self.nC), IND)), shape=(self.dim*self.nC, self.ntF)).tocsr() - Rf = self._deflationMatrix('F', withHanging=True) + Rf = self._deflationMatrix('F', withHanging=True, asOnes=True) return PXXX * Rf @@ -1169,9 +1277,9 @@ class Tree(InnerProducts): p = self._pointer(ind) w = self._levelWidth(p[-1]) - posX = [0,0] if xEdge == 'eX0' else [1, 0] if xEdge == 'eX1' else [0,1] if xEdge == 'eX2' else [1,1] - posY = [0,0] if yEdge == 'eY0' else [1, 0] if yEdge == 'eY1' else [0,1] if yEdge == 'eY2' else [1,1] - posZ = [0,0] if zEdge == 'eZ0' else [1, 0] if zEdge == 'eZ1' else [0,1] if zEdge == 'eZ2' else [1,1] + posX = [0,0] if xEdge == 'eX0' else [w, 0] if xEdge == 'eX1' else [0,w] if xEdge == 'eX2' else [w,w] + posY = [0,0] if yEdge == 'eY0' else [w, 0] if yEdge == 'eY1' else [0,w] if yEdge == 'eY2' else [w,w] + posZ = [0,0] if zEdge == 'eZ0' else [w, 0] if zEdge == 'eZ1' else [0,w] if zEdge == 'eZ2' else [w,w] ind1.append( self._ex2i[self._index([ p[0] , p[1] + posX[0], p[2] + posX[1], p[3]])] ) ind2.append( self._ey2i[self._index([ p[0] + posY[0], p[1] , p[2] + posY[1], p[3]])] + self.ntEx ) @@ -1181,9 +1289,9 @@ class Tree(InnerProducts): PXXX = sp.coo_matrix((np.ones(self.dim*self.nC), (range(self.dim*self.nC), IND)), shape=(self.dim*self.nC, self.ntE)).tocsr() - Rf = self._deflationMatrix('E', withHanging=False) + Re = self._deflationMatrix('E') - return PXXX * Rf + return PXXX * Re def _getEdgePxx(self): raise Exception('Not implemented') # this should be a reordering of the face inner product? @@ -1314,6 +1422,24 @@ class Tree(InnerProducts): if showIt:plt.show() + def plotImage(self, I, ax=None, showIt=True): + if self.dim == 3: raise Exception() + + if ax is None: ax = plt.subplot(111) + jet = cm = plt.get_cmap('jet') + cNorm = colors.Normalize(vmin=I.min(), vmax=I.max()) + scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet) + ax.set_xlim((self.x0[0], self.h[0].sum())) + ax.set_ylim((self.x0[1], self.h[1].sum())) + for ii, node in enumerate(self._sortedCells): + x0, sz = self._cellN(node), self._cellH(node) + ax.add_patch(plt.Rectangle((x0[0], x0[1]), sz[0], sz[1], facecolor=scalarMap.to_rgba(I[ii]), edgecolor='k')) + # if text: ax.text(self.center[0],self.center[1],self.num) + scalarMap._A = [] # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots + plt.colorbar(scalarMap) + if showIt: plt.show() + + if __name__ == '__main__': @@ -1332,9 +1458,9 @@ if __name__ == '__main__': else: return 0 - T = Tree([[(1,128)],[(1,128)],[(1,128)]],levels=7) + # T = Tree([[(1,128)],[(1,128)],[(1,128)]],levels=7) # T = Tree([128,128,128],levels=7) - # T = Tree([[(1,16)],[(1,16)]],levels=4) + T = Tree([[(1,16)],[(1,16)]],levels=4) # T = Tree([[(1,128)],[(1,128)]],levels=7) # T.refine(lambda xc:1, balance=False) # T._index([0,0,0]) @@ -1346,6 +1472,8 @@ if __name__ == '__main__': print time.time() - tic print T.nC + T.plotImage(np.random.rand(T.nC),showIt=True) + print T.getFaceInnerProduct() # print T.gridFz diff --git a/SimPEG/Tests.py b/SimPEG/Tests.py index c6aa5bc4..722d302b 100644 --- a/SimPEG/Tests.py +++ b/SimPEG/Tests.py @@ -4,6 +4,7 @@ from numpy.linalg import norm from SimPEG.Utils import mkvc, sdiag, diagEst from SimPEG import Utils from SimPEG.Mesh import TensorMesh, CurvilinearMesh, CylMesh +from SimPEG.Mesh.PointerTree import Tree import numpy as np import scipy.sparse as sp import unittest @@ -132,6 +133,32 @@ class OrderTest(unittest.TestCase): self.M = CurvilinearMesh([X, Y, Z]) return 1./nc + elif 'Tree' in self._meshType: + nc *= 2 + if 'uniform' in self._meshType: + h = [nc, nc, nc] + elif 'random' in self._meshType: + h1 = np.random.rand(nc)*nc*0.5 + nc*0.5 + h2 = np.random.rand(nc)*nc*0.5 + nc*0.5 + h3 = np.random.rand(nc)*nc*0.5 + nc*0.5 + h = [hi/np.sum(hi) for hi in [h1, h2, h3]] # normalize + else: + raise Exception('Unexpected meshType') + + levels = int(np.log(nc)/np.log(2)) + self.M = Tree(h[:self.meshDimension], levels=levels) + def function(xc): + r = xc - np.array([0.5]*len(xc)) + dist = np.sqrt(r.dot(r)) + if dist < 0.3: + return levels + return levels-1 + self.M.refine(function,balance=False) + self.M.number(balance=False) + # self.M.plotGrid(showIt=True) + max_h = max([np.max(hi) for hi in self.M.h]) + return max_h + def getError(self): """For given h, generate A[h], f and A(f) and return norm of error.""" return 1. diff --git a/tests/mesh/test_TreeOperators.py b/tests/mesh/test_TreeOperators.py new file mode 100644 index 00000000..9babe1a2 --- /dev/null +++ b/tests/mesh/test_TreeOperators.py @@ -0,0 +1,103 @@ +import numpy as np +import unittest +from SimPEG.Tests import OrderTest +import matplotlib.pyplot as plt + +#TODO: 'randomTensorMesh' +MESHTYPES = ['uniformTree'] #['randomTree', 'uniformTree'] +call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) +call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) +cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] +cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)] +cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy))) +cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey))) +cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz))) +cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez))) + + +class TestFaceDiv2D(OrderTest): + name = "Face Divergence 2D" + meshTypes = MESHTYPES + meshDimension = 2 + meshSizes = [16, 32] + + def getError(self): + #Test function + fx = lambda x, y: np.sin(2*np.pi*x) + fy = lambda x, y: np.sin(2*np.pi*y) + sol = lambda x, y: 2*np.pi*(np.cos(2*np.pi*x)+np.cos(2*np.pi*y)) + + Fc = cartF2(self.M, fx, fy) + F = self.M.projectFaceVector(Fc) + + divF = self.M.faceDiv.dot(F) + divF_ana = call2(sol, self.M.gridCC) + + err = np.linalg.norm((divF-divF_ana), np.inf) + + # self.M.plotImage(divF-divF_ana, showIt=True) + + return err + + def test_order(self): + self.orderTest() + +class TestFaceDiv3D(OrderTest): + name = "Face Divergence 3D" + meshTypes = MESHTYPES + meshSizes = [8, 16] + + def getError(self): + #Test function + fx = lambda x, y, z: np.sin(2*np.pi*x) + fy = lambda x, y, z: np.sin(2*np.pi*y) + fz = lambda x, y, z: np.sin(2*np.pi*z) + sol = lambda x, y, z: (2*np.pi*np.cos(2*np.pi*x)+2*np.pi*np.cos(2*np.pi*y)+2*np.pi*np.cos(2*np.pi*z)) + + Fc = cartF3(self.M, fx, fy, fz) + F = self.M.projectFaceVector(Fc) + + divF = self.M.faceDiv.dot(F) + divF_ana = call3(sol, self.M.gridCC) + + return np.linalg.norm((divF-divF_ana), np.inf) + + + def test_order(self): + self.orderTest() + + +class TestCurl(OrderTest): + name = "Curl" + meshTypes = MESHTYPES + meshSizes = [4, 8, 16, 32] + + def getError(self): + # fun: i (cos(y)) + j (cos(z)) + k (cos(x)) + # sol: i (sin(z)) + j (sin(x)) + k (sin(y)) + + funX = lambda x, y, z: np.cos(2*np.pi*y) + funY = lambda x, y, z: np.cos(2*np.pi*z) + funZ = lambda x, y, z: np.cos(2*np.pi*x) + + solX = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*z) + solY = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x) + solZ = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*y) + + Ec = cartE3(self.M, funX, funY, funZ) + E = self.M.projectEdgeVector(Ec) + + Fc = cartF3(self.M, solX, solY, solZ) + curlE_ana = self.M.projectFaceVector(Fc) + + curlE = self.M.edgeCurl.dot(E) + + err = np.linalg.norm((curlE - curlE_ana), np.inf) + + return err + + def test_order(self): + self.orderTest() + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mesh/test_pointerMesh.py b/tests/mesh/test_pointerMesh.py index c609f372..83672fe7 100644 --- a/tests/mesh/test_pointerMesh.py +++ b/tests/mesh/test_pointerMesh.py @@ -11,16 +11,21 @@ TOL = 1e-10 class TestSimpleQuadTree(unittest.TestCase): def test_counts(self): - - M = Tree([8,8]) + nc = 8 + h1 = np.random.rand(nc)*nc*0.5 + nc*0.5 + h2 = np.random.rand(nc)*nc*0.5 + nc*0.5 + h = [hi/np.sum(hi) for hi in [h1, h2]] # normalize + M = Tree(h) M._refineCell([0,0,0]) M._refineCell([0,0,1]) M.number() # M.plotGrid(showIt=True) - # assert sorted(M._cells) == [2, 34, 66, 99, 107, 115, 123, 129, 257, 386, 418, 450, 482] assert M.nhFx == 2 assert M.nFx == 9 - assert M.vol.sum() == 1.0 + + assert np.allclose(M.vol.sum(), 1.0) + + assert np.allclose(np.r_[M._areaFxFull, M._areaFyFull], M._deflationMatrix('F') * M.area) # def test_connectivity(self): @@ -48,10 +53,7 @@ class TestSimpleQuadTree(unittest.TestCase): # assert T._getNextCell([0,2,2], direction=1) == T._index([0,4,1]) # assert T._getNextCell([0,4,1], direction=1, positive=False) == [T._index([0,2,2]), [T._index([2,3,3]), T._index([3,3,3])]] - -class TestOperatorsQuadTree(unittest.TestCase): - - def test_counts(self): + def test_faceDiv(self): hx, hy = np.r_[1.,2,3,4], np.r_[5.,6,7,8] T = Tree([hx, hy], levels=2) @@ -76,7 +78,31 @@ class TestOperatorsQuadTree(unittest.TestCase): assert (M.faceDiv - T.permuteCC*T.faceDiv*T.permuteF.T).nnz == 0 -class TestOperatorsOcTree(unittest.TestCase): +class TestOcTree(unittest.TestCase): + + def test_counts(self): + nc = 8 + h1 = np.random.rand(nc)*nc*0.5 + nc*0.5 + h2 = np.random.rand(nc)*nc*0.5 + nc*0.5 + h3 = np.random.rand(nc)*nc*0.5 + nc*0.5 + h = [hi/np.sum(hi) for hi in [h1, h2, h3]] # normalize + M = Tree(h, levels=3) + M._refineCell([0,0,0,0]) + M._refineCell([0,0,0,1]) + M.number() + # M.plotGrid(showIt=True) + # assert M.nhFx == 2 + # assert M.nFx == 9 + + assert np.allclose(M.vol.sum(), 1.0) + + # assert np.allclose(M._areaFxFull, (M._deflationMatrix('F') * M.area)[:M.ntFx]) + # assert np.allclose(M._areaFyFull, (M._deflationMatrix('F') * M.area)[M.ntFx:(M.ntFx+M.ntFy)]) + # assert np.allclose(M._areaFzFull, (M._deflationMatrix('F') * M.area)[(M.ntFx+M.ntFy):]) + + # assert np.allclose(M._edgeExFull, (M._deflationMatrix('E') * M.edge)[:M.ntEx]) + # assert np.allclose(M._edgeEyFull, (M._deflationMatrix('E') * M.edge)[M.ntEx:(M.ntEx+M.ntEy)]) + # assert np.allclose(M._edgeEzFull, (M._deflationMatrix('E') * M.edge)[(M.ntEx+M.ntEy):]) def test_faceDiv(self): @@ -136,6 +162,23 @@ class TestOperatorsOcTree(unittest.TestCase): assert np.allclose(Mr.getFaceInnerProduct().todense(), (M.permuteF * M.getFaceInnerProduct() * M.permuteF.T).todense()) assert np.allclose(Mr.getEdgeInnerProduct().todense(), (M.permuteE * M.getEdgeInnerProduct() * M.permuteE.T).todense()) + def test_VectorIdenties(self): + hx, hy, hz = [[(1,4)], [(1,4)], [(1,4)]] + + M = Tree([hx, hy, hz], levels=2) + Mr = Mesh.TensorMesh([hx, hy, hz]) + + assert (M.faceDiv * M.edgeCurl).nnz == 0 + assert (Mr.faceDiv * Mr.edgeCurl).nnz == 0 + + hx, hy, hz = np.r_[1.,2,3,4], np.r_[5.,6,7,8], np.r_[9.,10,11,12] + + M = Tree([hx, hy, hz], levels=2) + Mr = Mesh.TensorMesh([hx, hy, hz]) + + assert np.max(np.abs((M.faceDiv * M.edgeCurl).todense().flatten())) < TOL + assert np.max(np.abs((Mr.faceDiv * Mr.edgeCurl).todense().flatten())) < TOL + if __name__ == '__main__':