blog‎ > ‎

An Example of Finite Element Method Simulation in Pure Java and HTML5 with Paper.js Library

posted Nov 26, 2012, 7:28 AM by Markus Sprunck   [ updated Apr 9, 2013, 12:51 AM ]
By Markus Sprunck; Revision: 1.2; Status: final; Last Content Change: Apr 9, 2013;

This sample application demonstrates the simplicity of Finite Element Method (FEM) for 2D structural mechanics. All the code is written in pure Java without additional libraries. 

Just the code for visualization is written in JavaScript and it draws directly to Canvas of HTML5 with support of Paper.js library.   

Expected Result

The following screen shot shows the result of a simple beam. A cantilever beam can be solved analytically and serves as a test case for the FEM simulation.



Figure 1: Display of the simple default FEM model

Run the main class of the sample and open the view in your html5 capable browser (Chrome, Firefox, Safari, or IE from version 10) and you should see the screen in Figure 1.

Finite Element Method Basics

If you are familiar with the basics you may skip this chapter. In principle FEM is quite simple. All is about creating a large linear equation system that is an approximation of the real physical behavior. 

The following four steps show a simplified view of what happens in structural mechanics - starting with a simple spring and ending with a solid triangular element for FEM.

Figure 2: Simple spring follows one linear equation

A simple spring has just one equation. The resulting force is linear dependent from the displacement. This is true as long the displacement is small compared with the size of the spring. 

Wikipedia says "In 1676 British physicist Robert Hooke discovered the principle behind springs' action, that the force it exerts is proportional to its extension, now called Hooke's law." [1] So this is not really new.

Figure 3: The 2D spring has two linear equation

The same spring in two dimensions needs two linear equations. The important point here is that the equations can also expressed with matrices. This makes it easier to because you avoid to write a lot equations.

Figure 4: Three springs have six linear equations 
and the stiffness matrix has just a filled diagonal.

Putting now three springs together to a triangle leads to six linear equations, because you have now tree points with two degrees of freedom. The stiffness matrix of this system is very simple, each node sums up the stiffness components from two springs. 

Figure 4: A triangle has six linear equations 
and  
the stiffness matrix is more complex.

A solid triangular has the same equations just the element stiffness matrix looks different. It is dependent on the material and geometry. 

Now there a just two additional things to do. 
  • Collect all the single element stiffness matrices of each single triangle to one large system stiffness matrix. In most of the cases this will be a banded matrix. Also determine the system with some displacement constraints and forces.
  • Rearrange the equations so that the known forces can be expressed as displacements. In FEM you will have always mixed linear equations.
A lot of real world problems can easily be reduced to two dimensions, but it is not difficult for you to extend this code also to three dimensional numerical tasks. Even if you are not familiar with physics - it makes fun to create own models and just playing with it. 

Source Code

The code was developed with Eclipse Juno (see also Download) and needs no additional Java libraries. All really interesting happens in the following class Solver.

// File #1: Solver.java 
package com.sprunck.fem.core;

import com.sprunck.fem.io.ModelUtil;

public class Solver extends Model {

    // Same for all elements
    private static final Matrix MATRIX_W = createElasticityMatrix();

    // Thickness of 2D structure in mm
    private static final double THICKNESS = 10.0f;

    // Poisson's Ratio of material, e.g. steel=0.27–0.30
    private static final double POISSON_RATIO = 0.2f;

    // Young's Modulus of material in N/mm^2
    private static final double YOUNGS_MODULUS = 1.6E+05f;

    public void run(String input) {
        // The bandwidth is the maximal nodeId difference in the model
        final int bandWidth = ModelUtil.parseModelFromString(this, input);

        // Stiffness matrix of all elements
        final double[][] stiffnesMatrix = createSystemStiffnessMatrix(bandWidth);
        final MatrixBanded stiffnessOriginal = new MatrixBanded(stiffnesMatrix);

        // Stiffness matrix with replacement of known forces
        rearangeGlobalStiffnesMatrix(stiffnesMatrix);
        final MatrixBanded stiffnessRearanged = new MatrixBanded(stiffnesMatrix);

        // Solve the linear equations
        setDeltaOut(MatrixBanded.solve(stiffnessRearanged, getForcesIn()));
        setForcesOut(stiffnessOriginal.times(getDeltasOut().transpose()).transpose());
    }

    private double[][] createSystemStiffnessMatrix(int bandWidth) {
        // will be done as band matrix to save memory and improve speed
        final double[][] result = new double[getNumberOfNodes() * 2][bandWidth];
        for (int elementID = 1; elementID <= getNumberOfElements(); elementID++) {
            final Matrix Ke = createElementStiffnessMatrix(elementID);
            // insert each element stiffness matrix into the system matrix
            for (int i = 1; i <= 3; i++) {
                for (int j = 1; j <= 3; j++) {
                    final int col = getNodeId(elementID, i) * 2;
                    final int row = getNodeId(elementID, j) * 2;
                    // put just the right side of symmetric matrix
                    if (col - row >= 0) {
                        result[row-2][col-row]  += Ke.getValue(i * 2 - 2, j * 2 - 2);
                        result[row-2][col-row+1]+= Ke.getValue(i * 2 - 1, j * 2 - 2);
                        result[row-1][col-row]  += Ke.getValue(i * 2 - 1, j * 2 - 1);
                    }
                    if (col - row >= 1) {
                        result[row-1][col-row-1]+= Ke.getValue(i * 2 - 2, j * 2 - 1);
                    }
                }
            }
        }
        return result;
    }

    private Matrix createElementStiffnessMatrix(int elementId) {
        final double area = calculateElementArea(elementId);
        final Matrix matrixB = createDeltaDifferentiationMatrix(elementId, area);
        return matrixB.times(MATRIX_W).times(matrixB.transpose()).mult(area*THICKNESS);
    }

    private double calculateElementArea(int elementId) {
        return 0.5f * ((getX(elementId, 3) - getX(elementId, 2)) * (getY(elementId, 1) 
            - getY(elementId, 2)) + (getX(elementId, 1) - getX(elementId, 2)) 
            * (getY(elementId, 2) - getY(elementId, 3)));
    }

    private Matrix createDeltaDifferentiationMatrix(int elementId, double area) {
        final double factor = 1.0 / (area * 2);
        final Matrix result = new Matrix(6, 3);
        result.setValue(0, 0, factor * (getY(elementId, 2) - getY(elementId, 3)));
        result.setValue(1, 0, 0);
        result.setValue(2, 0, factor * (getY(elementId, 3) - getY(elementId, 1)));
        result.setValue(3, 0, 0);
        result.setValue(4, 0, factor * (getY(elementId, 1) - getY(elementId, 2)));
        result.setValue(5, 0, 0);

        result.setValue(0, 1, 0);
        result.setValue(1, 1, factor * (getX(elementId, 3) - getX(elementId, 2)));
        result.setValue(2, 1, 0);
        result.setValue(3, 1, factor * (getX(elementId, 1) - getX(elementId, 3)));
        result.setValue(4, 1, 0);
        result.setValue(5, 1, factor * (getX(elementId, 2) - getX(elementId, 1)));

        result.setValue(0, 2, factor * (getX(elementId, 3) - getX(elementId, 2)));
        result.setValue(1, 2, factor * (getY(elementId, 2) - getY(elementId, 3)));
        result.setValue(2, 2, factor * (getX(elementId, 1) - getX(elementId, 3)));
        result.setValue(3, 2, factor * (getY(elementId, 3) - getY(elementId, 1)));
        result.setValue(4, 2, factor * (getX(elementId, 2) - getX(elementId, 1)));
        result.setValue(5, 2, factor * (getY(elementId, 1) - getY(elementId, 2)));
        return result;
    }

    private void rearangeGlobalStiffnesMatrix(double[][] K) {
        final int bandbreite = K[0].length;
        for (int nodeID = 1; nodeID <= getNumberOfNodes(); nodeID++) {
            if (isFixedX(nodeID)) {
                final int currentRow = nodeID * 2 - 1;
                for (int index = 1; index <= getNumberOfNodes(); index++) {
                    final int row = index * 2 - 1;
                    // upper part of equations
                    if (currentRow > row) {
                        setForceInX(index, getForceInX(index) 
                            - getDeltaInX(nodeID) * K[row - 1][currentRow - row]);
                        setForceInY(index, getForceInY(index) 
                            - getDeltaInX(nodeID) * K[row][currentRow - row - 1]);
                        K[row - 1][currentRow - row] = 0;
                        K[row][currentRow - row - 1] = 0;
                    }
                    // lower part of equations
                    if (currentRow <= row) {
                        if (row - currentRow + 1 <= bandbreite) {
                            setForceInX(index, getForceInX(index) - getDeltaInX(nodeID)
                                    * K[currentRow - 1][row - currentRow]);
                            setForceInY(index, getForceInY(index) - getDeltaInX(nodeID)
                                    * K[currentRow - 1][row - currentRow + 1]);
                            K[currentRow - 1][row - currentRow] = 0;
                            K[currentRow - 1][row - currentRow + 1] = 0;
                        }
                    }
                }
                K[currentRow - 1][0] = 1;
                setForceInX(nodeID, getDeltaInX(nodeID));
            }
            if (isFixedY(nodeID)) {
                final int currentRow = nodeID * 2;
                for (int index = 1; index <= getNumberOfNodes(); index++) {
                    final int row = index * 2;
                    // upper part of equations
                    if (currentRow >= row) {
                        setForceInX(index, getForceInX(index) 
                            - getDeltaInY(nodeID) * K[row - 2][currentRow - row + 1]);
                        setForceInY(index, getForceInY(index) 
                            - getDeltaInY(nodeID) * K[row - 1][currentRow - row + 1]);
                        K[row - 2][currentRow - row + 1] = 0;
                        K[row - 1][currentRow - row] = 0;
                    }
                    if (currentRow == row) {
                        K[currentRow - 1][1] = 0;
                    }
                    // lower part of equations
                    if (currentRow < row) {
                        if (row - currentRow + 1 <= bandbreite) {
                            setForceInX(index, getForceInX(index) - getDeltaInY(nodeID)
                                    * K[currentRow - 1][row - currentRow]);
                            setForceInY(index, getForceInY(index) - getDeltaInY(nodeID)
                                    * K[currentRow - 1][row - currentRow + 1]);
                            K[currentRow - 1][row - currentRow] = 0;
                            K[currentRow - 1][row - currentRow + 1] = 0;
                        }
                    }
                }
                setForceInY(nodeID, getDeltaInY(nodeID));
                K[currentRow - 1][0] = 1;
            }
        }

        // Now all forces are known and will be set to zero
        for (int nodeId = 0; nodeId < getNumberOfNodes() * 2; nodeId++) {
            if (Double.isNaN(getForcesIn().getValue(nodeId, 0))) {
                getForcesIn().setValue(nodeId, 0, 0.0f);
            }
        }
    }

    private static Matrix createElasticityMatrix() {
        final double factor = YOUNGS_MODULUS / (1 - POISSON_RATIO * POISSON_RATIO);
        final Matrix result = new Matrix(3, 3);
        result.setValue(0, 0, factor);
        result.setValue(1, 0, factor * POISSON_RATIO);
        result.setValue(2, 0, 0);
        result.setValue(0, 1, factor * POISSON_RATIO);
        result.setValue(1, 1, factor);
        result.setValue(2, 1, 0);
        result.setValue(0, 2, 0);
        result.setValue(1, 2, 0);
        result.setValue(2, 2, factor * POISSON_RATIO * (1 - POISSON_RATIO) / 2);
        return result;
    }
}
// File #2: Model.java - just for storing the model input and results 
package com.sprunck.fem.core;

public class Model {

    // Simple data structure
    private class Node {
        double x = 0.0f; // in mm
        double y = 0.0f; // in mm
        int nodeID = 0;
    }

    // Model with triangles and corners [elementId][cornerId]
    private Node[][] model;

    // Vector of input forces in N
    private Matrix forcesIn;

    // Vector of known displacements in mm
    private Matrix deltasIn;

    // Vector of resulting forces in N
    private Matrix forcesOut;

    // Vector of resulting displacements in mm
    private Matrix deltasOut;

    // Number of nodes
    private int numberOfNodes = 0;

    // Number of elements
    private int numberOfElements = 0;

    // ///////////////////////////////////////////////////////////////
    // getters, setters and helper for model

    public void initNodesByElement() {
        model = new Node[1 + numberOfElements][1 + numberOfNodes];
        for (int i = 0; i < 1 + numberOfElements; i++) {
            for (int k = 0; k < 1 + numberOfNodes; k++) {
                model[i][k] = new Node();
            }
        }
    }

    public int getNodeId(int elementId, int cornerId) {
        return model[elementId][cornerId].nodeID;
    }

    public void setX(int elementId, int cornerId, double value) {
        model[elementId][cornerId].x = value;
    }

    public void setY(int elementId, int cornerId, double value) {
        model[elementId][cornerId].y = value;
    }

    public double getX(int elementId, int cornerId) {
        return model[elementId][cornerId].x;
    }

    public double getY(int elementId, int cornerId) {
        return model[elementId][cornerId].y;
    }

    public void setNodeId(int elementId, int cornerId, int value) {
        model[elementId][cornerId].nodeID = value;
    }

    // ///////////////////////////////////////////////////////////////
    // getters, setters and helper for forcesIn

    public void initForcesIn() {
        final double[] newInputForces = new double[numberOfNodes * 2];
        for (int i = 0; i < numberOfNodes * 2; i++) {
            newInputForces[i] = Double.NaN;
        }
        forcesIn = new Matrix(newInputForces).transpose();
    }

    public Matrix getForcesIn() {
        return forcesIn;
    }

    public double getForceInY(int nodeId) {
        return forcesIn.getValue(nodeId * 2 - 1, 0);
    }

    public double getForceInX(int nodeId) {
        return forcesIn.getValue(nodeId * 2 - 2, 0);
    }

    public void setForceInY(int nodeId, double value) {
        forcesIn.setValue(nodeId * 2 - 1, 0, value);
    }

    public void setForceInX(int nodeId, double value) {
        forcesIn.setValue(nodeId * 2 - 2, 0, value);
    }

    // ///////////////////////////////////////////////////////////////
    // getters, setters and helper for deltasIn

    public void initDeltasIn() {
        final double[] newInputdeltas = new double[numberOfNodes * 2];
        for (int i = 0; i < numberOfNodes * 2; i++) {
            newInputdeltas[i] = Double.NaN;
        }
        deltasIn = new Matrix(newInputdeltas).transpose();
    }

    public Matrix getDeltasIn() {
        return deltasIn;
    }

    protected double getDeltaInY(int nodeId) {
        return deltasIn.getValue(nodeId * 2 - 1, 0);
    }

    protected double getDeltaInX(int nodeId) {
        return deltasIn.getValue(nodeId * 2 - 2, 0);
    }

    public void setDeltaInX(int nodeId, double value) {
        deltasIn.setValue(nodeId * 2 - 2, 0, value);
    }

    public void setSeltaInY(int nodeId, double value) {
        deltasIn.setValue(nodeId * 2 - 1, 0, value);
    }

    public boolean isFixedY(int nodeId) {
        return !Double.isNaN(deltasIn.getValue(nodeId * 2 - 1, 0));
    }

    public boolean isFixedX(int nodeId) {
        return !Double.isNaN(deltasIn.getValue(nodeId * 2 - 2, 0));
    }

    // ///////////////////////////////////////////////////////////////
    // getters, setters and helper for forcesOut

    public void setForcesOut(Matrix value) {
        forcesOut = value;
    }

    public Matrix getForcesOut() {
        return forcesOut;
    }

    public double getForceOutY(int nodeId) {
        return forcesOut.getValue(0, nodeId * 2 - 1);
    }

    public double getForceOutX(int nodeId) {
        return forcesOut.getValue(0, nodeId * 2 - 2);
    }

    // ///////////////////////////////////////////////////////////////
    // getters, setters and helper for deltasOut

    public void setDeltaOut(Matrix value) {
        deltasOut = value;
    }

    public Matrix getDeltasOut() {
        return deltasOut;
    }

    public double getDeltaOutY(int nodeId) {
        return deltasOut.getValue(0, nodeId * 2 - 1);
    }

    public double getDeltaOutX(int nodeId) {
        return deltasOut.getValue(0, nodeId * 2 - 2);
    }

    public double getDeltaYMean(int elementId) {
        return (getDeltaOutY(getNodeId(elementId, 3)) 
              + getDeltaOutY(getNodeId(elementId, 2)) 
              + getDeltaOutY(getNodeId(elementId, 1))) / 3.0f;
    }

    // ///////////////////////////////////////////////////////////////
    // getter and helper for numberOfNodes

    public void incementNumberOfNodes() {
        numberOfNodes++;
    }

    public int getNumberOfNodes() {
        return numberOfNodes;
    }

    // ///////////////////////////////////////////////////////////////
    // getter and helper for numberOfElements

    public void incementNumberOfElements() {
        numberOfElements++;
    }

    public int getNumberOfElements() {
        return numberOfElements;
    }

}
// File #3: MatrixBanded.java - stores the linear equations and is able to solve them 
package com.sprunck.fem.core;

final public class MatrixBanded {

    private static final int MAX_NUMBER_OF_ITTERATONS = 1000;

    private final int row;

    private final int col;

    private final double[][] matrix;

    public MatrixBanded(double[][] matrix) {

        if (matrix == null || matrix.length <= 0 || matrix[0].length <= 0) {
            throw new RuntimeException("Illegal parameter.");
        }

        row = matrix.length;
        col = matrix[0].length;
        this.matrix = new double[row][col];
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                this.matrix[i][j] = matrix[i][j];
            }
        }
    }

    public Matrix times(Matrix B) {
        final MatrixBanded A = this;

        if (A.row != B.row || B.col != 1) {
            throw new RuntimeException("Illegal matrix dimensions.");
        }

        final Matrix C = new Matrix(A.row, 1);
        for (int i = 0; i < A.row; i++) {
            for (int k = i - A.col; k - i < A.col; k++) {
                if (k >= 0 && k < B.row) {
                    C.matrix[i][0] += getValue(A, i, k) * B.matrix[k][0];
                }
            }
        }
        return C;
    }

    public static double getValue(final MatrixBanded bandmatrix, int row, int col) {
        if (col == row) {
            return bandmatrix.matrix[row][0];
        } else if (row > col) {
            return getValue(bandmatrix, col, row);
        } else if (col - row < bandmatrix.col) {
            return bandmatrix.matrix[row][col - row];
        }
        return 0.0f;
    }

    public static void setValue(final MatrixBanded bandmatrix, int row, int col, 
            double value) {
        if (col == row) {
            bandmatrix.matrix[row][0] = value;
        } else if (row > col) {
            setValue(bandmatrix, col, row, value);
        } else if (col - row < bandmatrix.col) {
            bandmatrix.matrix[row][col - row] = value;
        }
    }

    public static Matrix solve(final MatrixBanded A, final Matrix b) {
        Matrix x = new Matrix(b.row, b.col);
        Matrix r = b.minus(A.times(x));
        Matrix p = new Matrix(r);
        double rsold = r.transpose().times(r).getData()[0][0];
        int i;
        for (i = 1; i < MAX_NUMBER_OF_ITTERATONS; i++) {
            final Matrix Ap = A.times(p);
            final double alpha = rsold / p.transpose().times(Ap).getData()[0][0];
            x = x.plus(p.mult(alpha));
            r = r.minus(Ap.mult(alpha));
            final double rsnew = r.transpose().times(r).getData()[0][0];
            final double beta = rsnew / rsold;
            if (rsnew < 1e-5) {
                break;
            }
            p = r.plus(p.mult(beta));
            rsold = rsnew;
        }
        return x.transpose();
    }
}
// File #4: Matrix.java - helper class for matrix manipulations 
package com.sprunck.fem.core;

/**
 * This class is based on Matrix.java from Robert Sedgewick and Kevin Wayne
 * (http://introcs.cs.princeton.edu/java/95linear/Matrix.java.html). 
 * Thank you for sharing it.
 * 
 */
final public class Matrix {

    final int row;

    final int col;

    final double[][] matrix;

    public Matrix(int row, int col) {
        this.row = row;
        this.col = col;
        matrix = new double[row][col];
    }

    public Matrix(Matrix A) {
        this(A.matrix);
    }

    public Matrix(double[][] matrix) {
        row = matrix.length;
        col = matrix[0].length;
        this.matrix = new double[row][col];
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                this.matrix[i][j] = matrix[i][j];
            }
        }
    }

    public Matrix(double[] vector) {
        row = 1;
        col = vector.length;
        this.matrix = new double[row][col];
        for (int j = 0; j < col; j++) {
            this.matrix[0][j] = vector[j];
        }
    }

    // return C = A + B
    public Matrix plus(Matrix B) {
        final Matrix A = this;
        if (B.row != A.row || B.col != A.col) {
            throw new RuntimeException("Illegal matrix dimensions.");
        }
        final Matrix C = new Matrix(row, col);
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                C.matrix[i][j] = A.matrix[i][j] + B.matrix[i][j];
            }
        }
        return C;
    }

    // return C = A - B
    public Matrix minus(Matrix B) {
        final Matrix A = this;
        if (B.row != A.row || B.col != A.col) {
            throw new RuntimeException("Illegal matrix dimensions.");
        }
        final Matrix C = new Matrix(row, col);
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                C.matrix[i][j] = A.matrix[i][j] - B.matrix[i][j];
            }
        }
        return C;
    }

    public double getValue(int row, int col) {
        return matrix[row][col];
    }

    public double[][] getData() {
        return matrix;
    }

    // swap rows i and j
    void swap(int i, int j) {
        final double[] temp = matrix[i];
        matrix[i] = matrix[j];
        matrix[j] = temp;
    }

    // create and return the transpose of the invoking matrix
    public Matrix transpose() {
        final Matrix A = new Matrix(col, row);
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                A.matrix[j][i] = this.matrix[i][j];
            }
        }
        return A;
    }

    // does A = B exactly?
    public boolean isEqual(Matrix B) {
        final Matrix A = this;
        if (B.row != A.row || B.col != A.col) {
            throw new RuntimeException("Illegal matrix dimensions.");
        }
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                if (A.matrix[i][j] != B.matrix[i][j]) {
                    return false;
                }
            }
        }
        return true;
    }

    // return C = A * B
    public Matrix times(Matrix B) {
        final Matrix A = this;
        if (A.col != B.row) {
            throw new RuntimeException("Illegal matrix dimensions.");
        }
        final Matrix C = new Matrix(A.row, B.col);
        for (int i = 0; i < C.row; i++) {
            for (int j = 0; j < C.col; j++) {
                for (int k = 0; k < A.col; k++) {
                    C.matrix[i][j] += A.matrix[i][k] * B.matrix[k][j];
                }
            }
        }
        return C;
    }

    // return x = A^-1 b, assuming A is square and has full rank
    public Matrix solve(Matrix rhs) {
        if (row != col || rhs.row != col || rhs.col != 1) {
            throw new RuntimeException("Illegal matrix dimensions.");
        }

        // create copies of the data
        final Matrix A = new Matrix(this);
        final Matrix b = new Matrix(rhs);

        // Gaussian elimination with partial pivoting
        for (int i = 0; i < col; i++) {

            // find pivot row and swap
            int max = i;
            for (int j = i + 1; j < col; j++) {
                if (Math.abs(A.matrix[j][i]) > Math.abs(A.matrix[max][i])) {
                    max = j;
                }
            }
            A.swap(i, max);
            b.swap(i, max);

            // singular
            if (A.matrix[i][i] == 0.0) {
                throw new RuntimeException("Matrix is singular.");
            }

            // pivot within b
            for (int j = i + 1; j < col; j++) {
                b.matrix[j][0] -= b.matrix[i][0] * A.matrix[j][i] / A.matrix[i][i];
            }

            // pivot within A
            for (int j = i + 1; j < col; j++) {
                final double m = A.matrix[j][i] / A.matrix[i][i];
                for (int k = i + 1; k < col; k++) {
                    A.matrix[j][k] -= A.matrix[i][k] * m;
                }
                A.matrix[j][i] = 0.0f;
            }
        }

        // back substitution
        final Matrix x = new Matrix(col, 1);
        for (int j = col - 1; j >= 0; j--) {
            double t = 0.0f;
            for (int k = j + 1; k < col; k++) {
                t += A.matrix[j][k] * x.matrix[k][0];
            }
            x.matrix[j][0] = (b.matrix[j][0] - t) / A.matrix[j][j];
        }
        return x;

    }

    public void setValue(int row, int col, double value) {
        matrix[row][col] = value;
    }

    public Matrix mult(double alpha) {
        final Matrix A = this;
        final Matrix C = new Matrix(row, col);
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                C.matrix[i][j] = A.matrix[i][j] * alpha;
            }
        }
        return C;
    }
}
// File #5: ModelUtil.java - creates the default model and supports IO

package com.sprunck.fem.io;

import java.util.HashMap;
import java.util.LinkedList;
import java.util.List;

import com.sprunck.fem.core.Model;

public class ModelUtil {

    public static int parseModelFromString(Model femCore, String input) {

        class Point {
            public double x = 0.0f;
            public double y = 0.0f;
        }
        final List<Point> tempNodes = new LinkedList<Point>();
        tempNodes.add(new Point());

        int bandwidthExpected = 0;

        final List<Integer[]> tempElements = new LinkedList<Integer[]>();
        tempElements.add(new Integer[5]);

        final String[] lines = input.toString().split("\\n");
        for (final String line : lines) {
            if (!line.trim().isEmpty()) {
                final String[] args = line.split(",");
                if (0 == args[0].trim().compareToIgnoreCase("N")) {
                    femCore.incementNumberOfNodes();
                    final int number = Integer.valueOf(args[1].trim());
                    for (int index = tempNodes.size(); index <= number; index++) {
                        tempNodes.add(new Point());
                    }

                    final Integer first = Integer.valueOf(args[2].trim());
                    final Integer second = Integer.valueOf(args[3].trim());

                    tempNodes.get(number).x = first;
                    tempNodes.get(number).y = second;
                }
                if (0 == args[0].trim().compareToIgnoreCase("E")) {
                    femCore.incementNumberOfElements();
                    final int number = Integer.valueOf(args[1].trim());
                    for (int index = tempElements.size(); index <= number; index++) {
                        tempElements.add(new Integer[5]);
                    }

                    final Integer first = Integer.valueOf(args[2].trim());
                    final Integer second = Integer.valueOf(args[3].trim());
                    final Integer third = Integer.valueOf(args[4].trim());

                    tempElements.get(femCore.getNumberOfElements())[1] = number;
                    tempElements.get(femCore.getNumberOfElements())[2] = first;
                    tempElements.get(femCore.getNumberOfElements())[3] = second;
                    tempElements.get(femCore.getNumberOfElements())[4] = third;

                    final int max = Math.max(Math.max(first, second), third);
                    final int min = Math.min(Math.min(first, second), third);
                    final int bandwidthOfElement = (1 + max - min) * 2;
                    bandwidthExpected = Math.max(bandwidthExpected, bandwidthOfElement);
                }
                if (0 == args[0].trim().compareToIgnoreCase("D")) {

                    if (femCore.getDeltasIn() == null) {
                        femCore.initDeltasIn();
                    }

                    final int number = Integer.valueOf(args[1].trim());
                    if (0 == args[2].trim().compareToIgnoreCase("x")) {
                        femCore.setDeltaInX(number, Double.valueOf(args[3].trim()));
                    }
                    if (0 == args[2].trim().compareToIgnoreCase("y")) {
                        femCore.setSeltaInY(number, Double.valueOf(args[3].trim()));
                    }
                }
                if (0 == args[0].trim().compareToIgnoreCase("F")) {

                    if (femCore.getForcesIn() == null) {
                        femCore.initForcesIn();
                    }

                    final int number = Integer.valueOf(args[1].trim());
                    if (0 == args[2].trim().compareToIgnoreCase("x")) {
                        femCore.setForceInX(number, Double.valueOf(args[3].trim()));
                    }
                    if (0 == args[2].trim().compareToIgnoreCase("y")) {
                        femCore.setForceInY(number, Double.valueOf(args[3].trim()));
                    }
                }
            }
        }

        femCore.initNodesByElement();

        for (int i = 1; i <= femCore.getNumberOfElements(); i++) {
            femCore.setX(tempElements.get(i)[1], 1, 
                tempNodes.get(tempElements.get(i)[2]).x);
            femCore.setX(tempElements.get(i)[1], 2, 
                tempNodes.get(tempElements.get(i)[3]).x);
            femCore.setX(tempElements.get(i)[1], 3, 
                tempNodes.get(tempElements.get(i)[4]).x);

            femCore.setY(tempElements.get(i)[1], 1, 
                tempNodes.get(tempElements.get(i)[2]).y);
            femCore.setY(tempElements.get(i)[1], 2, 
                tempNodes.get(tempElements.get(i)[3]).y);
            femCore.setY(tempElements.get(i)[1], 3, 
                tempNodes.get(tempElements.get(i)[4]).y);

            femCore.setNodeId(tempElements.get(i)[1], 1, tempElements.get(i)[2]);
            femCore.setNodeId(tempElements.get(i)[1], 2, tempElements.get(i)[3]);
            femCore.setNodeId(tempElements.get(i)[1], 3, tempElements.get(i)[4]);
        }
        return bandwidthExpected;
    }

    public static String getModelAsJSON(Model model) {
        final HashMap<Integer, Boolean> nodeIds = new HashMap<Integer, Boolean>();
        final StringBuilder pre = new StringBuilder("var elements = [");
        final int numberOfElements = model.getNumberOfElements();
        for (int elementId = 1; elementId <= numberOfElements; elementId++) {
            pre.append("[");
            for (int cornerId = 1; cornerId < 4; cornerId++) {
                final int nodeId = model.getNodeId(elementId, cornerId);
                pre.append("\n{\"id\": " + nodeId //
                        + ", \"x_force\" : " + model.getForceOutX(nodeId) //
                        + ", \"y_force\" : " + model.getForceOutY(nodeId) //
                        + ", \"x_delta\" : " + model.getDeltaOutX(nodeId) //
                        + ", \"y_delta\" : " + model.getDeltaOutY(nodeId) //
                        + ", \"x_fixed\" : " + model.isFixedX(nodeId) //
                        + ", \"y_fixed\" : " + model.isFixedY(nodeId) //
                        + ", \"first\" : " + !nodeIds.containsKey(nodeId) //
                        + ", \"x\" : " + model.getX(elementId, cornerId) //
                        + ", \"y\" : " + model.getY(elementId, cornerId) //
                        + ", \"y_delta_mean\" : " + model.getDeltaYMean(elementId) //
                        + "  }");
                if (cornerId <= 3) {
                    pre.append(',');
                }
                nodeIds.put(nodeId, true);
            }
            pre.append("]");
            if (elementId < numberOfElements + 1) {
                pre.append(',');
            }
        }
        pre.append("];");
        return pre.toString();
    }

    public static String createDefaultModel(Model femCore) {
        final StringBuffer nodeText = new StringBuffer();
        final int maxCols = 20;
        final int maxRows = 5;
        final int scaleFactorX = 10;
        final int scaleFactorY = 10;
        for (int col = 1; col <= maxCols; col++) {
            for (int row = 1; row <= maxRows; row++) {
                final int nodeId = row + maxRows * (col - 1);
                nodeText.append("N, ").append(nodeId).append(", ")
                    .append(col * scaleFactorX).append(", ")
                    .append(row * scaleFactorY).append(",\n");
            }
        }

        for (int col = 1; col < maxCols; col++) {
            for (int row = 1; row < maxRows; row++) {
                final int firstElementId = row * 2 - 1 + (maxRows - 1) * 2 * (col - 1);
                final int secondElementId = row * 2 + (maxRows - 1) * 2 * (col - 1);
                final int node1Id = row + maxRows * (col - 1);
                final int node2Id = row + maxRows * col;
                final int node3Id = row + 1 + maxRows * (col - 1);
                final int node4Id = row + 1 + maxRows * (col + 1 - 1);
                nodeText.append("E, ").append(firstElementId).append(", ")
                    .append(node1Id).append(", ").append(node2Id)
                    .append(", ").append(node3Id).append(",\n");
                nodeText.append("E, ").append(secondElementId).append(", ")
                    .append(node2Id).append(", ")
                    .append(node4Id).append(", ").append(node3Id).append(",\n");
            }
        }

        nodeText.append("D, ").append(1).append(", y, ").append(0).append(",\n");
        for (int row = 1; row <= maxRows; row++) {
            nodeText.append("D, ").append(row).append(", x, ").append(0).append(",\n");
        }

        nodeText.append("F, ").append(maxRows * maxCols).append(", y, ")
                .append(30000.0).append(",\n");

        return nodeText.toString();
    }
}

// File #6: ResultTemplate.html - for visualization

<!DOCTYPE html>
<html>
<head>
    <script type="text/javascript" src="rendermodel.js"></script>
    <script type="text/javascript">
        // Next line will be replaced by the JSON model, don't change it ...
        XX_MODEL_PLACE_HOLDER;
    
        // Only executed our code once the DOM is ready.
        window.onload = function() {
            var canvas = document.getElementById('myCanvas');
            paper.setup(canvas);
            main(elements);
            paper.view.draw();
        }
    </script>
</head>
<body>
    <canvas id="myCanvas" resize></canvas>
</body>
</html>

// File #7: renderModel.js - renders the model with PaperScipt library direct to the HTML5 canvas

includeJavaScript('paper.js');

// Dimensions for display
var height = 400;
var width = 750;

// Scaling factors for stretch or compress during painting
var factorX = 2.5;
var factorY = 2.5;
var factorForce = 0.0005;
var factordelta = 10;

// Offset for model
var offset_x = 40;
var offset_y = 40;

// Offset and size for scale boxs
var offset_x_scale = width - 150; 
var scale_size_x = 25;
var scale_size_y = 250;

// Fixed scale
var minColor = 10;
var maxColor = 0;

function main(elements) {

    // draw all elements
    for ( var ele = 0; ele < elements.length; ele++) {
        var path = new paper.Path();
        path.closed = true;
        path.strokeColor = '#FFFFFF';
        path.strokeWidth = 0.2;
        path.selected = false;
        path.opacity = 0.9;
        for ( var nodeId = 0; nodeId < elements[ele].length; nodeId++) {
            var element = elements[ele][nodeId];
            path.fillColor = getColor(element.y_delta_mean);
            var point = new paper.Point(element.x * factorX + offset_x
                    + element.x_delta * factordelta, element.y * factorY
                    + element.y_delta * factordelta + offset_y);
            path.add(point);
        }
        ele.first = true;
    }

    for ( var ele = 0; ele < elements.length; ele++) {
        var path = new paper.Path();
        for ( var nodeId = 0; nodeId < elements[ele].length; nodeId++) {
            var element = elements[ele][nodeId];
            var point = new paper.Point(element.x * factorX + offset_x
                    + element.x_delta * factordelta, element.y * factorY
                    + element.y_delta * factordelta + offset_y);

            // ensure that all is just drawn once
            if (element.first) {

                var text = new paper.PointText(point);
                text.content = ' ' + element.id;
                text.fontSize = 6;
                text.justification = 'left';

                if (10000.0 < Math.abs(element.x_force)) {
                    drawVector(point, point.add(new paper.Point(element.x_force
                            * factorForce, 0)), true, (element.x_force > 0.0));
                }

                if (10000.0 < Math.abs(element.y_force)) {
                    drawVector(point, point.add(new paper.Point(0,
                            element.y_force * factorForce)), false,
                            (element.y_force > 0.0));
                }

                drawFixed(point, element.y_fixed, element.x_fixed);
            }
        }
    }

    drawBoxes();
    drawScale();
}

function drawVector(vectorStart, vectorEnd, horizontal, positive) {
    var arrow = new paper.Path();
    arrow.strokeWidth = 1.5;
    arrow.strokeColor = '#FF0000';
    arrow.add(vectorStart);
    arrow.add(vectorEnd);
    arrow.opacity = 1;

    var length = 5;
    var head = new paper.Path();
    head.strokeWidth = 1.0;
    head.strokeColor = '#FF0000';
    head.add(vectorEnd);
    head.opacity = 1;
    if (horizontal && !positive) {
        head.add(vectorEnd.add(new paper.Point(length, -length)));
        head.add(vectorEnd);
        head.add(vectorEnd.add(new paper.Point(length, length)));
    } else if (horizontal && positive) {
        head.add(vectorEnd.add(new paper.Point(-length, -length)));
        head.add(vectorEnd);
        head.add(vectorEnd.add(new paper.Point(-length, length)));
    } else if (!horizontal && positive) {
        head.add(vectorEnd.add(new paper.Point(length, -length)));
        head.add(vectorEnd);
        head.add(vectorEnd.add(new paper.Point(-length, -length)));
    } else if (!horizontal && !positive) {
        head.add(vectorEnd.add(new paper.Point(-length, length)));
        head.add(vectorEnd);
        head.add(vectorEnd.add(new paper.Point(length, length)));
    }
}

function drawFixed(vectorEnd, horizontal, vertical) {
    var length = 7;
    if (horizontal) {
        var head = new paper.Path();
        head.strokeWidth = 0.5;
        head.strokeColor = '#000000';
        head.add(vectorEnd);
        head.add(vectorEnd.add(new paper.Point(length * 0.75, length)));
        head.add(vectorEnd.add(new paper.Point(-length * 0.75, length)));
        head.add(vectorEnd);
    }

    if (vertical) {
        var head = new paper.Path();
        head.strokeWidth = 0.5;
        head.strokeColor = '#000000';
        head.add(vectorEnd);
        head.add(vectorEnd.add(new paper.Point(-length, -length * 0.75)));
        head.add(vectorEnd.add(new paper.Point(-length, length * 0.75)));
        head.add(vectorEnd);
    }
}

function drawBoxes() {
    var path = new paper.Path();
    path.closed = true;
    path.strokeWidth = 0.25;
    path.strokeColor = '#000000';
    path.selected = false;
    path.add(new paper.Point(0, 0));
    path.add(new paper.Point(0, height));
    path.add(new paper.Point(offset_x_scale - 5, height));
    path.add(new paper.Point(offset_x_scale - 5, 0));

    var path2 = new paper.Path();
    path2.closed = true;
    path2.strokeWidth = 0.25;
    path2.strokeColor = '#000000';
    path2.selected = false;
    path2.add(new paper.Point(width, 0));
    path2.add(new paper.Point(width, height));
    path2.add(new paper.Point(offset_x_scale, height));
    path2.add(new paper.Point(offset_x_scale, 0));
}

function drawScale() {
    var path = new paper.Path();
    var text = new paper.PointText(new paper.Point(10+offset_x_scale, 20));
    text.content = 'Average delta in y: ';
    text.justification = 'left';
    var delta = 0.10;
    for ( var index = 0.0; index <= 1.0; index += delta) {
        var path = new paper.Path();
        path.add(new paper.Point(10+offset_x_scale,                 
                40 + index * scale_size_y));
        path.add(new paper.Point(10+offset_x_scale,                 
                40 + (index + delta) * scale_size_y));
        path.add(new paper.Point(10+offset_x_scale + scale_size_x, 
                40 + (index + delta) * scale_size_y));
        path.add(new paper.Point(10+offset_x_scale + scale_size_x, 
                40 + index* scale_size_y));
        var text = new paper.PointText(
                new paper.Point(10+offset_x_scale+ scale_size_x * 1.2, 
                        40 + (index + delta / 2) * scale_size_y));

        var value = (minColor + (maxColor - minColor) * index);
        text.content = ' ' + (Math.round(100 * value) / 100.0) + ' mm';
        text.fontSize = 8;
        text.justification = 'left';
        path.closed = true;
        path.strokeWidth = 0.5;
        path.fillColor = getColor(value);
        path.opacity = 0.8;
        path.selected = false;
    }
}

// Color code helper /////////////////////////////////////////////////////////

var frequency = 0.5;

function getColor(value) {
    red = Math.sin(2 - frequency * value) * 127 + 128;
    green = Math.sin(1 - frequency * value) * 127 + 128;
    blue = Math.sin(4 - frequency * value) * 127 + 128;
    return '#' + integerToHex(red) + integerToHex(green) + integerToHex(blue);
}

function integerToHex(n) {
    n = Math.max(0, Math.min(parseInt(n, 10), 255));
    charFirst = "0123456789ABCDEF".charAt((n - n % 16) / 16);
    charSecond = "0123456789ABCDEF".charAt(n % 16);
    return charFirst + charSecond;
}

// Include helper ////////////////////////////////////////////////////////////

function includeJavaScript(filename) {
    // helper for lazy load of additional JavaScript libraries
    var head = document.getElementsByTagName('head')[0];
    var script = document.createElement('script');
    script.src = filename;
    script.type = 'text/javascript';
    head.appendChild(script)
}

// File #8: WebServer.java - just a helper for test (with main method of the application) 

package com.sprunck.fem.io;

import java.io.BufferedInputStream;
import java.io.BufferedOutputStream;
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.io.OutputStream;
import java.net.InetAddress;
import java.net.ServerSocket;
import java.net.Socket;
import java.net.UnknownHostException;

import com.sprunck.fem.core.Model;
import com.sprunck.fem.core.Solver;

public class WebServer extends Thread {

    private static final String NL = System.getProperty("line.separator");

    private final int port = 1234;

    private final Model sf;

    public WebServer(Model sf) {
        this.sf = sf;
        String hostname = "localhost";

        try {
            try {
                final InetAddress addr = InetAddress.getLocalHost();
                hostname = addr.getCanonicalHostName();
            } catch (final UnknownHostException e) {
                System.out.println("Error - " + e.getMessage() + NL);
            }

            final InetAddress addr = InetAddress.getLocalHost();
            hostname = addr.getCanonicalHostName();

            addr.getCanonicalHostName();
        } catch (final UnknownHostException e) {
            System.out.println(e.getMessage());
        }
        // System.out.println("listen at url: http://localhost:1234/index.html");
        System.out.println("webserver started        [http://" + hostname + ':' 
                + port + "/index.html]");

        final Thread serverThread = this;
        serverThread.start();
    }

    @Override
    public void run() {
        Socket connection = null;
        while (true) {
            try {
                final ServerSocket server = new ServerSocket(port);
                connection = server.accept();
                final OutputStream out = 
                    new BufferedOutputStream(connection.getOutputStream());
                final InputStream in = 
                    new BufferedInputStream(connection.getInputStream());
                final String request = readFirstLineOfRequest(in).toString();
                // System.out.println("get request " + request.toString());

                if (request.toLowerCase().startsWith("get /index.html")) {
                    // Create content of response
                    final String contentText = 
                        getPage("com/sprunck/fem/io/ResultTemplate.html").toString();
                    final byte[] content = contentText.getBytes();
                    // For HTTP/1.0 or later send a MIME header
                    if (request.indexOf("HTTP/") != -1) {
                        final String headerString = "HTTP/1.0 200 OK" + NL 
                                + "Server: FEM 1.0" + NL
                                + "Content-length: " + content.length + NL 
                                + "Content-type: text/html" + NL + NL;
                        final byte[] header = headerString.getBytes("ASCII");
                        out.write(header);
                    }
                    out.write(content);
                    out.flush();
                } else if (request.toLowerCase().startsWith("get /paper.js")) {
                    // Create content of response
                    final String contentText = 
                        getPage("com/sprunck/fem/io/paper.js").toString();
                    final byte[] content = contentText.getBytes();
                    // For HTTP/1.0 or later send a MIME header
                    if (request.indexOf("HTTP/") != -1) {
                        final String headerString = "HTTP/1.0 200 OK" + NL 
                                + "Server: FEM 1.0" + NL
                                + "Content-length: " + content.length + NL 
                                + "Content-type: text/javascript" + NL + NL;
                        final byte[] header = headerString.getBytes("ASCII");
                        out.write(header);
                    }
                    out.write(content);
                    out.flush();
                } else if (request.toLowerCase().startsWith("get /rendermodel.js")) {
                    // Create content of response
                    final String contentText = 
                        getPage("com/sprunck/fem/io/renderModel.js").toString();
                    final byte[] content = contentText.getBytes();
                    // For HTTP/1.0 or later send a MIME header
                    if (request.indexOf("HTTP/") != -1) {
                        final String headerString = "HTTP/1.0 200 OK" + NL 
                            + "Server: FEM 1.0" + NL
                            + "Content-length: " + content.length + NL 
                            + "Content-type: text/javascript" + NL + NL;
                        final byte[] header = headerString.getBytes("ASCII");
                        out.write(header);
                    }
                    out.write(content);
                    out.flush();
                } else if (request.startsWith("GET /terminate")) {
                    server.close();
                    throw new RuntimeException("terminate");
                }
                // Close the socket
                connection.close();
                in.close();
                out.close();
                server.close();
            } catch (final IOException e) {
                System.out.println(e.getMessage());
            }
        }
    }

    private StringBuffer readFirstLineOfRequest(final InputStream in) 
            throws IOException {
        final StringBuffer request;
        request = new StringBuffer(100);
        while (true) {
            final int character = in.read();
            if (character == '\n' || character == '\r' || character == -1) {
                break;
            }
            request.append((char) character);
        }
        return request;
    }

    public StringBuffer getPage(String filepath) {
        final StringBuffer fw = new StringBuffer(1000);
        try {
            final InputStream inputstream = 
                    this.getClass().getClassLoader().getResourceAsStream(filepath);
            final InputStreamReader is = new InputStreamReader(inputstream);
            final BufferedReader br = new BufferedReader(is);
            for (String s = br.readLine(); s != null; s = br.readLine()) {
                if (s.contains("XX_MODEL_PLACE_HOLDER")) {
                    fw.append(ModelUtil.getModelAsJSON(sf));
                } else {
                    fw.append(s).append('\n');
                }
            }
            inputstream.close();
        } catch (final Exception xc) {
            xc.printStackTrace();
        }
        return fw;
    }

    public static void main(String[] args) {
        final Solver fem = new Solver();
        fem.run(ModelUtil.createDefaultModel(fem));
        new WebServer(fem);
    }
}

Please, do not hesitate to contact me if you have any ideas for improvement and/or you find a bug in the sample code.  

Refererences

[1]    Wikipedia - Spring (device);
        http://en.wikipedia.org/wiki/Spring_(device)

[2]    Wikipedia - Matrix (mathematics)
        http://en.wikipedia.org/wiki/Matrix_(mathematics)

Change History

 Revision Date  Author  Description
 1.0  Nov 26, 2012  Markus Sprunck   first version
 1.1 Jan 21, 2013 Markus Sprunck smaller improvements
 1.2 Apr 9, 2013 Markus Sprunck source code now on GitHub

Top +1 Rated