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

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

File Solver.java:

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185``` ```package com.sprunck.fem.core; import com.sw_engineering_candies.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 Model.java - just for storing the model input and results

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205``` ```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 MatrixBanded.java - stores the linear equations and is able to solve them

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90``` ```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 Matrix.java - helper class for matrix manipulations

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206``` ```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 ModelUtil.java - creates the default model and supports IO

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194``` ```package com.sprunck.fem.io; import java.util.HashMap; import java.util.LinkedList; import java.util.List; import com.sw_engineering_candies.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 tempNodes = new LinkedList(); tempNodes.add(new Point()); int bandwidthExpected = 0; final List tempElements = new LinkedList(); 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 nodeIds = new HashMap(); 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 ResultTemplate.html - for visualization

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21``` ``` ```

File renderModel.js - renders the model with PaperScipt library direct to the HTML5 canvas

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220``` ```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 WebServer.java - just a helper for test (with main method of the application)

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170``` ```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.sw_engineering_candies.fem.core.Model; import com.sw_engineering_candies.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); } } ```

### Refererences

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

### Find Code on GitHub

#### 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