I developed an FEM toolkit in Java: FuturEye
FuturEye is a Java based Finite Element Method (FEM) Toolkit, providing concise, natural and easy understanding programming interfaces for users who wish to develop researching and/or engineering FEM algorithms for Forward and/or Inverse Problems.
The essential components of FEM are abstracted out, such as nodes, elements, meshes, degrees of freedom and shape function etc. The data and operations of these classes are encapsulated together properly. The classes that different from other existing objectoriented FEM softwares or libraries are function classes. The behavior of the function classes in Futureye is very similar to that in mathematical context. For example algebra of functions, function derivatives and composition of functions. Especially in FEM environment, shape functions, Jacobin of coordinate transforms and numerical integration are all based on the function classes. This feature leads to a more close integration between theory and computer implementation. FuturEye is designed to solve 1D,2D and 3D partial differential equations(PDE) of scalar and/or vector valued unknown functions. The start point of development of this toolkit is solving inverse problems of PDE. In order to solve inverse problems, usually some forward problems must be solved first and many exterior data manipulations should be performed during the solving processes. There are many classes defined for those data operations. However, the data processes are complicated in actual applications, we can not write down all the tools for manipulating data. The design of the basic classes allows operations to all aspect of data structure directly or indirectly in an easily understanding way. This is important for users who need write their own operations or algorithms on the bases of data structure in FuturEye. Some other existing FEM softwares or libraries may over encapsulate the data and operations for such inverse applications. The toolkit can be used for various purposes:

update to version 2.1
How to use WeakFormBuilder:
http://code.google.com/p/futureye/wi...self_weakforms package edu.uta.futureye.tutorial; import java.util.HashMap; import edu.uta.futureye.algebra.SolverJBLAS; import edu.uta.futureye.algebra.intf.Matrix; import edu.uta.futureye.algebra.intf.Vector; import edu.uta.futureye.core.Element; import edu.uta.futureye.core.Mesh; import edu.uta.futureye.core.NodeType; import edu.uta.futureye.function.AbstractFunction; import edu.uta.futureye.function.Variable; import edu.uta.futureye.function.basic.FC; import edu.uta.futureye.function.basic.FX; import edu.uta.futureye.function.intf.Function; import edu.uta.futureye.function.intf.ScalarShapeFunction; import edu.uta.futureye.io.MeshReader; import edu.uta.futureye.io.MeshWriter; import edu.uta.futureye.lib.assembler.AssemblerScalar; import edu.uta.futureye.lib.element.FEBilinearRectangle; import edu.uta.futureye.lib.element.FEBilinearRectangleRegular; import edu.uta.futureye.lib.element.FELinearTriangle; import edu.uta.futureye.lib.weakform.WeakFormBuilder; import edu.uta.futureye.util.container.ElementList; public class UseWeakFormBuilder { /** * <blockquote><pre> * Solve * k*\Delta{u} = f in \Omega * u(x,y) = 0, on boundary x=3.0 of \Omega * u_n + u = 0.01, on other boundary of \Omega * where * \Omega = [3,3]*[3,3] * k = 2 * f = 4*(x^2+y^2)+72 * u_n = \frac{\partial{u}}{\partial{n}} * n: outer unit normal of \Omega * </blockquote></pre> */ public static void rectangleTest() { //1.Read in a rectangle mesh from an input file with // format ASCII UCD generated by Gridgen MeshReader reader = new MeshReader("rectangle.grd"); Mesh mesh = reader.read2DMesh(); mesh.computeNodeBelongsToElements(); //2.Mark border types HashMap<NodeType, Function> mapNTF = new HashMap<NodeType, Function>(); //Robin type on boundary x=3.0 of \Omega mapNTF.put(NodeType.Robin, new AbstractFunction("x","y"){ @Override public double value(Variable v) { if(3.0v.get("x") < 0.01) return 1.0; //this is Robin condition else return 1.0; } }); //Dirichlet type on other boundary of \Omega mapNTF.put(NodeType.Dirichlet, null); mesh.markBorderNode(mapNTF); //3.Use element library to assign degrees of // freedom (DOF) to element ElementList eList = mesh.getElementList(); //FEBilinearRectangle bilinearRectangle = new FEBilinearRectangle(); //If the boundary of element parallel with coordinate use this one instead. //It will be faster than the old one. FEBilinearRectangleRegular bilinearRectangle = new FEBilinearRectangleRegular(); for(int i=1;i<=eList.size();i++) bilinearRectangle.assignTo(eList.at(i)); //4.Weak form. We use WeakFormBuilder to define weak form WeakFormBuilder wfb = new WeakFormBuilder() { /** * Override this function to define weak form */ @Override public Function makeExpression(Element e, Type type) { ScalarShapeFunction u = getScalarTrial(); ScalarShapeFunction v = getScalarTest(); //Call param() to get parameters, do NOT define functions here //except for constant functions (or class FC). Because functions //will be transformed to local coordinate system by param() Function fk = param(e,"k"); Function ff = param(e,"f"); switch(type) { case LHS_Domain: // k*(u_x*v_x + u_y*v_y) in \Omega return fk.M( u._d("x").M(v._d("x")) .A (u._d("y").M(v._d("y"))) ); case LHS_Border: // k*u*v on Robin boundary return fk.M(u.M(v)); case RHS_Domain: return ff.M(v); case RHS_Border: return v.M(0.01); default: return null; } } }; Function fx = FX.fx; Function fy = FX.fy; wfb.addParamters(FC.c(2.0), "k"); //Right hand side(RHS): f = 4*(x^2+y^2)+72 wfb.addParamters(fx.M(fx).A(fy.M(fy)).M(4.0).A(72.0),"f"); //5.Assembly process AssemblerScalar assembler = new AssemblerScalar(mesh, wfb.getScalarWeakForm()); System.out.println("Begin Assemble..."); assembler.assemble(); Matrix stiff = assembler.getStiffnessMatrix(); Vector load = assembler.getLoadVector(); //Boundary condition assembler.imposeDirichletCondition(FC.c0); System.out.println("Assemble done!"); //6.Solve linear system SolverJBLAS solver = new SolverJBLAS(); Vector u = solver.solveDGESV(stiff, load); System.out.println("u="); for(int i=1;i<=u.getDim();i++) System.out.println(String.format("%.3f", u.get(i))); //7.Output results to an Techplot format file MeshWriter writer = new MeshWriter(mesh); writer.writeTechplot("./tutorial/UseWeakFormBuilder2.dat", u); } public static void main(String[] args) { rectangleTest(); } } 
Update

Quote:
just a quick remark: In your plot, it should be Lagrange multiplier, not language multiplier :) 
thanks
many thanks!

Does anyone interested in developing this project?
I have a plan to extend and optimize the toolkit, FuturEye, if you are interested in developing FEM code, please contact me:)

So it is OpenSource (MIT) right? At ETH Zürich we developed once a Software for simulating soil salinisation (SimSalin) in Java including groundwater flow, evapotranspiration and transport, and we came to the point to see that Java is too slow. Since Java is so close related to C++, would it be worth the effort to rewrite part of your code in C++ and call the methods from Java? The discrete element software FARO for Rockfall protection nets for which I implement new Elements has a Java3D Interface and calls C++ methods and is among the fastest codes of its kind.

Since google code has been closed I have moved the project to Github: https://github.com/yuemingl/Futureye_v2

All times are GMT 4. The time now is 13:56. 