diff --git a/docs/api/JuliaFEM.elasticity_solver.md b/docs/api/JuliaFEM.elasticity_solver.md
index 4581c44..baafae9 100644
--- a/docs/api/JuliaFEM.elasticity_solver.md
+++ b/docs/api/JuliaFEM.elasticity_solver.md
@@ -1,2 +1,154 @@
# JuliaFEM.elasticity_solver
+## Internal
+
+---
+
+
+#### assemble!(fe, eldofs_, I, V) [¶](#method__assemble.1)
+Assemble global RHS to I,V ready for sparse format
+
+Parameters
+----------
+fe : local vector
+eldofs_ : Array
+ degrees of freedom
+I,V : Arrays for sparse matrix
+
+Notes
+-----
+eldofs can also be node ids for convenience. In that case dimension
+is calculated and eldofs are "extended" to problem dimension.
+
+
+*source:*
+[JuliaFEM/src/elasticity_solver.jl:171](https://github.com/JuliaFEM/JuliaFEM.jl/tree/e4544d2b8d034469617e908134cbae74008ba189/src/elasticity_solver.jl#L171)
+
+---
+
+
+#### assemble!(ke, eldofs_, I, J, V) [¶](#method__assemble.2)
+Assemble global stiffness matrix to I,J,V ready for sparse format
+
+Parameters
+----------
+ke : local matrix
+eldofs_ : Array
+ degrees of freedom
+I,J,V : Arrays for sparse matrix
+
+Notes
+-----
+eldofs can also be node ids for convenience. In that case dimension
+is calculated and eldofs are "extended" to problem dimension.
+
+
+*source:*
+[JuliaFEM/src/elasticity_solver.jl:130](https://github.com/JuliaFEM/JuliaFEM.jl/tree/e4544d2b8d034469617e908134cbae74008ba189/src/elasticity_solver.jl#L130)
+
+---
+
+
+#### calc_local_matrices!(X, u, R, Kt, N, dNdchi, lambda_, mu_, ipoints, iweights) [¶](#method__calc_local_matrices.1)
+Calculate local tangent stiffness matrix and residual force vector R = T - F
+
+
+*source:*
+[JuliaFEM/src/elasticity_solver.jl:68](https://github.com/JuliaFEM/JuliaFEM.jl/tree/e4544d2b8d034469617e908134cbae74008ba189/src/elasticity_solver.jl#L68)
+
+---
+
+
+#### eliminate_boundary_conditions(dirichletbc, I, J, V) [¶](#method__eliminate_boundary_conditions.1)
+Eliminate Dirichlet boundary conditions from matrix
+
+Parameters
+----------
+dirichletbc : array [dim x nnodes]
+I, J, V : sparse matrix arrays
+
+Returns
+-------
+I, J, V : boundary conditions removed
+
+Notes
+-----
+pros:
+- matrix assembly remains positive definite
+cons:
+- maybe inefficient because of extra sparse matrix operations. (It's hard to remove stuff from sparse matrix.)
+- if u != 0 in dirichlet boundary requires extra care
+
+Raises
+------
+Exception, if displacement boundary conditions given, i.e.
+DX=2 for some node, for example.
+
+
+
+*source:*
+[JuliaFEM/src/elasticity_solver.jl:218](https://github.com/JuliaFEM/JuliaFEM.jl/tree/e4544d2b8d034469617e908134cbae74008ba189/src/elasticity_solver.jl#L218)
+
+---
+
+
+#### eliminate_boundary_conditions(dirichletbc, I, V) [¶](#method__eliminate_boundary_conditions.2)
+Eliminate Dirichlet boundary conditions from vector
+
+Parameters
+----------
+dirichletbc : array [dim x nnodes]
+I, V : sparse vector arrays
+
+Returns
+-------
+I, V : boundary conditions removed
+
+Notes
+-----
+pros:
+- matrix assembly remains positive definite
+cons:
+- maybe inefficient because of extra sparse matrix operations. (It's hard to remove stuff from sparse matrix.)
+- if u != 0 in dirichlet boundary requires extra care
+
+Raises
+------
+Exception, if displacement boundary conditions given, i.e.
+DX=2 for some node, for example.
+
+
+*source:*
+[JuliaFEM/src/elasticity_solver.jl:257](https://github.com/JuliaFEM/JuliaFEM.jl/tree/e4544d2b8d034469617e908134cbae74008ba189/src/elasticity_solver.jl#L257)
+
+---
+
+
+#### interpolate{T<:Real}(field::Array{T<:Real, 1}, basis::Function, ip) [¶](#method__interpolate.1)
+Interpolate field variable using basis functions f for point ip.
+This function tries to be as general as possible and allows interpolating
+lot of different fields.
+
+Parameters
+----------
+field :: Array{Number, dim}
+ Field variable
+basis :: Function
+ Basis functions
+ip :: Array{Number, 1}
+ Point to interpolate
+
+
+*source:*
+[JuliaFEM/src/elasticity_solver.jl:30](https://github.com/JuliaFEM/JuliaFEM.jl/tree/e4544d2b8d034469617e908134cbae74008ba189/src/elasticity_solver.jl#L30)
+
+---
+
+
+#### solve_elasticity_increment!(X, u, du, elmap, nodalloads, dirichletbc, lambda, mu, N, dNdchi, ipoints, iweights) [¶](#method__solve_elasticity_increment.1)
+Solve one increment of elasticity problem
+
+
+*source:*
+[JuliaFEM/src/elasticity_solver.jl:278](https://github.com/JuliaFEM/JuliaFEM.jl/tree/e4544d2b8d034469617e908134cbae74008ba189/src/elasticity_solver.jl#L278)
+
diff --git a/docs/api/JuliaFEM.md b/docs/api/JuliaFEM.md
index fb09699..442b178 100644
--- a/docs/api/JuliaFEM.md
+++ b/docs/api/JuliaFEM.md
@@ -1,2 +1,104 @@
# JuliaFEM
+## Exported
+
+---
+
+
+#### add_elements(model, elements) [¶](#method__add_elements.1)
+Add new elements to model.
+Parameters
+----------
+list of dicts, dict = {element_type => elcode, elids => [node ids..]}
+
+Examples
+--------
+Create two tet4 element and add them:
+
+>>> m = new_model()
+>>> const TET4 = 0x6
+>>> el1 = Dict("element_type" => TET4, "node_ids" => [1, 2, 3, 4])
+>>> el2 = Dict("element_type" => TET4, "node_ids" => [4, 3, 2, 1])
+
+In dict key means element id
+
+>>> elements = Dict(1 => el1, 2 => el2)
+>>> add_elements(m, elements)
+
+
+
+*source:*
+[JuliaFEM/src/JuliaFEM.jl:102](https://github.com/JuliaFEM/JuliaFEM.jl/tree/e4544d2b8d034469617e908134cbae74008ba189/src/JuliaFEM.jl#L102)
+
+---
+
+
+#### get_elements(model, element_ids) [¶](#method__get_elements.1)
+Get subset of elements from model.
+Parameters
+----------
+element_ids : list of ints
+ Element id numbers
+
+Returns
+-------
+Dict
+{element_type = XXX, node_ids = [a, b, c, d, e, ..., n]}
+
+
+*source:*
+[JuliaFEM/src/JuliaFEM.jl:126](https://github.com/JuliaFEM/JuliaFEM.jl/tree/e4544d2b8d034469617e908134cbae74008ba189/src/JuliaFEM.jl#L126)
+
+---
+
+
+#### get_field(field_type, field_name) [¶](#method__get_field.1)
+Get field from model.
+
+Parameters
+----------
+field_type : Dict()
+ Target topology (model.model, model.nodes, model.elements,
+ model.element_nodes, model.element_gauss
+field_name : str
+ Field name
+
+create_if_doesnt_exist : bool, optional
+ If field doesn't exists, create one and return empty field
+
+Raises
+------
+Error, if field not found and create_if_doesnt_exist == false
+
+
+*source:*
+[JuliaFEM/src/JuliaFEM.jl:67](https://github.com/JuliaFEM/JuliaFEM.jl/tree/e4544d2b8d034469617e908134cbae74008ba189/src/JuliaFEM.jl#L67)
+
+---
+
+
+#### new_model() [¶](#method__new_model.1)
+Initialize empty model.
+
+Parameters
+----------
+None
+
+Returns
+-------
+New model struct
+
+
+*source:*
+[JuliaFEM/src/JuliaFEM.jl:41](https://github.com/JuliaFEM/JuliaFEM.jl/tree/e4544d2b8d034469617e908134cbae74008ba189/src/JuliaFEM.jl#L41)
+
+---
+
+
+#### JuliaFEM.Model [¶](#type__model.1)
+Basic model
+
+
+*source:*
+[JuliaFEM/src/JuliaFEM.jl:19](https://github.com/JuliaFEM/JuliaFEM.jl/tree/e4544d2b8d034469617e908134cbae74008ba189/src/JuliaFEM.jl#L19)
+
diff --git a/docs/api/index.md b/docs/api/index.md
index a8f397f..bcc500b 100644
--- a/docs/api/index.md
+++ b/docs/api/index.md
@@ -1,3 +1,43 @@
# API-INDEX
+## MODULE: JuliaFEM
+
+---
+
+## Methods [Exported]
+
+[add_elements(model, elements)](JuliaFEM.md#method__add_elements.1) Add new elements to model.
+
+[get_elements(model, element_ids)](JuliaFEM.md#method__get_elements.1) Get subset of elements from model.
+
+[get_field(field_type, field_name)](JuliaFEM.md#method__get_field.1) Get field from model.
+
+[new_model()](JuliaFEM.md#method__new_model.1) Initialize empty model.
+
+---
+
+## Types [Exported]
+
+[JuliaFEM.Model](JuliaFEM.md#type__model.1) Basic model
+
+## MODULE: JuliaFEM.elasticity_solver
+
+---
+
+## Methods [Internal]
+
+[assemble!(fe, eldofs_, I, V)](JuliaFEM.elasticity_solver.md#method__assemble.1) Assemble global RHS to I,V ready for sparse format
+
+[assemble!(ke, eldofs_, I, J, V)](JuliaFEM.elasticity_solver.md#method__assemble.2) Assemble global stiffness matrix to I,J,V ready for sparse format
+
+[calc_local_matrices!(X, u, R, Kt, N, dNdchi, lambda_, mu_, ipoints, iweights)](JuliaFEM.elasticity_solver.md#method__calc_local_matrices.1) Calculate local tangent stiffness matrix and residual force vector R = T - F
+
+[eliminate_boundary_conditions(dirichletbc, I, J, V)](JuliaFEM.elasticity_solver.md#method__eliminate_boundary_conditions.1) Eliminate Dirichlet boundary conditions from matrix
+
+[eliminate_boundary_conditions(dirichletbc, I, V)](JuliaFEM.elasticity_solver.md#method__eliminate_boundary_conditions.2) Eliminate Dirichlet boundary conditions from vector
+
+[interpolate{T<:Real}(field::Array{T<:Real, 1}, basis::Function, ip)](JuliaFEM.elasticity_solver.md#method__interpolate.1) Interpolate field variable using basis functions f for point ip.
+
+[solve_elasticity_increment!(X, u, du, elmap, nodalloads, dirichletbc, lambda, mu, N, dNdchi, ipoints, iweights)](JuliaFEM.elasticity_solver.md#method__solve_elasticity_increment.1) Solve one increment of elasticity problem
+