merging

FME/dsi_helper.pyx

@@ -208,10 +208,12 @@ def cg(double [:,:,:] EG, long [:,:] neighbours, long [:,:] elements,double [:,:
     cdef double [:,:] e2
     cdef long [:] idl  = np.zeros(4,dtype=np.int64) 
     cdef long [:] idr = np.zeros(4,dtype=np.int64) 
+    #loop over elements
     for e in range(ne):
         idl = elements[e,:]
         e1 = EG[e,:,:]
         flag[e] = 1
+        #loop over neighbours
         for n in range(4):
 
             neigh = neighbours[e,n]

FME/tet_mesh.py

@@ -454,64 +454,99 @@ def plot_mesh(self,propertyname,cmap=None):
             interpolate_before_map=True
         )
         p.show()
-    def lv_plot_isosurface(self,propertyname,isovalue,
-            colour='green',
-            name='IsoSurface',
-            interactive=False,
-            lv=None,
-            draw=True,
-            region=None
-            ):
+    def lv_plot_isosurface(self,propertyname,lv,**kwargs):
         #import lavavu in case its not imported
         try:
             import lavavu  #visualisation library   
         except ImportError:
             print("Cannot import Lavavu")
             return 
-        ##run the marching tetra algorithm        
-        reg = np.zeros(self.properties[propertyname].shape).astype(bool)
-        reg[:] = True
-        if region is not None:
-            reg = self.regions[region]
-        tri, ntri = marching_tetra(isovalue,self.elements,self.nodes,reg,self.properties[propertyname])
+        #if no isovalue is specified just use the average
+        property = self.properties[propertyname]
+        slices = [np.mean(property)]
+        colour='red'
+        if 'isovalue' in kwargs:
+            slices = [kwargs['isovalue']]
+        if 'slices' in kwargs:
+            slices = kwargs['slices']
+        if 'nslices' in kwargs:
+            slices = np.linspace(np.min(property),np.max(property),kwargs['nslices'])
+        if 'colour' in kwargs:
+            colour=kwargs['colour']
+        for isovalue in slices:
+            if isovalue < np.min(property) or isovalue > np.max(property):
+                print("No surface to create for isovalue")
+                isovalue=kwargs['isovalue']
+            reg = np.zeros(self.properties[propertyname].shape).astype(bool)
+            reg[:] = True
+            if 'region' in kwargs:
+                reg = kwargs['region']
+            name = propertyname+'_%f'%isovalue
+            if 'name' in kwargs:
+                name = kwargs['name']
+
+            tri, ntri = marching_tetra(isovalue,self.elements,self.nodes,reg,self.properties[propertyname])
+
+            ##convert from memoryview to np array
+            tri = np.array(tri)
+            ntri = np.array(ntri)[0]
+            ##create a triangle indices array and initialise to -1
+            tris = np.zeros((ntri,3)).astype(int)
+            tris[:,:] = -1
+            ##create a dict for storing nodes index where the key is the node as as a tuple. 
+            #A dict is preferable because it is very quick to check if a key exists
+            #assemble arrays for unique vertex and triangles defined by vertex indices
+            nodes = {}
+            n = 0 #counter
+            for i in range(ntri):
+                for j in range(3):
+                    if tuple(tri[i,j,:]) in nodes:
+                        tris[i,j] = nodes[tuple(tri[i,j,:])]
+                    else:
+                        nodes[tuple(tri[i,j,:])] = n
+                        tris[i,j] = n
+                        n+=1
+            #find the normal vector to the faces using the vertex order
+            a = tri[:ntri, 0, :] - tri[:ntri, 1, :]
+            b = tri[:ntri, 0, :] - tri[:ntri, 2, :]
+
+            crosses = np.cross(a, b)
+            crosses = crosses / (np.sum(crosses ** 2, axis=1) ** (0.5))[:, np.newaxis]
+
+            #get barycentre of faces and findproperty gradient from scalar field
+            tribc = np.mean(tri[:ntri,:,:],axis=1)
 
-        ##convert from memoryview to np array
-        tri = np.array(tri)
-        ntri = np.array(ntri)[0]
-        ##create a triangle indices array and initialise to -1
-        tris = np.zeros((ntri,3)).astype(int)
-        tris[:,:] = -1
-        ##create a dict for storing nodes index where the key is the node as as a tuple. 
-        #A dict is preferable because it is very quick to check if a key exists
-        #assemble arrays for unique vertex and triangles defined by vertex indices
-        nodes = {}
-        n = 0 #counter
-        for i in range(ntri):
-            for j in range(3):
-                if tuple(tri[i,j,:]) in nodes:
-                    tris[i,j] = nodes[tuple(tri[i,j,:])]
-                else:
-                    nodes[tuple(tri[i,j,:])] = n
-                    tris[i,j] = n
-                    n+=1
-        nodes_np = np.zeros((n,3))
-        for v in nodes.keys():
-            nodes_np[nodes[v],:] = np.array(v)
-        #if lv==None:
-        #    lv = lavavu.Viewer(border=True,quality=2)
-        surf = lv.triangles(name)
-        surf.vertices(nodes_np)
-        surf.indices(tris)
-        surf.colours(colour)
-        if interactive:
-            lv.interactive()
-        if draw:
-            lv.display()
+            propertygrad=self.eval_gradient(tribc,prop='strati')
+            propertygrad/=np.linalg.norm(propertygrad,axis=1)[:,None]
+
+            #dot product between gradient and normal indicates if faces are incorrectly ordered
+            dotproducts = (propertygrad * crosses).sum(axis=1)
+
+            #if dot product >0 then adjust triangle indexing
+            indices = (dotproducts>0).nonzero()[0]
+            tris[indices] = tris[indices,::-1]
+            nodes_np = np.zeros((n,3))
+            for v in nodes.keys():
+                nodes_np[nodes[v],:] = np.array(v)
+            surf = lv.triangles(name)
+            surf.vertices(nodes_np)
+            surf.indices(tris)
+            surf.colours(colour)
         return lv
-    def lv_plot_vector_field(self,propertyname,**kwargs):
+    def lv_plot_vector_field(self,propertyname,lv,**kwargs):
         try:
             import lavavu
         except ImportError:
             print("Cannot import Lavavu")
-            return
-
+        if 'colour' not in kwargs:
+            kwargs['colour'] = 'black'
+        vectorslicing = 100
+        if 'vectorslicing' in kwargs:
+            vectorslicing = kwargs['vectorslicing']
+        vector = self.property_gradients[propertyname]
+        #normalise
+        vector/=np.linalg.norm(vector,axis=1)[:,None]
+        vectorfield = lv.vectors(propertyname+"_grad",**kwargs)
+        vectorfield.vertices(self.barycentre[::vector_slicing,:])
+        vectorfield.vectors(vectors)
+        return