commit 954eb18d8d6b42b5846f2f0b534b6a15f8584877
Author: Andrej <andrej@al-x390>
Date:   Wed Feb 10 20:14:23 2021 +0100

    first commit

diff --git a/mesh.py b/mesh.py
new file mode 100644
index 0000000..9f2c170
--- /dev/null
+++ b/mesh.py
@@ -0,0 +1,174 @@
+#!/usr/bin/python3
+
+import numpy as np
+from scipy.spatial import Delaunay
+import matplotlib.pyplot as plt
+import shapely.geometry as geo
+
+#import sys
+#np.set_printoptions(threshold=sys.maxsize)
+
+def get_mesh(fd,fh,h0,bbox,pfix):
+    '''
+    A simple mesh generator based on:
+    Persson, G. Strang, A Simple Mesh Generator in MATLAB.
+    SIAM Review, Volume 46 (2), pp. 329-345, June 2004
+    '''
+    geps = 0.001*h0
+    ttol = 0.1
+    Fscale = 1.2
+    deltat = 0.2
+    deps = np.sqrt(np.finfo(float).eps)*h0
+    dptol = 0.001
+    maxIter = 500
+
+    X, Y = np.meshgrid(np.arange(bbox[0,0],bbox[1,0]+h0,h0),\
+            np.arange(bbox[0,1],bbox[1,1]+h0*0.5*np.sqrt(3),h0*0.5*np.sqrt(3)))
+    X[1:-1:2,:] += (0.5*h0)
+    p = np.column_stack((X.flatten('F'),Y.flatten('F')))
+
+    p = p[np.nonzero(np.apply_along_axis(fd,1,p) < geps)].squeeze()
+    r0 = 1/np.apply_along_axis(fh,1,p)**2
+    R0 = np.max(r0)
+    p = p[np.random.rand(p.shape[0]) < r0/R0]
+    if len(pfix):
+        p = np.vstack((pfix, p))
+    N = p.shape[0]
+
+    pold = np.Inf
+    counter = 0
+    while True:
+        if np.max(np.sqrt(np.sum((p-pold)**2,1))/h0) > ttol:
+            pold = np.copy(p)
+            T = Delaunay(p)
+            tri = T.simplices
+            pmid = (p[tri[:,0],:]+p[tri[:,1],:]+p[tri[:,2],:])/3
+            inner = np.apply_along_axis(fd,1,pmid) < -geps
+            tri = tri[inner,:]
+            bars = np.vstack((tri[:,[0,1]],tri[:,[1,2]],tri[:,[0,2]]))
+            bars.sort(1)
+            bars = remove_duplicate_row(bars)
+
+
+        bvec = p[bars[:,0],:] - p[bars[:,1],:]
+        L = np.sqrt(np.sum(bvec**2,1))
+        hbars = np.apply_along_axis(fh, 1,(p[bars[:,0],:]+p[bars[:,1],:])/2)
+        L0 = hbars*Fscale*np.sqrt(np.sum(L**2)/np.sum(hbars**2))
+        F = [np.max(L0 - [L], 0)]
+        Fvec = (F/L).T * bvec 
+
+        #sparse and back to full does not work the same in python as in matlab
+        #forced to do the loop, figure out better way
+        row = bars[:,[0,0,1,1]].flatten('F')
+        column = (np.ones(np.shape(F),dtype=np.int32).T @ [[0,1,0,1]]).flatten('F')
+        Ffrac = np.hstack((Fvec,-Fvec)).flatten('F')
+        Ftot = np.zeros([N,2])
+        for i in range(len(row)):
+            Ftot[row[i],column[i]] += Ffrac[i]
+        if len(pfix):
+            Ftot[0:pfix.shape[0]-1,:] = 0
+        p += (deltat*Ftot)
+
+        #bring back to boundary
+        d = np.apply_along_axis(fd,1,p)
+        ix = d>0 #index of points outside boundary
+        dx = np.column_stack((p[ix,0]+deps, p[ix,1]))
+        dy = np.column_stack((p[ix,0], p[ix,1]+deps))
+
+        gradx = (np.apply_along_axis(fd,1,dx)-d[ix])/deps
+        grady = (np.apply_along_axis(fd,1,dy)-d[ix])/deps
+        p[ix,:] -= np.column_stack((d[ix]*gradx,d[ix]*grady))
+
+        #print(np.max(np.sqrt(np.sum(deltat*Ftot[d<-geps,:]**2,1))/h0))
+
+        counter += 1
+        #print(counter)
+
+        if Ftot[d<-geps,:].size == 0 or counter > maxIter:
+            print(counter)
+            return p, tri
+        if np.max(np.sqrt(np.sum(deltat*Ftot[d<-geps,:]**2,1))/h0) < dptol:
+            print(counter)
+            return p, tri
+
+def remove_duplicate_row(M):
+    #remove dupplicate rows of the array
+    A = np.ascontiguousarray(M).view(np.dtype((np.void,M.dtype.itemsize*M.shape[1])))
+    _, idx = np.unique(A, return_index=True)
+    return M[idx]
+
+def print_mesh(p, tri, title='', output=''):
+    plt.figure()
+    plt.triplot(p[:,0], p[:,1], tri, color='black')
+    plt.axis('equal')
+    plt.axis(False)
+    if title:
+        plt.title(title)
+    if output:
+        plt.savefig(output)
+        plt.close()
+        return
+    plt.show()
+    return
+
+def ddiff(d1, d2):
+    return max(d1, -d2)
+
+def dunion(d1, d2):
+    return min(d1, d2)
+
+def dintersect(d1, d2):
+    return max(d1, d2)
+
+def dcircle(p, xc, yc, r):
+    return np.sqrt((p[0]-xc)**2+(p[1]-yc)**2)-r
+
+def drectangle(p, x1, x2, y1, y2):
+    return -min(min(min(-y1+p[1],y2-p[1]),-x1+p[0]), x2-p[0])
+
+def protate(p, phi):
+    A = np.array([[np.cos(phi),-np.sin(phi)],[np.sin(phi),np.cos(phi)]])
+    return (A@p.T).T
+
+def _dpoly(p, poly, ring):
+    ds = ring.distance(p)
+    inside = ring.contains(p) or poly.contains(p)
+    return (-1)**inside*ds
+
+def dpoly(p, poly):
+    return _dpoly(geo.asPoint(p), geo.asPolygon(poly),geo.asLinearRing(poly))
+
+if __name__ == "__main__":
+    #print("test 1")
+    #fd = lambda x: ddiff(dcircle(x, 0, 0, 1), dcircle(x, 0.5, -0.2, 0.2))
+    #fh = lambda x: 0.8
+    #bbox = np.array([[-1, -1],[1, 1]])
+    #p, tri = get_mesh(fd, fh, 0.1, bbox, [])
+    #print_mesh(p, tri)
+
+    print("test 2")
+    import time
+    poly = np.array([[0,0],[0,1],[1,1],[1,2],[2,2],[2,1],[3,1],[3,0]])
+    polyG = geo.asPolygon(poly)
+    ringG = geo.asLinearRing(poly)
+    fd = lambda x: dpoly(x,poly)
+    fh = lambda x: 0.15-0.22*fd(x)
+    bbox = np.array([[0,0],[3,3]])
+    t = time.time()
+    p, tri = get_mesh(fd, fh, 0.1, bbox, poly)
+    elap = time.time()-t
+    print("prvi = {}s".format(elap))
+    fd = lambda x: _dpoly(geo.asPoint(x), polyG, ringG)
+    fh = lambda x: 0.15-0.22*fd(x)
+    t = time.time()
+    p, tri = get_mesh(fd, fh, 0.1, bbox, poly)
+    elap = time.time()-t
+    print("drugi = {}s".format(elap))
+    print_mesh(p, tri)
+    
+
+    #print("test 3")
+    #fd = lambda x: drectangle(protate(x,np.pi/4),0,1,0,1)
+    #bbox = np.array([[-2,-2],[2,2]])
+    #p, tri = get_mesh(fd, fh, 0.08, bbox, [])
+    #print_mesh(p, tri)