diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 51e4a29..b8edd56 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -111,13 +111,13 @@ using Tensors  # For basis functions (Vec type)
 import Calculus  # For symbolic differentiation in basis generation
 
 import FEMSparse
+import FEMQuad  # Still using vendor FEMQuad for now
 
-@reexport using FEMBase
-import FEMBase: get_unknown_field_name, get_unknown_field_dimension,
-                assemble!, update!, initialize!
-using FEMBase: get_problems
+# Note: Consolidating FEMBase and FEMBasis into JuliaFEM
+# Previously: @reexport using FEMBase
+# Now: Include files directly below
 
-# Consolidate FEMBasis.jl into src/basis/ (Phase 1)
+# Consolidate jl into src/basis/ (Phase 1)
 include("basis/abstract.jl")
 include("basis/subs.jl")
 include("basis/vandermonde.jl")
@@ -135,13 +135,31 @@ include("basis/nurbs_surface.jl")
 include("basis/nurbs_solid.jl")
 include("basis/math.jl")
 
+# Consolidate FEMBase.jl into src/ (Phase 1 continued)
+# Order matters: fields → types → sparse → elements → integrate → problems → assembly
+include("fields/fields.jl")          # Field system (DCTI, DVTI, etc.)
+include("core_types.jl")             # Node, IP, IntegrationPoint
+
+# Compatibility shim: Create FEMBase module for vendor packages EARLY
+# This must come before preprocess.jl or any code that uses FEMBase.something
+include("fembase_compat.jl")
+
+include("sparse/sparse.jl")          # SparseMatrixCOO, SparseVectorCOO
+include("elements/elements.jl")      # Element type and interface
+include("elements/elements_lagrange.jl")  # Lagrange element specifics
+include("elements/integrate.jl")     # Integration utilities
+include("assembly/problems.jl")      # Problem types
+include("assembly/assembly.jl")      # Assembly framework
+include("solvers/solvers_base.jl")   # Base solver types
+include("analysis.jl")               # Analysis and AbstractResultsWriter
+
 using TimerOutputs
 export @timeit, print_timer
 
 import Base: getindex, setindex!, convert, length, size, isapprox,
-             similar, first, last, vec,
-             ==, +, -, *, /, haskey, copy, push!, isempty, empty!,
-             append!, read, copy
+   similar, first, last, vec,
+   ==, +, -, *, /, haskey, copy, push!, isempty, empty!,
+   append!, read, copy
 
 using AbaqusReader
 using AsterReader
@@ -167,17 +185,17 @@ export assemble!, postprocess!
 include("problems_mortar.jl")
 include("problems_mortar_3d.jl")
 export calculate_normals, calculate_normals!, project_from_slave_to_master,
-       project_from_master_to_slave, Mortar, get_slave_elements,
-       get_polygon_clip
+   project_from_master_to_slave, Mortar, get_slave_elements,
+   get_polygon_clip
 include("io.jl")
 export Xdmf, h5file, xmffile, xdmf_filter, new_dataitem, update_xdmf!, save!
 include("solvers.jl")
 export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver,
-       get_unknown_field_name, get_formulation_type, get_problems,
-       get_field_problems, get_boundary_problems,
-       get_field_assembly, get_boundary_assembly,
-       initialize!, create_projection, eliminate_interior_dofs,
-       is_field_problem, is_boundary_problem
+   get_unknown_field_name, get_formulation_type, get_problems,
+   get_field_problems, get_boundary_problems,
+   get_field_assembly, get_boundary_assembly,
+   initialize!, create_projection, eliminate_interior_dofs,
+   is_field_problem, is_boundary_problem
 include("solvers_modal.jl")
 export Modal
 include("problems_contact.jl")
@@ -188,12 +206,12 @@ export Contact
 module Preprocess
 end
 
-using FEMBase, SparseArrays, LinearAlgebra
+using SparseArrays, LinearAlgebra
 include("preprocess.jl")
 export create_elements, Mesh, add_node!, add_nodes!,
-       add_element_to_element_set!, add_node_to_node_set!,
-       find_nearest_nodes, find_nearest_node, reorder_element_connectivity!,
-       create_node_set_from_element_set!, filter_by_element_set
+   add_element_to_element_set!, add_node_to_node_set!,
+   find_nearest_nodes, find_nearest_node, reorder_element_connectivity!,
+   create_node_set_from_element_set!, filter_by_element_set
 include("preprocess_abaqus_reader.jl")
 export abaqus_read_mesh, create_surface_elements, create_nodal_elements
 include("preprocess_aster_reader.jl")
@@ -206,8 +224,8 @@ end
 
 include("postprocess_utils.jl")
 export calc_nodal_values!, get_nodal_vector, get_nodal_dict, copy_field!,
-       calculate_area, calculate_center_of_mass, calculate_second_moment_of_mass,
-       extract
+   calculate_area, calculate_center_of_mass, calculate_second_moment_of_mass,
+   extract
 
 include("deprecations.jl")
 
@@ -215,13 +233,13 @@ export SparseMatrixCOO, SparseVectorCOO, optimize!, resize_sparse
 export DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment
 export FieldProblem, BoundaryProblem, Problem, Node, Element, Assembly
 export Poi1, Seg2, Seg3, Tri3, Tri6, Tri7, Quad4, Quad8, Quad9,
-       Tet4, Tet10, Pyr5, Wedge6, Wedge15, Hex8, Hex20, Hex27
+   Tet4, Tet10, Pyr5, Wedge6, Wedge15, Hex8, Hex20, Hex27
 export update!, add_elements!, get_unknown_field_name, add!,
-       is_field_problem, is_boundary_problem, get_gdofs,
-       initialize!, get_integration_points, group_by_element_type,
-       get_unknown_field_dimension, get_connectivity
+   is_field_problem, is_boundary_problem, get_gdofs,
+   initialize!, get_integration_points, group_by_element_type,
+   get_unknown_field_dimension, get_connectivity
 export get_nonzero_rows, get_local_coordinates, inside, IP, get_element_type,
-       get_elements, AbstractProblem, IntegrationPoint, filter_by_element_type,
-       get_element_id, get_nonzero_columns, resize_sparse, resize_sparsevec
+   get_elements, AbstractProblem, IntegrationPoint, filter_by_element_type,
+   get_element_id, get_nonzero_columns, resize_sparse, resize_sparsevec
 
 end
diff --git a/src/analysis.jl b/src/analysis.jl
new file mode 100644
index 0000000..ae501e9
--- /dev/null
+++ b/src/analysis.jl
@@ -0,0 +1,105 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+abstract type AbstractAnalysis end
+abstract type AbstractResultsWriter end
+
+mutable struct Analysis{A<:AbstractAnalysis}
+    name :: String
+    problems :: Vector{Problem}
+    fields :: Dict{String, AbstractField}
+    results_writers :: Vector{AbstractResultsWriter}
+    properties :: A
+end
+
+"""
+    Analysis(A, analysis_name)
+
+Create a new analysis of type `A`, where `A` is a subtype of `AbstractAnalysis`.
+Analysis can be e.g. `Linear` for linear quasistatic analysis, `Nonlinear` for
+nonlinear quasistatic analysis or `Modal` for natural frequency analysis.
+
+# Examples
+
+```julia
+analysis = Analysis(Linear, "linear quasistatic analysis of beam structure")
+```
+"""
+function Analysis(::Type{A}, name::String="$A Analysis") where A<:AbstractAnalysis
+    analysis = Analysis{A}(name, [], Dict(), [], A())
+    @info("Creating a new analysis of type $A with name `$name`.")
+    return analysis
+end
+
+function add_problem!(analysis::Analysis, problems::Problem...)
+    for problem in problems
+        @info("Adding problem `$(problem.name)` to analysis `$(analysis.name)`.")
+        push!(analysis.problems, problem)
+    end
+    return nothing
+end
+
+add_problems!(analysis::Analysis, problems::Union{Vector, Tuple}) = add_problem!(analysis, problems...)
+
+"""
+    add_problems!(analysis, problem...)
+
+Add problem(s) to analysis.
+
+# Examples
+
+Add two problems, `beam` and `bc`, to linear quasistatic analysis:
+```julia
+beam = Problem(Beam, "beam structure", 6)
+bc = Problem(Dirichlet, "fix", 6, "displacement")
+analysis = Analysis(Linear, "linear quasistatic analysis")
+add_problems!(analysis, beam, bc)
+```
+"""
+add_problems!(analysis::Analysis, problems...) = add_problem!(analysis, problems...)
+
+function get_problems(analysis::Analysis)
+    return analysis.problems
+end
+
+function get_problem(analysis::Analysis, problem_name)
+    problems = filter(p -> p.name == problem_name, get_problems(analysis))
+    if length(problems) != 1
+        error("Several problem with name $problem_name found from analysis $(analysis.name).")
+    end
+    return first(problems)
+end
+
+function add_results_writer!(analysis::Analysis, writer::W) where W<:AbstractResultsWriter
+    push!(analysis.results_writers, writer)
+    return nothing
+end
+
+function get_results_writers(analysis::Analysis)
+    return analysis.results_writers
+end
+
+function run!(::Analysis{A}) where A<:AbstractAnalysis
+    @info("This is a placeholder function for running an analysis $A for a set of problems.")
+end
+
+function write_results!(::Analysis{A}, ::W) where {A<:AbstractAnalysis, W<:AbstractResultsWriter}
+    @info("Writing the results of analysis $A is not supported by a results writer $W")
+    return nothing
+end
+
+function write_results!(analysis)
+    results_writers = get_results_writers(analysis)
+    if isempty(results_writers)
+        @info("No result writers attached to the analysis $(analysis.name). " *
+              "In order to get results of the analysis stored to the disk, one " *
+              "must attach some results writer to the analysis using " *
+              "add_results_writer!, e.g. xdmf_writer = Xdmf(\"results\"); " *
+              "add_results_writer!(analysis, xdmf_writer)")
+        return nothing
+    end
+    for results_writer in results_writers
+        write_results!(analysis, results_writer)
+    end
+    return nothing
+end
diff --git a/src/assembly/assembly.jl b/src/assembly/assembly.jl
new file mode 100644
index 0000000..501c0ed
--- /dev/null
+++ b/src/assembly/assembly.jl
@@ -0,0 +1,196 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+function isapprox(a1::Assembly, a2::Assembly)
+    T = isapprox(a1.K, a2.K)
+    T &= isapprox(a1.C1, a2.C1)
+    T &= isapprox(a1.C2, a2.C2)
+    T &= isapprox(a1.D, a2.D)
+    T &= isapprox(a1.f, a2.f)
+    T &= isapprox(a1.g, a2.g)
+    return T
+end
+
+function assemble_prehook!(::Problem, ::T) where T<:Number end
+
+function assemble_posthook!(::Problem, ::T) where T<:Number end
+
+"""
+    assemble_elements!(problem, assembly, elements, time)
+
+Assemble elements for problem.
+
+This should be overridden with own `assemble_elements!`-implementation.
+"""
+function assemble_elements!(problem::Problem, assembly::Assembly,
+                            elements::Vector{T}, time) where T<:AbstractElement{E} where E
+    elements2 = convert(Vector{Element}, elements)
+    assemble!(assembly, problem, elements2, time)
+end
+
+function assemble!(problem::Problem, time)
+
+    assemble_prehook!(problem, time)
+    elements = get_elements(problem)
+    assembly = get_assembly(problem)
+
+    if !isempty(assembly)
+        @warn("Problem assembly is not empty before assembling. This is probably " *
+              "causing unexpected results. To remove old assembly, use " *
+              "`empty!(problem.assembly)`", typeof(problem), problem.name)
+        assemble_posthook!(problem, time)
+        return nothing
+    end
+
+    if isempty(elements)
+        @warn("There is no elements defined in problem. Before assembling a " *
+              "problem, elements must be added using " *
+              "`add_elements!(problem, elements)`.", typeof(problem), problem.name)
+        assemble_posthook!(problem, time)
+        return nothing
+    end
+
+    first_element = first(elements)
+    unknown_field_name = get_unknown_field_name(problem)
+    if !haskey(first_element, unknown_field_name)
+        #=
+        warn("Assembling elements for problem $(problem.name): seems that ",
+             "problem is uninitialized. To initialize problem, use ",
+             "`initialize!(problem, time)`.")
+        info("Initializing problem $(problem.name) at time $time automatically.")
+        =#
+        initialize!(problem, time)
+    end
+
+    for (element_type, elements) in group_by_element_type(elements)
+        assemble_elements!(problem, assembly, elements, time)
+    end
+    assemble_posthook!(problem, time)
+    return nothing
+end
+
+function assemble!(problem::Problem)
+    @warn("assemble!(problem) will be deprecated. Use assemble!(problem, time)")
+    assemble!(problem, 0.0)
+end
+
+function assemble_mass_matrix!(problem::Problem, time::Float64)
+    if !isempty(problem.assembly.M)
+        @info("Mass matrix for is already assembled, not assembling.",
+              problem.name)
+        return nothing
+    end
+    elements = get_elements(problem)
+    for (element_type, elements) in group_by_element_type(get_elements(problem))
+        assemble_mass_matrix!(problem::Problem, elements, time)
+    end
+    return nothing
+end
+
+function assemble_mass_matrix!(problem::Problem, elements::Vector{E}, time) where E<:AbstractElement{M_,B} where {M_,B}
+    nnodes = length(first(elements))
+    dim = get_unknown_field_dimension(problem)
+    M = zeros(nnodes, nnodes)
+    N = zeros(1, nnodes)
+    NtN = zeros(nnodes, nnodes)
+    ldofs = zeros(Int, nnodes)
+    for element in elements
+        fill!(M, 0.0)
+        for ip in get_integration_points(element, 2)
+            detJ = element(ip, time, Val{:detJ})
+            rho = element("density", ip, time)
+            w = ip.weight*rho*detJ
+            eval_basis!(B, N, ip)
+            N = element(ip, time)
+            mul!(NtN, transpose(N), N)
+            rmul!(NtN, w)
+            for i=1:nnodes^2
+                M[i] += NtN[i]
+            end
+        end
+        for (i, j) in enumerate(get_connectivity(element))
+            @inbounds ldofs[i] = (j-1)*dim
+        end
+        for i=1:dim
+            add!(problem.assembly.M, ldofs.+i, ldofs.+i, M)
+        end
+    end
+    return
+end
+
+"""
+    assemble_mass_matrix!(problem, elements::Vector{Element{Tet10}}, time)
+
+Assemble Tet10 mass matrices using special method. If Tet10 has constant metric
+if can be integrated analytically to gain performance.
+"""
+function assemble_mass_matrix!(problem::Problem, elements::Vector{E}, time) where E<:AbstractElement{M_, Tet10} where M_
+    nnodes = length(Tet10)
+    dim = get_unknown_field_dimension(problem)
+    M = zeros(nnodes, nnodes)
+    N = zeros(1, nnodes)
+    NtN = zeros(nnodes, nnodes)
+    ldofs = zeros(Int, nnodes)
+
+    M_CM = 1.0/2520.0 * [
+         6  1  1  1 -4 -6 -4 -4 -6 -6
+         1  6  1  1 -4 -4 -6 -6 -4 -6
+         1  1  6  1 -6 -4 -4 -6 -6 -4
+         1  1  1  6 -6 -6 -6 -4 -4 -4
+        -4 -4 -6 -6 32 16 16 16 16  8
+        -6 -4 -4 -6 16 32 16  8 16 16
+        -4 -6 -4 -6 16 16 32 16  8 16
+        -4 -6 -6 -4 16  8 16 32 16 16
+        -6 -4 -6 -4 16 16  8 16 32 16
+        -6 -6 -4 -4  8 16 16 16 16 32]
+
+    function is_CM(::AbstractElement{M, Tet10}, X; rtol=1.0e-6) where M
+        isapprox(X[5],  1/2*(X[1]+X[2]); rtol=rtol) || return false
+        isapprox(X[6],  1/2*(X[2]+X[3]); rtol=rtol) || return false
+        isapprox(X[7],  1/2*(X[3]+X[1]); rtol=rtol) || return false
+        isapprox(X[8],  1/2*(X[1]+X[4]); rtol=rtol) || return false
+        isapprox(X[9],  1/2*(X[2]+X[4]); rtol=rtol) || return false
+        isapprox(X[10], 1/2*(X[3]+X[4]); rtol=rtol) || return false
+        return true
+    end
+
+
+    n_CM = 0
+    for element in elements
+        for (i, j) in enumerate(get_connectivity(element))
+            @inbounds ldofs[i] = (j-1)*dim
+        end
+
+        X = element("geometry", time)
+        rho = element("density", time)
+        if is_CM(element, X) && length(rho) == 1
+            ip = (1.0/3.0, 1.0/3.0, 1.0/3.0)
+            detJ = element(ip, time, Val{:detJ})
+            rho = element("density", ip, time)
+            CM_s = detJ*rho
+            n_CM += 1
+            for i=1:dim
+                add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, CM_s * M_CM)
+            end
+        else
+            fill!(M, 0.0)
+            for ip in get_integration_points(element, 2)
+                detJ = element(ip, time, Val{:detJ})
+                rho = element("density", ip, time)
+                w = ip.weight*rho*detJ
+                eval_basis!(Tet10, N, ip)
+                N = element(ip, time)
+                mul!(NtN, transpose(N), N)
+                rmul!(NtN, w)
+                for i=1:nnodes^2
+                    M[i] += NtN[i]
+                end
+            end
+            for i=1:dim
+                add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, M)
+            end
+        end
+    end
+    @info("$n_CM of $(length(elements)) was constant metric.")
+    return
+end
diff --git a/src/assembly/problems.jl b/src/assembly/problems.jl
new file mode 100644
index 0000000..90e8ecb
--- /dev/null
+++ b/src/assembly/problems.jl
@@ -0,0 +1,484 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+abstract type AbstractProblem end
+abstract type FieldProblem<:AbstractProblem end
+abstract type BoundaryProblem<:AbstractProblem end
+abstract type MixedProblem<:AbstractProblem end
+
+"""
+General linearized problem to solve
+    (K₁+K₂)Δu  +   C1'*Δλ = f₁+f₂
+         C2Δu  +     D*Δλ = g
+"""
+mutable struct Assembly
+
+    M :: SparseMatrixCOO  # mass matrix
+
+    # for field assembly
+    K :: SparseMatrixCOO   # stiffness matrix
+    Kg :: SparseMatrixCOO  # geometric stiffness matrix
+    f :: SparseMatrixCOO   # force vector
+    fg :: SparseMatrixCOO  #
+
+    # for boundary assembly
+    C1 :: SparseMatrixCOO
+    C2 :: SparseMatrixCOO
+    D :: SparseMatrixCOO
+    g :: SparseMatrixCOO
+    c :: SparseMatrixCOO
+
+    u :: Vector{Float64}  # solution vector u
+    u_prev :: Vector{Float64}  # previous solution vector u
+    u_norm_change :: Real  # change of norm in u
+
+    la :: Vector{Float64}  # solution vector la
+    la_prev :: Vector{Float64}  # previous solution vector u
+    la_norm_change :: Real # change of norm in la
+
+    removed_dofs :: Vector{Int} # manually remove dofs from assembly
+end
+
+function Assembly()
+    return Assembly(
+        SparseMatrixCOO(),
+        SparseMatrixCOO(),
+        SparseMatrixCOO(),
+        SparseMatrixCOO(),
+        SparseMatrixCOO(),
+        SparseMatrixCOO(),
+        SparseMatrixCOO(),
+        SparseMatrixCOO(),
+        SparseMatrixCOO(),
+        SparseMatrixCOO(),
+        [], [], Inf,
+        [], [], Inf,
+        [])
+end
+
+function empty!(assembly::Assembly)
+    empty!(assembly.M)
+    empty!(assembly.K)
+    empty!(assembly.Kg)
+    empty!(assembly.f)
+    empty!(assembly.fg)
+    empty!(assembly.C1)
+    empty!(assembly.C2)
+    empty!(assembly.D)
+    empty!(assembly.g)
+    empty!(assembly.c)
+end
+
+function isempty(assembly::Assembly)
+    T = isempty(assembly.M)
+    T &= isempty(assembly.K)
+    T &= isempty(assembly.Kg)
+    T &= isempty(assembly.f)
+    T &= isempty(assembly.fg)
+    T &= isempty(assembly.C1)
+    T &= isempty(assembly.C2)
+    T &= isempty(assembly.D)
+    T &= isempty(assembly.g)
+    T &= isempty(assembly.c)
+    return T
+end
+
+"""
+    Problem{P<:AbstractProblem}
+
+Defines a new problem of type `P`, where `P` characterizes the physics of the
+problem. `P` can be for example `Elasticity`, if the physics of the system is
+described by Cauchy's stress equation ∇⋅σ + b = ̈ρu, or `Heat`, if the physics
+of the problem is described by heat equation -∇⋅(k∇u) = f.
+
+"""
+mutable struct Problem{P<:AbstractProblem}
+    name :: AbstractString           # descriptive name for the problem
+    dimension :: Int                 # degrees of freedom per node
+    parent_field_name :: AbstractString # (optional) name of the parent field e.g. "displacement"
+    elements :: Vector{Element}
+    dofmap :: Dict{Element, Vector{Int}} # connects the element local dofs to the global dofs
+    assembly :: Assembly
+    fields :: Dict{String, AbstractField}
+    postprocess_fields :: Vector{String}
+    properties :: P
+end
+
+"""
+    Problem(problem_type, problem_name, problem_dimension)
+
+Construct a new field problem.
+
+`problem_type` must be a subtype of `FieldProblem` (`Elasticity`, `Heat`, etc..).
+`problem_dimensions` is the number of degrees of freedom each node is containing.
+
+# Examples
+
+To create vector-valued elasticity problem, having 3 dofs / node:
+```julia
+problem1 = Problem(Elasticity, "test problem", 3)
+```
+
+To create scalar-valued Poisson problem:
+```julia
+problem2 = Problem(Heat, "test problem 2", 1)
+```
+
+"""
+function Problem(::Type{P}, name::AbstractString, dimension::Int) where P<:FieldProblem
+    parent_field_name = "none"
+    elements = []
+    dofmap = Dict()
+    assembly = Assembly()
+    fields = Dict()
+    postprocess_fields = Vector()
+    properties = P()
+    problem = Problem{P}(name, dimension, parent_field_name, elements, dofmap,
+                         assembly, fields, postprocess_fields, properties)
+    @info("Creating a new problem of type $P, having name `$name` and " *
+          "dimension $dimension dofs/node.")
+    return problem
+end
+
+"""
+    Problem(problem_type, problem_name, problem_dimension, parent_field_name)
+
+Construct a new boundary problem.
+
+`problem_type` must be a subtype of `BoundaryProblem` (`Dirichlet`, `Contact`,
+etc..). `problem_dimensions` is the number of degrees of freedom each node is
+containing. `parent_field_name` is describing the field, where the boundary
+problem is affecting.
+
+# Examples
+
+To create a Dirichlet boundary condition for a vector-valued elasticity problem,
+having 3 dofs / node:
+```julia
+bc1 = Problem(Dirichlet, "fix displacement on support", 3, "displacement")
+```
+
+To create a Dirichlet boundary condition for scalar-valued Poisson problem:
+```julia
+bc2 = Problem(Dirichlet, "fix surface temperature", 1, "temperature")
+```
+"""
+function Problem(::Type{P}, name, dimension, parent_field_name) where P<:BoundaryProblem
+    elements = []
+    dofmap = Dict()
+    assembly = Assembly()
+    fields = Dict()
+    postprocess_fields = Vector()
+    properties = P()
+    problem = Problem{P}(name, dimension, parent_field_name, elements, dofmap,
+                         assembly, fields, postprocess_fields, properties)
+    @info("Creating a new boundary problem of type $P, having name `$name` and " *
+          "dimension $dimension dofs/node. This boundary problems fixes field " *
+          "`$parent_field_name`.")
+    return problem
+end
+
+function get_formulation_type(::Problem)
+    return :incremental
+end
+
+"""
+    get_unknown_field_dimension(problem)
+
+Return the dimension of the unknown field of this problem.
+"""
+function get_unknown_field_dimension(problem::Problem)
+    return problem.dimension
+end
+
+"""
+    get_unknown_field_name(problem)
+
+Default function if unknown field name is not defined for some problem.
+"""
+function get_unknown_field_name(::P) where P<:AbstractProblem
+    @warn("The name of unknown field (e.g. displacement, temperature, ...) of the " *
+          "problem type must be given by defining a function " *
+          "`get_unknown_field_name(::$P)`")
+    return "N/A"
+end
+
+""" Return the name of the unknown field of this problem. """
+function get_unknown_field_name(problem::Problem{P}) where P
+    return get_unknown_field_name(problem.properties)
+end
+
+""" Return the name of the parent field of this (boundary) problem. """
+function get_parent_field_name(problem::Problem{P}) where P<:BoundaryProblem
+    return problem.parent_field_name
+end
+
+function get_unknown_field_name(::P) where P<:BoundaryProblem
+    return "lambda"
+end
+
+is_field_problem(::Problem) = false
+is_field_problem(::Problem{P}) where {P<:FieldProblem} = true
+is_boundary_problem(::Problem) = false
+is_boundary_problem(::Problem{P}) where {P<:BoundaryProblem} = true
+
+function get_elements(problem::Problem)
+    return problem.elements
+end
+
+function update!(problem::P, attr::Pair{String, String}...) where P<:AbstractProblem
+    for (name, value) in attr
+        setfield!(problem, Meta.parse(name), Meta.parse(value))
+    end
+end
+
+"""
+    function initialize!(problem_type, element_name, time)
+
+Initialize the element ready for calculation, where `problem_type` is the type
+of the problem (Elasticity, Dirichlet, etc.), `element_name` is the name of a
+constructed element (see Element(element_type, connectivity_vector)) and `time`
+is the starting time of the initializing process.
+"""
+function initialize!(problem::Problem, element::AbstractElement, time::Float64)
+    field_name = get_unknown_field_name(problem)
+    field_dim = get_unknown_field_dimension(problem)
+    nnodes = length(element)
+    if field_dim == 1  # scalar field
+        empty_field = tuple(zeros(nnodes)...)
+    else  # vector field
+        # FIXME: the most effective way to do
+        # ([0.0,0.0], [0.0,0.0], ..., [0.0,0.0]) ?
+        empty_field = tuple(map((x)->zeros(field_dim)*x, 1:nnodes)...)
+    end
+
+    # initialize primary field
+    if !haskey(element, field_name)
+        update!(element, field_name, time => empty_field)
+    end
+
+    # if a boundary problem, initialize also a field for the main problem
+    is_boundary_problem(problem) || return
+    field_name = get_parent_field_name(problem)
+    if !haskey(element, field_name)
+        update!(element, field_name, time => empty_field)
+    end
+end
+
+function initialize!(problem::Problem, time::Float64=0.0)
+    for element in get_elements(problem)
+        initialize!(problem, element, time)
+    end
+end
+
+function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector)
+
+    # resize & fill with zeros vectors if length mismatch with current solution
+
+    if length(u) != length(assembly.u)
+        resize!(assembly.u, length(u))
+        fill!(assembly.u, 0.0)
+    end
+
+    if length(la) != length(assembly.la)
+        resize!(assembly.la, length(la))
+        fill!(assembly.la, 0.0)
+    end
+
+    # copy current solutions to previous ones and add/replace new solution
+    # TODO: here we have couple of options and they need to be clarified
+    # for total formulation we are solving total quantity Ku = f while in
+    # incremental formulation we solve KΔu = f and u = u + Δu
+    assembly.u_prev = copy(assembly.u)
+    assembly.la_prev = copy(assembly.la)
+
+    if get_formulation_type(problem) == :total
+        assembly.u = u
+        assembly.la = la
+    elseif get_formulation_type(problem) == :incremental
+        assembly.u += u
+        assembly.la = la
+    elseif get_formulation_type(problem) == :forwarddiff
+        assembly.u += u
+        assembly.la += la
+    else
+        @info("$(problem.name): unknown formulation type, don't know what to do with results")
+        error("serious failure with problem formulation: $(get_formulation_type(problem))")
+    end
+
+    # calculate change of norm
+    assembly.u_norm_change = norm(assembly.u - assembly.u_prev)
+    assembly.la_norm_change = norm(assembly.la - assembly.la_prev)
+    return assembly.u, assembly.la
+end
+
+"""
+    get_global_solution(problem, assembly)
+
+Return a global solution (u, la) for a problem.
+
+Notes
+-----
+If the length of solution vector != number of nodes, i.e. the field dimension is
+something else than 1, reshape vectors so that their length matches to the
+number of nodes. This helps to get nodal results easily.
+"""
+function get_global_solution(problem::Problem, assembly::Assembly)
+    u = assembly.u
+    la = assembly.la
+    field_dim = get_unknown_field_dimension(problem)
+    if field_dim == 1
+        return u, la
+    else
+        nnodes = round(Int, length(u)/field_dim)
+        u = reshape(u, field_dim, nnodes)
+        u = Vector{Float64}[u[:,i] for i in 1:nnodes]
+        la = reshape(la, field_dim, nnodes)
+        la = Vector{Float64}[la[:,i] for i in 1:nnodes]
+        return u, la
+    end
+end
+
+function update!(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64) where P<:FieldProblem
+    u, la = get_global_solution(problem, assembly)
+    field_name = get_unknown_field_name(problem)
+    # update solution u for elements
+    for element in elements
+        connectivity = get_connectivity(element)
+        update!(element, field_name, time => tuple(u[connectivity]...))
+    end
+end
+
+function update!(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64) where P<:BoundaryProblem
+    u, la = get_global_solution(problem, assembly)
+    parent_field_name = get_parent_field_name(problem) # displacement
+    field_name = get_unknown_field_name(problem) # lambda
+    # update solution and lagrange multipliers for boundary elements
+    for element in elements
+        connectivity = get_connectivity(element)
+        update!(element, parent_field_name, time => tuple(u[connectivity]...))
+        update!(element, field_name, time => tuple(la[connectivity]...))
+    end
+end
+
+"""
+    add_element!(problem, element1, element2, ...)
+
+Add element(s) to the problem.
+"""
+function add_element!(problem, elements...)
+    for element in elements
+        push!(problem.elements, element)
+    end
+    return nothing
+end
+
+"""
+    add_elements!(problem, element_set_1, element_set_2, ...)
+
+Add vectors/tuples of element(s) to the problem.
+"""
+function add_elements!(problem, element_sets::Union{Vector,Tuple}...)
+    for elements in element_sets
+        nelements = length(elements)
+        @info("Adding $nelements elements to problem `$(problem.name)`")
+        add_element!(problem, elements...)
+    end
+    return nothing
+end
+
+add_elements!(problem, elements::Element...) = add_element!(problem, elements...)
+
+function add_elements!(problem, elements_or_lists_of_elements...)
+    for item in elements_or_lists_of_elements
+        add_elements!(problem, item)
+    end
+end
+
+get_assembly(problem::Problem) = problem.assembly
+Base.length(problem::Problem) = length(problem.elements)
+
+function update!(problem::Problem, field_name::AbstractString, data)
+    #if haskey(problem.fields, field_name)
+    #    update!(problem.fields[field_name], field_name::AbstractString, data)
+    #else
+    #    problem.fields[field_name] = Field(data)
+    #end
+    update!(problem.elements, field_name::AbstractString, data)
+end
+
+function haskey(problem::Problem, field_name::AbstractString)
+    return haskey(problem.fields, field_name)
+end
+
+function getindex(problem::Problem, field_name::String)
+    return problem.fields[field_name]
+end
+
+#""" Return field calculated to nodal points for elements in problem p. """
+function (problem::Problem)(field_name::String, time::Float64)
+    #if haskey(problem, field_name)
+    #    return problem[field_name](time)
+    #end
+    f = Dict{Int, Any}()
+    for element in get_elements(problem)
+        haskey(element, field_name) || continue
+        for (c, v) in zip(get_connectivity(element), element(field_name, time))
+            if haskey(f, c)
+                if !isapprox(f[c], v)
+                    @info("several values for single node when returning field $field_name")
+                    @info("already have: $(f[c]), and trying to set $v")
+                end
+            else
+                f[c] = v
+            end
+        end
+    end
+    #f == nothing && return f
+    #update!(problem, field_name, time => f)
+    return f
+end
+
+function push!(problem::Problem, elements...)
+    push!(problem.elements, elements...)
+end
+
+function push!(problem::Problem, elements_::Vector...)
+    for elements in elements_
+        push!(problem.elements, elements...)
+    end
+end
+
+"""
+    set_gdofs!(problem, element)
+
+Set element global degrees of freedom.
+"""
+function set_gdofs!(problem, element, dofs)
+    problem.dofmap[element] = dofs
+end
+
+"""
+    get_gdofs(problem, element)
+
+Return the global degrees of freedom for element.
+
+First make lookup from problem dofmap. If not defined there, make implicit
+assumption that dofs follow formula `gdofs = [dim*(nid-1)+j for j=1:dim]`,
+where `nid` is node id and `dim` is the dimension of problem. This formula
+arranges dofs so that first comes all dofs of node 1, then node 2 and so on:
+(u11, u12, u13, u21, u22, u23, ..., un1, un2, un3) for 3 dofs/node setting.
+"""
+function get_gdofs(problem::Problem, element::AbstractElement)
+    if haskey(problem.dofmap, element)
+        return problem.dofmap[element]
+    end
+    conn = get_connectivity(element)
+    if length(conn) == 0
+        error("element connectivity not defined, cannot determine global ",
+              "degrees of freedom for element #: $(element.id)")
+    end
+    dim = get_unknown_field_dimension(problem)
+    gdofs = [dim*(i-1)+j for i in conn for j=1:dim]
+    return gdofs
+end
diff --git a/src/basis/abstract.jl b/src/basis/abstract.jl
index 0aac09d..aff1233 100644
--- a/src/basis/abstract.jl
+++ b/src/basis/abstract.jl
@@ -2,7 +2,7 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE
 
 # AbstractBasis type and interface
-# Consolidated from FEMBasis.jl package
+# Consolidated from jl package
 
 using Tensors
 using LinearAlgebra
diff --git a/src/basis/create_basis.jl b/src/basis/create_basis.jl
index be6cfea..59173e9 100644
--- a/src/basis/create_basis.jl
+++ b/src/basis/create_basis.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 __precompile__(false)
 
@@ -22,7 +22,7 @@ function calculate_interpolation_polynomials(p, V)
         N = Expr(:call, :+)
         for (ai, bi) in zip(solution, args)
             isapprox(ai, 0.0) && continue
-            push!(N.args, Calculus.simplify( :( $ai * $bi ) ))
+            push!(N.args, Calculus.simplify(:($ai * $bi)))
         end
         push!(basis, N)
     end
@@ -55,17 +55,17 @@ function create_basis(name, description, X::Vector{<:Vecish{D}}, basis::Vector)
     return create_basis(name, description, Vec.(X), basis, dbasis)
 end
 
-function create_basis(name, description, X::Vector{<:Vecish{D, T}}, basis, dbasis) where {D, T}
+function create_basis(name, description, X::Vector{<:Vecish{D,T}}, basis, dbasis) where {D,T}
     N = length(X)
     @debug "create basis given basis functions and derivatives" name description X basis dbasis
 
     Q = Expr(:block)
-    for i=1:N
+    for i = 1:N
         push!(Q.args, :(N[$i] = $(basis[i])))
     end
 
     V = Expr(:block)
-    for i=1:N
+    for i = 1:N
         push!(V.args, :(dN[$i] = Vec(float.(tuple($(dbasis[:, i]...))))))
     end
 
diff --git a/src/basis/lagrange_hexahedrons.jl b/src/basis/lagrange_hexahedrons.jl
index 053f11e..783bb54 100644
--- a/src/basis/lagrange_hexahedrons.jl
+++ b/src/basis/lagrange_hexahedrons.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 code = create_basis_and_eval(
     :Hex8,
diff --git a/src/basis/lagrange_pyramids.jl b/src/basis/lagrange_pyramids.jl
index 516e066..ef45615 100644
--- a/src/basis/lagrange_pyramids.jl
+++ b/src/basis/lagrange_pyramids.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 # Kaltenbacher, Manfred. Numerical simulation of mechatronic sensors and actuators: finite elements for computational multiphysics. Springer, 2015.
 code = create_basis_and_eval(
diff --git a/src/basis/lagrange_quadrangles.jl b/src/basis/lagrange_quadrangles.jl
index 927f8cc..eb5f863 100644
--- a/src/basis/lagrange_quadrangles.jl
+++ b/src/basis/lagrange_quadrangles.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 code = create_basis_and_eval(
     :Quad4,
diff --git a/src/basis/lagrange_segments.jl b/src/basis/lagrange_segments.jl
index 6750d06..c221ae7 100644
--- a/src/basis/lagrange_segments.jl
+++ b/src/basis/lagrange_segments.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 code = create_basis_and_eval(
     :Seg2,
diff --git a/src/basis/lagrange_tetrahedrons.jl b/src/basis/lagrange_tetrahedrons.jl
index efe803d..abf51d2 100644
--- a/src/basis/lagrange_tetrahedrons.jl
+++ b/src/basis/lagrange_tetrahedrons.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 code = create_basis_and_eval(
     :Tet4,
diff --git a/src/basis/lagrange_triangles.jl b/src/basis/lagrange_triangles.jl
index 7a443dc..ea7f5d8 100644
--- a/src/basis/lagrange_triangles.jl
+++ b/src/basis/lagrange_triangles.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 code = create_basis_and_eval(
     :Tri3,
diff --git a/src/basis/lagrange_wedges.jl b/src/basis/lagrange_wedges.jl
index c78cafc..295eb78 100644
--- a/src/basis/lagrange_wedges.jl
+++ b/src/basis/lagrange_wedges.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 # Kaltenbacher, Manfred. Numerical simulation of mechatronic sensors and actuators: finite elements for computational multiphysics. Springer, 2015.
 create_basis_and_eval(
diff --git a/src/basis/math.jl b/src/basis/math.jl
index 88bac36..cd67bc0 100644
--- a/src/basis/math.jl
+++ b/src/basis/math.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 """
     interpolate(B, T, xi)
@@ -148,7 +148,7 @@ BasisInfo(Tri3)
 
 # output
 
-FEMBasis.BasisInfo{FEMBasis.Tri3,Float64}([0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0; 0.0 0.0], [0.0 0.0; 0.0 0.0], 0.0)
+BasisInfo{Tri3,Float64}([0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0; 0.0 0.0], [0.0 0.0; 0.0 0.0], 0.0)
 
 ```
 
diff --git a/src/basis/nurbs.jl b/src/basis/nurbs.jl
index f104f53..a0455e0 100644
--- a/src/basis/nurbs.jl
+++ b/src/basis/nurbs.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 import Base: size, length
 
diff --git a/src/basis/nurbs_segment.jl b/src/basis/nurbs_segment.jl
index 06ac906..1e5bf05 100644
--- a/src/basis/nurbs_segment.jl
+++ b/src/basis/nurbs_segment.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 """ NURBS segment. """
 mutable struct NSeg <: AbstractBasis{1}
diff --git a/src/basis/nurbs_solid.jl b/src/basis/nurbs_solid.jl
index 6a4730d..3046ef2 100644
--- a/src/basis/nurbs_solid.jl
+++ b/src/basis/nurbs_solid.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 mutable struct NSolid <: AbstractBasis{3}
     order_u :: Int
diff --git a/src/basis/nurbs_surface.jl b/src/basis/nurbs_surface.jl
index bf20abb..f83d600 100644
--- a/src/basis/nurbs_surface.jl
+++ b/src/basis/nurbs_surface.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 mutable struct NSurf <: AbstractBasis{2}
     order_u :: Int
diff --git a/src/basis/subs.jl b/src/basis/subs.jl
index a4aa6e2..cda2cc2 100644
--- a/src/basis/subs.jl
+++ b/src/basis/subs.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 function subs(p::Number, ::Any)
     return p
diff --git a/src/basis/vandermonde.jl b/src/basis/vandermonde.jl
index 254d1b2..d24702e 100644
--- a/src/basis/vandermonde.jl
+++ b/src/basis/vandermonde.jl
@@ -1,5 +1,5 @@
 # This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE
 
 """
     vandermonde_matrix(polynomial, coordinates)
diff --git a/src/core_types.jl b/src/core_types.jl
new file mode 100644
index 0000000..6d92054
--- /dev/null
+++ b/src/core_types.jl
@@ -0,0 +1,72 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+const Node = Vector{Float64}
+
+abstract type AbstractPoint end
+
+mutable struct Point{P<:AbstractPoint}
+    id :: Int
+    weight :: Float64
+    coords :: Tuple{Vararg{Float64}}
+    fields :: Dict{String, AbstractField}
+    properties :: P
+end
+
+function setindex!(point::Point, val::Pair{Float64, T}, field_name) where T
+    point.fields[field_name] = field(val)
+end
+
+function getindex(point::Point, field_name)
+    return point.fields[field_name]
+end
+
+function getindex(point::Point, idx::Int)
+    return point.coords[idx]
+end
+
+function haskey(point::Point, field_name)
+    return haskey(point.fields, field_name)
+end
+
+function (point::Point)(field_name, time)
+    interpolate(point.fields[field_name], time)
+end
+
+function Base.iterate(point::Point)
+    return Base.iterate(point.coords)
+end
+
+function Base.iterate(point::Point, i::Int)
+    return Base.iterate(point.coords, i)
+end
+
+function update!(point::Point, field_name, val::Pair{Float64, T}) where T
+    if haskey(point, field_name)
+        update!(point[field_name], val)
+    else
+        point[field_name] = val
+    end
+end
+
+#= TODO: in future
+type Node <: AbstractPoint
+end
+
+type MaterialPoint <: AbstractPoint
+end
+=#
+
+struct IntegrationPoint <: AbstractPoint
+end
+
+const IP = Point{IntegrationPoint}
+
+function IP(id, weight, coords::Tuple)
+    return IP(id, weight, coords, Dict(), IntegrationPoint())
+end
+
+function IP(id, weight, coords::Vector)
+    @warn "Consider giving coordinates as tuple."
+    return IP(id, weight, tuple(coords...), Dict(), IntegrationPoint())
+end
diff --git a/src/elements/elements.jl b/src/elements/elements.jl
new file mode 100644
index 0000000..cb1addf
--- /dev/null
+++ b/src/elements/elements.jl
@@ -0,0 +1,527 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+"""
+    AbstractFieldSet{N<:Int}
+
+Abstract supertype for all field sets, where `N` is the length of the discrete
+fields (typically is the number of the nodes in element).
+"""
+abstract type AbstractFieldSet{N} end
+
+"""
+    EmptyFieldSet{N} <: AbstractFieldSet{N}
+
+Empty field set used as a default for all elements.
+"""
+struct EmptyFieldSet{N} <: AbstractFieldSet{N}
+end
+
+const DefaultFieldSet = EmptyFieldSet
+
+"""
+    AbstractElement{M<:AbstractFieldSet, B<:AbstractBasis}
+
+Abstract supertype for all elements.
+"""
+abstract type AbstractElement{M<:AbstractFieldSet, B<:AbstractBasis} end
+
+mutable struct Element{M,B} <: AbstractElement{M,B}
+    id :: Int
+    connectivity :: Vector{Int}
+    integration_points :: Vector{IP}
+    dfields :: Dict{Symbol, AbstractField}
+    sfields :: M
+    properties :: B
+end
+
+"""
+    Element(topology, connectivity)
+
+Construct a new element where `topology` is the topological type of the element
+and connectivity contains node numbers where element is connected.
+
+# Topological types
+
+## 1d elements
+- `Seg2`
+- `Seg3`
+
+## 2d elements
+- `Tri3`
+- `Tri6`
+- `Tri7`
+- `Quad4`
+- `Quad8`
+- `Quad9`
+
+## 3d elements
+- `Tet4`
+- `Tet10`
+- `Hex8`
+- `Hex20`
+- `Hex27`
+- `Pyr5`
+- `Wedge6`
+- `Wedge15`
+
+# Examples
+
+```julia
+element = Element(Tri3, (1, 2, 3))
+```
+"""
+function Element(::Type{T}, connectivity::NTuple{N, Int}) where {N, T<:AbstractBasis}
+    return Element(T, DefaultFieldSet, connectivity)
+end
+
+function Element(::Type{T}, ::Type{M}, connectivity::NTuple{N, Int}) where {N, M<:AbstractFieldSet, T<:AbstractBasis}
+    element_id = -1
+    topology = T()
+    integration_points = Point{IntegrationPoint}[]
+    dfields = Dict{Symbol,AbstractField}()
+    sfields = M{N}()
+    element = Element(element_id, collect(connectivity), integration_points,
+                      dfields, sfields, topology)
+    return element
+end
+
+function Element(::Type{T}, connectivity::Vector{Int}) where T<:AbstractBasis
+    return Element(T, (connectivity...,))
+end
+
+function get_element_id(element::AbstractElement)
+    return element.id
+end
+
+function get_element_type(::AbstractElement{M,T}) where {M,T}
+    return T
+end
+
+function is_element_type(::AbstractElement{M,T}, element_type) where {M,T}
+    return T === element_type
+end
+
+function filter_by_element_type(element_type, elements)
+    return Iterators.filter(element -> is_element_type(element, element_type), elements)
+end
+
+function get_connectivity(element::AbstractElement)
+    return element.connectivity
+end
+
+"""
+    group_by_element_type(elements)
+
+Given a vector of elements, group elements by element type to several vectors.
+Returns a dictionary, where key is the element type and value is a vector
+containing all elements of type `element_type`.
+"""
+function group_by_element_type(elements)
+    eltypes = map(T -> typeof(T), elements)
+    elgroups = Dict(T => T[] for T in eltypes)
+    for element in elements
+        T = typeof(element)
+        push!(elgroups[T], element)
+    end
+    return elgroups
+end
+
+### dfields - dynamically defined fields
+
+# This is the "old" field system, where fields are defined to dictionary.
+# It is known that this approach is having a performance issue caused by
+# type instability.
+
+function has_dfield(element, field_name)
+    return haskey(element.dfields, field_name)
+end
+
+function get_dfield(element, field_name)
+    return getindex(element.dfields, field_name)
+end
+
+function create_dfield!(element, field_name, field_::AbstractField)
+    T = typeof(field_)
+    if has_dfield(element, field_name)
+        @debug("Replacing the content of a field $field_name with a new field of type $T.")
+    else
+        @debug("Creating a new dfield $field_name of type $T")
+    end
+    element.dfields[field_name] = field_
+    return
+end
+
+function create_dfield!(element, field_name, field_data)
+    create_dfield!(element, field_name, field(field_data))
+end
+
+function update_dfield!(element, field_name, field_data)
+    if has_dfield(element, field_name)
+        field = get_dfield(element, field_name)
+        @debug("Update $field_name with data $field_data")
+        update_field!(field, field_data)
+    else
+        create_dfield!(element, field_name, field_data)
+    end
+end
+
+# A helper function to pick element data from dictionary
+function pick_data_(element, field_data)
+    connectivity = get_connectivity(element)
+    N = length(connectivity)
+    picked_data = ntuple(i -> getindex(field_data, connectivity[i]), N)
+    return picked_data
+end
+
+function update_dfield!(element, field_name, (time, field_data)::Pair{Float64, Dict{Int,V}}) where V
+    update_dfield!(element, field_name, time => pick_data_(element, field_data))
+end
+
+function update_dfield!(element, field_name, field_data::Dict{Int,V}) where V
+    update_dfield!(element, field_name, pick_data_(element, field_data))
+end
+
+function update_dfield!(element, field_name, field_data::Function)
+    if hasmethod(field_data, Tuple{Element, Any, Any})
+        element.dfields[field_name] = field((ip, time) -> field_data(element, ip, time))
+    else
+        element.dfields[field_name] = field(field_data)
+    end
+end
+
+function interpolate_dfield(element, field_name, time)
+    field = get_dfield(element, field_name)
+    return interpolate(field, time)
+end
+
+### sfields statically defined fields
+
+# A new-style field system, where fields are defined in sfields <: AbstractFieldSet
+# during the initialization of element.
+
+function has_sfield(element, field_name)
+    return isdefined(element.sfields, field_name)
+end
+
+function get_sfield(element, field_name)
+    return getfield(element.sfields, field_name)
+end
+
+function update_sfield!(element, field_name, field_data)
+    field = get_sfield(element, field_name)
+    update!(field, field_data)
+end
+
+function interpolate_sfield(element, field_name, time)
+    field = get_sfield(element, field_name)
+    return interpolate(field, time)
+end
+
+### dfield & sfield -- common routines
+
+function has_field(element, field_name)
+    return has_sfield(element, field_name) || has_dfield(element, field_name)
+end
+
+function get_field(element, field_name)
+    if has_sfield(element, field_name)
+        return get_sfield(element, field_name)
+    else
+        return get_dfield(element, field_name)
+    end
+end
+
+function update_field!(element, field_name, field_data)
+    if has_sfield(element, field_name)
+        update_sfield!(element, field_name, field_data)
+    else
+        update_dfield!(element, field_name, field_data)
+    end
+end
+
+function interpolate_field(element, field_name::Symbol, time)
+    if has_sfield(element, field_name)
+        return interpolate_sfield(element, field_name, time)
+    elseif has_dfield(element, field_name)
+        return interpolate_dfield(element, field_name, time)
+    else
+        error("Cannot interpolate from field $field_name: no such field.")
+    end
+end
+
+function interpolate(element::AbstractElement, field_name, time)
+    return interpolate_field(element, field_name, time)
+end
+
+function update_field!(elements::Vector{Element}, field_name, field_data)
+    for element in elements
+        update_field!(element, field_name, field_data)
+    end
+end
+
+# Update fields when given a dictionary or time => dictionary:
+# pick data from dictionary diven by the connectivity information of element
+#=
+function update_field!(element::AbstractElement, field::F,
+                       data::Dict{T,V}) where {F<:DVTI,T,V}
+    connectivity = get_connectivity(element)
+    N = length(connectivity)
+    picked_data = ntuple(i -> data[connectivity[i]], N)
+    update_field!(field, picked_data)
+end
+
+function update_field!(element::AbstractElement, field::F,
+                       ddata::Pair{Float64, Dict{T,V}}) where {F<:DVTV,T,V}
+    time, data = ddata
+    connectivity = get_connectivity(element)
+    N = length(connectivity)
+    picked_data = ntuple(i -> data[connectivity[i]], N)
+    update_field!(field, time => picked_data)
+end
+=#
+
+"""
+    interpolate(element, field_name, time)
+
+Interpolate field `field_name` from element at given `time`.
+
+# Example
+```
+element = Element(Seg2, [1, 2])
+data1 = Dict(1 => 1.0, 2 => 2.0)
+data2 = Dict(1 => 2.0, 2 => 3.0)
+update!(element, "my field", 0.0 => data1)
+update!(element, "my field", 1.0 => data2)
+interpolate(element, "my field", 0.5)
+
+# output
+
+(1.5, 2.5)
+
+```
+"""
+function interpolate(element::AbstractElement, field_name::String, time::Float64)
+    field = element[field_name]
+    result = interpolate(field, time)
+    if isa(result, Dict)
+        connectivity = get_connectivity(element)
+        return tuple((result[i] for i in connectivity)...)
+    else
+        return result
+    end
+end
+
+function info_update_field(elements, field_name, data)
+    nelements = length(elements)
+    @info("Updating field `$field_name` for $nelements elements.")
+end
+
+function info_update_field(elements, field_name, data::Float64)
+    nelements = length(elements)
+    @info("Updating field `$field_name` => $data for $nelements elements.")
+end
+
+"""
+    update!(elements, field_name, data)
+
+Given a list of elements, field name and data, update field to elements. Data
+is passed directly to the `field`-function.
+
+# Examples
+
+Create two elements with topology `Seg2`, one is connecting to nodes (1, 2) and
+the other is connecting to (2, 3). Some examples of updating fields:
+
+```julia
+elements = [Element(Seg2, [1, 2]), Element(Seg2, [2, 3])]
+X = Dict(1 => 0.0, 2 => 1.0, 3 => 2.0)
+u = Dict(1 => 0.0, 2 => 0.0, 3 => 0.0)
+update!(elements, "geometry", X)
+update!(elements, "displacement", 0.0 => u)
+update!(elements, "youngs modulus", 210.0e9)
+update!(elements, "time-dependent force", 0.0 => 0.0)
+update!(elements, "time-dependent force", 1.0 => 100.0)
+```
+
+When using dictionaries in definition of fields, key of dictionary corresponds
+to node id, that is, updating field `geometry` in the example above is updating
+values `(0.0, 1.0)` for the first elements and values `(1.0, 2.0)` to the second
+element. For time dependent field, syntax `time => data` is used. If field is
+initialized without time-dependency, it cannot be changed to be time-dependent
+afterwards. If unsure, it's better to initialize field with time dependency.
+
+"""
+function update!(elements, field_name, data)
+    info_update_field(elements, field_name, data)
+    for element in elements
+        update!(element, field_name, data)
+    end
+end
+
+
+## Interpolate fields in spatial direction
+
+const ConstantField = Union{DCTI, DCTV}
+const VariableFields = Union{DVTV, DVTI}
+const DictionaryFields = Union{DVTVd, DVTId}
+
+function interpolate_field(::AbstractElement, field::ConstantField, ip, time)
+    return interpolate_field(field, time)
+end
+
+function interpolate_field(element::AbstractElement, field::VariableFields, ip, time)
+    data = interpolate_field(field, time)
+    basis = get_basis(element, ip, time)
+    N = length(basis)
+    return sum(data[i]*basis[i] for i=1:N)
+end
+
+function interpolate_field(element::AbstractElement, field::DictionaryFields, ip, time)
+    data = interpolate_field(field, time)
+    basis = element(ip, time)
+    N = length(element)
+    c = get_connectivity(element)
+    return sum(data[c[i]]*basis[i] for i=1:N)
+end
+
+function interpolate_field(::AbstractElement, field::CVTV, ip, time)
+    return field(ip, time)
+end
+
+function interpolate(element::AbstractElement, field_name, ip, time)
+    field = get_field(element, field_name)
+    interpolate_field(element, field, ip, time)
+end
+
+
+## Other stuff
+
+function get_basis(element::AbstractElement{M,B}, ip, ::Any) where {M,B}
+    T = typeof(first(ip))
+    N = zeros(T, 1, length(element))
+    eval_basis!(B, N, tuple(ip...))
+    return N
+end
+
+function get_dbasis(element::AbstractElement{M,B}, ip, ::Any) where {M,B}
+    T = typeof(first(ip))
+    dN = zeros(T, size(element)...)
+    eval_dbasis!(B, dN, tuple(ip...))
+    return dN
+end
+
+function (element::Element)(ip, time::Float64=0.0)
+    return get_basis(element, ip, time)
+end
+
+#"""
+#Examples
+#julia> el = Element(Quad4, [1, 2, 3, 4]);
+#julia> el([0.0, 0.0], 0.0, 1)
+#1x4 Array{Float64,2}:
+# 0.25  0.25  0.25  0.25
+#julia> el([0.0, 0.0], 0.0, 2)
+#2x8 Array{Float64,2}:
+# 0.25  0.0   0.25  0.0   0.25  0.0   0.25  0.0
+# 0.0   0.25  0.0   0.25  0.0   0.25  0.0   0.25
+#"""
+function (element::Element)(ip, time::Float64, dim::Int)
+    dim == 1 && return get_basis(element, ip, time)
+    Ni = vec(get_basis(element, ip, time))
+    N = zeros(dim, length(element)*dim)
+    for i=1:dim
+        N[i,i:dim:end] += Ni
+    end
+    return N
+end
+
+function (element::Element)(ip, time, ::Type{Val{:Jacobian}})
+    X = element("geometry", time)
+    J = jacobian(element.properties, X, ip)
+    return J
+end
+
+function (element::Element)(ip, time::Float64, ::Type{Val{:detJ}})
+    J = element(ip, time, Val{:Jacobian})
+    n, m = size(J)
+    if n == m  # volume element
+        return det(J)
+    end
+    JT = transpose(J)
+    if size(JT, 2) == 1  # boundary of 2d problem, || ∂X/∂ξ ||
+        return norm(JT)
+    else # manifold on 3d problem, || ∂X/∂ξ₁ × ∂X/∂ξ₂ ||
+        return norm(cross(JT[:,1], JT[:,2]))
+    end
+end
+
+function (element::Element)(ip, time::Float64, ::Type{Val{:Grad}})
+    J = element(ip, time, Val{:Jacobian})
+    return inv(J)*get_dbasis(element, ip, time)
+end
+
+function (element::Element)(field_name::String, ip, time::Float64, ::Type{Val{:Grad}})
+    X = element("geometry", time)
+    u = element(field_name, time)
+    return grad(element.properties, u, X, ip)
+end
+
+
+function get_integration_points(element::AbstractElement{E}) where E
+    # first time initialize default integration points
+    if length(element.integration_points) == 0
+        ips = get_integration_points(element.properties)
+        element.integration_points = [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
+    end
+    return element.integration_points
+end
+
+""" This is a special case, temporarily change order
+of integration scheme mainly for mass matrix.
+"""
+function get_integration_points(element::AbstractElement{E}, change_order::Int) where E
+    ips = get_integration_points(element.properties, Val{change_order})
+    return [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
+end
+
+""" Find inverse isoparametric mapping of element. """
+function get_local_coordinates(element::AbstractElement, X::Vector, time::Float64; max_iterations=10, tolerance=1.0e-6)
+    haskey(element, "geometry") || error("element geometry not defined, cannot calculate inverse isoparametric mapping")
+    dim = size(element, 1)
+    dim == length(X) || error("manifolds not supported.")
+    xi = zeros(dim)
+    dX = element("geometry", xi, time) - X
+    for i=1:max_iterations
+        J = element(xi, time, Val{:Jacobian})'
+        xi -= J \ dX
+        dX = element("geometry", xi, time) - X
+        norm(dX) < tolerance && return xi
+    end
+    debug("get_local_coordinates", X, dX, xi)
+    error("Unable to find inverse isoparametric mapping for element $element for X = $X")
+end
+
+""" Test is X inside element. """
+function inside(element::AbstractElement{M,B}, X, time) where {M,B}
+    xi = get_local_coordinates(element, X, time)
+    return inside(B, xi)
+end
+
+## Convenience functions
+
+# element("displacement", 0.0)
+function (element::Element)(field_name::String, time::Float64)
+    return interpolate(element, field_name, time)
+end
+
+# element("displacement", (0.0, 0.0), 0.0)
+function (element::Element)(field_name::String, ip, time::Float64)
+    return interpolate(element, field_name, ip, time)
+end
+
+function element_info!(bi::BasisInfo{T}, element::AbstractElement{M,T}, ip, time) where {M,T}
+    X = interpolate(element, "geometry", time)
+    eval_basis!(bi, X, ip)
+    return bi.J, bi.detJ, bi.N, bi.grad
+end
diff --git a/src/elements/elements_lagrange.jl b/src/elements/elements_lagrange.jl
new file mode 100644
index 0000000..2a0e718
--- /dev/null
+++ b/src/elements/elements_lagrange.jl
@@ -0,0 +1,53 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+struct Poi1 <: AbstractBasis{0} end
+
+function get_basis(::E, ::Any, ::Any) where E<:AbstractElement{M,Poi1} where M
+    return [1]
+end
+
+function get_dbasis(::E, ::Any, ::Any) where E<:AbstractElement{M,Poi1} where M
+    return [0]
+end
+
+function (::Element{M,Poi1})(::Any, ::Float64, ::Type{Val{:detJ}}) where M
+    return 1.0
+end
+
+function get_integration_order(::Poi1)
+    return 1
+end
+
+function get_integration_points(::Poi1, ::Int)
+    return [(1.0, (0.0,))]
+end
+
+function size(::Type{Poi1})
+    return (0, 1)
+end
+
+function length(::Type{Poi1})
+    return 1
+end
+
+function get_reference_element_coordinates(::Type{Poi1})
+    Vector{Float64}[[0.0]]
+end
+
+function inside(::Union{Type{Seg2},Type{Seg3},Type{Quad4},
+        Type{Quad8},Type{Quad9},Type{Pyr5},
+        Type{Hex8},Type{Hex20},
+        Type{Hex27}}, xi)
+    return all(-1.0 .<= xi .<= 1.0)
+end
+
+function inside(::Union{Type{Tri3},Type{Tri6},Type{Tri7},
+        Type{Tet4},Type{Tet10}}, xi)
+    return all(xi .>= 0.0) && (sum(xi) <= 1.0)
+end
+
+function get_reference_coordinates(::E) where E<:AbstractElement{M,B} where {M,B}
+    return get_reference_element_coordinates(B)
+end
+
diff --git a/src/elements/integrate.jl b/src/elements/integrate.jl
new file mode 100644
index 0000000..571da66
--- /dev/null
+++ b/src/elements/integrate.jl
@@ -0,0 +1,57 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+# Default number of integration points for each element. First rule is the
+# default integration rule returned by `get_integration_points(element)`.
+# Sometimes we want to increase integration order, e.g. when integrating mass
+# matrix or boundary conditions. For that reason, additional rules are provied
+# in list, so e.g. `get_integration_points(element, 1)` returns the second rule,
+# `get_integration_points(element, 2)` third rule and so on. Rules should be
+# ordered so that picking next one integrates more accurately.
+integration_rule_mapping = (
+    :Seg2 => (:GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5),
+    :Seg3 => (:GLSEG3, :GLSEG4, :GLSEG5),
+    :NSeg => (:GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5),
+    :Quad4 => (:GLQUAD4, :GLQUAD9, :GLQUAD16, :GLQUAD25),
+    :Quad8 => (:GLQUAD9, :GLQUAD16, :GLQUAD25),
+    :Quad9 => (:GLQUAD9, :GLQUAD16, :GLQUAD25),
+    :NSurf => (:GLQUAD9, :GLQUAD16, :GLQUAD25),
+    :Hex8 => (:GLHEX8, :GLHEX27, :GLHEX64, :GLHEX125),
+    :Hex20 => (:GLHEX27, :GLHEX64, :GLHEX125),
+    :Hex27 => (:GLHEX27, :GLHEX64, :GLHEX125),
+    :NSolid => (:GLHEX27, :GLHEX64, :GLHEX125),
+    :Tri3 => (:GLTRI1, :GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12),
+    :Tri6 => (:GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12),
+    :Tri7 => (:GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12),
+    :Tet4 => (:GLTET1, :GLTET4, :GLTET5, :GLTET15),
+    :Tet10 => (:GLTET4, :GLTET5, :GLTET15),
+    :Pyr5 => (:GLPYR5,),
+    :Wedge6 => (:GLWED6, :GLWED21),
+    :Wedge15 => (:GLWED21,))
+
+for (E, R) in integration_rule_mapping
+    for i in 1:length(R)
+        P = Val{R[i]}
+        order = Val{i - 1}
+        local code  # Explicitly declare as local to avoid warning
+        if isequal(i, 1)
+            code = quote
+                function get_integration_points(element::$E)
+                    return FEMQuad.get_quadrature_points($P)
+                end
+            end
+        else
+            code = quote
+                function get_integration_points(element::$E, ::Type{$order})
+                    return FEMQuad.get_quadrature_points($P)
+                end
+            end
+        end
+        eval(code)
+    end
+end
+
+# All good codes needs a special case. Here we have it: Poi1
+function get_integration_points(::Poi1)
+    [(1.0, (0.0,))]
+end
diff --git a/src/fembase_compat.jl b/src/fembase_compat.jl
new file mode 100644
index 0000000..8d48aaa
--- /dev/null
+++ b/src/fembase_compat.jl
@@ -0,0 +1,22 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE
+
+"""
+Compatibility shim for vendor packages that expect FEMBase types.
+
+Since we've consolidated FEMBase into JuliaFEM, we need to provide
+the FEMBase module namespace for backward compatibility with code
+that uses FEMBase.function_name().
+
+This creates a minimal FEMBase module with function forwarding.
+Type aliases are added after all types are defined.
+"""
+module FEMBase
+
+# Note: We can only create function aliases here, not type aliases,
+# because not all types have been defined yet when this module is included.
+
+# Forward declarations for functions that exist at this point
+# We'll add more after types are defined
+
+end # module FEMBase
diff --git a/src/fields/fields.jl b/src/fields/fields.jl
new file mode 100644
index 0000000..24146ce
--- /dev/null
+++ b/src/fields/fields.jl
@@ -0,0 +1,376 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+"""
+    AbstractField
+
+Abstract supertype for all fields in JuliaFEM.
+"""
+abstract type AbstractField end
+
+function length(f::F) where F<:AbstractField
+    return length(f.data)
+end
+
+function size(f::F) where F<:AbstractField
+    return size(f.data)
+end
+
+function ==(x::F, y) where F<:AbstractField
+    return ==(x.data, y)
+end
+
+function ==(x, y::F) where F<:AbstractField
+    return ==(x, y.data)
+end
+
+function ==(x::F, y::F) where F<:AbstractField
+    return ==(x.data, y.data)
+end
+
+function getindex(f::F, i::Int) where F<:AbstractField
+    return getindex(f.data, i)
+end
+
+function interpolate_field(field::AbstractField, ::Any)
+    return field.data
+end
+
+function update_field!(field::AbstractField, data)
+   field.data = data
+end
+
+"""
+    DCTI{T} <: AbstractField
+
+Discrete, constant, time-invariant field.
+
+This field is constant in both spatial direction and time direction,
+i.e. df/dX = 0 and df/dt = 0.
+
+# Example
+
+```jldoctest
+julia> DCTI(1)
+FEMBase.DCTI{Int64}(1)
+```
+"""
+mutable struct DCTI{T} <: AbstractField
+    data :: T
+end
+
+function getindex(field::DCTI, ::Int)
+    return field.data
+end
+
+"""
+    DVTI{N,T} <: AbstractField
+
+Discrete, variable, time-invariant field.
+
+This is constant in time direction, but not in spatial direction, i.e. df/dt = 0
+but df/dX != 0. The basic structure of data is `Tuple`, and it is implicitly
+assumed that length of field matches to the number of shape functions, so that
+interpolation in spatial direction works.
+
+# Example
+
+```jldoctest
+julia> DVTI(1, 2, 3)
+FEMBase.DVTI{3,Int64}((1, 2, 3))
+```
+"""
+mutable struct DVTI{N,T} <: AbstractField
+    data :: NTuple{N,T}
+end
+
+function DVTI(data...)
+    return DVTI(data)
+end
+
+"""
+    DCTV{T} <: AbstractField
+
+Discrete, constant, time variant field. This type of field can change in time
+direction but not in spatial direction.
+
+# Example
+
+Field having value 5 at time 0.0 and value 10 at time 1.0:
+
+```jldoctest
+julia> DCTV(0.0 => 5, 1.0 => 10)
+FEMBase.DCTV{Int64}(Pair{Float64,Int64}[0.0=>5, 1.0=>10])
+```
+
+"""
+mutable struct DCTV{T} <: AbstractField
+    data :: Vector{Pair{Float64,T}}
+end
+
+function DCTV(data::Pair{Float64,T}...) where T
+    return DCTV(collect(data))
+end
+
+function update_field!(f::DCTV, data::Pair{Float64, T}) where T
+    if isapprox(last(f.data).first, data.first)
+        f.data[end] = data
+    else
+        push!(f.data, data)
+    end
+end
+
+function interpolate_field(field::DCTV, time)
+    time < first(field.data).first && return first(field.data).second
+    time > last(field.data).first && return last(field.data).second
+    for i=reverse(1:length(field))
+        isapprox(field.data[i].first, time) && return field.data[i].second
+    end
+    for i=length(field.data):-1:2
+        t0 = field.data[i-1].first
+        t1 = field.data[i].first
+        if t0 < time < t1
+            y0 = field.data[i-1].second
+            y1 = field.data[i].second
+            dy = y1-y0
+            dt = t1-t0
+            return y0 + (time-t0)*dy/dt
+        end
+    end
+end
+
+"""
+    DVTV{N,T} <: AbstractField
+
+Discrete, variable, time variant field. The most general discrete field can
+change in both temporal and spatial direction.
+
+# Example
+
+```jldoctest
+julia> DVTV(0.0 => (1, 2), 1.0 => (2, 3))
+FEMBase.DVTV{2,Int64}(Pair{Float64,Tuple{Int64,Int64}}[0.0=>(1, 2), 1.0=>(2, 3)])
+```
+"""
+mutable struct DVTV{N,T} <: AbstractField
+    data :: Vector{Pair{Float64,NTuple{N,T}}}
+end
+
+function DVTV(data::Pair{Float64,NTuple{N,T}}...) where {N,T}
+    return DVTV(collect(data))
+end
+
+function update_field!(f::DVTV, data::Pair{Float64, NTuple{N,T}}) where {N,T}
+    if isapprox(last(f.data).first, data.first)
+        f.data[end] = data
+    else
+        push!(f.data, data)
+    end
+end
+
+function interpolate_field(field::DVTV{N,T}, time) where {N,T}
+    time < first(field.data).first && return first(field.data).second
+    time > last(field.data).first && return last(field.data).second
+    for i=reverse(1:length(field))
+        isapprox(field.data[i].first, time) && return field.data[i].second
+    end
+    for i=length(field.data):-1:2
+        t0 = field.data[i-1].first
+        t1 = field.data[i].first
+        if t0 < time < t1
+            y0 = field.data[i-1].second
+            y1 = field.data[i].second
+            dt = t1-t0
+            return map((a,b) -> a + (time-t0)*(b-a)/dt, y0, y1)
+        end
+    end
+end
+
+"""
+    CVTV <: AbstractField
+
+Continuous, variable, time variant field.
+
+# Example
+
+```jldoctest
+julia> f = CVTV((xi,t) -> xi*t)
+FEMBase.CVTV(#1)
+```
+"""
+mutable struct CVTV <: AbstractField
+    data :: Function
+end
+
+function (f::CVTV)(xi, time)
+    return f.data(xi, time)
+end
+
+"""
+    DVTId(X::Dict)
+
+Discrete, variable, time invariant dictionary field.
+"""
+mutable struct DVTId{T} <: AbstractField
+    data :: Dict{Int, T}
+end
+
+function update_field!(field::DVTId{T}, data::Dict{Int, T}) where T
+    merge!(field.data, data)
+end
+
+"""
+    DVTVd(time => data::Dict)
+
+Discrete, variable, time variant dictionary field.
+"""
+mutable struct DVTVd{T} <: AbstractField
+    data :: Vector{Pair{Float64,Dict{Int,T}}}
+end
+
+function DVTVd(data::Pair{Float64,Dict{Int,T}}...) where T
+    return DVTVd(collect(data))
+end
+
+function interpolate_field(field::DVTVd{T}, time) where T
+    time >= last(field.data).first && return last(field.data).second
+    time <= first(field.data).first && return first(field.data).second
+    for i=reverse(1:length(field))
+        isapprox(field.data[i].first, time) && return field.data[i].second
+    end
+    for i=length(field.data):-1:2
+        t0 = field.data[i-1].first
+        t1 = field.data[i].first
+        if t0 < time < t1
+            y0 = field.data[i-1].second
+            y1 = field.data[i].second
+            f = (time-t0)/(t1-t0)
+            new_data = empty(y0)
+            for i in keys(y0)
+                new_data[i] = f*y0[i] + (1-f)*y1[i]
+            end
+            return new_data
+        end
+    end
+end
+
+function update_field!(f::DVTVd, data::Pair{Float64,Dict{Int,T}}) where T
+    if isapprox(last(f.data).first, data.first)
+        f.data[end] = data
+    else
+        push!(f.data, data)
+    end
+end
+
+function new_field(data)
+    return DCTI(data)
+end
+
+function new_field(data...)
+    return DVTI(data)
+end
+
+function new_field(data::NTuple{N,T}) where {N,T}
+    return DVTI(data)
+end
+
+function new_field(data::Pair{Float64,T}...) where T
+    return DCTV(collect(data))
+end
+
+function new_field(data::Pair{Float64,NTuple{N,T}}...) where {N,T}
+    return DVTV(collect(data))
+end
+
+function new_field(data::Function)
+    return CVTV(data)
+end
+
+function new_field(data::Pair{Int, T}...) where T
+    return DVTId(Dict(data))
+end
+
+function new_field(data::Pair{Float64, NTuple{N, Pair{Int, T}}}...) where {N,T}
+    return DVTVd(collect(t => Dict(d) for (t, d) in data))
+end
+
+function new_field(data::Dict{Int,T}) where T
+    return DVTId(data)
+end
+
+function new_field(data::Pair{Float64, Dict{Int, T}}...) where T
+    return DVTVd(collect(data))
+end
+
+"""
+    field(x)
+
+Create new field. Field type is deduced from data type.
+"""
+function field(data...)
+    return new_field(data...)
+end
+
+"""
+    interpolate(field, time)
+
+Interpolate field in time direction.
+
+# Examples
+
+For time invariant fields [`DCTI`](@ref), [`DVTI`](@ref), [`DVTId`](@ref)
+solution is trivially the data inside field as fields does not depend from
+the time:
+
+```jldoctest
+julia> a = field(1.0)
+FEMBase.DCTI{Float64}(1.0)
+
+julia> interpolate(a, 0.0)
+1.0
+```
+
+```jldoctest
+julia> a = field((1.0, 2.0))
+FEMBase.DVTI{2,Float64}((1.0, 2.0))
+
+julia> interpolate(a, 0.0)
+(1.0, 2.0)
+```
+
+```jldoctest
+julia> a = field(1=>1.0, 2=>2.0)
+FEMBase.DVTId{Float64}(Dict(2=>2.0,1=>1.0))
+
+julia> interpolate(a, 0.0)
+Dict{Int64,Float64} with 2 entries:
+  2 => 2.0
+  1 => 1.0
+```
+
+DVTId trivial solution is returned. For time variant fields DCTV, DVTV, DVTVd
+linear interpolation is performed.
+
+# Other notes
+
+First algorithm checks that is time out of range, i.e. time is smaller than
+time of first frame or larger than last frame. If that is the case, return
+first or last frame. Secondly algorithm finds is given time exact match to
+time of some frame and return that frame. At last, we find correct bin so
+that t0 < time < t1 and use linear interpolation.
+
+"""
+function interpolate(field::AbstractField, time)
+    return interpolate_field(field, time)
+end
+
+"""
+    interpolate(a, b)
+
+A helper function for interpolate routines. Given iterables `a` and `b`,
+calculate c = aᵢbᵢ. Length of `a` can be less than `b`, but not vice versa.
+"""
+function interpolate(a, b)
+    @assert length(a) <= length(b)
+    return sum(a[i]*b[i] for i=1:length(a))
+end
diff --git a/src/preprocess.jl b/src/preprocess.jl
index 57e1bea..2b07dc1 100644
--- a/src/preprocess.jl
+++ b/src/preprocess.jl
@@ -122,7 +122,7 @@ end
 Add an element into the mesh. ´elid´ is the element id, ´eltype´ is the type of
 the element and ´connectivity´ is the connectivity of the element.
 """
-function FEMBase.add_element!(mesh::Mesh, elid, eltype, connectivity)
+function add_element!(mesh::Mesh, elid, eltype, connectivity)
     mesh.elements[elid] = connectivity
     mesh.element_types[elid] = eltype
     return nothing
@@ -133,7 +133,7 @@ end
 
 Add elements into the mesh.
 """
-function FEMBase.add_elements!(mesh::Mesh, elements::Dict{Int, Tuple{Symbol, Vector{Int}}})
+function add_elements!(mesh::Mesh, elements::Dict{Int, Tuple{Symbol, Vector{Int}}})
     for (elid, (eltype, elcon)) in elements
         add_element!(mesh, elid, eltype, elcon)
     end
diff --git a/src/solvers.jl b/src/solvers.jl
index e3992ea..1abd2f2 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -572,7 +572,7 @@ function has_converged(solver::Solver{Nonlinear})
 end
 
 """ Default solver for quasistatic nonlinear problems. """
-function FEMBase.run!(solver::Solver{Nonlinear})
+function run!(solver::Solver{Nonlinear})
 
     time = solver.properties.time
     problems = get_problems(solver)
@@ -638,7 +638,7 @@ function Linear()
     return Linear(0.0)
 end
 
-function FEMBase.run!(analysis::Analysis{Linear})
+function run!(analysis::Analysis{Linear})
     time = analysis.properties.time
     @info("Running linear quasistatic analysis `$(analysis.name)` at time $time.")
     problems = get_problems(analysis)
diff --git a/src/solvers/solvers_base.jl b/src/solvers/solvers_base.jl
new file mode 100644
index 0000000..99f1ef4
--- /dev/null
+++ b/src/solvers/solvers_base.jl
@@ -0,0 +1,53 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+mutable struct LinearSystem{Tv, Ti<:Integer}
+    M :: SparseMatrixCSC{Tv, Ti}
+    K :: SparseMatrixCSC{Tv, Ti}
+    Kg :: SparseMatrixCSC{Tv, Ti}
+    C1 :: SparseMatrixCSC{Tv, Ti}
+    C2 :: SparseMatrixCSC{Tv, Ti}
+    D :: SparseMatrixCSC{Tv, Ti}
+    f :: SparseVector{Tv, Ti}
+    fg :: SparseVector{Tv, Ti}
+    g :: SparseVector{Tv, Ti}
+    u :: SparseVector{Tv, Ti}
+    la :: SparseVector{Tv, Ti}
+    dim :: Int
+end 
+
+function LinearSystem(dim::Int)
+    return LinearSystem(spzeros(dim, dim), spzeros(dim, dim),
+                        spzeros(dim, dim), spzeros(dim, dim),
+                        spzeros(dim, dim), spzeros(dim, dim),
+                        spzeros(dim), spzeros(dim), spzeros(dim),
+                        spzeros(dim), spzeros(dim), dim)
+end
+
+abstract type AbstractLinearSystemSolver end
+
+function solve!(::LinearSystem, ::Solver) where Solver<:AbstractLinearSystemSolver
+    @info("This is a placeholder function for solving linear systems. To solve " *
+          "linear systems, you must define a function " *
+          "solve!(system::LinearSystem, solver::$Solver)")
+end
+
+function can_solve(::LinearSystem, ::Solver) where Solver<:AbstractLinearSystemSolver
+    return (true, "OK")
+end
+
+function solve!(ls::LinearSystem, solvers::Vector{S}) where S<:AbstractLinearSystemSolver
+    for solver in solvers
+        Solver = typeof(solver)
+        cansolve, msg = can_solve(ls, solver)
+        if !cansolve
+            @info("Solver $Solver cannot solve linear system: $msg")
+            continue
+        end
+        timeit("solve linear system using solver $Solver") do
+            solve!(ls, solver)
+        end
+        return
+    end
+    error("Failed to solve linear system.")
+end
diff --git a/src/solvers_modal.jl b/src/solvers_modal.jl
index ca30cc4..7dd6037 100644
--- a/src/solvers_modal.jl
+++ b/src/solvers_modal.jl
@@ -29,6 +29,9 @@ end
 
 A helper function to calculate P = D^-1*M
 """
+# TEMPORARILY DISABLED: Vendor package Mortar2D expects old FEMBase.AbstractProblem
+# TODO: Re-enable after vendor packages are consolidated or updated
+#= 
 function calc_projection(problem::T) where
         {T<:Union{Problem{Mortar}, Problem{Mortar2D}}}
 
@@ -53,8 +56,9 @@ function calc_projection(problem::T) where
 
     return s, m, P
 end
+=#
 
-function FEMBase.eliminate_boundary_conditions!(problem::P, K, M, f) where {P}
+function eliminate_boundary_conditions!(problem::P, K, M, f) where {P}
     isempty(problem.assembly.C2) && return nothing
     C1 = sparse(problem.assembly.C1)
     C2 = sparse(problem.assembly.C2)
@@ -76,7 +80,10 @@ end
 
 Eliminate Mortar boundary condition from matrices K, M and force vector f.
 """
-function FEMBase.eliminate_boundary_conditions!(problem::T, K, M, f) where
+# TEMPORARILY DISABLED: Vendor package Mortar2D expects old FEMBase.AbstractProblem
+# TODO: Re-enable after vendor packages are consolidated or updated
+#=
+function eliminate_boundary_conditions!(problem::T, K, M, f) where
                             {T <: Union{Problem{Mortar}, Problem{Mortar2D}}}
     @info("Eliminating mesh tie constraint $(problem.name) using static condensation")
     s, m, P = calc_projection(problem)
@@ -89,8 +96,9 @@ function FEMBase.eliminate_boundary_conditions!(problem::T, K, M, f) where
     M[:,:] .= Q'*M*Q
     return nothing
 end
+=#
 
-function FEMBase.run!(solver::Solver{Modal})
+function run!(solver::Solver{Modal})
     time = solver.properties.time
     problems = get_problems(solver)
     properties = solver.properties
diff --git a/src/sparse/sparse.jl b/src/sparse/sparse.jl
new file mode 100644
index 0000000..148fc8c
--- /dev/null
+++ b/src/sparse/sparse.jl
@@ -0,0 +1,202 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+using SparseArrays
+import SparseArrays: sparse, sparsevec
+
+mutable struct SparseMatrixCOO{T<:Real}
+    I :: Vector{Int}
+    J :: Vector{Int}
+    V :: Vector{T}
+end
+
+const SparseVectorCOO = SparseMatrixCOO
+
+function SparseMatrixCOO()
+    return SparseMatrixCOO{Float64}([], [], [])
+end
+
+function SparseVectorCOO(I::Vector, V::Vector)
+    return SparseVectorCOO(I, fill!(similar(I), 1), V)
+end
+
+function convert(::Type{SparseMatrixCOO}, A::SparseMatrixCSC)
+    return SparseMatrixCOO(findnz(A)...)
+end
+
+function convert(::Type{SparseVectorCOO}, A::SparseVector)
+    return SparseVectorCOO(findnz(A)...)
+end
+
+function convert(::Type{SparseMatrixCOO}, A::Matrix)
+    idx = findall(!iszero, A)
+    I = getindex.(idx, 1)
+    J = getindex.(idx, 2)
+    V = [A[i] for i in idx]
+    return SparseMatrixCOO(I, J, V)
+end
+
+function convert(::Type{SparseMatrixCOO}, b::Vector)
+    I = findall(!iszero, b)
+    J = fill(1, size(I))
+    V = b[I]
+    return SparseMatrixCOO(I, J, V)
+end
+
+SparseArrays.sparse(A::SparseMatrixCOO) = sparse(A.I, A.J, A.V)
+SparseArrays.sparse(A::SparseMatrixCOO, n::Int, m::Int) = sparse(A.I, A.J, A.V, n, m)
+SparseArrays.sparse(A::SparseMatrixCOO, n::Int, m::Int, f::Function) = sparse(A.I, A.J, A.V, n, m, f)
+Base.Matrix(A::SparseMatrixCOO) = Matrix(sparse(A))
+Base.Matrix(A::SparseMatrixCOO, n::Int, m::Int) = Matrix(sparse(A, n, m))
+
+SparseArrays.sparsevec(b::SparseVectorCOO) = sparsevec(b.I, b.V)
+SparseArrays.sparsevec(b::SparseVectorCOO, n::Int) = sparsevec(b.I, b.V, n)
+Base.Vector(b::SparseVectorCOO) = Vector(sparsevec(b))
+Base.Vector(b::SparseVectorCOO, n::Int) = Vector(sparsevec(b, n))
+
+function add!(A::SparseMatrixCOO, I::Int, J::Int, V::Float64)
+    push!(A.I, I)
+    push!(A.J, J)
+    push!(A.V, V)
+    return nothing
+end
+
+function add!(A::SparseMatrixCOO, I::Int, V::Float64)
+    push!(A.I, I)
+    push!(A.J, 1)
+    push!(A.V, V)
+    return nothing
+end
+
+function empty!(A::SparseMatrixCOO)
+    empty!(A.I)
+    empty!(A.J)
+    empty!(A.V)
+    return nothing
+end
+
+function append!(A::SparseMatrixCOO, B::SparseMatrixCOO)
+    append!(A.I, B.I)
+    append!(A.J, B.J)
+    append!(A.V, B.V)
+    return nothing
+end
+
+function isempty(A::SparseMatrixCOO)
+    return isempty(A.I) && isempty(A.J) && isempty(A.V)
+end
+
+"""
+    add!(K, dofs1, dofs2, ke)
+
+Add local element matrix `ke` to sparse matrix `K` for indices defined by `dofs1`
+and `dofs2`. This basically does `A[dofs1, dofs2] = A[dofs1, dofs2] + data`.
+
+# Examples
+
+```julia
+S = [3, 4]
+M = [6, 7, 8]
+ke = [5 6 7; 8 9 10]
+K = SparseMatrixCOO()
+add!(K, S, M, ke)
+Matrix(A)
+
+# output
+
+4x8 Array{Float64,2}:
+ 0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
+ 0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
+ 0.0  0.0  0.0  0.0  0.0  5.0  6.0   7.0
+ 0.0  0.0  0.0  0.0  0.0  8.0  9.0  10.0
+```
+"""
+function add!(A::SparseMatrixCOO, dofs1::AbstractVector{Int}, dofs2::AbstractVector{Int}, data)
+    n, m = length(dofs1), length(dofs2)
+    @assert length(data) == n*m
+    k = 1
+    for j=1:m
+        for i=1:n
+            add!(A, dofs1[i], dofs2[j], data[k])
+            k += 1
+        end
+    end
+    return nothing
+end
+
+""" Add sparse matrix of CSC to COO. """
+function add!(A::SparseMatrixCOO, B::SparseMatrixCSC)
+    i, j, v = findnz(B)
+    C = SparseMatrixCOO(i, j, v)
+    append!(A, C)
+end
+
+""" Add new data to COO Sparse vector. """
+function add!(A::SparseMatrixCOO, dofs::Vector{Int}, data::Array{Float64}, dim::Int=1)
+    if length(dofs) != length(data)
+        @error("Dimension mismatch when adding data to sparse vector!", dofs, data)
+        error("Simulation stopped.")
+    end
+    append!(A.I, dofs)
+    append!(A.J, dim*ones(Int, length(dofs)))
+    append!(A.V, vec(data))
+end
+
+""" Add SparseVector to SparseVectorCOO. """
+function add!(a::SparseVectorCOO, b::SparseVector)
+    i, v = findnz(b)
+    c = SparseVectorCOO(i, v)
+    append!(a, c)
+    return
+end
+
+"""
+    get_nonzero_rows(A)
+
+Returns indices of all nonzero rows from a sparse matrix `A`.
+"""
+function get_nonzero_rows(A)
+    return sort(unique(A.rowval))
+end
+
+"""
+    get_nonzero_columns(A)
+
+Returns indices of all nonzero columns from a sparse matrix `A`.
+"""
+function get_nonzero_columns(A)
+    return get_nonzero_rows(copy(transpose(A)))
+end
+
+function size(A::SparseMatrixCOO)
+    isempty(A) && return (0, 0)
+    return maximum(A.I), maximum(A.J)
+end
+
+function size(A::SparseMatrixCOO, idx::Int)
+    return size(A)[idx]
+end
+
+""" Resize sparse matrix A to (higher) dimension n x m. """
+function resize_sparse(A, n, m)
+    idx = findall(!iszero, A)
+    I = getindex.(idx, 1)
+    J = getindex.(idx, 2)
+    V = [A[i] for i in idx]
+    return sparse(I, J, V, n, m)
+end
+
+""" Resize sparse vector b to (higher) dimension n. """
+function resize_sparsevec(b, n)
+    return sparsevec(b.nzind, b.nzval, n)
+end
+
+""" Approximative comparison of two matrices A and B. """
+function isapprox(A::SparseMatrixCOO, B::SparseMatrixCOO)
+    A2 = sparse(A)
+    B2 = sparse(B, size(A2)...)
+    return isapprox(A2, B2)
+end
+
+isapprox(A::SparseMatrixCOO, B) = isapprox(Matrix(A), B)
+isapprox(A, B::SparseMatrixCOO) = isapprox(A, Matrix(B))