FREEFEM: a PDE solver
Summary
The project freefem has two branches: the 2D solver (freefem++) is mature and well documented and used by many while the 3D solver (freefem3D) is still under construction but one is welcome to use it.
The code is open source and can be downloaded from the freefem web site
To solve a Partial Differential Equation (PDE) or set of equations with freefem one must first write an algorithm which uses only linear systems (one or more) of PDE as intermediate steps.
Then each system must be well posed, i.e. provided with admissible boundary conditions and written in variational form. For instance Dirichlet's problem
defines in an open set D of the plane a scalar function v(x,y) at every point (x,y) of D; v would be the electrostatic potential in the vacuum with charge density f(x,y) around objects of boundary əD at potential g(x,y).
The automatic triangulation, the finite element solver and the graphics are integrated in freefem and driven by a freefem script.

Example
The Dirichlet problem above can be solved with freefem++ by the following script
border C(t=0,2*pi){x=5*cos(t); y=5*sin(t);}
mesh Th = buildmesh (C(50));
fespace Vh(Th,P1); // says that we will use linear triangular elements
Vh u,v; // defines u and v as piecewise-P1 continuous functions
func f= x*y; // definition of a function called f
real cpu=clock(); // optional
solve Poisson(u,v)
= int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) // bilinear part
- int2d(Th)( f*v) // right hand side of the variational form
+ on(C,u=0) ; // Dirichlet boundary condition (g=0)
plot(u);
cout << " CPU time = " << clock()-cpu << endl;