diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 0a27ee5..b1ce349 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -6,14 +6,31 @@ This is JuliaFEM -- Finite Element Package
 """
 module JuliaFEM
 
-#using Logging
-#@Logging.configure(level=DEBUG)
-#using Lexicon
+module API
+include("api.jl")
+
+# export ....
+
+end
+
+module Preprocess
+# include("abaqus_reader.jl") <-- ERROR: LoadError: LoadError: UndefVarError: Model not defined
+include("aster_reader.jl")
+end
+
+module Postprocess
+include("xdmf.jl")
+end
+
+""" JuliaFEM testing routines. """
+module Test
+include("test.jl")
+end
+
+module Core
 
 import Base: +, -, /, *, push!, convert, getindex, setindex!, length, similar, call, vec, endof
 
-#importall Base
-
 """
 A very simple debugging macro. It prints debug message if environment variable
 JULIAFEM_DEBUG is found.
@@ -27,9 +44,11 @@ macro debug(msg)
     end
     return :( println("DEBUG: ", $msg) )
 end
+
 function set_debug_on!()
     ENV["JULIAFEM_DEBUG"] = 1;
 end
+
 function set_debug_off!()
     pop!(ENV, "JULIAFEM_DEBUG");
 end
@@ -38,7 +57,7 @@ export @debug, set_debug_on!, set_debug_off!
 
 using ForwardDiff
 autodiffcache = ForwardDiffCache()
-export derivative, jacobian, hessian
+# export derivative, jacobian, hessian
 
 """ Simple linspace extension to arrays.
 
@@ -64,6 +83,7 @@ function Base.resize!(A::SparseMatrixCSC, m::Int64, n::Int64)
     A.m = m
 end
 
+
 # fields, see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
 include("fields.jl")
 #include("basis.jl")  # interpolation of discrete fields
@@ -92,19 +112,12 @@ include("assembly.jl")
 include("solvers.jl")
 include("directsolver.jl") # parallel sparse direct solver for non-linear problems
 
-### API ###
-include("api.jl")
-
 ### MORTAR STUFF ###
 include("mortar.jl")  # mortar projection
+end
 
-# PRE AND POSTPROCESS
-include("xdmf.jl")
-include("abaqus_reader.jl")
-
-include("test.jl")
+module Interfaces
+include("interfaces.jl")
+end
 
 end # module
-
-FEM = JuliaFEM
-
diff --git a/src/api.jl b/src/api.jl
index c045ac8..3a9c551 100644
--- a/src/api.jl
+++ b/src/api.jl
@@ -1,6 +1,10 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+# define these
+abstract Problem
+abstract Element
+
 abstract BoundaryCondition
 
 type NeumannBC{ S <: AbstractString} <: BoundaryCondition
diff --git a/src/aster_reader.jl b/src/aster_reader.jl
new file mode 100644
index 0000000..7a15fe2
--- /dev/null
+++ b/src/aster_reader.jl
@@ -0,0 +1,25 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+
+
+function parse(mesh::ASCIIString, ::Type{Val{:CODE_ASTER_MAIL}})
+    model = Dict{ASCIIString, Any}()
+    header = nothing
+    data = ASCIIString[]
+    for line in split(mesh, '\n')
+        length(line) != 0 || continue
+        info("line: $line")
+        if is_aster_mail_keyword(strip(line))
+            header = parse_aster_header(line)
+            empty!(data)
+            continue
+        end
+        if line == "FINSF"
+            info(data)
+            header = nothing
+            process_aster_section!(model, join(data, ""), header, Val{header[1]})
+        end
+    end
+    return model
+end
diff --git a/src/interfaces.jl b/src/interfaces.jl
index 6fccba3..7d546b6 100644
--- a/src/interfaces.jl
+++ b/src/interfaces.jl
@@ -1,24 +1,6 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-module interfaces
-
-using Logging
-@Logging.configure(level=DEBUG)
-
-#using JuliaFEM.elasticity_solver
-export solve_elasticity_interface!
-
-"""
-This is generic interface that reads data from data model, solves elasticity
-problem and updates model.
-
-Parameters
-----------
-model : to be defined
-"""
-function solve_elasticity_interface!()
-  return 0
-end
-
+function foo()
+    return "bar"
 end
diff --git a/src/test.jl b/src/test.jl
index ed62fba..b69a347 100644
--- a/src/test.jl
+++ b/src/test.jl
@@ -1,9 +1,6 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-""" JuliaFEM testing routines. """
-module Test
-
 using Base.Test
 
 abstract TestResult
@@ -159,4 +156,3 @@ end
 
 export @test, run_test, print_test_statistics
 
-end
diff --git a/src/xdmf.jl b/src/xdmf.jl
index 446175b..c117a63 100644
--- a/src/xdmf.jl
+++ b/src/xdmf.jl
@@ -155,9 +155,7 @@ function xdmf_new_field(grid, name, source, data)
 
     typ = string(typeof(data))
     datatype = "unknown"
-    @debug("typeof: $typ")
     for j in ["Int", "Float"]
-        @debug(j)
         if contains(typ, j)
             datatype = j
         end
diff --git a/test/test_assembly.jl b/test/test_assembly.jl
index 297dccb..21859d9 100644
--- a/test/test_assembly.jl
+++ b/test/test_assembly.jl
@@ -3,10 +3,9 @@
 
 module AssemblyTests
 
-using JuliaFEM
 using JuliaFEM.Test
-using JuliaFEM: Seg2, Quad4, HeatProblem, DirichletProblem, assemble
-using JuliaFEM: condensate, reconstruct!
+using JuliaFEM.Core: Seg2, Quad4, HeatProblem, DirichletProblem, assemble
+using JuliaFEM.Core: condensate, reconstruct!
 
 function test_static_condensation()
     nodes = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
diff --git a/test/test_aster_reader.jl b/test/test_aster_reader.jl
new file mode 100644
index 0000000..ec2830f
--- /dev/null
+++ b/test/test_aster_reader.jl
@@ -0,0 +1,50 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+module AsterReaderTests
+
+using JuliaFEM
+using JuliaFEM.Test
+
+using JuliaFEM: parse
+
+function test_read_mesh()
+mesh = """
+    COOR_2D
+    N1          0.0 0.0
+    N2          1.0 0.0
+    N3          1.0 1.0
+    N4          0.0 1.0
+    FINSF
+
+    QUAD4
+    E1         N1 N2 N3 N4
+    FINSF
+
+    SEG2
+    E2         N3 N4
+    FINSF
+
+    GROUP_NO NOM=NALL
+    N1 N2
+    FINSF
+
+    GROUP_MA NOM=BODY1
+    E1 E2
+    FINSF
+
+    FIN
+    """
+
+    m = parse(mesh, Val{:CODE_ASTER_MAIL})
+    @test m["nodes"]["N1"] == [0.0, 0.0]
+    @test m["elements"]["E1"] == ["QUAD4", ["N1", "N2", "N3", "N4"]]
+    @test m["elements"]["E2"] == ["SEG2", ["N3", "N4"]]
+    @test m["elsets"]["BODY1"] == ["E1", "E2"]
+    @test m["nsets"]["NALL"] == ["N1", "N2"]
+
+end
+
+test_read_mesh()
+
+end
diff --git a/test/test_basis.jl b/test/test_basis.jl
index a5fa952..0f248f0 100644
--- a/test/test_basis.jl
+++ b/test/test_basis.jl
@@ -4,11 +4,8 @@
 module BasisTests
 
 using JuliaFEM.Test
-
-using JuliaFEM
-using JuliaFEM: AbstractElement, Element
-
-import JuliaFEM: get_basis, get_dbasis
+using JuliaFEM.Core: AbstractElement, Element
+import JuliaFEM.Core: get_basis, get_dbasis
 
 abstract TestElement <: AbstractElement
 
diff --git a/test/test_directsolver.jl b/test/test_directsolver.jl
index 6245b1e..b87e2a0 100644
--- a/test/test_directsolver.jl
+++ b/test/test_directsolver.jl
@@ -4,11 +4,9 @@
 module DirectSolverTests
 
 using JuliaFEM.Test
-using JuliaFEM
-
-using JuliaFEM: Seg2, Quad4
-using JuliaFEM: PlaneStressElasticityProblem, DirichletProblem
-using JuliaFEM: DirectSolver
+using JuliaFEM.Core: Seg2, Quad4
+using JuliaFEM.Core: PlaneStressElasticityProblem, DirichletProblem
+using JuliaFEM.Core: DirectSolver
 
 function test_solver_multiple_dirichlet_bc()
 
diff --git a/test/test_dirichlet.jl b/test/test_dirichlet.jl
index faba641..b5a881a 100644
--- a/test/test_dirichlet.jl
+++ b/test/test_dirichlet.jl
@@ -4,8 +4,7 @@
 module TestDirichletBoundaryCondition
 
 using JuliaFEM.Test
-using JuliaFEM
-using JuliaFEM: Seg2, DirichletProblem, Assembly, assemble
+using JuliaFEM.Core: Seg2, DirichletProblem, Assembly, assemble
 
 function test_dirichlet_problem_1_dim()
     element = Seg2([1, 2])
diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl
index 4a6410c..5c4427d 100644
--- a/test/test_elasticity.jl
+++ b/test/test_elasticity.jl
@@ -4,7 +4,7 @@
 module ElasticityTests
 
 using JuliaFEM.Test
-using JuliaFEM: Seg2, Quad4, PlaneStressElasticityProblem, solve!
+using JuliaFEM.Core: Seg2, Quad4, PlaneStressElasticityProblem, solve!
 
 function test_elasticity_volume_load()
     element = Quad4([1, 2, 3, 4])
diff --git a/test/test_elements.jl b/test/test_elements.jl
index a3bbb5b..2c8d2f9 100644
--- a/test/test_elements.jl
+++ b/test/test_elements.jl
@@ -5,8 +5,8 @@ module ElementTests
 
 using JuliaFEM.Test
 
-using JuliaFEM: AbstractElement, Element, Field, FieldSet, test_element
-import JuliaFEM: get_basis, get_dbasis
+using JuliaFEM.Core: AbstractElement, Element, Field, FieldSet, test_element
+import JuliaFEM.Core: get_basis, get_dbasis
 import Base: size
 
 """ Prototype element
diff --git a/test/test_fields.jl b/test/test_fields.jl
index f9cca2b..e0b93ee 100644
--- a/test/test_fields.jl
+++ b/test/test_fields.jl
@@ -4,10 +4,9 @@
 
 module FieldTests
 
-using JuliaFEM
 using JuliaFEM.Test
 
-using JuliaFEM: Field
+using JuliaFEM.Core: Field
 
 function test_create_field()
     f = Field(1.0)
diff --git a/test/test_global_assembly.jl b/test/test_global_assembly.jl
index c8f5349..9518a43 100644
--- a/test/test_global_assembly.jl
+++ b/test/test_global_assembly.jl
@@ -4,8 +4,8 @@
 module AssemblyTests
 
 using JuliaFEM.Test
-using JuliaFEM: Quad4, Seg2, FieldSet, Field, HeatProblem
-using JuliaFEM: Assembly, assemble!
+using JuliaFEM.Core: Quad4, Seg2, FieldSet, Field, HeatProblem
+using JuliaFEM.Core: Assembly, assemble!
 
 """assemble a simple two element problem and solve"""
 function test_assembly()
diff --git a/test/test_heat.jl b/test/test_heat.jl
index e9a1cab..57171f4 100644
--- a/test/test_heat.jl
+++ b/test/test_heat.jl
@@ -7,7 +7,7 @@ module HeatTests  # always wrap tests to module ending with "Tests"
 
 using JuliaFEM.Test  # always use JuliaFEM.Test, not Base.Test
 
-using JuliaFEM: Seg2, Quad4, HeatProblem, assemble
+using JuliaFEM.Core: Seg2, Quad4, HeatProblem, assemble
 
 function test_one_element()  # always start test function with name test_
 
diff --git a/test/test_interfaces.jl b/test/test_interfaces.jl
index ddb5133..d5ce5c7 100644
--- a/test/test_interfaces.jl
+++ b/test/test_interfaces.jl
@@ -1,7 +1,12 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-using FactCheck
-using Logging
-@Logging.configure(level=DEBUG)
+module InterfaceTests
 
+using JuliaFEM.Test
+
+function test_foo()
+    @test 1+2 == 3
+end
+
+end
diff --git a/test/test_math.jl b/test/test_math.jl
deleted file mode 100644
index 1d97002..0000000
--- a/test/test_math.jl
+++ /dev/null
@@ -1,41 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-module MathTests
-using JuliaFEM.Test
-
-function test_math()
-    @test 1+1 == 2
-end
-
-#using JuliaFEM: interpolate
-#=
-facts("test interpolation of different field variables") do
-    N(xi) = [
-      (1-xi[1])*(1-xi[2])/4
-      (1+xi[1])*(1-xi[2])/4
-      (1+xi[1])*(1+xi[2])/4
-      (1-xi[1])*(1+xi[2])/4
-    ]
-    dNdξ(ξ) = [-(1-ξ[2])/4.0    -(1-ξ[1])/4.0
-                (1-ξ[2])/4.0    -(1+ξ[1])/4.0
-                (1+ξ[2])/4.0     (1+ξ[1])/4.0
-               -(1+ξ[2])/4.0     (1-ξ[1])/4.0]
-    F1 = [36.0, 36.0, 36.0, 36.0]
-    F2 = [36.0 36.0 36.0 36.0]
-    F3 = F2'
-    F4 = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'
-    F5 = F4'
-    F6 = [36, 36, 36, 36]
-
-    @fact interpolate(F1, N, [0.0, 0.0]) --> 36.0
-    @fact interpolate(F2, N, [0.0, 0.0]) --> 36.0
-    @fact interpolate(F3, N, [0.0, 0.0]) --> 36.0
-    @fact interpolate(F4, N, [0.0, 0.0]) --> [5.0; 0.5]
-    @fact interpolate(F5, N, [0.0, 0.0]) --> [5.0; 0.5]
-    @fact interpolate(F5, dNdξ, [0.0, 0.0]) --> [5.0 0.0; 0.0 0.5]
-    @fact interpolate(F6, N, [0.0, 0.0]) --> 36
-end
-=#
-
-end
diff --git a/test/test_mortar.jl b/test/test_mortar.jl
index 0075fd0..7089d62 100644
--- a/test/test_mortar.jl
+++ b/test/test_mortar.jl
@@ -3,12 +3,11 @@
 
 module MortarTests
 
-using JuliaFEM
 using JuliaFEM.Test
 
-using JuliaFEM: Element, Seg2, Quad4, MortarProblem, Assembly, assemble!
-using JuliaFEM: PlaneStressElasticityProblem, DirichletProblem, DirectSolver
-using JuliaFEM: project_from_slave_to_master, project_from_master_to_slave
+using JuliaFEM.Core: Element, Seg2, Quad4, MortarProblem, Assembly, assemble!
+using JuliaFEM.Core: PlaneStressElasticityProblem, DirichletProblem, DirectSolver
+using JuliaFEM.Core: project_from_slave_to_master, project_from_master_to_slave
 
 function get_test_2d_model()
     # this is hand calculated and given as an example in my thesis
diff --git a/test/test_potential_energy.jl b/test/test_potential_energy.jl
index 8144ff0..0dbbe37 100644
--- a/test/test_potential_energy.jl
+++ b/test/test_potential_energy.jl
@@ -3,12 +3,13 @@
 
 module ElementTests
 
-using JuliaFEM
 using JuliaFEM.Test
 
-using JuliaFEM: AbstractProblem, Problem
-using JuliaFEM: Element, Seg2, Quad4
-using JuliaFEM: IntegrationPoint, solve!
+using JuliaFEM.Core: AbstractProblem, Problem
+using JuliaFEM.Core: Element, Seg2, Quad4
+using JuliaFEM.Core: IntegrationPoint, solve!
+
+import JuliaFEM.Core: get_unknown_field_name, get_unknown_field_type, get_potential_energy
 
 abstract HeatProblem <: AbstractProblem
 
@@ -16,16 +17,16 @@ function HeatProblem(dim::Int=1, elements=[])
     return Problem{HeatProblem}(dim, elements)
 end
 
-function JuliaFEM.get_unknown_field_name{P<:HeatProblem}(::Type{P})
+function get_unknown_field_name{P<:HeatProblem}(::Type{P})
     return "temperature"
 end
 
-function JuliaFEM.get_unknown_field_type{P<:HeatProblem}(::Type{P})
+function get_unknown_field_type{P<:HeatProblem}(::Type{P})
     return Float64
 end
 
 """ Calculate a potential Π = Wint - Wext of system. """
-function JuliaFEM.get_potential_energy(problem::Problem{HeatProblem}, element::Element{Quad4}, ip::IntegrationPoint, time::Number; variation=nothing)
+function get_potential_energy(problem::Problem{HeatProblem}, element::Element{Quad4}, ip::IntegrationPoint, time::Number; variation=nothing)
     k = element("temperature thermal conductivity", ip, time)
     f = element("temperature load", ip, time)
     T = element("temperature", ip, time, variation)
@@ -36,7 +37,7 @@ function JuliaFEM.get_potential_energy(problem::Problem{HeatProblem}, element::E
     return Wint - Wext
 end
 
-function JuliaFEM.get_potential_energy(problem::Problem{HeatProblem}, element::Element{Seg2}, ip::IntegrationPoint, time::Number; variation=nothing)
+function get_potential_energy(problem::Problem{HeatProblem}, element::Element{Seg2}, ip::IntegrationPoint, time::Number; variation=nothing)
     T = element("temperature", ip, time, variation)[1]
     T_ext = element("temperature external", ip, time)[1]
     coeff = element("temperature coefficient", ip, time)[1]
diff --git a/test/test_solver.jl b/test/test_solver.jl
index f19a68e..bbd66f2 100644
--- a/test/test_solver.jl
+++ b/test/test_solver.jl
@@ -4,11 +4,10 @@
 module SolverTests
 
 using JuliaFEM.Test
-using JuliaFEM
 
-using JuliaFEM: Seg2, Quad4
-using JuliaFEM: DirichletProblem, HeatProblem
-using JuliaFEM: LinearSolver
+using JuliaFEM.Core: Seg2, Quad4
+using JuliaFEM.Core: DirichletProblem, HeatProblem
+using JuliaFEM.Core: LinearSolver
 
 function test_linearsolver()
     el1 = Quad4([1, 2, 3, 4])
diff --git a/test/test_symbolic.jl b/test/test_symbolic.jl
index dd8d752..74aafb8 100644
--- a/test/test_symbolic.jl
+++ b/test/test_symbolic.jl
@@ -5,10 +5,9 @@
 
 module SymbolicFieldTests
 
-using JuliaFEM
-using JuliaFEM: Basis, Field, FieldSet, Expression, diff, grad
-
 using JuliaFEM.Test
+using JuliaFEM.Core: Basis, Field, FieldSet, Expression, diff, grad
+
 
 function get_basis()
     basis(xi) = 1/4*[
diff --git a/test/test_types.jl b/test/test_types.jl
index 7f620c2..7651f90 100644
--- a/test/test_types.jl
+++ b/test/test_types.jl
@@ -3,10 +3,9 @@
 
 module TypesTests
 
-using JuliaFEM
 using JuliaFEM.Test
 
-using JuliaFEM: Field, FieldSet
+using JuliaFEM.Core: Field, FieldSet
 
 function test_foo()
     @test 1+1 == 2
diff --git a/test/test_virtual_work.jl b/test/test_virtual_work.jl
index e6dd270..9a65a2f 100644
--- a/test/test_virtual_work.jl
+++ b/test/test_virtual_work.jl
@@ -4,9 +4,9 @@
 module TestAutoDiffWeakForm
 
 using JuliaFEM.Test
-using JuliaFEM
 
-using JuliaFEM: Problem, AbstractProblem, CG, Element, IntegrationPoint, Quad4, solve!
+using JuliaFEM.Core: Problem, AbstractProblem, CG, Element, IntegrationPoint, Quad4, solve!
+import JuliaFEM.Core: get_unknown_field_name, get_unknown_field_type, get_residual_vector
 
 abstract PlaneStressElasticityProblem <: AbstractProblem
 
@@ -14,15 +14,15 @@ function PlaneStressElasticityProblem(dim::Int=2, elements=[])
     return Problem{PlaneStressElasticityProblem}(dim, elements)
 end
 
-function JuliaFEM.get_unknown_field_name{P<:PlaneStressElasticityProblem}(::Type{P})
+function get_unknown_field_name{P<:PlaneStressElasticityProblem}(::Type{P})
     return "displacement"
 end
 
-function JuliaFEM.get_unknown_field_type{P<:PlaneStressElasticityProblem}(::Type{P})
+function get_unknown_field_type{P<:PlaneStressElasticityProblem}(::Type{P})
     return Vector{Float64}
 end
 
-function JuliaFEM.get_residual_vector{EL<:CG}(problem::Problem{PlaneStressElasticityProblem}, element::Element{EL}, ip::IntegrationPoint, time::Number; variation=nothing)
+function get_residual_vector{EL<:CG}(problem::Problem{PlaneStressElasticityProblem}, element::Element{EL}, ip::IntegrationPoint, time::Number; variation=nothing)
 
 
     basis = element(ip, time)