Deckelnik/Eliott method implementation (2)

I was able to significantly improve and speed up my code, which was needed to run more tests with finer meshes:

- Added ILU preconditioner and changed iterative solver from GMRES to CGS (conjugate gradient squared)

- Switched mesh generation to 'dumb' gridding with squares or hexagons. Since we don't actually need a mesh for any specific domain, I generate a very simple grid that contains required region, and restrict to that grid. Another benefit is that all triangles come out exactly same, so I simplified a number of matrix assembly steps by computing stiffness and mass sub-matrices once.

- Profiled my code and removed several bottlenecks. Most important of those was fixing the ways I handled the sparse matrix

This allowed me to run a couple more refinements for last week's results. Unfortunately it looks like I was wrong, I do not have O(h^2) convergence:

Update: I found an error in my reference triangle calculations (item 2 in list above).
Update 2: Above error with indexing only affected square meshes. Equilateral triangles in hex meshes are not affected by this error.

I have also checked again if my overall FEM implementation was ok, using a simple example u=cos(pi x)cos(pi y) on a [-1 1]x[-1 1] square. It was correct:

 

 So, I am still unclear why I did not get the Deckelnik unfitted method to work.