diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..5c624b4b86fcdb92f8b04bd56ab533f6996d4191
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,95 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+env/
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# PyInstaller
+# Usually these files are written by a python script from a template
+# before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*,cover
+.hypothesis/
+
+# pyc files:
+*.pyc
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+
+# Flask instance folder
+instance/
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# IPython Notebook
+.ipynb_checkpoints
+
+# pyenv
+.python-version
+
+# celery beat schedule file
+celerybeat-schedule
+
+# dotenv
+.env
+
+# virtualenv
+venv/
+ENV/
+
+# Spyder project settings
+.spyderproject
+
+# vim files
+*.sw*
+
+# folder
+tmp/
+.idea
diff --git a/tfields/mesh3D.py b/tfields/mesh3D.py
new file mode 100644
index 0000000000000000000000000000000000000000..39112538c1e53ee2a64200075b3cc341367ba03b
--- /dev/null
+++ b/tfields/mesh3D.py
@@ -0,0 +1,1708 @@
+#!/usr/bin/env
+# encoding: utf-8
+"""
+Author: Daniel Boeckenhoff
+Mail: daniel.boeckenhoff@ipp.mpg.de
+
+part of tfields library
+"""
+import numpy as np
+import os
+import sympy
+import warnings
+import tfields
+import ioTools
+import mplTools
+import decoTools
+import pyTools
+import symTools
+from sympy.abc import y, z
+from scipy.spatial import ConvexHull
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+import loggingTools
+import cuttingTree
+
+logger = loggingTools.Logger(__name__)
+
+
+class Mesh3D(tfields.TensorMaps):
+ # pylint: disable=R0904
+ """
+ Points3D child used as vertices combined with faces to build a geometrical mesh of triangles
+ Examples:
+ >>> m = tfields.Mesh3D([[1,2,3], [3,3,3], [0,0,0], [5,6,7]], faces=[[0, 1, 2], [1, 2, 3]])
+ >>> m
+ Mesh3D([[ 1., 2., 3.],
+ [ 3., 3., 3.],
+ [ 0., 0., 0.],
+ [ 5., 6., 7.]])
+ >>> m.faces
+ array([[0, 1, 2],
+ [1, 2, 3]])
+
+ conversion to points only
+ >>> tfields.Points3D(m)
+ Points3D([[ 1., 2., 3.],
+ [ 3., 3., 3.],
+ [ 0., 0., 0.],
+ [ 5., 6., 7.]])
+
+ Empty instances
+ >>> m = tfields.Mesh3D([]);
+
+ going from Mesh3D to Triangles3D instance is easy and will be cached.
+ >>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], faces=[[0, 1, 2]]);
+ >>> m.triangles
+ Triangles3D([[ 1., 0., 0.],
+ [ 0., 1., 0.],
+ [ 0., 0., 0.]])
+
+ a list of scalars is assigned to each face
+ >>> m.faceScalars
+ array([], shape=(1, 0), dtype=float64)
+ >>> mScalar = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], [[0, 1, 2]], faceScalars=[.5]);
+ >>> mScalar.faceScalars
+ array([[ 0.5]])
+
+ adding together two meshes:
+ >>> m2 = tfields.Mesh3D([[1,0,0],[2,0,0],[0,3,0]], [[0,1,2]], faceScalars=[.7])
+ >>> msum = tfields.Mesh3D([mScalar, m2])
+ >>> msum
+ Mesh3D([[ 1., 0., 0.],
+ [ 0., 1., 0.],
+ [ 0., 0., 0.],
+ [ 1., 0., 0.],
+ [ 2., 0., 0.],
+ [ 0., 3., 0.]])
+ >>> msum.faces
+ array([[0, 1, 2],
+ [3, 4, 5]])
+
+ Saving and reading
+ >>> from tempfile import NamedTemporaryFile
+ >>> outFile = NamedTemporaryFile(suffix='.npz')
+ >>> m.savez(outFile.name)
+ >>> _ = outFile.seek(0)
+ >>> m1 = tfields.Mesh3D.createFromFile(outFile.name)
+ >>> bool(np.all(m == m1))
+ True
+ >>> m1.faces
+ array([[0, 1, 2]])
+
+ """
+
+ __slots__ = ['faces', 'faceScalars', 'coordSys', '_cache']
+
+ def __new__(cls, tensors, **kwargs):
+ if not issubclass(type(tensors), Mesh3D):
+ kwargs['dim'] = 3
+ return super(Mesh3D, cls).__new__(cls, tensors, **kwargs)
+
+ @classmethod
+ def _updateSlotKwargs(cls, kwargs, skipCache=True):
+ faces = kwargs.pop('faces', None)
+ faceScalars = kwargs.pop('faceScalars', None)
+
+ # Add faces
+ if faces is None or len(faces) == 0:
+ faces = np.empty((0, 3), dtype=int)
+ faces = np.array(faces, dtype=int)
+ kwargs['faces'] = faces
+
+ # Add faceScalars
+ if faceScalars is None or len(faceScalars) == 0:
+ faceScalars = np.empty((len(faces), 0), dtype=float)
+ if not len(faceScalars) == len(faces):
+ raise ValueError("Length of faceScalars has to be the same as faces lenght "
+ "({0}) but is {1}.".format(len(faces), len(faceScalars)))
+ # demand 2d structure of faceScalars
+ faceScalars = np.array(faceScalars, dtype=float)
+ if len(faceScalars.shape) == 1:
+ faceScalars = faceScalars.reshape(faces.shape[0], -1)
+ kwargs['faceScalars'] = faceScalars
+
+ super(Mesh3D, cls)._updateSlotKwargs(kwargs, skipCache=skipCache)
+
+ @classmethod
+ def createFromObjFile(cls, filePath, *groupNames):
+ """
+ Factory method
+ Given a filePath to a obj/wavefront file, construct the object
+ """
+ ioCls = ioTools.ObjFile
+ with ioCls(filePath, 'r') as f:
+ f.process()
+ vertices, faces = f.getVerticesFaces(*groupNames, firstFace=0)
+
+ log = logger.new()
+ if len(vertices) == 0:
+ return cls([])
+ faceLenghts = [len(face) for face in faces]
+ for i in reversed(range(len(faceLenghts))):
+ length = faceLenghts[i]
+ if length == 3:
+ continue
+ if length == 4:
+ log.warning("Given a Rectangle. I will split it but "
+ "sometimes the order is different.")
+ faces.insert(i + 1, faces[i][2:] + faces[i][:1])
+ faces[i] = faces[i][:3]
+ else:
+ raise NotImplementedError()
+ mesh = cls(vertices, faces=faces)
+ if groupNames:
+ mesh = mesh.clean()
+ return mesh
+
+ @classmethod
+ def createFromInpFile(cls, filePath, **kwargs):
+ """
+ Factory method
+ Given a filePath to a inp file, construct the object
+ """
+ import transcoding as tc
+ transcoding = tc.getTranscoding('inp')
+ content = transcoding.read(filePath)
+ part = content['parts'][0]
+ vertices = np.array([part['x'], part['y'], part['z']]).T / 1000
+ indices = np.array(part['nodeIndex']) - 1
+ if not list(indices) == range(len(indices)):
+ raise ValueError("node index skipped")
+ faces = np.array([part['nodeIndex{i}'.format(i=i)] for i in range(3)]).T - 1
+ return cls(vertices, faces=faces)
+
+ @classmethod
+ def createMeshGrid(cls, *baseVectors, **kwargs):
+ if not len(baseVectors) == 3:
+ raise AttributeError("3 baseVectors vectors required")
+
+ indices = [0, -1]
+ coords = range(3)
+ baseLengthsAbove1 = [len(b) > 1 for b in baseVectors]
+ # if one plane is given: rearrange indices and coords
+ if not all(baseLengthsAbove1):
+ indices = [0]
+ for i, b in enumerate(baseLengthsAbove1):
+ if not b:
+ coords = [i]
+ break
+
+ baseVectors = list(baseVectors)
+ mParts = []
+ for ind in indices:
+ for coord in coords:
+ basePart = baseVectors[:]
+ basePart[coord] = np.array([baseVectors[coord][ind]],
+ dtype=float)
+
+ mParts.append(cls.createMeshPlane(*basePart))
+ inst = cls.__new__(cls, mParts, **kwargs)
+ return inst
+
+ @classmethod
+ def createMeshPlane(cls, *baseVectors, **kwargs):
+ points = tfields.Points3D.createMeshGrid(*baseVectors)
+ fixCoord = None
+ for coord in range(3):
+ if len(baseVectors[coord]) > 1:
+ continue
+ if len(baseVectors[coord]) == 0:
+ continue
+ fixCoord = coord
+
+ variableCoords = list(range(3))
+ variableCoords.pop(variableCoords.index(fixCoord))
+
+ faces = []
+ base0, base1 = baseVectors[variableCoords[0]], baseVectors[variableCoords[1]]
+ for i1 in range(len(base1) - 1):
+ for i0 in range(len(base0) - 1):
+ pointIdxTopLeft = len(base1) * (i0 + 0) + (i1 + 0)
+ pointIdxTopRight = len(base1) * (i0 + 0) + (i1 + 1)
+ pointIdxBotLeft = len(base1) * (i0 + 1) + (i1 + 0)
+ pointIdxBotRight = len(base1) * (i0 + 1) + (i1 + 1)
+ faces.append([pointIdxTopLeft, pointIdxTopRight, pointIdxBotLeft])
+ faces.append([pointIdxTopRight, pointIdxBotLeft, pointIdxBotRight])
+ inst = cls.__new__(cls, points, faces=faces, **kwargs)
+ return inst
+
+ @decoTools.cached_property()
+ def triangles(self):
+ """
+ with the decorator, this should be handled like an attribute though it is a function
+
+ """
+ if self.faces.size == 0:
+ return tfields.Triangles3D([])
+ return tfields.Triangles3D(self[self.faces.flatten()], triangleScalars=self.faceScalars)
+
+ @decoTools.cached_property()
+ def planes(self):
+ if self.faces.size == 0:
+ return tfields.Planes3D([])
+ return tfields.Planes3D(self.getCentroids(), self.getNormVectors())
+
+ def saveTxt(self, filePath):
+ if self.coordSys != self.CARTESIAN:
+ cpy = self.copy()
+ cpy.coordinateTransform(self.CARTESIAN)
+ else:
+ cpy = self
+ with ioTools.TextFile(filePath, 'w') as f:
+ matrix = []
+ for i, face in enumerate(self.faces):
+ matrix.append(self.faceScalars[i, :])
+ matrix.extend(self[face])
+ f.writeMatrix(matrix, seperator=' ', lineBreak='\n')
+
+ @classmethod
+ def createFromTxtFile(cls, filePath):
+ return tfields.Triangles3D.createFromTxtFile(filePath).getMesh3D()
+
+ def __getattr__(self, name):
+ """
+ getter methods are forwarded to self.triangles
+ Examples:
+ >>> m = tfields.Mesh3D([]);
+ >>> m.getAreas
+ <bound method Triangles3D.getAreas of Triangles3D([], shape=(0, 3), dtype=float64)>
+
+ """
+ if name.startswith('get'):
+ if not hasattr(tfields.Triangles3D, name) or name == 'getMask':
+ raise AttributeError("Could not forward attribute {0}".format(name))
+ else:
+ return getattr(self.triangles, name)
+ else:
+ raise AttributeError("No attribute with name {0}".format(name))
+
+ def getNFaces(self):
+ return self.faces.shape[0]
+
+ def getMean(self, *args, **kwargs):
+ """
+ Forward this manually since getMean is derived already.
+ """
+ return self.triangles.getMean(*args, **kwargs)
+
+ def getStd(self, *args, **kwargs):
+ """
+ Forward this manually since getStd is derived already.
+ """
+ return self.triangles.getStd(*args, **kwargs)
+
+ def getScalars(self):
+ return self.faceScalars
+
+ def getScalarArrays(self):
+ return self.faceScalars.T
+
+ def getScalarDepth(self):
+ return self.faceScalars.shape[1]
+
+ def setScalarArray(self, scalarIndex, scalarArray):
+ """
+ >>> m = tfields.Mesh3D([[1,2,3], [3,3,3], [0,0,0], [5,6,7]], [[0, 1, 2], [1, 2, 3]],
+ ... faceScalars=[[1,2,3,4,5], [6,7,8,9,0]])
+ >>> m.setScalarArray(1, [42, 84])
+ >>> m.faceScalars
+ array([[ 1., 42., 3., 4., 5.],
+ [ 6., 84., 8., 9., 0.]])
+
+ """
+ if scalarIndex == self.getScalarDepth():
+ self.appendScalars(scalarArray)
+ else:
+ self.faceScalars[:, scalarIndex] = scalarArray
+ # try:
+ # self.faceScalars[:, scalarIndex] = scalarArray
+ # except:
+ # tfields.saveNested("~/tmp/bug1e-6.nest.npz", (self, scalarArray,
+ # scalarIndex))
+ # raise
+
+ def removeScalars(self, scalarIndex=None):
+ if scalarIndex is not None:
+ raise NotImplementedError()
+ self.faceScalars = np.empty(self.faceScalars.shape[0],
+ dtype=self.faceScalars.dtype)
+
+ def appendScalars(self, scalars):
+ """
+ Similar to list.append but in axis 1
+ """
+ self.faceScalars = np.concatenate([self.faceScalars, np.array([scalars]).T], 1)
+
+ def stackScalars(self, *stack, **kwargs):
+ """
+ add all faceScalars of stack meshes to self.faceScalars.
+ Args:
+ *stack (Mesh3D): input meshes are required to have the
+ same or simiar faces. This is not tested for time reasons
+ though.
+ **kwargs:
+ mapping (str):
+ 'centroid': Use centroids for check
+ which face belongs to which.
+ 'order': Assume all faces are in the same order.
+ Examples:
+ >>> m = tfields.Mesh3D([[1,0,0],[0,1,0],[0,0,1], [0,0,0]],
+ ... [[0,1,2], [2,3,0], [2,3,1]],
+ ... faceScalars=[1,2,3])
+ >>> c = m.copy()
+
+ If you exactly know, the two meshes are the same:
+ >>> m.stackScalars(c, mapping='order')
+ >>> all(m.getScalarArrays()[0, :] == [2.,4.,6.])
+ True
+
+ Using the default compares centroids
+ >>> m = c.copy()
+ >>> m.stackScalars(c, mapping='centroid')
+ >>> all(m.getScalarArrays()[0, :] == [2.,4.,6.])
+ True
+
+ """
+ mapping = kwargs.pop('mapping', 'centroid')
+ for stackObj in stack:
+ """
+ faceMap is a list of tuples: first meshFaceIndex second
+ selfFaceIndex
+ """
+ if mapping == 'order':
+ if not isinstance(stackObj, Mesh3D):
+ raise NotImplementedError()
+ faceRange = range(stackObj.getNFaces())
+ faceMap = [(i, i) for i in faceRange]
+ objScalars = stackObj.faceScalars
+ elif mapping == 'centroid':
+ if isinstance(stackObj, Mesh3D):
+ faceRange = range(stackObj.getNFaces())
+ centroids = stackObj.getCentroids()
+ objScalars = stackObj.faceScalars
+ elif isinstance(stackObj, tfields.ScalarField3D):
+ faceRange = range(stackObj.getNPoints())
+ centroids = stackObj
+ if mapping != 'centroid':
+ raise NotImplementedError()
+ objScalars = stackObj.scalars
+ else:
+ raise NotImplementedError()
+ closestCentroidIndices = \
+ centroids.closestPoints(self.getCentroids())
+ faceMap = [(i, j) for j, i in zip(faceRange, closestCentroidIndices)]
+ else:
+ raise NotImplementedError()
+ faceMap = np.array(faceMap)
+
+ """
+ Rearrange scalars according to faceMap
+ This is the part that takes longest
+ """
+ scalars = np.full(self.faceScalars.shape,
+ 0.,
+ dtype=objScalars.dtype)
+ for i in set(faceMap[:, 0]):
+ scalars[i] = objScalars[faceMap[faceMap[:, 0] == i,
+ 1]].sum(axis=0)
+
+ """
+ Add all scalars
+ """
+ self.faceScalars = self.faceScalars + scalars
+
+ def toOneSegment(self, mirrorZ=True):
+ """
+ Map the points to the first segment and mirror to positive z
+ if mirrorZOption is True. Special w7x method
+ Examples:
+ Build a mesh in three segments
+ >>> scalars = np.array([[1,2], [3,4], [5,6], [7,8]], dtype=float)
+ >>> m = tfields.Mesh3D([[6,-1,1], [6,0,1], [6,1,1],
+ ... [6,-1,0.5], [6,0,0.5], [6,1,0.5]],
+ ... [[0, 3, 4], [0, 1, 4], [1,4,5], [1,2,5]],
+ ... faceScalars=scalars)
+ >>> c = m.copy()
+ >>> c.coordinateTransform(c.CYLINDER)
+ >>> c[:, 1] *= -1
+ >>> c[:, 2] *= -1
+
+ >>> m = tfields.Mesh3D([m, c])
+ >>> m2 = m.copy()
+ >>> m2.toSegment(1)
+ >>> m3 = m.copy()
+ >>> m3.toSegment(2)
+ >>> mAll = tfields.Mesh3D([m, m2, m3])
+ >>> mAll.toOneSegment()
+ >>> bool((mAll.faceScalars == scalars * 6).all())
+ True
+
+ """
+ with self.tempCoordSys(self.CYLINDER):
+ dropCut = (-2 * np.pi / 10 < y) & (y < 2 * np.pi / 10)
+ if mirrorZ:
+ dropCut = dropCut & (z > 0)
+ dropCutMask = self.getMask(dropCut)
+ faceKeepMask = self.getFaceMask(dropCutMask)
+ excludedMesh = self.copy()
+ self.removeFaces(~faceKeepMask)
+ excludedMesh.removeFaces(faceKeepMask)
+ # remove 0 faces to be faster
+ from sympy.abc import s
+ zeroMask = tfields.getMask(excludedMesh.faceScalars,
+ cutExpression=(s == 0),
+ coords=[s] * excludedMesh.getScalarDepth())
+ excludedMesh.removeFaces(zeroMask)
+
+ # to one segment
+ super(Mesh3D, self).toOneSegment(mirrorZ=mirrorZ)
+ centroidSF = tfields.ScalarField3D(excludedMesh.getCentroids(),
+ excludedMesh.faceScalars)
+ centroidSF.toOneSegment(mirrorZ=mirrorZ)
+
+ self.stackScalars(centroidSF)
+
+ def pickScalars(self, *scalarIndices):
+ """
+ Reduces the faceScalars to only the indices given.
+ Examples:
+ >>> m = tfields.Mesh3D([[1,2,3], [3,3,3], [0,0,0], [5,6,7]],
+ ... [[0, 1, 2], [1, 2, 3]],
+ ... faceScalars=[[1,2,3,4,5], [6,7,8,9,0]])
+ >>> m.pickScalars(1, 3, 4)
+ >>> m.faceScalars
+ array([[ 2., 4., 5.],
+ [ 7., 9., 0.]])
+
+ """
+ self.faceScalars = self.faceScalars[:, list(scalarIndices)]
+
+ def pointsInMesh(self, points, delta, method='baryc', assignMultiple=False):
+ """
+ Check whether points lie within triangles with Barycentric Technique
+ """
+ masks = self.triangles.pointsInTriangles(points, delta, method='baryc',
+ assignMultiple=assignMultiple)
+ return masks
+
+ def convertNaN(self, value=0.):
+ super(Mesh3D, self).convertNaN(value)
+ nanIndicesScalars = np.isnan(self.faceScalars)
+ self.faceScalars[nanIndicesScalars] = value
+
+ def cutScalars(self, cutExpression, coords=None,
+ replaceValue=np.nan, scalarIndex=None, inplace=False):
+ """
+ Set a threshold to the scalars.
+ Args:
+ cutExpression (sympy cut expression or list of those):
+ threshold(sympy cut expression): cut scalars globaly
+ threshold(list of sympy cut expressions): set on threshold for every scalar array
+ Examples:
+ >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [0,1,0], [0,0,1]],
+ ... faces=[[0,1,2], [0,1,3]],
+ ... faceScalars=[[1, 1], [2, 2]])
+
+ Cuting all scalars at once
+ >>> from sympy.abc import s
+ >>> m.cutScalars(s <= 1., replaceValue=0.).faceScalars
+ array([[ 0., 0.],
+ [ 2., 2.]])
+
+ Cutting scalars different:
+ >>> m.cutScalars([s <= 1, s >= 2], replaceValue=0.).faceScalars
+ array([[ 0., 1.],
+ [ 2., 0.]])
+
+ Cuttin one special scalar Array only
+ >>> m.cutScalars(s <= 1, replaceValue=0., scalarIndex=1).faceScalars
+ array([[ 1., 0.],
+ [ 2., 2.]])
+
+ Using a list of cut expressions to cut every scalar index different
+
+ """
+ if inplace:
+ inst = self
+ else:
+ inst = self.copy()
+
+ if isinstance(cutExpression, list):
+ if scalarIndex is not None:
+ raise ValueError("scalarIndex must be None, "
+ "if cutExpression is list of cutExpressions")
+ if not len(cutExpression) == inst.getScalarDepth():
+ raise ValueError("lenght of cutExpression must meet scalar depth")
+ for si, ce in enumerate(cutExpression):
+ inst.cutScalars(ce, coords=coords,
+ replaceValue=replaceValue,
+ scalarIndex=si, inplace=True)
+ else:
+ if coords is None:
+ freeSymbols = cutExpression.free_symbols
+ if len(freeSymbols) > 1:
+ raise ValueError('coords must be given if multiple variables are given')
+ elif len(freeSymbols) == 0:
+ raise NotImplementedError("Expressiongs like {cutExpression} "
+ "are not understood for coords".format(**locals()))
+ coords = list(freeSymbols) * inst.getScalarDepth()
+ scalarArrays = inst.getScalars()
+ if scalarIndex is not None:
+ scalarArrays = scalarArrays[:, scalarIndex:scalarIndex + 1]
+
+ maskBelow = tfields.getMask(scalarArrays,
+ cutExpression=cutExpression,
+ coords=[coords[scalarIndex]])
+ scalarArrays[maskBelow] = replaceValue
+ inst.faceScalars[:, scalarIndex:scalarIndex + 1] = scalarArrays
+ else:
+ maskBelow = tfields.getMask(scalarArrays,
+ cutExpression=cutExpression,
+ coords=coords)
+ scalarArrays[maskBelow] = replaceValue
+ inst.faceScalars = scalarArrays
+ if not inplace:
+ return inst
+
+ def getFaceMask(self, mask):
+ """
+ Examples:
+ >>> m = tfields.Mesh3D([[1,2,3], [3,3,3], [0,0,0], [5,6,7]],
+ ... [[0, 1, 2], [1, 2, 3]],
+ ... faceScalars=[[1,2,3,4,5], [6,7,8,9,0]])
+ >>> from sympy.abc import x,y,z
+ >>> vertexMask = m.getMask(z < 6)
+ >>> faceMask = m.getFaceMask(vertexMask)
+ >>> faceMask
+ array([ True, False], dtype=bool)
+
+ Returns:
+ mask of faces with all vertices in mask
+ """
+ faceDeleteMask = np.full((self.faces.shape[0]), False, dtype=bool)
+ indices = np.array(range(self.getNPoints()))
+ deleteIndices = set(indices[~mask]) # set speeds up everything
+ for i, face in enumerate(self.faces):
+ for index in face:
+ if index in deleteIndices:
+ faceDeleteMask[i] = True
+ break
+
+ return ~faceDeleteMask
+
+ def getRemovedVertices(self, vertexDeleteMask):
+ """
+ Return copy of self without vertices where vertexDeleteMask is True
+ Copy because self is immutable
+
+ Examples:
+ >>> m = tfields.Mesh3D([[0,0,0], [1,1,1], [2,2,2], [0,0,0],
+ ... [3,3,3], [4,4,4], [5,5,5]],
+ ... [[0, 1, 2], [0, 1, 3], [3, 4, 5], [3, 4, 1],
+ ... [3, 4, 6]],
+ ... faceScalars=[[1,2], [3,4], [5,6], [7,8], [9,0]])
+ >>> c = m.getRemovedVertices([True, True, True, False, False,
+ ... False, False])
+ >>> c
+ Mesh3D([[ 0., 0., 0.],
+ [ 3., 3., 3.],
+ [ 4., 4., 4.],
+ [ 5., 5., 5.]])
+ >>> c.faces
+ array([[0, 1, 2],
+ [0, 1, 3]])
+ >>> c.faceScalars
+ array([[ 5., 6.],
+ [ 9., 0.]])
+
+ """
+ log = logger.new()
+ vertexDeleteMask = np.array(vertexDeleteMask)
+ log.verbose("Remove {0} vertices.".format(vertexDeleteMask.sum()))
+ # built instance that only contains the vaild points
+ inst = self[~vertexDeleteMask].copy()
+
+ moveUpCounter = np.zeros(self.faces.shape, dtype=int)
+
+ # correct faces:
+ deleteIndices = np.arange(self.getNPoints())[vertexDeleteMask]
+ for p in deleteIndices:
+ moveUpCounter[self.faces > p] -= 1
+
+ faceKeepMask = self.getFaceMask(~vertexDeleteMask)
+ inst.faces = (self.faces + moveUpCounter)[faceKeepMask]
+ inst.faceScalars = self.faceScalars[faceKeepMask]
+ return inst
+
+ def cleaned(self, stale=True, duplicates=True):
+ """
+ Args:
+ stale (bool): remove stale vertices
+ duplicates (bool): replace duplicate vertices by originals
+ Examples:
+ >>> m = tfields.Mesh3D([[0,0,0], [1,1,1], [2,2,2], [0,0,0],
+ ... [3,3,3], [4,4,4], [5,6,7]],
+ ... [[0, 1, 2], [3, 4, 5]],
+ ... faceScalars=[[1,2,3,4,5], [6,7,8,9,0]])
+ >>> c = m.clean()
+ >>> c
+ Mesh3D([[ 0., 0., 0.],
+ [ 1., 1., 1.],
+ [ 2., 2., 2.],
+ [ 3., 3., 3.],
+ [ 4., 4., 4.]])
+ >>> c.faces
+ array([[0, 1, 2],
+ [0, 3, 4]])
+
+ Returns:
+ copy of self without stale vertices and duplicat points
+ """
+ log = logger.new()
+ log.verbose("Cleaning up.")
+ # remove stale vertices
+ if stale:
+ vertexDeleteMask = self.staleVertices()
+ else:
+ vertexDeleteMask = np.full(self.shape[0], False, dtype=bool)
+ # remove duplicates in order to not have any artificial separations
+ inst = self
+ if duplicates:
+ inst = self.copy()
+ log.verbose("Finding Duplicates")
+ dups = tfields.duplicates(self, axis=0)
+ for i, dupi in zip(range(self.shape[0]), dups):
+ if dupi != i:
+ log.verbose("Found Duplicate point @ index {0}".format(i))
+ vertexDeleteMask[i] = True
+ # redirect faces
+ log.verbose("Run trough all faces to let it point to the"
+ "original")
+ for f in range(self.getNFaces()):
+ if i in self.faces[f]:
+ index = tfields.index(self.faces[f], i)
+ inst.faces[f][index] = dupi
+
+ return inst.getRemovedVertices(vertexDeleteMask)
+
+ def clean(self, *args, **kwargs):
+ """
+ Deprecated
+ """
+ warnings.warn("Name clean is deprecated. take cleaned instead", DeprecationWarning)
+ return self.cleaned(*args, **kwargs)
+
+ def getScalarMap(self, mask):
+ """
+ Return copy of self without vertices where mask is True
+ Copy because self is immutable
+ """
+ # built instance that only contains the vaild points
+ faceKeepMask = self.getFaceMask(mask)
+ scalarMap = np.arange(self.getNFaces())[faceKeepMask]
+ return scalarMap
+
+ def removeFaces(self, faceDeleteMask):
+ """
+ Remove faces where faceDeleteMask is True
+ Examples:
+ >>> m = tfields.Mesh3D([[1,2,3], [3,3,3], [0,0,0], [5,6,7]],
+ ... [[0, 1, 2], [1, 2, 3]],
+ ... faceScalars=[[1,2], [6,7]])
+ >>> m.removeFaces([True, False])
+ >>> m.faces
+ array([[1, 2, 3]])
+
+ """
+ faceDeleteMask = np.array(faceDeleteMask, dtype=bool)
+ self.faces = self.faces[~faceDeleteMask]
+ self.faceScalars = self.faceScalars[~faceDeleteMask]
+
+ def keepFaces(self, faceMask=None, faces=None, faceIndices=None):
+ """
+ Inverse method like removeFaces
+ Args:
+ faceMask (np.array):
+ faces (list of list of int)
+ faceIndices (list of int)
+ """
+ if faces is None:
+ faces = []
+ if faceIndices is None:
+ faceIndices = []
+ if faceMask is None:
+ faceMask = np.full(self.faces.shape[0], False, dtype=bool)
+
+ for i, face in enumerate(self.faces):
+ # np. version of if face in faces:
+ if any((face == f).all() for f in faces):
+ faceIndices.append(i)
+
+ for ind in faceIndices:
+ faceMask[ind] = True
+
+ self.removeFaces(~faceMask)
+
+ def staleVertices(self):
+ """
+ Returns:
+ Mask for all vertices that are stale i.e. are not refered by faces
+ """
+ staleMask = np.full(self.getNPoints(), False, dtype=bool)
+ used = set(self.faces.flatten())
+ for i in range(self.getNPoints()):
+ if i not in used:
+ staleMask[i] = True
+ return staleMask
+
+ def getFaces(self, vertex=None):
+ """
+ Args:
+ vertex (None / int / array of length 3)
+ """
+ if vertex is None:
+ return self.faces
+ if isinstance(vertex, int):
+ vertex = self[vertex]
+ if not (isinstance(vertex, list) or isinstance(vertex, np.ndarray)):
+ raise TypeError("Vertex has wrong type {0}".format(type(vertex)))
+ index = tfields.index(self, vertex, axis=0)
+ faces = []
+ for face in self.faces:
+ if index in face:
+ faces.append(face)
+ return faces
+
+ def _inputToFaceIndices(self, arg):
+ """
+ convert an input to a faceIndices list
+ Returns:
+ list
+ """
+ arg = np.array(arg)
+ if arg.dtype == bool:
+ # mask
+ return np.arange(self.faces.shape[0])[arg]
+ if len(arg.shape) > 1:
+ # face
+ raise NotImplementedError()
+ else:
+ return arg
+
+ def _inputToFaceMask(self, arg):
+ """
+ convert an input to a face mask
+ Returns:
+ np.array, dtype=bool
+ """
+ arg = np.array(arg)
+ if arg.dtype == bool:
+ # mask
+ return arg
+ if len(arg.shape) > 1:
+ # face
+ raise NotImplementedError()
+ else:
+ # faceIndices
+ tmp = np.full(self.faces.shape[0], False)
+ tmp[arg] = True
+ return tmp
+
+ def getParts(self, faceGroupIndicesList):
+ """
+ Args:
+ faceGroupIndicesList (list of int)
+ """
+ log = logger.new()
+ faceIndices = range(len(self.faces))
+ parts = []
+ log.verbose("Run through all {0} groups and partition mesh"
+ .format(len(faceGroupIndicesList)))
+ for f, faceGroupIndices in enumerate(faceGroupIndicesList):
+ log.verbose("Group {0} / {1}".format(f, len(faceGroupIndicesList)))
+ mesh = self.copy()
+ # for speed up:
+ faceGroupIndices = set(faceGroupIndices)
+ faceDeleteMask = [True
+ if i not in faceGroupIndices
+ else False
+ for i in faceIndices]
+ mesh.removeFaces(faceDeleteMask)
+ mesh = mesh.getRemovedVertices(mesh.staleVertices())
+ parts.append(mesh)
+ return parts
+
+ def getLinkedFaces(self, skipFaces=None):
+ """
+ Retrieve the faceIndices that are connected grouped together
+ Args:
+ skipFaces: faceSelector (mask, faces, faceIndices)
+ Returns:
+ list of list of int: groups of face indices that are linked
+
+ Examples:
+ >>> import tfields
+ >>> a = tfields.Mesh3D([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]],
+ ... faces=[[0, 1, 2], [0, 2, 3]])
+ >>> b = a.copy()
+
+ >>> b[:, 0] += 2
+ >>> m = tfields.Mesh3D([a, b])
+ >>> groupIndices = m.getLinkedFaces()
+ >>> parts = m.getParts(groupIndices)
+ >>> aa, ba = parts
+ >>> bool((aa.faces == a.faces).all())
+ True
+ >>> bool((ba.faces == b.faces).all())
+ True
+ >>> bool((aa == a).all())
+ True
+ >>> bool((ba == b).all())
+ True
+
+ """
+ faces = self.faces
+ if skipFaces is not None:
+ mask = ~self._inputToFaceMask(skipFaces)
+ faces = faces[mask]
+ faceGroupIndicesList = pyTools.setTools.disjointGroupIndices(faces)
+ if skipFaces is not None:
+ faceIndices = np.arange(self.faces.shape[0])
+ faceGroupIndicesList = [faceIndices[mask][group]
+ for group in faceGroupIndicesList]
+ return faceGroupIndicesList
+
+ def getRegion(self, seedFace, **kwargs):
+ """
+ Grow a region from the seedFace until breaking criterion is reached
+ Breaking criterion is specified in kwargs
+ Args:
+ seedFace (faceMask or faces or faceIndices):
+ **kwargs: keys:
+ maxAngle: breaking criterion specified for the normal
+ vectors not to deviate from neighbours more than maxAngle
+ Examples:
+ Get only one side of a cube:
+ >>> import tfields
+ >>> import numpy as np
+ >>> base = [np.linspace(0, 1, 10),
+ ... np.linspace(0, 1, 10),
+ ... np.linspace(0, 1, 10)]
+ >>> mesh = tfields.Mesh3D.createMeshGrid(*base).cleaned()
+
+ Some small mistake occured in the test. Check that.
+ # Select the first face as a seedFace
+ # >>> faceGroups = mesh.getRegion([0], maxAngle=np.pi * 2 / 8)
+ # >>> parts = mesh.getParts(faceGroups)
+
+ # Should only return one group. does not yet -> TODO!
+ # >>> len(parts) == 1
+
+ """
+ log = logger.new()
+ if not kwargs:
+ log.warning("No boundaries specified")
+ return np.arange(self.faces.shape[0])
+
+ faceIndices = list(self._inputToFaceIndices(seedFace))
+
+ # get break condition from kwargs
+ maxAngle = kwargs.pop('maxAngle', None)
+
+ norms = self.triangles.getNormVectors()
+ meanVector = np.mean(norms[faceIndices], axis=0)
+
+ excludedFaceIndices = set()
+ length = 0
+ while len(faceIndices) > length:
+ length = len(faceIndices)
+ for f, face in enumerate(self.faces):
+ vertexIndices = list(set(pyTools.flatten(self.faces[faceIndices])))
+ for index in vertexIndices:
+ if index not in face:
+ continue
+ if f in faceIndices:
+ continue
+ if f in excludedFaceIndices:
+ continue
+ norm = norms[f]
+ angle = np.arccos(np.einsum("...j,...j", meanVector, norm))
+ if abs(angle) > maxAngle:
+ excludedFaceIndices.add(f)
+ continue
+ log.verbose("Found new neighbour at face index "
+ "{f}".format(**locals()))
+ faceIndices.append(f)
+ if not len(faceIndices) > length:
+ log.info("Found no neighbours")
+ return faceIndices
+
+ def getSides(self, mainAxes=None, deviation=2 * np.pi / 8):
+ """
+ Grouping together face indices that have normal vectors in the
+ limits of +- deviation or +- pi + deviation.
+ Examples:
+ Get only one side of a cube:
+ >>> import tfields
+ >>> import numpy as np
+ >>> base = [np.linspace(0, 1, 2),
+ ... np.linspace(0, 1, 4),
+ ... np.linspace(0, 1, 4)]
+ >>> mesh = tfields.Mesh3D.createMeshGrid(*base).cleaned()
+
+ Select the first face as a seedFace
+ >>> faceGroups = mesh.getSides([[1,0,0],[0,1,0],[0,0,1]])
+ >>> parts = mesh.getParts(faceGroups)
+ >>> len(parts) == 6
+ True
+
+ Faces that have inconsistant norm vector direction are no problem
+ To show that, we invert the normal vector of one
+ face in the middle of the cube
+ >>> mesh.faces[8] = [5, 9, 6]
+ >>> faceGroups2 = mesh.getSides([[1,0,0],[0,1,0],[0,0,1]])
+ >>> parts2 = mesh.getParts(faceGroups2)
+ >>> len(parts2) == 6
+ True
+
+ """
+ if mainAxes is None:
+ axes = self.getMainAxes()
+ else:
+ axes = tfields.Points3D(mainAxes)
+ n = np.apply_along_axis(np.linalg.norm, 0, axes.T).reshape(-1, 1)
+ axes = axes / n
+
+ norms = self.triangles.getNormVectors()
+ norms = tfields.Points3D(norms)
+
+ faceGroupIndices = []
+ for vector in axes:
+ angles = np.arccos(np.einsum("...ij,...j", norms, vector))
+ mask = np.logical_or(abs(angles) < deviation,
+ abs(angles - np.pi) < deviation)
+ tmp = self.getLinkedFaces(skipFaces=~mask)
+ faceGroupIndices += tmp
+ return faceGroupIndices
+
+ def getMeshMap(self, subMesh, delta=1e-9):
+ """
+ Returns:
+ Mesh3D: meshMap (see mapToCut), can be used as meshMap to retrieve
+ subMesh from self instance
+ Examples:
+ >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,2,0], [1,2,0]],
+ ... faces=[[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
+ ... faceScalars=[[1],[2],[3],[4]])
+ >>> from sympy.abc import y
+ >>> mCut, mapMesh = m.mapToCut(y < 1.5, atIntersection='split')
+ >>> mm = m.getMeshMap(mCut)
+ >>> bool((mm == mapMesh).all())
+ True
+ >>> bool((mm.faceScalars == mapMesh.faceScalars).all())
+ True
+ >>> bool((mm.faces == mapMesh.faces).all())
+ True
+
+ """
+ # log = logger.new()
+ faceIndices = np.arange(self.faces.shape[0])
+ cents = tfields.Points3D(subMesh.getCentroids())
+ scalars = []
+ # nCents = cents.shape[0]
+ # for i in range(nCents):
+ # faceMask = self.pointsInMesh(cents[i: i + 1], delta=delta)[0]
+ # if True not in faceMask:
+ # log.warning("No point face was assigned to this. I will retry.")
+ # faceMask = self.pointsInMesh(cents[i: i + 1], delta=delta * 100)[0]
+ # if True not in faceMask:
+ # raise ValueError()
+ # else:
+ # log.info("Centroid {i} / {nCents}".format(**locals()))
+ # scalars.append(faceIndices[faceMask])
+ mask = self.pointsInMesh(cents, delta=delta)
+ scalars = [faceIndices[faceMask] for faceMask in mask]
+ inst = subMesh.copy()
+ inst.setScalarArray(0, scalars)
+ return inst
+
+
+ def _distFromPlane(self, point, plane):
+ return plane['normal'].dot(point) + plane['d']
+
+ def _getSegmentPlaneIntersection(self, p0, p1, plane):
+ distance0 = self._distFromPlane(p0, plane)
+ distance1 = self._distFromPlane(p1, plane)
+ p0OnPlane = abs(distance0) < np.finfo(float).eps
+ p1OnPlane = abs(distance1) < np.finfo(float).eps
+ outPoints = []
+ direction = 0
+ if p0OnPlane:
+ outPoints.append(p0)
+
+ if p1OnPlane:
+ outPoints.append(p1)
+ # remove duplicate points
+ if len(outPoints) > 1:
+ outPoints = tfields.Points3D(np.unique(outPoints, axis=0))
+ if p0OnPlane and p1OnPlane:
+ return outPoints, direction
+
+ if distance0 * distance1 > np.finfo(float).eps:
+ return outPoints, direction
+
+ direction = np.sign(distance0)
+ if abs(distance0) < np.finfo(float).eps:
+ return outPoints, direction
+ elif abs(distance1) < np.finfo(float).eps:
+ return outPoints, direction
+ if abs(distance0 - distance1) > np.finfo(float).eps:
+ t = distance0 / (distance0 - distance1)
+ else:
+ return outPoints, direction
+
+ outPoints.append(p0 + t * (p1 - p0))
+ # remove duplicate points
+ if len(outPoints) > 1:
+ outPoints = tfields.Points3D(np.unique(outPoints, axis=0))
+ return outPoints, direction
+
+ def _intersect(self, triangle, plane, facePointsRejected):
+ nTrue = facePointsRejected.count(True)
+ lonelyBool = True if nTrue == 1 else False
+ index = facePointsRejected.index(lonelyBool)
+ s0, d0 = self._getSegmentPlaneIntersection(triangle[0], triangle[1], plane)
+ s1, d1 = self._getSegmentPlaneIntersection(triangle[1], triangle[2], plane)
+ s2, d2 = self._getSegmentPlaneIntersection(triangle[2], triangle[0], plane)
+
+ singlePoint = triangle[index]
+ couplePoints = [triangle[j] for j in range(3)
+ if not facePointsRejected[j] == lonelyBool]
+
+ # TODO handle special cases. For now triangles with at least two points on plane are excluded
+ newPoints = None
+ faces = None
+
+ if len(s0) == 2:
+ # both points on plane
+ return None, None
+ if len(s1) == 2:
+ # both points on plane
+ return None, None
+ if len(s2) == 2:
+ # both points on plane
+ return None, None
+ if lonelyBool:
+ # one new triangle
+ if len(s0) == 1 and len(s1) == 1:
+ newPoints = np.array([couplePoints[0], couplePoints[1], s0[0], s1[0]])
+ faces = np.array([[0, 2, 1], [1, 2, 3]])
+ elif len(s1) == 1 and len(s2) == 1:
+ newPoints = np.array([couplePoints[0], couplePoints[1], s1[0], s2[0]])
+ faces = np.array([[0, 1, 2], [0, 2, 3]])
+ elif len(s0) == 1 and len(s2) == 1:
+ newPoints = np.array([couplePoints[0], couplePoints[1], s0[0], s2[0]])
+ faces = np.array([[0, 1, 2], [1, 3, 2]])
+ else:
+ return None, None
+
+ else:
+ # two new triangles
+ if len(s0) == 1 and len(s1) == 1:
+ newPoints = np.array([singlePoint, s0[0], s1[0]])
+ faces = np.array([[0, 2, 1]])
+ elif len(s1) == 1 and len(s2) == 1:
+ newPoints = np.array([singlePoint, s1[0], s2[0]])
+ faces = np.array([[0, 2, 1]])
+ elif len(s0) == 1 and len(s2) == 1:
+ newPoints = np.array([singlePoint, s0[0], s2[0]])
+ faces = np.array([[0, 1, 2]])
+ else:
+ return None, None
+ return newPoints, faces
+
+ def mapToCut(self, cutOrMeshMap, coordSys=None, atIntersection="remove"):
+ """
+ Partition the mesh with the cuts given
+ Args:
+ cutOrMeshMaps (cutExpression or meshMap):
+ if cutOrMeshMap is a cutExpression: cut like in self.cut but also return
+ a map from new to old faces
+ if cutOrMeshMap is a meshMap: map without cutting according to the map given
+ coordSys (str): coordSys for cutExpression
+ atIntersection (str): option switch for cutExpression (see
+ tfields.Mesh3D.cut)
+ Examples:
+ Cut away the first face / select the last face only
+ >>> m = tfields.Mesh3D([[1,2,3], [3,3,3], [0,0,0], [5,6,7]],
+ ... faces=[[0, 1, 2], [1, 2, 3]],
+ ... faceScalars=[[1,2], [6,7]])
+ >>> mNew = tfields.Mesh3D([[3,3,3], [0,0,0], [5,6,7]], [[0, 1, 2]], faceScalars=[[1]])
+ >>> mCut, _ = m.mapToCut(mNew)
+ >>> mCut.faceScalars
+ array([[ 6., 7.]])
+
+ Cut with condition will return the manual/instruction on how to
+ conduct the cut fast (mapMesh)
+ >>> from sympy.abc import y
+ >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,2,0], [1,2,0]],
+ ... faces=[[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
+ ... faceScalars=[[1],[2],[3],[4]])
+ >>> mCut, mapMesh = m.mapToCut(y < 1.5, atIntersection='split')
+ >>> mCut.faceScalars.T
+ array([[ 1., 2., 3., 3., 4.]])
+
+ Applying the mapMesh results in the same result as applying the cut
+ with the condition
+ >>> mCut2, _ = m.mapToCut(mapMesh)
+ >>> bool((mCut == mCut2).all())
+ True
+ >>> bool((mCut.faceScalars == mCut2.faceScalars).all())
+ True
+
+ """
+ cutExpression = None
+ meshMap = None
+ # check what cutOrMeshMap is and set it
+ if isinstance(cutOrMeshMap, Mesh3D):
+ meshMap = cutOrMeshMap
+ else:
+ cutExpression = cutOrMeshMap
+
+ if cutExpression is not None:
+ eps = 0.000000001
+ # direct return if self is empty
+ if len(self) == 0:
+ return self.copy(), self.copy()
+
+ # mask for points that do not fulfill the cut expression
+ mask = self.getMask(cutExpression, coordSys=coordSys)
+ # remove the points
+ inst = self.getRemovedVertices(~mask)
+ scalarMap = self.getScalarMap(mask)
+ scalarMap = scalarMap.reshape((-1, 1)).astype(float)
+
+ if not any(~mask):
+ # no faces have been removed.
+ return inst, Mesh3D([inst], faceScalars=scalarMap)
+
+ if all(~mask):
+ # all faces have been removed)
+ return inst, inst.copy()
+
+ # add points and faces intersecting with the plane
+ if atIntersection == 'split' or atIntersection == 'splitRough':
+ cutExpressionParts = symTools.getExpressionParts(cutExpression)
+ if len(cutExpressionParts) > 1:
+ onEdgeMesh = self.copy()
+ if atIntersection == 'splitRough':
+ """
+ the following is, to speed up the process. Problem is, that
+ triangles can exist, where all points lie outside the cut,
+ but part of the area
+ still overlaps with the cut.
+ These are at the intersection line between two cuts.
+ """
+ faceIntersMask = np.full((self.faces.shape[0]), False, dtype=bool)
+ for i, face in enumerate(self.faces):
+ facePointsRejected = [-mask[f] for f in face]
+ faceOnEdge = any(facePointsRejected) and not all(facePointsRejected)
+ if faceOnEdge:
+ faceIntersMask[i] = True
+ onEdgeMesh.removeFaces(-faceIntersMask)
+
+ for exprPart in cutExpressionParts:
+ onEdgeMesh = onEdgeMesh.cut(exprPart, atIntersection='split')
+ newMesh = onEdgeMesh
+ elif len(cutExpressionParts) == 1:
+
+ points = [sympy.symbols('x0, y0, z0'), sympy.symbols('x1, y1, z1'), sympy.symbols('x2, y2, z2')]
+ fpr = sympy.symbols('fpr2, fpr1, fpr2')
+ planeSympy = symTools.getPlane(cutExpression)
+ n = np.array(planeSympy.normal_vector).astype(float)
+ d = -n.dot(np.array(planeSympy.p1).astype(float))
+ plane = {'normal': n, 'd': d}
+
+ normVectors = self.getNormVectors()
+ newPoints = np.empty((0, 3))
+ newFaces = np.empty((0, 3))
+ newFaceScalars = []
+ newNormVectors = []
+ newScalarMap = []
+ nNew = 0
+ for i, face in enumerate(self.faces):
+ """
+ facePointsRejected is a mask for each face that is True, where
+ a Point is on the rejected side of the plane
+ """
+ facePointsRejected = [~mask[f] for f in face]
+ faceOnEdge = any(facePointsRejected) and not all(facePointsRejected)
+ if faceOnEdge:
+ """
+ define the two lines that are intersecting the plane:
+ one point will be alone on one side of the cutting plane (singlePoint)
+ the other two will be on ther other side (couplePoints)
+ """
+ nTrue = facePointsRejected.count(True)
+ # the bool that occures on
+ lonelyBool = True if nTrue == 1 else False
+
+ triangle_points = [self[f] for f in face]
+ """
+ Add the intersection points and faces
+ """
+ # tick = time()
+ if lonelyBool:
+ # singlePoint on cut away side
+ newP, newF = self._intersect(triangle_points, plane, facePointsRejected)
+ if newP is not None:
+ newP = np.array(newP)
+ np.place(newP, np.logical_and(0 < newP, newP <= eps), 0.0)
+ newPoints = np.concatenate((newPoints, newP))
+ newFaces = np.concatenate((newFaces, newF + nNew))
+ newFaceScalars.extend([self.faceScalars[i]] * 2)
+ newNormVectors.extend([normVectors[i]] * 2)
+ newScalarMap.extend([i] * 2)
+ nNew += 4
+ else:
+ pass
+ else:
+ # singlePoint on keep side
+ newP, newF = self._intersect(triangle_points, plane, facePointsRejected)
+ if newP is not None:
+ newP = np.array(newP)
+ np.place(newP, np.logical_and(0 < newP, newP <= eps), 0.0)
+ newPoints = np.concatenate((newPoints, newP))
+ newFaces = np.concatenate((newFaces, newF + nNew))
+ newFaceScalars.extend([self.faceScalars[i]] * 1)
+ newNormVectors.extend([normVectors[i]] * 1)
+ newScalarMap.extend([i])
+ nNew += 3
+ else:
+ pass
+
+ newPoints = [[float(p[0]), float(p[1]), float(p[2])] for p in newPoints]
+ newMesh = self.__class__(newPoints, faces=newFaces,
+ coordSys=inst.getCoordSys(),
+ faceScalars=newFaceScalars)
+ newMesh.alignNormVectors(newNormVectors)
+ newScalarMap = np.array(newScalarMap).reshape((-1, 1)).astype(float)
+ scalarMap = np.concatenate((scalarMap, newScalarMap))
+ inst = self.__class__([inst, newMesh])
+ else:
+ raise ValueError("CutExpression is not splitable.")
+
+ elif atIntersection == 'remove':
+ pass
+ else:
+ raise AttributeError("No atIntersection method called {atIntersection} "
+ "implemented".format(**locals()))
+ meshMap = inst.copy()
+ meshMap.faceScalars = scalarMap
+ elif meshMap is not None:
+ if coordSys is not None or atIntersection != "remove":
+ log = logger.new()
+ log.warning("coordSys and atIntersection arguments are not default.")
+ if len(meshMap) > 0:
+ inst = Mesh3D(meshMap,
+ faces=meshMap.faces,
+ faceScalars=self.faceScalars[meshMap.getScalarArrays()[0].astype(int)],
+ coordSys=meshMap.coordSys)
+ else:
+ inst = Mesh3D([])
+ else:
+ raise TypeError("cutOrMeshMap is not understood.")
+ return inst, meshMap
+
+ def cut(self, cutExpression, coordSys=None, atIntersection="remove"):
+ """
+ cut method for Mesh3D.
+ If the same cut to the same mesh should be applied multiple times with different scalars:
+ use mapToCut!
+ Args:
+ see Points3D.cut
+ additionally:
+
+ atIntersection (str): instruction on what to do, when a cut will intersect a triangle.
+ Options: "remove" (Default)
+ "split" - Create new triangles that make up the old one.
+
+ Examples:
+ define the cut
+ >>> from sympy.abc import x,y,z
+ >>> cutExpr = x > 1.5
+
+ >>> m = tfields.Mesh3D.createMeshGrid(np.linspace(0,3,4),
+ ... np.linspace(0,3,4),
+ ... np.linspace(0, 0, 1), faceScalars=np.linspace(4, 8, 18))
+ >>> mNew = m.cut(cutExpr)
+ >>> mNew.getNPoints()
+ 8
+ >>> mNew.getNFaces()
+ 6
+ >>> float(mNew[:, 0].min())
+ 2.0
+
+ Cutting with the split option will create new triangles on the edge:
+ >>> mSplit = m.cut(cutExpr, atIntersection='split')
+ >>> float(mSplit[:, 0].min())
+ 1.5
+ >>> mSplit.getNPoints()
+ 29
+ >>> mSplit.getNFaces()
+ 15
+
+ Cuttin with a sympy BooleanFunction:
+ >>> cutExprBoolFun = (x > 1.5) & (y < 1.5) & (y >0.2) & (z > -0.5)
+ >>> mSplit = m.cut(cutExprBoolFun, atIntersection='split')
+
+ Returns:
+ copy of cut mesh
+
+ """
+ inst, _ = self.mapToCut(cutExpression, coordSys=coordSys, atIntersection=atIntersection)
+ return inst
+
+ def _clusterDbscan(self, maxClusterDist=15, minSamples=20, outType=None):
+ """
+ runs DBSCAN alcorithm to find seperate mesh clusters
+ KwArgs:
+ outType (type): type of cluster-points. E.g. np.array / Points3D / list...
+ Returns:
+ clusters (list of pointsLists), excludedPoints (pointsList)
+ """
+ raise NotImplementedError("")
+ clusters = [None]
+ excluded = None
+ return clusters, excluded
+
+ def _clusterBox(self,
+ min0=None, max0=None, nSections0=None,
+ min1=None, max1=None, nSections1=None,
+ min2=None, max2=None, nSections2=None,
+ clusterMeshes=None, outType=None, cutTree=False, withTransform=False, **cutKwargs):
+ log = logger.new()
+
+ if withTransform:
+ targetMesh, transformation, translate = self.transformDirections()
+ else:
+ targetMesh = self
+ transformation = np.eye(3)
+ translate = np.zeros(3, )
+
+ if cutTree and (clusterMeshes is None):
+ minmax = np.array([[min0, min1, min2], [max0, max1, max2]])
+ xRange = list(np.linspace(minmax[0,0], minmax[1,0], nSections0 + 1))[1:-1]
+ yRange = list(np.linspace(minmax[0,1], minmax[1,1], nSections1 + 1))[1:-1]
+ zRange = list(np.linspace(minmax[0,2], minmax[1,2], nSections2 + 1))[1:-1]
+ cuts = {'x': xRange, 'y': yRange, 'z': zRange}
+ box = {'x': minmax[:,0], 'y': minmax[:,1], 'z': minmax[:,2]}
+ log.info("Start creating cutting tree")
+ atIntersection = cutKwargs.get('atIntersection', 'remove')
+ cutTree = cuttingTree.Node(parent=None, mesh=targetMesh, remainingCuts=cuts,
+ atIntersection=atIntersection, box=box)
+ log.info("Done building the tree structure")
+
+ leaves = cutTree.findLeaves()
+ leaves = cuttingTree.Node.sortLeaves(leaves)
+
+ clusterMeshes = [leaf.getGlobalCutMesh() for leaf in leaves]
+ clusteredObjs = [leaf.mesh for leaf in leaves]
+
+ # TODO exclude Mesh
+ excludedObjs = Mesh3D([])
+ else:
+ clusteredObjs, excludedObjs, clusterMeshes, transformation = super(Mesh3D, self)._clusterBox(min0=min0, max0=max0,
+ nSections0=nSections0,
+ min1=min1, max1=max1,
+ nSections1=nSections1,
+ min2=min2, max2=max2,
+ nSections2=nSections2,
+ clusterMeshes=clusterMeshes,
+ outType=outType,
+ withTransform=withTransform,
+ **cutKwargs)
+ return clusteredObjs, excludedObjs, clusterMeshes, transformation
+
+ def _minimum_bounding_rectangle(self, points2D):
+ """
+ Find the smallest bounding rectangle for a set of points.
+ Returns a set of points representing the corners of the bounding box.
+
+ :param points: an nx2 matrix of coordinates
+ :rval: an nx2 matrix of coordinates
+ """
+ from scipy.ndimage.interpolation import rotate
+ pi2 = np.pi / 2.
+
+ # get the convex hull for the points
+ hull_points = points2D[ConvexHull(points2D).vertices]
+
+ # calculate edge angles
+ edges = np.zeros((len(hull_points) - 1, 2))
+ edges = hull_points[1:] - hull_points[:-1]
+
+ angles = np.zeros((len(edges)))
+ angles = np.arctan2(edges[:, 1], edges[:, 0])
+
+ angles = np.abs(np.mod(angles, pi2))
+ angles = np.unique(angles)
+
+ # find rotation matrices
+ # XXX both work
+ rotations = np.vstack([
+ np.cos(angles),
+ np.cos(angles - pi2),
+ np.cos(angles + pi2),
+ np.cos(angles)]).T
+ # rotations = np.vstack([
+ # np.cos(angles),
+ # -np.sin(angles),
+ # np.sin(angles),
+ # np.cos(angles)]).T
+ rotations = rotations.reshape((-1, 2, 2))
+
+ # apply rotations to the hull
+ rot_points = np.dot(rotations, hull_points.T)
+
+ # find the bounding points
+ min_x = np.nanmin(rot_points[:, 0], axis=1)
+ max_x = np.nanmax(rot_points[:, 0], axis=1)
+ min_y = np.nanmin(rot_points[:, 1], axis=1)
+ max_y = np.nanmax(rot_points[:, 1], axis=1)
+
+ # find the box with the best area
+ areas = (max_x - min_x) * (max_y - min_y)
+ best_idx = np.argmin(areas)
+
+ # return the best box
+ x1 = max_x[best_idx]
+ x2 = min_x[best_idx]
+ y1 = max_y[best_idx]
+ y2 = min_y[best_idx]
+ r = rotations[best_idx]
+
+ rval = np.zeros((4, 2))
+ rval[0] = np.dot([x1, y2], r)
+ rval[1] = np.dot([x2, y2], r)
+ rval[2] = np.dot([x2, y1], r)
+ rval[3] = np.dot([x1, y1], r)
+
+ return rval, r
+
+ def transformDirections(self):
+ """
+ Transforms the coordinate system such that it faces towards the mean of all triangle normals.
+ This only works if the triangle normals mean is a reasonable choice for the problem.
+ An improved version in the future should consider to weight the normals by triangle area
+ The minimal value in this coordinate system is 0, the maximal value is 1
+ :param inMesh: input mesh that the transformation is based on
+ :return: normalized mesh with transformed vertices, faces and faceScalars as in inMesh,
+ tranformationMatrix and translate transforms the normalizedMesh back to its original shape
+ """
+
+ zNormal = np.array([0.0, 0.0, 1.0])
+ meanNormal = np.mean(self.getNormVectors(), axis=0)
+ meanNormal = meanNormal / np.linalg.norm(meanNormal)
+
+ cosMeanNormalZ = meanNormal.dot(zNormal)
+ if cosMeanNormalZ == -1.0:
+ # mean normal points in opposite direction as zNormal
+ rotateMeanNormal = np.diag([1.0, 1.0, -1.0])
+ elif cosMeanNormalZ == 1.0:
+ # already in z direction
+ rotateMeanNormal = np.eye(3)
+ else:
+ # calculate transformation matrix to let meanNormal be transformed to zNormal
+ normalsCrossProduct = np.cross(zNormal, meanNormal)
+ crossProductNorm = np.linalg.norm(normalsCrossProduct)
+ cosMeanNormalZ = meanNormal.dot(zNormal)
+ # skew-symmetric cross product matrix of normalsCrossProduct
+ nux = np.array([[0, -normalsCrossProduct[2], normalsCrossProduct[1]],
+ [normalsCrossProduct[2], 0, -normalsCrossProduct[0]],
+ [-normalsCrossProduct[1], normalsCrossProduct[0], 0]])
+ # formula found at
+ # https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
+ rotateMeanNormal = np.eye(3) + nux + nux.dot(nux) * (1 - cosMeanNormalZ) / (crossProductNorm ** 2)
+ vertices3D = self.dot(rotateMeanNormal)
+ unrotatedMesh = Mesh3D(vertices3D, faces=self.faces, faceScalars=self.faceScalars)
+
+ vertices2D = vertices3D[:, :2]
+ boundingRectangleVertices2D, rot = self._minimum_bounding_rectangle(vertices2D)
+
+ rotatedRectangle = boundingRectangleVertices2D.dot(rot.T)
+
+ # longer side rotated parallel to x-axis
+ side0 = np.linalg.norm(rotatedRectangle[0, :] - rotatedRectangle[1, :])
+ side1 = np.linalg.norm(rotatedRectangle[0, :] - rotatedRectangle[3, :])
+ if side0 < side1:
+ rot = rot.dot(np.array([[0, -1], [1, 0]]))
+ rotMatrix3D = np.eye(3)
+ rotMatrix3D[:2, :2] = rot.T
+
+ vertices3D = vertices3D.dot(rotMatrix3D)
+
+ # find mins and max for scaling
+ mins = np.min(vertices3D, axis=0)
+ maxs = np.max(vertices3D, axis=0)
+ # scaling matrix
+ scale = np.diag([1.0 / (ma - mi) if (abs(ma - mi) > np.finfo(float).eps) else 1.0
+ for ma, mi in zip(maxs, mins)])
+ # translation vector
+ translate = np.array([mi / (mi - ma) if (abs(ma - mi) > np.finfo(float).eps) else mi
+ for ma, mi in zip(maxs, mins)])
+ normalizedVertices = vertices3D.dot(scale) + translate
+ normalizedMesh = Mesh3D(normalizedVertices, faces=self.faces, faceScalars=self.faceScalars)
+ tranformationMatrix = rotateMeanNormal.dot(rotMatrix3D.dot(scale))
+
+ return normalizedMesh, tranformationMatrix, translate
+
+ def alignNormVectors(self, normVectors):
+ """
+ Orientate the faces such, that their normVectors align to the normVectors given.
+ Examples
+ >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [-1,0,0], [0,1,0], [0,0,1]],
+ ... [[0, 1, 3], [1, 3, 4], [1, 3, 2]]);
+ >>> newNorms = m.getNormVectors() * -1
+ >>> m.alignNormVectors(newNorms)
+ >>> m.faces
+ array([[0, 3, 1],
+ [1, 4, 3],
+ [1, 2, 3]])
+
+ """
+ if not self.getNFaces() == 0:
+ # vector product < 0
+ mask = np.einsum('...i,...i', self.getNormVectors(), normVectors) < 0
+ """
+ the line:
+ " self.faces[:, [1, 2]][mask] = self.faces[:, [2, 1]][mask] "
+ would be a nice solution, but numpy does not mutate the [1, 2] but returns a copy
+
+ """
+ temp = np.copy(self.faces[mask, 1])
+ self.faces[mask, 1] = self.faces[mask, 2]
+ self.faces[mask, 2] = temp
+
+ def saveObj(self, filePath, **kwargs):
+ """
+ Save obj as wavefront/.obj file
+ """
+ obj = kwargs.pop('object', None)
+ group = kwargs.pop('group', None)
+
+ cmap = kwargs.pop('cmap', 'jet')
+ faceScalarIndex = kwargs.pop('faceScalarIndex', None)
+
+ filePath = filePath.replace('.obj', '')
+ fileDir, fileName = os.path.split(filePath)
+
+ if not (self.faceScalars.size == 0 or faceScalarIndex is None):
+ scalars = self.faceScalars[:, faceScalarIndex]
+ minScalar = scalars[-np.isnan(scalars)].min()
+ maxScalar = scalars[-np.isnan(scalars)].max()
+ vmin = kwargs.pop('vmin', minScalar)
+ vmax = kwargs.pop('vmax', maxScalar)
+ if vmin == vmax:
+ if vmin == 0.:
+ vmax = 1.
+ else:
+ vmin = 0.
+ norm = colors.Normalize(vmin, vmax)
+ colorMap = plt.get_cmap(cmap)
+ else:
+ # this is a switch for not coloring the triangles and thus not producing the materials
+ norm = None
+
+ if norm is not None:
+ fileNameMat = fileName + '_frame_{0}.mat'.format(faceScalarIndex)
+ scalars[np.isnan(scalars)] = minScalar - 1
+ sortedScalars = scalars[scalars.argsort()]
+ sortedScalars[sortedScalars == minScalar - 1] = np.nan
+ sortedFaces = self.faces[scalars.argsort()]
+ scalarSet = np.unique(sortedScalars)
+ scalarSet[scalarSet == minScalar - 1] = np.nan
+ with open(os.path.join(fileDir + fileNameMat), 'w') as mf:
+ for s in scalarSet:
+ if np.isnan(s):
+ mf.write("newmtl nan")
+ mf.write("Kd 0 0 0\n\n")
+ else:
+ mf.write("newmtl mtl_{0}\n".format(s))
+ mf.write("Kd {c[0]} {c[1]} {c[2]}\n\n".format(c=colorMap(norm(s))))
+ else:
+ sortedFaces = self.faces
+
+ # writing of the obj file
+ with open(filePath + '.obj', 'w') as f:
+ f.write("# File saved with tfields Mesh3D.saveObj method\n\n")
+ if norm is not None:
+ f.write("mtllib ./{0}\n\n".format(fileNameMat))
+ if obj is not None:
+ f.write("o {0}\n".format(obj))
+ if group is not None:
+ f.write("g {0}\n".format(group))
+ for vertex in self:
+ f.write("v {v[0]} {v[1]} {v[2]}\n".format(v=vertex))
+
+ lastScalar = None
+ for i, face in enumerate(sortedFaces + 1):
+ if norm is not None:
+ if not lastScalar == sortedScalars[i]:
+ lastScalar = sortedScalars[i]
+ f.write("usemtl mtl_{0}\n".format(lastScalar))
+ f.write("f {f[0]} {f[1]} {f[2]}\n".format(f=face))
+
+ def plot(self, **kwargs):
+ """
+ Frowarding to plotTools.plotMesh
+ """
+ # old version with plotTools
+ # import plotTools as pt
+ # return pt.plotMesh(self, **kwargs)
+
+ plotVertices = kwargs.pop('plotVertices', False)
+ faceScalarIndex = kwargs.pop('faceScalarIndex', None)
+ if faceScalarIndex is not None:
+ warnings.warn("faceScalarIndex is renamed scalarIndex. Also default is not 0 any more.",
+ DeprecationWarning)
+ kwargs['scalarIndex'] = kwargs.pop('scalarIndex', faceScalarIndex)
+
+ scalarsDemanded = any([v in kwargs for v in ['vmin', 'vmax', 'cmap']])
+ defaultScalarIndex = None if not scalarsDemanded else 0
+ # defaultScalarIndex = None
+ scalarIndex = kwargs.pop('scalarIndex', defaultScalarIndex)
+ if scalarIndex is not None:
+ if not self.faceScalars.shape[0] == 0:
+ kwargs['color'] = self.faceScalars[:, scalarIndex]
+ if plotVertices:
+ return mplTools.plotArray(self, **kwargs)
+ # deprecated dim=3 default. I want to go for 2 as everywhere
+ dimDefined = False
+ if 'axis' in kwargs:
+ dimDefined = True
+ if 'zAxis' in kwargs:
+ if kwargs['zAxis'] is not None:
+ kwargs['dim'] = 3
+ else:
+ kwargs['dim'] = 2
+ dimDefined = True
+ if 'dim' in kwargs:
+ dimDefined = True
+
+ if not dimDefined:
+ warnings.warn("Dimension is not defined. I want to change"
+ " default dimension to 2. Keep 3 for now. "
+ "If you want 3 specify keyword 'dim'.",
+ DeprecationWarning)
+ kwargs['dim'] = 3
+
+ return mplTools.plotMesh(self, self.faces, **kwargs)
+
+
+if __name__ == '__main__':
+ import doctest
+
+ # doctest.run_docstring_examples(Mesh3D.clean, globals())
+ doctest.testmod()
diff --git a/tfields/planes3D.py b/tfields/planes3D.py
new file mode 100644
index 0000000000000000000000000000000000000000..8d7dab3a9628fa0e142bd1253d3a3961d474ead2
--- /dev/null
+++ b/tfields/planes3D.py
@@ -0,0 +1,48 @@
+#!/usr/bin/env
+# encoding: utf-8
+"""
+Author: Daniel Boeckenhoff
+Mail: daniel.boeckenhoff@ipp.mpg.de
+
+part of npTools library
+"""
+import npTools
+import sympy
+import mplTools as mpt
+import loggingTools
+logger = loggingTools.Logger(__name__)
+
+
+class Planes3D(npTools.VectorField3D):
+ """
+ Point-NormVector representaion of planes
+ """
+
+ def symbolic(self):
+ """
+ Returns:
+ list: list with sympy.Plane objects
+ """
+ return [sympy.Plane(point, normal_vector=vector)
+ for point, vector in zip(self.points, self.vectors)]
+
+ def plot(self, **kwargs):
+ """
+ forward to Mesh3D plotting
+ """
+ artists = []
+ for i in range(self.points.shape[0]):
+ artists.append(mpt.plotPlane(self.points[0],
+ self.vectors[0],
+ **kwargs))
+ # symbolic = self.symbolic()
+ # planeMeshes = [npTools.Mesh3D([pl.arbitrary_point(t=(i + 1) * 1. / 2 * np.pi)
+ # for i in range(4)],
+ # faces=[[0, 1, 2], [1, 2, 3], [2, 3, 0]]) for pl in symbolic]
+ # artists = [m.plot(**kwargs) for m in planeMeshes]
+ return artists
+
+
+if __name__ == '__main__':
+ import doctest
+ doctest.testmod()
diff --git a/tfields/scalarField3D.py b/tfields/scalarField3D.py
new file mode 100644
index 0000000000000000000000000000000000000000..8926a77f091d2d74011e6094e5fb217765edbdc0
--- /dev/null
+++ b/tfields/scalarField3D.py
@@ -0,0 +1,346 @@
+#!/usr/bin/env
+# encoding: utf-8
+"""
+Author: Daniel Boeckenhoff
+Mail: daniel.boeckenhoff@ipp.mpg.de
+
+part of npTools library
+"""
+import numpy as np
+import npTools
+import loggingTools
+logger = loggingTools.Logger(__name__)
+
+
+class ScalarField3D(npTools.Points3D):
+ """Class to Store scalar Fields in 3D.
+ Multiple scalars can be saved. Thus scalars is a multidimensional array.
+
+ Note:
+
+ Args:
+ inputArray: see Points3D
+ scalars (list like object)
+
+ Attributes:
+ x1 (list): coordinate x
+ x2 (list): coordinate y
+ x3 (list): coordinate z
+ scalars (list, can be list of lists.)
+
+ Examples:
+ No need to fille a scalar to the field
+ >>> empty = npTools.ScalarField3D([[1,2,3], [4,5,6], [7,8,9]])
+ >>> empty.scalars
+ array([], shape=(3, 0), dtype=float64)
+
+ Multiple scalars also ok.
+ >>> s = npTools.ScalarField3D.createTest()
+ >>> s
+ points: [[ 3. 0. 6.]
+ [ 3. 0. 9.]]
+ scalars: [[ 0. 2. 0.]
+ [ 1. 2. 3.]]
+
+ Saving and reading works
+ >>> from tempfile import NamedTemporaryFile
+ >>> outFile = NamedTemporaryFile(suffix='.npz')
+ >>> s.savez(outFile.name)
+ >>> _ = outFile.seek(0)
+ >>> s1 = npTools.ScalarField3D.createFromFile(outFile.name)
+ >>> bool(np.all(s == s1))
+ True
+ >>> s1.scalars
+ array([[ 0., 2., 0.],
+ [ 1., 2., 3.]])
+
+ Adding scalar fields
+ >>> s1[:, 2] = 42
+ >>> s1.scalars[:, 0] = 24
+ >>> s2 = npTools.ScalarField3D([s, s1, s])
+ >>> s2
+ points: [[ 3. 0. 6.]
+ [ 3. 0. 9.]
+ [ 3. 0. 42.]
+ [ 3. 0. 42.]
+ [ 3. 0. 6.]
+ [ 3. 0. 9.]]
+ scalars: [[ 0. 2. 0.]
+ [ 1. 2. 3.]
+ [ 24. 2. 0.]
+ [ 24. 2. 3.]
+ [ 0. 2. 0.]
+ [ 1. 2. 3.]]
+
+ >>> s3 = npTools.ScalarField3D([s, s1, s], scalars = [[0], [1], [2], [3], [4], [5]])
+ >>> s3
+ points: [[ 3. 0. 6.]
+ [ 3. 0. 9.]
+ [ 3. 0. 42.]
+ [ 3. 0. 42.]
+ [ 3. 0. 6.]
+ [ 3. 0. 9.]]
+ scalars: [[0]
+ [1]
+ [2]
+ [3]
+ [4]
+ [5]]
+
+ """
+ __slots__ = ['coordSys', 'scalars']
+
+ def __new__(cls, inputArray, scalars=None, **kwargs):
+
+ if any([isinstance(a, ScalarField3D) for a in inputArray]) and scalars is None:
+ if not all([isinstance(a, ScalarField3D) for a in inputArray]):
+ # TODO: could allow if all scalars are none
+ raise TypeError("Copy constructor of ScalarField3D only accepts"
+ "ScalarField3D instances.")
+ if not len(set([o.scalars.shape for o in inputArray])) == 1:
+ raise TypeError("Shape of scalars that should be concatenated must be the same.")
+ scalars = np.concatenate([o.scalars for o in inputArray], axis=0)
+
+ obj = npTools.Points3D.__new__(cls, inputArray, **kwargs)
+
+ # add scalars
+ if scalars is None:
+ scalars = np.empty((len(obj), 0))
+ scalars = np.array(scalars)
+ if not len(scalars) == len(obj):
+ raise ValueError("Length of scalars ({0}) has to be the same as "
+ "points length ({1}) but is not.".format(len(scalars), len(obj)))
+ obj.scalars = scalars
+ return obj
+
+ def __repr__(self):
+ return "points: %s\nscalars: %s" % (self, self.scalars)
+
+ def cut(self, cutExpression, coordSys=None):
+ """
+ Examples:
+ >>> import npTools
+ >>> x, y, z = sympy.symbols('x y z')
+ >>> p = npTools.ScalarField3D([[1., 2., 3.], [4., 5., 6.], [1, 2, -6],
+ ... [-5, -5, -5], [1,0,-1], [0,1,-1]], scalars = [1,2,3,4,5,6])
+ >>> p.cut(x > 0)
+ points: [[ 1. 2. 3.]
+ [ 4. 5. 6.]
+ [ 1. 2. -6.]
+ [ 1. 0. -1.]]
+ scalars: [1 2 3 5]
+
+ cut combination
+ >>> p.cut((x > 0) & (z < 0))
+ points: [[ 1. 2. -6.]
+ [ 1. 0. -1.]]
+ scalars: [3 5]
+
+ """
+ mask = self.getMask(cutExpression, coordSys=coordSys)
+
+ inst = self[mask].copy()
+ inst.scalars = inst.scalars[mask]
+
+ return inst
+
+ @classmethod
+ def createTest(cls):
+ p = npTools.MeshGrid3D(3, 4, 1, 0, 2, 1, 6, 9, 2, iterOrder=[1, 0, 2], coordSys=cls.CYLINDER)
+ scalars = [[i * 1., 2, 3 * i] for i in range(len(p))]
+ return cls(p, scalars)
+
+ def evalAtPoint(self, point):
+ """
+ Examples:
+ >>> s = npTools.ScalarField3D.createTest()
+ >>> s.evalAtPoint(npTools.Points3D([[3, 0, 6]]))
+ array([[ 0., 2., 0.]])
+ >>> s.evalAtPoint([3, 0, 9])
+ array([[ 1., 2., 3.]])
+ >>> s.evalAtPoint([[0, 0, 0]])
+ array([], shape=(0, 3), dtype=float64)
+
+ """
+ return self.scalars[np.equal(self, point).all(axis=1)]
+
+ def getScalarDepth(self):
+ return self.scalars.shape[1]
+
+ def getScalarArrays(self):
+ """
+ Examples:
+ >>> s = npTools.ScalarField3D.createTest()
+ >>> s.getScalarArrays()
+ array([[ 0., 1.],
+ [ 2., 2.],
+ [ 0., 3.]])
+
+ """
+ return self.scalars.T
+
+ def setScalarArray(self, scalarIndex, scalarArray):
+ """
+ Examples:
+ >>> s = npTools.ScalarField3D.createTest()
+ >>> s.setScalarArray(1, [42, 84])
+ >>> s.getScalarArrays()
+ array([[ 0., 1.],
+ [ 42., 84.],
+ [ 0., 3.]])
+
+ >>> s.setScalarArray(3, [-1, -1])
+ >>> s.getScalarArrays()
+ array([[ 0., 1.],
+ [ 42., 84.],
+ [ 0., 3.],
+ [ -1., -1.]])
+
+ """
+ if scalarIndex == self.getScalarDepth():
+ self.appendScalars(scalarArray)
+ self.scalars[:, scalarIndex] = np.array(scalarArray)
+
+ def appendScalars(self, scalars):
+ """
+ Examples:
+ >>> s = npTools.ScalarField3D.createTest()
+ >>> s.appendScalars([42, 84])
+ >>> s.getScalarArrays()
+ array([[ 0., 1.],
+ [ 2., 2.],
+ [ 0., 3.],
+ [ 42., 84.]])
+
+ """
+ self.scalars = np.concatenate([self.scalars, np.array([scalars]).T], 1)
+
+ def dropScalars(self, indices=None):
+ """
+ Examples:
+ >>> s = npTools.ScalarField3D.createTest()
+ >>> s.appendScalars([42, 84])
+
+ give index
+ >>> s.dropScalars(-2)
+ >>> s.getScalarArrays()
+ array([[ 0., 1.],
+ [ 2., 2.],
+ [ 42., 84.]])
+
+ or indices
+ >>> s.dropScalars([0, 1])
+ >>> s.getScalarArrays()
+ array([[ 42., 84.]])
+
+ """
+ if indices:
+ self.scalars = np.delete(self.scalars, indices, 1)
+
+ def triangulate(self, method="areaConserving", **kwargs):
+ """
+ Args:
+ method (str): triangulation method
+
+ valid values:
+ - "areaConserving"
+ kwargs:
+ vicinity (list of list of int): map of point
+ neigbourship relations.
+
+ Examples:
+ Imagine you have a camera with some pixels.
+ Each pixel has a center point which can be mapped to the area viewed in 3d.
+ >>> points = np.array([[0,0,0],[1,0,0],[2,0,0],[0,1,0],[1,1,0],[2,1,0],[3,1,0],
+ ... [0,2,0],[1,2,0],[1.7,1.7,0],[3,2,1],[1,3,1],[2,3,1],[3,2,1]])
+
+ Every pixel and thus every point has a value belonging to it.
+ >>> scalars = points * 42.1
+ >>> sf = npTools.ScalarField3D(points, scalars=scalars)
+
+ It is known which point belongs to which pixel
+ >>> pixels = np.array([[0,0],[1,0],[2,0],[0,1],[1,1],[2,1],[3,1],
+ ... [0,2],[1,2],[2,2],[3,2],[1,3],[2,3],[3,3]])
+
+ Triangulation:
+ >>> t = sf.triangulate(method="areaConserving", vicinity=pixels)
+ >>> t.getNTriangles()
+ 84
+ >>> t[9:12] # 3rd triangle
+ Triangles3D([[ 1. , 1. , 0. ],
+ [ 1. , 0.5 , 0. ],
+ [ 0.66666667, 0.66666667, 0. ]])
+ >>> t.triangleScalars[3]
+ array([ 42.1, 42.1, 0. ])
+
+ """
+ if method == "areaConserving":
+ vicinity = kwargs.get('vicinity')
+ if vicinity is None:
+ raise AttributeError("kwarg vicinity is mandatory for areaConserving method.")
+ pointsVicinityMap = ScalarField3D(self, scalars=vicinity)
+ pixelTriangleCornerPoints = []
+ if self.scalars.size == 0:
+ pixelTriangleScalars = None
+ else:
+ pixelTriangleScalars = []
+ relPixels = [[1, 0], [1, -1], [0, -1]]
+ for i, pixel in enumerate(pointsVicinityMap.scalars):
+ locNeighbours = []
+ for relP in relPixels:
+ mask = (pointsVicinityMap.scalars == pixel + relP).all(axis=1)
+ locNeighbours.append(mask)
+
+ pointI = pointsVicinityMap[i]
+ # locNeighbours.append(locNeighbours[0])
+ locNeighbours = np.array(locNeighbours)
+ # if not locNeighbours.sum() == 9:
+ # continue
+ for n in range(2):
+ if not locNeighbours[n: n + 2].any(axis=1).all():
+ # not 2 neighbour points for triangle n
+ continue
+
+ j = np.where(locNeighbours[n])[0][0]
+ k = np.where(locNeighbours[n + 1])[0][0]
+
+ pointJ = pointsVicinityMap[j]
+ pointK = pointsVicinityMap[k]
+
+ # points bisecting pointI and closePoints
+ halfPointIJ = 0.5 * (pointI + pointJ)
+ halfPointIK = 0.5 * (pointI + pointK)
+ halfPointJK = 0.5 * (pointJ + pointK)
+
+ triangle = npTools.Triangles3D([pointI, pointJ, pointK])
+ centroidPoint = triangle.getCentroids()
+
+ pixelTriangleCornerPoints.extend([pointI, halfPointIJ, centroidPoint])
+ pixelTriangleCornerPoints.extend([pointI, halfPointIK, centroidPoint])
+ pixelTriangleCornerPoints.extend([pointJ, halfPointIJ, centroidPoint])
+ pixelTriangleCornerPoints.extend([pointJ, halfPointJK, centroidPoint])
+ pixelTriangleCornerPoints.extend([pointK, halfPointIK, centroidPoint])
+ pixelTriangleCornerPoints.extend([pointK, halfPointJK, centroidPoint])
+
+ if pixelTriangleScalars is not None:
+ pixelTriangleScalars.extend([self.scalars[i]] * 2)
+ pixelTriangleScalars.extend([self.scalars[j]] * 2)
+ pixelTriangleScalars.extend([self.scalars[k]] * 2)
+
+ return npTools.Triangles3D(pixelTriangleCornerPoints,
+ triangleScalars=pixelTriangleScalars)
+
+ else:
+ raise NotImplementedError("Method {0} not implemented.".format(method))
+
+ def plot(self, **kwargs):
+ scalarIndex = kwargs.pop('scalarIndex', 0)
+ if not (self.scalars.size == 0 or scalarIndex is None):
+ kwargs['c'] = kwargs.pop('c', self.scalars.T[scalarIndex])
+ return super(ScalarField3D, self).plot(**kwargs)
+
+
+if __name__ == '__main__':
+ import doctest
+ import sympy # NOQA: F401
+ doctest.testmod()
diff --git a/tfields/triangles3D.py b/tfields/triangles3D.py
new file mode 100644
index 0000000000000000000000000000000000000000..452c6c27ab0c63203239c5d1b05c7851df495e4d
--- /dev/null
+++ b/tfields/triangles3D.py
@@ -0,0 +1,831 @@
+#!/usr/bin/env
+# encoding: utf-8
+"""
+Author: Daniel Boeckenhoff
+Mail: daniel.boeckenhoff@ipp.mpg.de
+
+part of npTools library
+"""
+import numpy as np
+import warnings
+import npTools
+import ioTools
+import decoTools
+import pyTools
+import loggingTools
+
+logger = loggingTools.Logger(__name__)
+
+
+class Triangles3D(npTools.Points3D):
+ # pylint: disable=R0904
+ """
+ Points3D child restricted to n * 3 Points. 3 Points always group together to one triangle.
+ Examples:
+ >>> t = npTools.Triangles3D([[1,2,3], [3,3,3], [0,0,0]]); t
+ Triangles3D([[ 1., 2., 3.],
+ [ 3., 3., 3.],
+ [ 0., 0., 0.]])
+ >>> t.triangleScalars
+ array([], shape=(1, 0), dtype=float64)
+
+ >>> _ = t.coordinateTransform(t.CYLINDER)
+ >>> a = [npTools.Points3D([[1],[2],[3]], coordSys=npTools.Points3D.CYLINDER),
+ ... npTools.Points3D([[3],[3],[3]]), npTools.Points3D([[5],[1],[3]])]
+ >>> npTools.Triangles3D(a, coordSys=npTools.Points3D.CARTESIAN)
+ Triangles3D([[-0.41614684, 0.90929743, 3. ],
+ [ 3. , 3. , 3. ],
+ [ 5. , 1. , 3. ]])
+
+ You can add triangleScalars to each triangle
+ >>> t.appendScalars([1])
+
+ """
+ __slots__ = ['coordSys', 'triangleScalars', '_cache']
+
+ def __new__(cls, inputArray, triangleScalars=None, **kwargs):
+
+ if any([isinstance(a, Triangles3D) for a in inputArray]) and triangleScalars is None:
+ if not all([isinstance(a, Triangles3D) for a in inputArray]):
+ # TODO: could allow if all triangleScalars are none
+ raise TypeError("Copy constructor of Triangles3D only"
+ "accepts Triangles3D instances.")
+ if not len(set([o.triangleScalars.shape for o in inputArray])) == 1:
+ raise TypeError("Shape of triangleScalars that "
+ "should be concatenated must be the same.")
+ triangleScalars = np.concatenate([o.triangleScalars for o in inputArray], axis=0)
+
+ obj = npTools.Points3D.__new__(cls, inputArray, **kwargs)
+ if not obj.getNPoints() % 3 == 0:
+ raise TypeError("Points 3D does not describe triangles. "
+ "Size of obj is {0}, has no divider 3".format(len(obj)))
+
+ # Add triangleScalars
+ if triangleScalars is None:
+ triangleScalars = np.empty((len(obj) // 3, 0))
+ triangleScalars = np.array(triangleScalars)
+ if not len(triangleScalars) == obj.getNPoints() / 3:
+ raise ValueError("Length of triangleScalars has to be the same as"
+ "number of triangles ({0}) but is {1}."
+ .format(obj.getNPoints() / 3, len(triangleScalars)))
+ if len(triangleScalars.shape) == 1:
+ triangleScalars = triangleScalars.reshape(-1, 1)
+ obj.triangleScalars = triangleScalars
+
+ return obj
+
+ def saveTxt(self, filePath):
+ self.getMesh3D().save(filePath)
+
+ @classmethod
+ def createFromTxtFile(cls, filePath):
+ with ioTools.TextFile(filePath, 'r') as f:
+ matrix = f.process(seperator=' ', lineBreak='\n')
+ parsedMatrix = []
+ for line in matrix:
+ if line[0].startswith('#'):
+ continue
+ if '' in line:
+ pyTools.array.removeValues(line, '')
+ if not line:
+ continue
+ elif len(line) == 3:
+ parsedMatrix.append(map(float, line))
+ elif len(line) == 1:
+ parsedMatrix.append([eval(line[0])])
+ else:
+ raise ValueError("Length of line is not 1 or 3 but %s" % len(line))
+
+ infoPeriodicity = 4
+ scalars = [x[0] for i, x in enumerate(parsedMatrix) if i % infoPeriodicity == 0]
+ pyTools.array.removePositions(parsedMatrix,
+ [i for i in range(len(parsedMatrix))
+ if i % infoPeriodicity == 0])
+
+ return cls(parsedMatrix, triangleScalars=scalars)
+
+ def convertNaN(self, value=0.):
+ super(Triangles3D, self).convertNaN(value)
+ nanIndicesScalars = np.isnan(self.triangleScalars)
+ self.triangleScalars[nanIndicesScalars] = value
+
+ def appendScalars(self, scalars):
+ """
+ Examples:
+ >>> t = npTools.Triangles3D([[1,2,3], [3,3,3], [0,0,0]], [[1001]])
+ >>> t.appendScalars([42])
+ >>> t.triangleScalars
+ array([[1001, 42]])
+
+ """
+ self.triangleScalars = np.concatenate([self.triangleScalars, np.array([scalars]).T], 1)
+
+ def getScalarArrays(self):
+ """
+ Returns:
+ np.array: transposed scalar arrays
+ """
+ return self.triangleScalars.T
+
+ def getNTriangles(self):
+ """
+ Returns:
+ int: number of triangles
+ """
+ return self.getNPoints() // 3
+
+ def getMask(self, cutExpression, coordSys=None):
+ """
+ Triangle3D implementation
+ """
+ mask = super(Triangles3D, self).getMask(cutExpression, coordSys=coordSys)
+ mask = mask.reshape((self.getNTriangles(), 3))
+ mask = mask.all(axis=1)
+ mask = np.array([mask] * 3).T.reshape((self.getNPoints()))
+ return mask
+
+ def cut(self, cutExpression, coordSys=None):
+ """
+ Default cut method for Triangles3D
+ Examples:
+ >>> x, y, z = sympy.symbols('x y z')
+ >>> t = npTools.Triangles3D([[1., 2., 3.], [-4., 5., 6.], [1, 2, -6],
+ ... [5, -5, -5], [1,0,-1], [0,1,-1], [-5, -5, -5],
+ ... [1,0,-1], [0,1,-1]])
+ >>> tc = t.cut(x >= 0)
+ >>> tc
+ Triangles3D([[ 5., -5., -5.],
+ [ 1., 0., -1.],
+ [ 0., 1., -1.]])
+ >>> tc.triangleScalars
+ array([], shape=(1, 0), dtype=float64)
+ >>> t.appendScalars([1,2,3])
+ >>> tc2 = t.cut(x >= 0)
+ >>> tc2.triangleScalars
+ array([[ 2.]])
+
+ """
+ mask = self.getMask(cutExpression, coordSys=coordSys)
+
+ inst = self[mask].copy()
+ inst.triangleScalars = inst.triangleScalars[mask.reshape((self.getNTriangles(), 3)).T[0]]
+
+ return inst
+
+ @decoTools.cached_property()
+ def areas(self):
+ """
+ Cached method to retrieve areas of triangles
+ """
+ indices = range(self.getNTriangles())
+ aIndices = [i * 3 for i in indices]
+ bIndices = [i * 3 + 1 for i in indices]
+ cIndices = [i * 3 + 2 for i in indices]
+
+ # define 3 vectors building the triangle and take their norm
+ a = np.linalg.norm(self[aIndices, :] - self[bIndices, :], axis=1)
+ b = np.linalg.norm(self[aIndices, :] - self[cIndices, :], axis=1)
+ c = np.linalg.norm(self[bIndices, :] - self[cIndices, :], axis=1)
+ # sort by length for numerical stability
+ lengths = np.concatenate((a.reshape(-1, 1), b.reshape(-1, 1), c.reshape(-1, 1)), axis=1)
+ lengths.sort()
+ a, b, c = lengths.T
+
+ return 0.25 * np.sqrt((a + (b + c)) * (c - (a - b)) * (c + (a - b)) * (a + (b - c)))
+
+ def getAreas(self):
+ """
+ Calculate area with "heron's formula"
+ Examples:
+ >>> m = npTools.Mesh3D([[1,0,0], [0,0,1], [0,0,0]], [[0, 1, 2]]);
+ >>> m.triangles.getAreas()
+ array([ 0.5])
+
+ >>> m = npTools.Mesh3D([[1,0,0], [0,1,0], [0,0,0], [0,0,1]], [[0, 1, 2], [1, 2, 3]]);
+ >>> m.triangles.getAreas()
+ array([ 0.5, 0.5])
+
+ >>> m = npTools.Mesh3D([[1,0,0], [0,1,0], [1,1,0], [0,0,1], [1,0,1]],
+ ... [[0, 1, 2], [0, 3, 4]]);
+ >>> m.triangles.getAreas()
+ array([ 0.5, 0.5])
+
+ """
+ return self.areas
+
+ def getAreasWithTransformation(self, transformationMatrix):
+ """
+ Calculate area with "heron's formula".
+ The triangle points are transformed with transformationMatrix
+ """
+ indices = range(self.getNTriangles())
+ aIndices = [i * 3 for i in indices]
+ bIndices = [i * 3 + 1 for i in indices]
+ cIndices = [i * 3 + 2 for i in indices]
+
+ # define 3 vectors building the triangle, transform it back and take their norm
+
+ a = np.linalg.norm(np.linalg.solve(transformationMatrix.T, (self[aIndices, :] - self[bIndices, :]).T), axis=0)
+ b = np.linalg.norm(np.linalg.solve(transformationMatrix.T, (self[aIndices, :] - self[cIndices, :]).T), axis=0)
+ c = np.linalg.norm(np.linalg.solve(transformationMatrix.T, (self[bIndices, :] - self[cIndices, :]).T), axis=0)
+
+ # sort by length for numerical stability
+ lengths = np.concatenate((a.reshape(-1, 1), b.reshape(-1, 1), c.reshape(-1, 1)), axis=1)
+ lengths.sort()
+ a, b, c = lengths.T
+
+ return 0.25 * np.sqrt((a + (b + c)) * (c - (a - b)) * (c + (a - b)) * (a + (b - c)))
+
+ def getCorners(self):
+ """
+ Returns:
+ three np.arrays with corner points of triangles
+ """
+ indices = range(self.getNTriangles())
+ aIndices = [i * 3 for i in indices]
+ bIndices = [i * 3 + 1 for i in indices]
+ cIndices = [i * 3 + 2 for i in indices]
+
+ a = self[aIndices, :]
+ b = self[bIndices, :]
+ c = self[cIndices, :]
+ return a, b, c
+
+ def getCircumcenters(self):
+ """
+ Semi baricentric method to calculate circumcenter points of the triangles
+
+ Examples:
+ >>> m = npTools.Mesh3D([[0,0,0], [1,0,0], [-1,0,0], [0,1,0], [0,0,1]],
+ ... [[0, 1, 3],[0, 2, 3],[1,2,4], [1, 3, 4]]);
+ >>> m.triangles.getCircumcenters()
+ Points3D([[ 5.00000000e-01, 5.00000000e-01, 0.00000000e+00],
+ [ -5.00000000e-01, 5.00000000e-01, 0.00000000e+00],
+ [ 0.00000000e+00, 0.00000000e+00, 2.22044605e-16],
+ [ 3.33333333e-01, 3.33333333e-01, 3.33333333e-01]])
+
+ """
+ pointA, pointB, pointC = self.getCorners()
+ a = np.linalg.norm(pointC - pointB, axis=1) # side length of side opposite to pointA
+ b = np.linalg.norm(pointC - pointA, axis=1)
+ c = np.linalg.norm(pointB - pointA, axis=1)
+ bary1 = a ** 2 * (b ** 2 + c ** 2 - a ** 2)
+ bary2 = b ** 2 * (a ** 2 + c ** 2 - b ** 2)
+ bary3 = c ** 2 * (a ** 2 + b ** 2 - c ** 2)
+ matrices = np.concatenate((pointA, pointB, pointC), axis=1).reshape(pointA.shape + (3,))
+ # transpose the inner matrix
+ matrices = np.einsum('...ji', matrices)
+ vectors = np.array((bary1, bary2, bary3)).T
+ # matrix vector product for matrices and vectors
+ P = np.einsum('...ji,...i', matrices, vectors)
+ P /= vectors.sum(axis=1).reshape((len(vectors), 1))
+ return npTools.Points3D(P)
+
+ @decoTools.cached_property()
+ def centroids(self):
+ """
+ this version is faster but takes much more ram also.
+ So much that i get memory error with a 32 GB RAM
+ """
+ nT = self.getNTriangles()
+ mat = np.ones((1, 3)) / 3.0
+ # matrix product calculatesq center of all triangles
+ return npTools.Points3D(np.dot(mat, self.reshape(nT, 3, 3))[0], coordSys=self.coordSys)
+
+ """
+ Old version:
+ pointA, pointB, pointC = self.getCorners()
+ return Points3D(1. / 3 * (pointA + pointB + pointC)), coordSys=self.coordSys
+ This versioin was slightly slower (110 % of new version)
+ Memory usage of new version is better for a factor of 4 or so.
+ Not really reliable method of measurement
+ """
+
+ def getCentroids(self):
+ """
+ Examples:
+ >>> m = npTools.Mesh3D([[0,0,0], [1,0,0], [-1,0,0], [0,1,0], [0,0,1]],
+ ... [[0, 1, 3],[0, 2, 3],[1,2,4], [1, 3, 4]]);
+ >>> m.triangles.getCentroids()
+ Points3D([[ 0.33333333, 0.33333333, 0. ],
+ [-0.33333333, 0.33333333, 0. ],
+ [ 0. , 0. , 0.33333333],
+ [ 0.33333333, 0.33333333, 0.33333333]])
+
+ """
+ return self.centroids
+
+ def getTriangleCenters(self):
+ """
+ Deprecated! use centroids
+ """
+ warnings.warn("deprecated", DeprecationWarning)
+ return self.centroids
+
+ def getEdgeVectors(self):
+ """
+ Return two vectors that define two triangle edges
+ """
+ a, b, c = self.getCorners()
+
+ ab = b - a
+ ac = c - a
+
+ return ab, ac
+
+ def getNormVectors(self):
+ """
+ Examples:
+ >>> m = npTools.Mesh3D([[0,0,0], [1,0,0], [-1,0,0], [0,1,0], [0,0,1]],
+ ... [[0, 1, 3],[0, 2, 3],[1,2,4], [1, 3, 4]]);
+ >>> m.triangles.getNormVectors()
+ array([[ 0. , 0. , 1. ],
+ [ 0. , 0. , -1. ],
+ [ 0. , 1. , 0. ],
+ [ 0.57735027, 0.57735027, 0.57735027]])
+
+ """
+ ab, ac = self.getEdgeVectors()
+ vectors = np.cross(ab, ac)
+ norms = np.apply_along_axis(np.linalg.norm, 0, vectors.T).reshape(-1, 1)
+ # cross product may be zero, so replace zero norms by ones to divide vectors by norms
+ np.place(norms, norms == 0.0, 1.0)
+ return vectors / norms
+
+ def getBaricentricCoordinates(self, point, delta=0.):
+ """
+ Determine baricentric coordinates like
+
+ [u,v,w] = [ab, ac, ab x ac]^-1 * ap
+ where ax is vector from point a to point x
+
+ return baricentric coordinates of triangle
+ last coordinate is prism
+ Examples:
+ empty Meshes return right formatted array
+ >>> m = npTools.Mesh3D([], faces=[])
+ >>> m.triangles.getBaricentricCoordinates(np.array([0.2, 0.2, 0]))
+ array([], dtype=float64)
+
+ >>> m2 = npTools.Mesh3D([[0,0,0], [2,0,0], [0,2,0], [0,0,2]], [[0, 1, 2], [0, 2, 3]]);
+ >>> m2.triangles.getBaricentricCoordinates(np.array([0.2, 0.2, 0]), delta=2.)
+ array([[ 0.1, 0.1, 0. ],
+ [ 0.1, 0. , 0.1]])
+
+ if point lies in the plane, return nan, else inf for w if delta is exactly 0.
+ >>> m2.triangles.getBaricentricCoordinates(np.array([0.2, 0.2, 0]), delta=0.)
+ array([[ 0.1, 0.1, nan],
+ [ 0.1, 0. , inf]])
+
+ If you define triangles that have colinear side vectors or in general lead to
+ not invertable matrices [ab, ac, ab x ac] a warning wil be thrown but no error
+ >>> m3 = npTools.Mesh3D([[0,0,0], [2,0,0], [4,0,0], [0,1,0]], [[0, 1, 2], [0, 1, 3]]);
+ >>> import pytest
+ >>> with pytest.warns(UserWarning) as record:
+ ... bc = m3.triangles.getBaricentricCoordinates(np.array([0.2, 0.2, 0]), delta=0.3)
+ >>> bc
+ array([[ nan, nan, nan],
+ [ 0.1, 0.2, 0. ]])
+ >>> len(record)
+ 1
+ >>> "Could not invert matrix" in str(record[0].message)
+ True
+
+ Returns:
+ arrays of u, v, w for each triangle
+ """
+ if self.getNTriangles() == 0:
+ return np.array([])
+
+ try:
+ point = np.reshape(point, 3)
+ except ValueError:
+ raise ValueError("point must be castable to shape 3 but is of shape {0}"
+ .format(point.shape))
+ except:
+ raise
+
+ a, b, c = self.getCorners()
+
+ ap = point - a
+ # matrix vector product for matrices and vectors
+ barCoords = np.einsum('...ji,...i', self.barycentricMatrix, ap)
+ with warnings.catch_warnings():
+ warnings.filterwarnings('ignore', message="invalid value encountered in divide")
+ warnings.filterwarnings('ignore', message="divide by zero encountered in divide")
+ barCoords[:, 2] /= delta # allow devide by 0.
+ return barCoords
+
+ @decoTools.cached_property()
+ def barycentricMatrix(self):
+ """
+ cached barycentric matrix for faster calculations
+ """
+ a, b, c = self.getCorners()
+ ab = b - a
+ ac = c - a
+
+ # get norm vector
+ norm = np.cross(ab, ac)
+ normLen = np.linalg.norm(norm, axis=1)
+ normLen = normLen.reshape((1,) + normLen.shape)
+ colinearMask = (normLen == 0)
+ with warnings.catch_warnings():
+ warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning)
+ # prevent divide by 0
+ colinSlice = np.where(~colinearMask.T)
+ norm[np.where(~colinearMask.T)] = \
+ norm[np.where(~colinearMask.T)] / normLen.T[np.where(~colinearMask.T)]
+
+ matrix = np.concatenate((ab, ac, norm), axis=1).reshape(ab.shape + (3,))
+ matrix = np.einsum('...ji', matrix) # transpose the inner matrix
+
+ # invert matrix if possible
+ # matrixI = np.linalg.inv(matrix) # one line variant without exception
+ matrixI = []
+ for mat in matrix:
+ try:
+ matrixI.append(np.linalg.inv(mat))
+ except np.linalg.linalg.LinAlgError as e:
+ if str(e) == "Singular matrix":
+ warnings.warn("Singular matrix: Could not invert matrix "
+ "{0}. Return nan matrix."
+ .format(str(mat).replace('\n', ',')))
+ matrixI.append(np.full((3, 3), np.nan))
+ return np.array(matrixI)
+
+ def pointsInTriangles(self, points3D, delta=0., method='baryc', assignMultiple=False):
+ """
+ Args:
+ points3D (Points3D instance)
+
+ optional:
+ delta (float): normal distance to a triangle, that the points
+ is concidered to be contained in the triangle.
+ method (str): choose, which method to use for the calculation.
+ baryc is recommended
+ assignMultiple (bool): if True, one point may belong to multiple
+ triangles at the same time. In the other case the first
+ occurence will be True the other False
+
+ Returns:
+ 2-d numpy mask which is True, where a point is in a face
+ face ind ->
+ points index
+ 1 0 0
+ | 1 0 0
+ V 0 0 1
+
+ if assignMultiple == False:
+ there is just one True for each points index
+ face indices can have multiple true values
+
+ """
+ log = logger.new()
+ log.verbose("Checking for points in Triangles")
+ try:
+ points3D = points3D.reshape((-1, 3))
+ except ValueError:
+ raise ValueError("points3D must be castable to shape (nPoints, 3) "
+ "but is of shape {0}".format(points3D.shape))
+ except:
+ raise
+
+ if self.getNTriangles() == 0:
+ return np.empty((points3D.shape[0], 0), dtype=bool)
+
+ masks = np.zeros((points3D.getNPoints(), self.getNTriangles()), dtype=bool)
+ for i, point in enumerate(points3D):
+ if i % 1000 == 0:
+ log.verbose("Status: {i} points processed. Yet {nP} are found "
+ "to be in triangles.".format(nP=masks.sum() - 1, **locals()))
+ masks[i] = self.pointInTriangles(point, delta, method=method)
+ log.verbose("Status: All points processed. Yet {nP} are found "
+ "to be in triangles.".format(nP=masks.sum() - 1, **locals()))
+
+ if len(masks) == 0:
+ return masks
+
+ if not assignMultiple:
+ """
+ index of first faces containing the point for each point. This gets the
+ first membership index always. ignoring multiple triangle memberships
+ """
+ faceMembershipIndices = np.argmax(masks, axis=1)
+ # True if point lies in any triangle
+ membershipBools = (masks.sum(axis=1) != 0)
+
+ masks = np.full(masks.shape, False, dtype=bool)
+ for j, valid in enumerate(membershipBools):
+ if valid:
+ masks[j, faceMembershipIndices[j]] = True
+ """
+ masks is now a matrix:
+ face ind ->
+ points index
+ 1 0 0
+ | 1 0 0
+ V 0 0 1
+
+ there is just one True for each points index
+ face indices can have multiple true values
+ """
+ return masks
+
+ def pointInTriangles(self, point, delta=0., method='baryc'):
+ """
+ Returns:
+ np.array: boolean mask, True where point in a triangle within delta
+ Args:
+ point (list of len 3)
+ delta (float): acceptance in +- norm vector direction
+ """
+ if method == 'baryc':
+ return self.pointInTrianglesBaryc(point, delta=delta)
+ elif method == 'baryc2':
+ return self.pointInTrianglesBarycentric2(point, delta=delta)
+ else:
+ raise NotImplementedError
+
+ def pointInTrianglesBaryc(self, point, delta=0.):
+ """
+ Examples:
+ >>> m = npTools.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], [[0, 1, 2]]);
+ >>> m.triangles.pointInTrianglesBaryc(npTools.Points3D([[0.2, 0.2, 0]]))
+ array([ True], dtype=bool)
+
+ wrong format of point will raise a ValueError
+ >>> m.triangles.pointInTrianglesBaryc(npTools.Points3D([[0.2, 0.2, 0], [0.2, 0.2, 0]]))
+ Traceback (most recent call last):
+ ...
+ ValueError: point must be castable to shape 3 but is of shape (2, 3)
+
+ All Triangles are tested
+ >>> m2 = npTools.Mesh3D([[1,0,0], [0,1,0], [0,0,0], [4,0,0], [4, 4, 0], [8, 0, 0]],
+ ... [[0, 1, 2], [3, 4, 5]]);
+
+ >>> m2.triangles.pointInTrianglesBaryc(np.array([0.2, 0.2, 0]))
+ array([ True, False], dtype=bool)
+ >>> m2.triangles.pointInTrianglesBaryc(np.array([5, 2, 0]))
+ array([False, True], dtype=bool)
+
+ delta allows to accept points that lie within delta orthogonal to the tringle plain
+ >>> m2.triangles.pointInTrianglesBaryc(np.array([0.2, 0.2, 9000]), 0.0)
+ array([False, False], dtype=bool)
+ >>> m2.triangles.pointInTrianglesBaryc(np.array([0.2, 0.2, 0.1]), 0.2)
+ array([ True, False], dtype=bool)
+
+ If you define triangles that have colinear side vectors or in general lead to
+ not invertable matrices the you will always get False
+ >>> m3 = npTools.Mesh3D([[0,0,0], [2,0,0], [4,0,0], [0,1,0]], [[0, 1, 2], [0, 1, 3]]);
+ >>> import pytest
+ >>> with pytest.warns(UserWarning) as record:
+ ... mask = m3.triangles.pointInTrianglesBaryc(np.array([0.2, 0.2, 0]), delta=0.3)
+ >>> mask
+ array([False, True], dtype=bool)
+
+ """
+ if self.getNTriangles() == 0:
+ return np.array([], dtype=bool)
+
+ u, v, w = self.getBaricentricCoordinates(point, delta=delta).T
+
+ if delta == 0.:
+ w[np.isnan(w)] = 0. # division by 0 in getBaricentricCoordinates makes w = 0 nan.
+
+ with warnings.catch_warnings():
+ warnings.filterwarnings('ignore', message="invalid value encountered in less_equal")
+ orthogonalAcceptance = (abs(w) <= 1)
+ barycentricBools = ((v <= 1.) & (v >= 0.)) & ((u <= 1.) & (u >= 0.)) & ((v + u <= 1.0))
+ return np.array(barycentricBools & orthogonalAcceptance)
+
+ def pointInTrianglesBarycentric2(self, point, delta=0.):
+ """
+ Check whether point lies within triangles with Barycentric Technique
+ It could lie within multiple triangles!
+ Examples:
+ >>> m = npTools.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], [[0, 1, 2]]);
+ >>> m.triangles.pointInTrianglesBarycentric2(npTools.Points3D([[0.2, 0.2, 0]]))
+ array([ True], dtype=bool)
+
+ wrong format of point will raise a ValueError
+ >>> m.triangles.pointInTrianglesBarycentric2(npTools.Points3D([[0.2, 0.2, 0],
+ ... [0.2, 0.2, 0]]))
+ Traceback (most recent call last):
+ ...
+ ValueError: point must be castable to shape 3 but is of shape (2, 3)
+
+ All Triangles are tested
+ >>> m2 = npTools.Mesh3D([[1,0,0], [0,1,0], [0,0,0], [4,0,0],
+ ... [4, 4, 0], [8, 0, 0]],
+ ... [[0, 1, 2], [3, 4, 5]]);
+ >>> m2.triangles.pointInTrianglesBarycentric2(np.array([0.2, 0.2, 0]))
+ array([ True, False], dtype=bool)
+ >>> m2.triangles.pointInTrianglesBarycentric2(np.array([5, 2, 0]))
+ array([False, True], dtype=bool)
+
+ delta allows to accept points that lie within delta orthogonal to the tringle plain
+ >>> m2.triangles.pointInTrianglesBarycentric2(np.array([0.2, 0.2, 9000]), 0.0)
+ array([False, False], dtype=bool)
+ >>> m2.triangles.pointInTrianglesBarycentric2(np.array([0.2, 0.2, 0.1]), 0.2)
+ array([ True, False], dtype=bool)
+
+ Returns:
+ array of bools for each triangle
+ """
+ try:
+ point = np.reshape(point, 3)
+ except ValueError:
+ raise ValueError("point must be castable to shape 3 but is of shape {0}"
+ .format(point.shape))
+ except:
+ raise
+
+ indices = range(self.getNTriangles())
+ aIndices = [i * 3 for i in indices]
+ bIndices = [i * 3 + 1 for i in indices]
+ cIndices = [i * 3 + 2 for i in indices]
+
+ # compute base vectors
+ AP = (self[aIndices, :] - point).T
+ AB = (self[aIndices, :] - self[bIndices, :]).T
+ AC = (self[aIndices, :] - self[cIndices, :]).T
+
+ # compute barycentric coordinates
+ np.seterr(divide='ignore', invalid='ignore')
+
+ denom = AB[1] * AC[0] - AB[0] * AC[1]
+ v = (AB[1] * AP[0] - AB[0] * AP[1]) / denom
+ u = -(AC[1] * AP[0] - AC[0] * AP[1]) / denom
+ # u = (AP[0] - v * AC[0]) / AB[0] # also true, but problematic with AB[0]
+
+ # calculate bools
+ orthogonalAcceptance = abs(AP[2] - u * AB[2] - v * AC[2]) <= delta
+ barycentricBools = ((v <= 1.) & (v >= 0.)) & ((u <= 1.) & (u >= 0.)) & ((v + u <= 1.0))
+ np.seterr(divide='warn', invalid='warn')
+ return np.array(barycentricBools & orthogonalAcceptance)
+
+ def pointOnTriangleSide(self, point):
+ """
+ Returns:
+ np.array: boolean mask which is true, if point is on one side ray
+ of a triangle
+ Examples:
+ >>> m = npTools.Mesh3D([[0,0,0], [1,0,0], [-1,0,0], [0,1,0], [0,0,1]],
+ ... [[0, 1, 3],[0, 2, 3],[1,2,4]]);
+
+ Corner points are found
+ >>> m.triangles.pointOnTriangleSide(npTools.Points3D([[0,1,0]]))
+ array([ True, True, False], dtype=bool)
+
+ Side points are found, too
+ >>> m.triangles.pointOnTriangleSide(npTools.Points3D([[0.5,0,0.5]]))
+ array([False, False, True], dtype=bool)
+
+ """
+ u, v, w = self.getBaricentricCoordinates(point, 1.).T
+
+ orthogonalAcceptance = (w == 0) # point should lie in triangle
+ barycentricBools = (((0. <= v) & (v <= 1.)) & (u == 0.)) | \
+ (((0. <= u) & (u <= 1.)) & (v == 0.)) | \
+ (v + u == 1.0)
+ return np.array(barycentricBools & orthogonalAcceptance)
+
+ def getWeightedAreas(self, weightScalarIndex=0):
+ """
+ Args:
+ weightScalarIndex (int): index to scalars to weight with
+ """
+ # set weights to 1.0 if weightScalarIndex is None
+ if weightScalarIndex is None:
+ weights = np.ones(self.getNTriangles())
+ elif self.getScalarArrays().shape[0] == 0:
+ weights = np.ones(self.getNTriangles())
+ else:
+ weights = self.getScalarArrays()[weightScalarIndex]
+ areas = self.getAreas()
+ return areas * weights
+
+ def getNormedWeightedAreas(self, weightScalarIndex=0, transform=None):
+ """
+ Args:
+ weightScalarIndex (int): index to scalars to weight with
+ Returns:
+ A_i * w_i / sum_{i=0}^{N}(A_i * w_i)
+ """
+ # set weights to 1.0 if weightScalarIndex is None
+ if transform is None:
+ transform = np.eye(3)
+ if weightScalarIndex is None:
+ weights = np.ones(self.getNTriangles())
+ elif self.getScalarArrays().shape[0] == 0:
+ weights = np.ones(self.getNTriangles())
+ else:
+ weights = self.getScalarArrays()[weightScalarIndex]
+ areas = self.getAreasWithTransformation(transform)
+ return areas * weights / np.dot(areas, weights)
+
+ def getMean(self, weightScalarIndex=0):
+ """
+ Args:
+ weightScalarIndex (int): index to scalars to weight with
+ Returns:
+ mean of the centroids weighted by scalars @ <weightScalarIndex>
+ """
+ centroids = self.getCentroids()
+ weightedAreas = self.getNormedWeightedAreas(weightScalarIndex=weightScalarIndex)
+
+ # sum up all triangle centers weighted by the area and divided by the total area
+ return np.dot(centroids.T, weightedAreas)
+
+ def getStd(self, mean=None, weightScalarIndex=0):
+ """
+ Args:
+ mean(Points3D)
+ weightScalarIndex (int): index to scalars to weight with
+ """
+ if mean is None:
+ mean = self.getMean(weightScalarIndex=weightScalarIndex)
+ centroids = self.getCentroids()
+ weightedAreas = self.getNormedWeightedAreas(weightScalarIndex=weightScalarIndex)
+ squaredDiff = (centroids - mean) ** 2
+
+ var = np.dot(squaredDiff.T, weightedAreas)
+ return np.sqrt(var)
+
+ def getDirection(self, weightScalarIndex=0, transform=None):
+ """
+ Returns: Eigenvector corresponding to greatest eigenvalue multiplied by greatest eigenvalue
+ """
+ # calculate covariance matrix
+ log = logger.new()
+ try:
+ cov = np.cov(self.getCentroids().T,
+ ddof=0,
+ aweights=self.getNormedWeightedAreas(weightScalarIndex=weightScalarIndex, transform=transform))
+ # calculate eigenvalues and eigenvectors of covariance
+
+ evals, evecs = np.linalg.eigh(cov)
+ except:
+ log.warning("np.linalg.eig() failed. use zeros instead!")
+ evals = np.zeros(3)
+ evecs = np.zeros((3, 3))
+ idx = np.argmax(evals)
+ # return eigenvector corresponding to greatest eigenvalue multiplied by greatest eigenvalue
+ # always positive in x-direction
+ return evals.max() * (evecs[:, idx] if evecs[0, idx] > 0 else -evecs[:, idx])
+
+ def getCovEig(self, weightScalarIndex=0):
+ """
+ Args:
+ weightScalarIndex (int): index to scalars to weight with
+ """
+ cov = np.cov(self.getCentroids().T,
+ ddof=0,
+ aweights=self.getNormedWeightedAreas(weightScalarIndex=weightScalarIndex))
+ # calculate eigenvalues and eigenvectors of covariance
+ evals, evecs = np.linalg.eigh(cov)
+ idx = evals.argsort()[::-1]
+ evals = evals[idx]
+ evecs = evecs[:, idx]
+ e = np.concatenate((evecs, evals.reshape(1, 3)))
+ return e.T.reshape(12, )
+
+ def getMainAxes(self):
+ """
+ Args:
+ weightScalarIndex (int): index to scalars to weight with
+ Returns:
+ Main Axes eigen-vectors
+ """
+ mean = self.getMean(None)
+ relativeCentroids = self.getCentroids() - mean
+ x = relativeCentroids
+ cov = np.cov(x.T,
+ ddof=0,
+ aweights=self.getNormedWeightedAreas(weightScalarIndex=None))
+ # calculate eigenvalues and eigenvectors of covariance
+ evals, evecs = np.linalg.eigh(cov)
+ return (evecs * evals.T).T
+
+ def getMesh3D(self):
+ """
+ Forward to an alternative representation of the triangeles as mesh
+ """
+ return npTools.Mesh3D(self,
+ [[3 * i, 3 * i + 1, 3 * i + 2]
+ for i in range(self.getNTriangles())],
+ faceScalars=self.triangleScalars)
+
+ def plot(self, **kwargs):
+ """
+ Forward to Mesh3D plot
+ """
+ return self.getMesh3D().plot(**kwargs)
+
+
+if __name__ == '__main__':
+ import doctest
+ import sympy # NOQA: F401
+
+ doctest.testmod()
diff --git a/tfields/vectorField3D.py b/tfields/vectorField3D.py
new file mode 100644
index 0000000000000000000000000000000000000000..468ee5278f6599313728e899fb1f436936fd5b05
--- /dev/null
+++ b/tfields/vectorField3D.py
@@ -0,0 +1,160 @@
+#!/usr/bin/env
+# encoding: utf-8
+"""
+Author: Daniel Boeckenhoff
+Mail: daniel.boeckenhoff@ipp.mpg.de
+
+part of npTools library
+"""
+import numpy as np
+import npTools
+import mplTools as mpt
+import loggingTools
+logger = loggingTools.Logger(__name__)
+
+
+class VectorField3D(object):
+ """
+ Container class for two points3D instances points and vectors
+ Attributes:
+ points
+ vectors
+ """
+ def __init__(self, points, vectors, **kwargs):
+ """
+ Examples:
+ >>> p = npTools.Points3D([[1, 1],[2, 2],[3, 3]], coordSys=npTools.Points3D.CYLINDER)
+ >>> v = npTools.Points3D([[4, 3],[5, 3],[6, 3]])
+ >>> vf = npTools.VectorField3D(p, v, coordSys=npTools.Points3D.CARTESIAN)
+ >>> vf
+ points: [[-0.41614684 0.90929743 3. ]
+ [-0.41614684 0.90929743 3. ]], vectors: [[ 4. 5. 6.]
+ [ 3. 3. 3.]]
+ >>> vf2 = npTools.VectorField3D(p, v, coordSys=npTools.Points3D.CYLINDER)
+ >>> vf2 == vf
+ True
+
+ """
+ coordSys = kwargs.pop('coordSys', None)
+ if not isinstance(points, npTools.Points3D):
+ self.points = npTools.Points3D(points, coordSys=coordSys)
+ else:
+ self.points = points.copy()
+ if coordSys is None:
+ coordSys = points.getCoordSys()
+ if not isinstance(vectors, npTools.Points3D):
+ self.vectors = npTools.Points3D(vectors, coordSys=coordSys)
+ else:
+ self.vectors = vectors.copy()
+ if coordSys is None:
+ coordSys = vectors.getCoordSys()
+ if coordSys is not None: # meaning a coordinate system has been specified
+ self.coordinateTransform(coordSys)
+ if not len(self.vectors) == len(self.points):
+ raise ValueError("Vectors and points must be of same dimension!")
+
+ @classmethod
+ def createTest(cls):
+ """
+ A test factory method
+ """
+ p = npTools.MeshGrid3D(3, 4, 2, 0, 2, 3, 6, 9, 4, iterOrder=[1, 0, 2],
+ coordSys=npTools.Points3D.CYLINDER)
+ v = npTools.MeshGrid3D(3, 4, 2, 0, 2, 3, 6, 9, 4, iterOrder=[1, 0, 2],
+ coordSys=npTools.Points3D.CYLINDER)
+ return cls(p, v)
+
+ def coordinateTransform(self, coordSys):
+ """
+ Args:
+ coordSys (str): see Points3D.
+ """
+ self.points.coordinateTransform(coordSys)
+ self.vectors.coordinateTransform(coordSys)
+
+ def mirrorCoordinate(self, coordinate, condition=None):
+ """
+ Args:
+ coordSys (str): see Points3D.
+ """
+ if condition is not None:
+ raise NotImplementedError("Mirror with condition not defined in vectorfield")
+ self.points.mirrorCoordinate(coordinate)
+ self.vectors.mirrorCoordinate(coordinate)
+
+ def changeIterOrder(self, iterOrder):
+ """
+ Set a new iterorder
+ """
+ if not issubclass(self.points, npTools.MeshGrid3D):
+ raise TypeError("Can only change iterOrder for MeshGrid3D")
+ self.points.changeIterOrder(iterOrder)
+ self.vectors.changeGridIterOrder(iterOrder, *self.points.baseLengths,
+ actualIterOrder=self.points.iterOrder)
+
+ def getCoordSys(self):
+ """
+ Returns:
+ str: coordSys
+ """
+ return self.points.getCoordSys()
+
+ def copy(self):
+ """
+ Returns:
+ copy of object
+ """
+ return self.__class__(self.points, self.vectors)
+
+ def __getitem__(self, index):
+ return self.points[index], self.vectors[index]
+
+ def __repr__(self):
+ return "points: %s, vectors: %s" % (self.points, self.vectors)
+
+ def __eq__(self, other):
+ if self.getCoordSys() != other.getCoordSys():
+ other = other.copy()
+ other.coordinateTransform(self.getCoordSys())
+ return bool(np.all(np.isclose(self.points, other.points)) and
+ np.all(np.isclose(self.vectors, other.vectors)))
+
+ def getVectors(self):
+ """
+ Returns:
+ vectors
+ """
+ return self.vectors.copy()
+
+ def getPoints(self):
+ """
+ Returns:
+ point bases
+ """
+ return self.points.copy()
+
+ def evalAtPoint(self, point):
+ """
+ Examples:
+ >>> vf = npTools.VectorField3D.createTest()
+ >>> vf.evalAtPoint(npTools.Points3D([[3, 0, 6]]))
+ MeshGrid3D([[ 3., 0., 6.]])
+ >>> vf.evalAtPoint([3, 1, 6])
+ MeshGrid3D([[ 3., 1., 6.]])
+ >>> vf.evalAtPoint([[0, 0, 0]])
+ MeshGrid3D([], shape=(0, 3), dtype=float64)
+
+ """
+ return self.vectors[np.equal(self.points, point).all(axis=1)]
+
+ def plot(self, **kwargs):
+ """
+ Frowarding to plotTools.plotVectorField
+ """
+ artist = mpt.plotVectorField(self.points, self.vectors, **kwargs)
+ return artist
+
+
+if __name__ == '__main__':
+ import doctest
+ doctest.testmod()