Poisson2D

Approximate a solution of the Poisson equation using finite element method.

HTML

<script src="http://www.numericjs.com/lib/numeric-1.2.6.min.js"></script>
<script src="http://ssrb.github.io/assets/schur2/js/CanvasMatrix.js"></script>
<script src="http://ssrb.github.io/assets/schur2/js/mesh.js"></script>
<script src="http://ssrb.github.io/assets/schur2/js/Controls.js"></script>
<body>
    <script id="shader-vs" type="x-shader/x-vertex">
        attribute vec2 aPos;
        attribute float aSol;
        varying float sol;

        uniform mat4 mvMatrix;
        uniform mat4 prMatrix;

        uniform float solMin, solMax;

        void main(void) {
            gl_Position = prMatrix * mvMatrix * vec4(aPos, 10. * aSol, 1.);
            sol = 1. - (aSol - solMin) / (solMax - solMin);
        }
    </script>
    <script id="shader-fs" type="x-shader/x-fragment">
        precision highp float;
        varying float sol;

        vec3 heatcolor() {

            float h = 4. * sol,
            s = 1.,
            v = 0.8,
            f = h - floor(h),
            p = v * (1. - s),
            q = v * (1. - s * f),
            t = v * (1. - s * (1. - f));

            if (h <= 1.) {
                return vec3(v, t, p);
            }

            if (h <= 2.) {
                return vec3(q, v, p);
            }

            if (h <= 3.) {
                return vec3(p, v, t);
            }

            return vec3(p, q, v);
        }

        void main(void) {
            gl_FragColor = vec4(heatcolor(), 1.);
        }
    </script>
    <canvas id="canvas" width="400" height="300"></canvas>
</body>

JavaScript

function main() {
    Ab = assembleStiffnessMatrixAndLoadVector();
    var solution = solveLinearSystem(Ab[0], Ab[1]);
    displaySolution(solution);
}

function assembleStiffnessMatrixAndLoadVector() {

    var vertex2Border = [];
    for (i = 0; i < border.length; ++i) {
        vertex2Border[border[i]] = 1;
    }

    var nbVertices = vertices.length / 2,
        nbTriangles = triangles.length / 3;

    var A = numeric.rep([nbVertices, nbVertices], 0);
    var b = numeric.rep([nbVertices], 0);

    for (var ti = 0; ti < nbTriangles; ++ti) {

        var q = [];
        for (i = 0; i < 3; ++i) {
            var si = 3 * ti + i,
                sj = 3 * ti + ((i + 1) % 3);
            q[i] = [vertices[2 * triangles[si]] - vertices[2 * triangles[sj]],
            vertices[2 * triangles[si] + 1] - vertices[2 * triangles[sj] + 1]];
        }

        var area = 0.5 * numeric.det([q[0], q[1]]);

        for (i = 0; i < 3; ++i) {
            var vi = triangles[3 * ti + i];
            if (vertex2Border[vi] != 1) {
                for (j = 0; j < 3; ++j) {
                    var vj = triangles[3 * ti + j];
                    if (vertex2Border[vj] != 1) {
                        var qi = (i + 1) % 3,
                            qj = (j + 1) % 3;
                        A[vi][vj] += numeric.dot(q[qi], q[qj]) / (4 * area);
                    }
                }
                b[vi] += -area / 3;
            }
        }
    }

    return [A, b];
}

function solveLinearSystem(A, b) {
    return numeric.ccsLUPSolve(numeric.ccsLUP(numeric.ccsSparse(A)), b);
}

var c_w, c_h, prMatrix, mvMat, mvMatLoc, rotMat;

function displaySolution(sol) {

    var solMin = sol[0],
        solMax = sol[0];

    for (i = 0; i < sol.length; ++i) {
        if (sol[i] < solMin) {
            solMin = sol[i];
        } else if (sol[i] > solMax) {
            solMax = sol[i];
        }
    }

    function anim() {
        drawScene();
        requestAnimationFrame(anim);
    }

    initGL();
...