AllKnittStuff.morpho

//Knit Stuff
import meshtools
import symmetry
import constants
import optimize
import povray
import plot
import functionals 

fn addSymmetryToKnit(m,h,k,N,M,fix_ends = false){
    /* List the x-y coordinates of the end-points.
    Makes NxM unit cells

    The end points are labeled like so:
          c   e
        _ |_ _| _ 
       |  |   |  |
    d---\-'   '-/---f
         \     /
         a     b

    Thus, for the periodicity, 
    a <-> c
    b <-> e
    d <-> f

    d a b are alway taken to be on the bottom/left of the knit
    c e f are ont the top/right
    */
    // loop over units in X direction adding bottom/top symmetry 
    // I.E. a <-> c, b <-> e
    var dx = (1 + 2*k)
    var dy = h

    var disp_x = Matrix([(1+2*k)*N, 0, 0]) // This is the discrete period of the knit in the x direction
    var disp_y = Matrix([0, h*M, 0]) // This is the discrete period of the knit in the y direction

    var ty = Translate(disp_y)

    var ends_yT = Selection(m)
    var ends_xT = Selection(m)
    for (i in 0...N){
        
        var xa = -0.5 + i * dx
        var ya = -0.5*h

        var xb = 0.5 + i * dx
        var yb = -0.5*h

        var xc = -0.5 + i * dx
        var yc = 0.5*h + dy * (M-1)

        var xe = 0.5 + i * dx
        var ye = 0.5*h + dy * (M-1)

        var a = Selection(m, fn (x,y,z) abs(x-xa)<0.01 && abs(y-ya)<0.01)
        a.addgrade(0)
        var b = Selection(m, fn (x,y,z) abs(x-xb)<0.01 && abs(y-yb)<0.01)
        b.addgrade(0)
        var c = Selection(m, fn (x,y,z) abs(x-xc)<0.01 && abs(y-yc)<0.01)
        c.addgrade(0)
        var e = Selection(m, fn (x,y,z) abs(x-xe)<0.01 && abs(y-ye)<0.01)
        e.addgrade(0)


        // Impose discrete translation symmetry in y direction for points a, c, b and e
        var ends_y = a.union(c.union(b.union(e)))
        ends_yT = ends_yT.union(ends_y)

        m.addsymmetry(ty, ends_y)
    }

    var tx = Translate(disp_x)

    for (j in 0...M){

        var xd = -(k+0.5)
        var yd = 0.5*h-1 + dy * j

        var xf = (k+0.5) + dx * (N-1)
        var yf = 0.5*h-1 + dy * j

        // Select the end-points on the mesh

        var d = Selection(m, fn (x,y,z) abs(x-xd)<0.01 && abs(y-yd)<0.01)
        d.addgrade(0)
        var f = Selection(m, fn (x,y,z) abs(x-xf)<0.01 && abs(y-yf)<0.01)
        f.addgrade(0)
        
        // Add symmetries       
        // Impose discrete translation symmetry in x direction for points d and f
        var ends_x = d.union(f) 
        ends_xT = ends_xT.union(ends_x)
        m.addsymmetry(tx, ends_x)
    }

    return [tx,ty,ends_xT,ends_yT]
}

fn knitStart(h = 1.2, k = 0.5, dx = 0.05, N = 10, M = 10){
    // N Number of unit cells in x dir
    // M Number of unit cells in y dir


    // var h = 1.2 // Decides the vertical spacing between the knits. This number is eye-balled right now
    // var k = 0.5 // ratio of the minor axis to the major axis of the ellipse 
    // var dx = 0.05*Pi
    dx = dx * Pi

    var m1 = LineMesh(fn (t) [k*sin(t)-(k+0.5), 0.5*h+cos(t), 0], 0.5*Pi..Pi:dx, closed=false) // Left thread of the ellipse
    var m2 = LineMesh(fn (t) [k*sin(t)+(k+0.5), 0.5*h+cos(t), 0], -0.5*Pi..-Pi:(-dx), closed=false) // Right thread of the ellipse
    var m3 = LineMesh(fn (t) [sin(t)*(0.5 + cos(t)^2), cos(t)-0.5*h, 0.5*cos(3*t)], -0.5*Pi..0.5*Pi:dx, closed=false) // Part of the dumbbell curve

    var merge = MeshMerge([m1])
    merge.addmesh([m2])
    merge.addmesh([m3])
    
    var merge = MeshMerge([m1,m2,m3])
    var m = merge.merge()

    // now we have the unit cell make an NxM grid of them
    var unitCellMerge = MeshMerge([])
    for (i in 0...N){
        for (j in 0...M){
            var newM = m.clone()
            for (mi in 0...newM.count()) {
                newM.setvertexposition(mi,newM.vertexposition(mi) + Matrix([2*(0.5+k)*(i),(h)*(j),0]))
            }
            unitCellMerge.addmesh([newM])
        }
    }
    var finalM = unitCellMerge.merge()

    finalM.addgrade(1)
    var sym = addSymmetryToKnit(finalM,h,k,N,M)
    return [finalM,sym[0],sym[1],sym[2],sym[3]]
}

var fix_ends = true
var MeshAndSymm = knitStart(N = 1, M = 1)
var m = MeshAndSymm[0]
var tx = MeshAndSymm[1]
var ty = MeshAndSymm[2]
var ends_x = MeshAndSymm[3]
var ends_y = MeshAndSymm[4]
var problem = OptimizationProblem(m)

// Add Hertzian sphere potential 

// define the cutoff for the Hertzian potential
var sigma = 0.05
// var lv = PairwisePotential(fn (r) cutoff(r, (1-r/sigma)^(5/2), sigma), fn (r) cutoff(r, -5/(2*sigma)*(1-r/sigma)^(3/2), sigma))
var lv = PairwisePotential(fn (r) (1-r/sigma)^(5/2), fn (r) -5/(2*sigma)*(1-r/sigma)^(3/2), cutoff=sigma)
problem.addenergy(lv, prefactor=100) // Using a large prefactor to ensure impenetrability

// Add line curvature and length energies
var lc = LineCurvatureSq()
var len = Length()

problem.addenergy(len)
problem.addenergy(lc)


var opt = ShapeOptimizer(problem, m)
    var yvals = []
    for (mi in 0...m.count()) {
        yvals.append(m.vertexposition(mi)[1])
    }
var actualYStretch = max(yvals)-min(yvals)
print("y-dist is ${actualYStretch} at start")

if (fix_ends){

  opt.fix(ends_x)
  opt.fix(ends_y)

}

opt.stepsize = 0.00001 // 0.00001 works
opt.conjugategradient(10000) // Conjugate gradient works ~150 times faster than relax!

var yvals = []
for (mi in 0...m.count()) {
    yvals.append(m.vertexposition(mi)[1])
}
var actualYStretch = max(yvals)-min(yvals)
print("y-dist is ${actualYStretch}")


var g2 = plotmesh(m, grade=1)
g2.title = "After"
Show(g2)
//print m.connectivitymatrix(0,0).count()
m.save("unit_cell.mesh")


// loop to try quasistatcially streching the knitt
/*
var h =1.2
var dy = 0.01

for (i in 0..100){
    // Matrix to multiply each vertex by, this will strectch the top and bottom by a 
    //constant factor dy while moving the other points by a proptional ammout
    var Factor = Matrix([[1,0,0],[0,dy/(h+i*dy)+1,0],[0,0,1])
    ty.v = ty.v *Factor
    for (mi in 0...m.count()) {
    m.setvertexposition(mi,m.vertexposition(mi)*Factor)
    opt.conjugategradient(10000) // Conjugate gradient works ~150 times faster than relax!
    plotmesh(m, grade=1)

}

*/