### More cloth

Verlet integration is used to simulate cloth.
The code is not optimized (and very slow) and self-intersection of the cloth is not prevented.

en.wikipedia.org/wiki/verlet_integration

#physics #cloth

```// More cloth. Created by Reinder Nijhoff 2020
// @reindernijhoff
//
// https://turtletoy.net/turtle/b3acf08303
//

const turtle = new Slowpoke();
const cameraPos = [-7,3,8];
const cameraLookAt = [-0.1,-1.3,0];
const viewProjectionMatrix = setupCamera(cameraPos, cameraLookAt);
const polygons = new Polygons();
const shape = 1; // min=0, max=1, step=1 (Sphere, Box)
const simulationSteps = 175; // min=1, max=250, step=1

// standard verlet integration
class point {
constructor(p, mass) {
this.p = [...p];
this.p_prev = [...p];
// The cloth has more weigth in the center to make the simulation more stable.
this.mass = mass / dot3(this.p, this.p);
this.acc = [0, -9.8, 0];
}
attach(p2, stiffness) {
}
update(dt) {
let vel = sub3(this.p, this.p_prev);
vel = scale3(vel, .9);
this.p_prev = this.p;
this.p = p;
}
solveCollisions() {
if (shape == 0) { // collision with sphere
if (len3(this.p) <= 1) {
this.p = normalize3(this.p_prev);
this.p_prev = [...this.p];
}
} else { // colision with block
const clamp = (a,b,t) => Math.max(Math.min(b,t),a);
if (Math.abs(this.p[0]) < 1 && Math.abs(this.p[2]) < 1 && this.p[1] <= 1 && this.p[1] >= .9) {
const p = this.p_prev;
this.p = [p[0], clamp(.9, 1, p[1]+.01), p[2]];
this.p_prev = [...this.p];
}
}
}
solve() {
this.solveCollisions();
}
}

constructor(p1, p2, stiffness) {
this.p1 = p1;
this.p2 = p2;
this.stiffness = stiffness;
this.distance = dist3(p1.p, p2.p);
}
solve() {
const v = sub3(this.p1.p, this.p2.p);
const d = len3(v);
const diff = (this.distance - d)/d;
const s = (1/this.p1.mass) / (1/this.p1.mass + 1/this.p2.mass);
}
}

class world {
constructor() {
this.points = [];
}
this.points.push(p);
return p;
}
update(dt, relax) {
this.points.forEach(p => p.update(dt));
for (let i=0; i<relax; i++) {
this.points.forEach(p => p.solve());
}
}
}

// cloth specific code
class cloth {
constructor(world, w, h, s, center) {
this.width = w;
this.height = h;
const ps = this.points = [];
const tr = this.triangles = [];
const trgrid = [];

const mass = 1;

for (let x=0; x<w; x++) {
this.points[x] = [];
for (let y=0; y<h; y++) {
const p = new point(add3([(x-(w-1)/2)*s, 0, (y-(h-1)/2)*s], center), mass);
}
}

const stiffness = .95;
const stiffness2 = .1;

for (let x=0; x<w; x++) {
for (let y=0; y<h; y++) {
const p = ps[x][y];
if (y < h-1) p.attach(ps[x][y+1], stiffness);
if (x < w-1) p.attach(ps[x+1][y], stiffness);
if (x < w-1 && y < h-1) p.attach(ps[x+1][y+1], stiffness);

if (y < h-2) p.attach(ps[x][y+2], stiffness2);
if (x < w-2) p.attach(ps[x+2][y], stiffness2);
if (x > 1  && y < h-2) p.attach(ps[x-2][y+2], stiffness2);
}
}

for (let x=0; x<w-1; x++) {
trgrid[x] = [];
for (let y=0; y<h-1; y++) {
const t1 = new triangle(0,x,y,ps[x][y],ps[x+1][y],ps[x+1][y+1]);
const t2 = new triangle(1,x,y,ps[x][y],ps[x][y+1],ps[x+1][y+1]);

trgrid[x][y*2+0] = t1; trgrid[x][y*2+1] = t2;
tr.push(t1); tr.push(t2);
}
}

const trg = (x, y) => (x >= 0 && x < w-1 && y >= 0 && y < h*2-1) ? trgrid[x][y] : false;
for (let x=0; x<w-1; x++) {
for (let y=0; y<h-1; y++) {
trgrid[x][y*2+0].setNeighbours(trg(x, y*2-1), trg(x+1, y*2+1), trg(x, y*2+1));
trgrid[x][y*2+1].setNeighbours(trg(x-1, y*2+0), trg(x, y*2+2), trg(x, y*2+0));
}
}
}
}

// triangles used for visualisation
class triangle {
constructor(o, x, y, p1, p2, p3) {
this.o = o, this.x = x, this.y = y;
this.p1 = p1, this.p2 = p2, this.p3 = p3;
}
setNeighbours(t1, t2, t3) {
this.t1 = t1, this.t2 = t2, this.t3 = t3;
}
draw(t) {
const p1 = transform4([...this.p1.p,1], viewProjectionMatrix);
const p2 = transform4([...this.p2.p,1], viewProjectionMatrix);
const p3 = transform4([...this.p3.p,1], viewProjectionMatrix);

const s = 50;
const p = polygons.create();
[p2[0]/p2[3]*s, -p2[1]/p2[3]*s],
[p3[0]/p3[3]*s, -p3[1]/p3[3]*s]);

const edge = (t1, t2) => t2 &&
Math.sign(dot3(t2.getNormal(), sub3(t2.p1.p, cameraPos))) ==
Math.sign(dot3(t1.getNormal(), sub3(t1.p1.p, cameraPos)));

if (this.o == 0) {
if (this.x % 2 == 1 || edge(this, this.t2)) p.addSegments(p.cp[1], p.cp[2]);
if (this.y % 2 == 0 || edge(this, this.t1)) p.addSegments(p.cp[0], p.cp[1]);
} else {
if (this.x % 2 == 0 || edge(this, this.t1)) p.addSegments(p.cp[0], p.cp[1]);
if (this.y % 2 == 1 || edge(this, this.t2)) p.addSegments(p.cp[1], p.cp[2]);
}

if (this.x % 4 < 2) {
const pl1 = this.o == 1 ? this.p1.p :  this.p3.p;
const pl2 = this.p2.p;
const pl3 = this.o == 1 ? this.p3.p :  this.p1.p;

const subdiv = 4;

const pl1m3 = sub3(pl1, pl3);
const pl1m2 = sub3(pl1, pl2);

for (let i=0; i<=subdiv; i++) {
const sp = add3(pl2, scale3(pl1m2, i/subdiv));
const ep = add3(pl3, scale3(pl1m3, i/subdiv));

const spp = transform4([...sp,1], viewProjectionMatrix);
const epp = transform4([...ep,1], viewProjectionMatrix);

}
}
polygons.draw(t, p);
}
getDist() {
return this.dist ? this.dist :
this.dist = Math.max(Math.max(dist3(this.p1.p, cameraPos),
dist3(this.p2.p, cameraPos)),
dist3(this.p3.p, cameraPos));
}
getNormal() {
return this.normal ? this.normal :
this.normal = cross3(sub3(this.p1.p, this.p3.p), sub3(this.p2.p, this.p3.p));
}
}

const w = new world();
const r = new cloth(w, 77, 53, .1, [0,1,0]);

w.points.sort((a,b) => Math.random()-.5);

// The walk function will be called until it returns false.
function walk(i) {
if (i < simulationSteps) {
w.update(1/30, 8);

turtle.jump(i*190/simulationSteps-95, 95);
turtle.goto((i+1)*190/simulationSteps-95, 95);

if (i == simulationSteps-1) { // sort front to back
r.triangles = r.triangles.sort( (a,b) => a.getDist() - b.getDist());
}
} else {
r.triangles[i - simulationSteps].draw(turtle);
}
return i < r.triangles.length - 1 + simulationSteps;
}

function setupCamera(camPos, camLookat) {
const viewMatrix = lookAt4m(camPos, camLookat, [0,1,0]);
const projectionMatrix = perspective4m(0.25, 1);
return multiply4m(projectionMatrix, viewMatrix);
}

// vec3 functions
function scale3(a,b) { return [a[0]*b,a[1]*b,a[2]*b]; }
function len3(a) { return Math.sqrt(dot3(a,a)); }
function dist3(a,b) { return Math.sqrt((a[0]-b[0])**2+(a[1]-b[1])**2+(a[2]-b[2])**2); }
function normalize3(a) { return scale3(a,1/len3(a)); }
function add3(a,b) { return [a[0]+b[0],a[1]+b[1],a[2]+b[2]]; }
function sub3(a,b) { return [a[0]-b[0],a[1]-b[1],a[2]-b[2]]; }
function dot3(a,b) { return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; }
function cross3(a,b) { return [a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0]]; }
// vec4 functions
function transform4(a,b) {
const d=new Float32Array(4);
for(let c=0;4>c;c++)d[c]=b[c]*a[0]+b[c+4]*a[1]+b[c+8]*a[2]+b[c+12];
return d;
}
// mat4 functions
function lookAt4m(a,b,d) { // pos, lookAt, up
const c=new Float32Array(16);
b=normalize3(sub3(a,b));
d=normalize3(cross3(d,b));
const e=normalize3(cross3(b,d));
c[0]=d[0];c[1]=e[0];c[2]=b[0];c[3]=0;
c[4]=d[1];c[5]=e[1];c[6]=b[1];c[7]=0;
c[8]=d[2];c[9]=e[2];c[10]=b[2];c[11]=0;
c[12]=-(d[0]*a[0]+d[1]*a[1]+d[2]*a[2]);
c[13]=-(e[0]*a[0]+e[1]*a[1]+e[2]*a[2]);
c[14]=-(b[0]*a[0]+b[1]*a[1]+b[2]*a[2]);
c[15]=1;
return c;
}
function multiply4m(a,b) {
const d=new Float32Array(16);
for(let c=0;16>c;c+=4)
for(let e=0;4>e;e++)
d[c+e]=b[c+0]*a[0+e]+b[c+1]*a[4+e]+b[c+2]*a[8+e]+b[c+3]*a[12+e];
return d;
}
function perspective4m(a,b) { // fovy, aspect
const c=(new Float32Array(16)).fill(0,0);
c[5]=1/Math.tan(a/2);
c[0]=c[5]/b;
c[10]=c[11]=-1;
return c;
}

////////////////////////////////////////////////////////////////
// Polygon Clipping utility code - Created by Reinder Nijhoff 2019
// https://turtletoy.net/turtle/a5befa1f8d
////////////////////////////////////////////////////////////////

function Polygons() {
const polygonList = [];
const Polygon = class {
constructor() {
this.cp = [];       // clip path: array of [x,y] pairs
this.dp = [];       // 2d lines [x0,y0],[x1,y1] to draw
this.aabb = [];     // AABB bounding box
}
// add point to clip path and update bounding box
let xmin = 1e5, xmax = -1e5, ymin = 1e5, ymax = -1e5;
(this.cp = [...this.cp, ...points]).forEach( p => {
xmin = Math.min(xmin, p[0]), xmax = Math.max(xmax, p[0]);
ymin = Math.min(ymin, p[1]), ymax = Math.max(ymax, p[1]);
});
this.aabb = [(xmin+xmax)/2, (ymin+ymax)/2, (xmax-xmin)/2, (ymax-ymin)/2];
}
// add segments (each a pair of points)
points.forEach(p => this.dp.push(p));
}
for (let i = 0, l = this.cp.length; i < l; i++) {
this.dp.push(this.cp[i], this.cp[(i + 1) % l]);
}
}
draw(t) {
for (let i = 0, l = this.dp.length; i < l; i+=2) {
t.jump(this.dp[i]), t.goto(this.dp[i + 1]);
}
}
const tp = new Polygon();
tp.cp.push([-1e5,-1e5],[1e5,-1e5],[1e5,1e5],[-1e5,1e5]);
const dx = Math.sin(a) * d,   dy = Math.cos(a) * d;
const cx = Math.sin(a) * 200, cy = Math.cos(a) * 200;
for (let i = 0.5; i < 150 / d; i++) {
tp.dp.push([dx * i + cy,   dy * i - cx], [dx * i - cy,   dy * i + cx]);
tp.dp.push([-dx * i + cy, -dy * i - cx], [-dx * i - cy, -dy * i + cx]);
}
tp.boolean(this, false);
this.dp = [...this.dp, ...tp.dp];
}
inside(p) {
let int = 0; // find number of i ntersection points from p to far away
for (let i = 0, l = this.cp.length; i < l; i++) {
if (this.segment_intersect(p, [0.1, -1000], this.cp[i], this.cp[(i + 1) % l])) {
int++;
}
}
return int & 1; // if even your outside
}
boolean(p, diff = true) {
// bouding box optimization by ge1doot.
if (Math.abs(this.aabb[0] - p.aabb[0]) - (p.aabb[2] + this.aabb[2]) >= 0 &&
Math.abs(this.aabb[1] - p.aabb[1]) - (p.aabb[3] + this.aabb[3]) >= 0) return this.dp.length > 0;

// polygon diff algorithm (narrow phase)
const ndp = [];
for (let i = 0, l = this.dp.length; i < l; i+=2) {
const ls0 = this.dp[i];
const ls1 = this.dp[i + 1];
// find all intersections with clip path
const int = [];
for (let j = 0, cl = p.cp.length; j < cl; j++) {
const pint = this.segment_intersect(ls0, ls1, p.cp[j], p.cp[(j + 1) % cl]);
if (pint !== false) {
int.push(pint);
}
}
if (int.length === 0) {
// 0 intersections, inside or outside?
if (diff === !p.inside(ls0)) {
ndp.push(ls0, ls1);
}
} else {
int.push(ls0, ls1);
// order intersection points on line ls.p1 to ls.p2
const cmpx = ls1[0] - ls0[0];
const cmpy = ls1[1] - ls0[1];
int.sort( (a,b) =>  (a[0] - ls0[0]) * cmpx + (a[1] - ls0[1]) * cmpy -
(b[0] - ls0[0]) * cmpx - (b[1] - ls0[1]) * cmpy);

for (let j = 0; j < int.length - 1; j++) {
if ((int[j][0] - int[j+1][0])**2 + (int[j][1] - int[j+1][1])**2 >= 0.001) {
if (diff === !p.inside([(int[j][0]+int[j+1][0])/2,(int[j][1]+int[j+1][1])/2])) {
ndp.push(int[j], int[j+1]);
}
}
}
}
}
return (this.dp = ndp).length > 0;
}
//port of http://paulbourke.net/geometry/pointlineplane/Helpers.cs
segment_intersect(l1p1, l1p2, l2p1, l2p2) {
const d   = (l2p2[1] - l2p1[1]) * (l1p2[0] - l1p1[0]) - (l2p2[0] - l2p1[0]) * (l1p2[1] - l1p1[1]);
if (d === 0) return false;
const n_a = (l2p2[0] - l2p1[0]) * (l1p1[1] - l2p1[1]) - (l2p2[1] - l2p1[1]) * (l1p1[0] - l2p1[0]);
const n_b = (l1p2[0] - l1p1[0]) * (l1p1[1] - l2p1[1]) - (l1p2[1] - l1p1[1]) * (l1p1[0] - l2p1[0]);
const ua = n_a / d;
const ub = n_b / d;
if (ua >= 0 && ua <= 1 && ub >= 0 && ub <= 1) {
return [l1p1[0] + ua * (l1p2[0] - l1p1[0]), l1p1[1] + ua * (l1p2[1] - l1p1[1])];
}
return false;
}
};
return {
list: () => polygonList,
create: () => new Polygon(),
draw: (turtle, p, addToVisList=true) => {
for (let j = 0; j < polygonList.length && p.boolean(polygonList[j]); j++);
p.draw(turtle);
}
};
}

////////////////////////////////////////////////////////////////
// Slowpoke utility code. Created by Reinder Nijhoff 2019
////////////////////////////////////////////////////////////////

function Slowpoke(x, y) {
const linesDrawn = {};
class Slowpoke extends Turtle {
goto(x, y) {
const p = Array.isArray(x) ? [...x] : [x, y];
if (this.isdown()) {
const o = [this.x(), this.y()];
const h1 = o[0].toFixed(2)+'_'+p[0].toFixed(2)+o[1].toFixed(2)+'_'+p[1].toFixed(2);
const h2 = p[0].toFixed(2)+'_'+o[0].toFixed(2)+p[1].toFixed(2)+'_'+o[1].toFixed(2);
if (linesDrawn[h1] || linesDrawn[h2]) {
super.up();
super.goto(p);
super.down();
return;
}
linesDrawn[h1] = linesDrawn[h2] = true;
}
super.goto(p);
}
}
return new Slowpoke(x,y);
}```