diff --git a/test/runtests.jl b/test/runtests.jl
index a4c349f..2e5c321 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,7 +4,7 @@
 using Base.Test
 using TimerOutputs
 
-pkg_dir = Pkg.dir("JuliaFEM")
+pkg_dir = dirname(@__DIR__)
 maybe_test_files = readdir(joinpath(pkg_dir, "test"))
 is_test_file(fn) = startswith(fn, "test_") & endswith(fn, ".jl")
 test_files = filter(is_test_file, maybe_test_files)
diff --git a/test/test_contact_2d_finite_sliding.jl b/test/test_contact_2d_finite_sliding.jl
index 624b30b..4155a5d 100644
--- a/test/test_contact_2d_finite_sliding.jl
+++ b/test/test_contact_2d_finite_sliding.jl
@@ -7,7 +7,7 @@ using JuliaFEM.Postprocess
 using JuliaFEM.Testing
 
 @testset "2d curved block with frictionless finite sliding contact using forwarddiff" begin
-    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d_curved.med"
+    meshfile = @__DIR__() * "/testdata/block_2d_curved.med"
     mesh = aster_read_mesh(meshfile)
 
     upper = Problem(Elasticity, "upper", 2)
diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl
index 354db5c..fcfc82e 100644
--- a/test/test_elasticity_2d_linear_with_surface_load.jl
+++ b/test/test_elasticity_2d_linear_with_surface_load.jl
@@ -13,7 +13,7 @@ using JuliaFEM.Testing
 =#
 @testset "test 2d linear elasticity with surface + volume load" begin
     meshfile = "/geometry/2d_block/BLOCK_1elem.med"
-    mesh = aster_read_mesh(Pkg.dir("JuliaFEM")*meshfile)
+    mesh = aster_read_mesh(dirname(@__DIR__)*meshfile)
 
     # field problem
     block = Problem(Elasticity, "BLOCK", 2)
diff --git a/test/test_elasticity_2d_nonlinear_with_surface_load.jl b/test/test_elasticity_2d_nonlinear_with_surface_load.jl
index 5dbf3ff..4e51d1f 100644
--- a/test/test_elasticity_2d_nonlinear_with_surface_load.jl
+++ b/test/test_elasticity_2d_nonlinear_with_surface_load.jl
@@ -7,7 +7,7 @@ using JuliaFEM.Testing
 
 @testset "test 2d nonlinear elasticity with surface load" begin
     meshfile = "/geometry/2d_block/BLOCK_1elem.med"
-    mesh = aster_read_mesh(Pkg.dir("JuliaFEM")*meshfile)
+    mesh = aster_read_mesh(dirname(@__DIR__)*meshfile)
 
     # field problem
     block = Problem(Elasticity, "BLOCK", 2)
diff --git a/test/test_elasticity_3d_linear_with_surface_load.jl b/test/test_elasticity_3d_linear_with_surface_load.jl
index 0ed6154..1499746 100644
--- a/test/test_elasticity_3d_linear_with_surface_load.jl
+++ b/test/test_elasticity_3d_linear_with_surface_load.jl
@@ -49,7 +49,7 @@ using JuliaFEM.Testing
 end
 
 function solve_rod_model_elasticity(eltype)
-    fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med"
+    fn = @__DIR__() * "/testdata/rod_short.med"
     mesh = aster_read_mesh(fn, eltype)
     element_sets = join(keys(mesh.element_sets), ", ")
     info("element sets: $element_sets")
diff --git a/test/test_elasticity_3d_unit_block.jl b/test/test_elasticity_3d_unit_block.jl
index f52874d..1e66fd6 100644
--- a/test/test_elasticity_3d_unit_block.jl
+++ b/test/test_elasticity_3d_unit_block.jl
@@ -7,7 +7,7 @@ using JuliaFEM.Postprocess
 using JuliaFEM.Testing
 
 function calc_model(mesh_name; with_volume_load=false, debug_print=false)
-    meshfile = Pkg.dir("JuliaFEM")*"/test/testdata/3d_block.med"
+    meshfile = @__DIR__()*"/testdata/3d_block.med"
     mesh = aster_read_mesh(meshfile, mesh_name)
 
     block = Problem(Elasticity, "BLOCK", 3)
diff --git a/test/test_elasticity_hollow_sphere_with_surface_pressure.jl b/test/test_elasticity_hollow_sphere_with_surface_pressure.jl
index 1d023e0..2d4fdd3 100644
--- a/test/test_elasticity_hollow_sphere_with_surface_pressure.jl
+++ b/test/test_elasticity_hollow_sphere_with_surface_pressure.jl
@@ -88,7 +88,7 @@ if a=0.9, b=1.0, ν=1/3, E = 24580 and p = 7317 equation yields
 http://mms2.ensmp.fr/emms_paris/plasticite3D/exercices/eSpherePress.pdf
 =#
 function test_wedge_sphere(model, u_CA, S_CA)
-    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh_file = @__DIR__() * "/testdata/primitives.med"
     mesh = aster_read_mesh(mesh_file, model)
     body = Problem(Elasticity, "hollow sphere 1/8 model", 3)
     body.elements = create_elements(mesh, "HOLLOWSPHERE8")
diff --git a/test/test_elasticity_med_pyr5_point_load.jl b/test/test_elasticity_med_pyr5_point_load.jl
index 6b8403c..ec46e9c 100644
--- a/test/test_elasticity_med_pyr5_point_load.jl
+++ b/test/test_elasticity_med_pyr5_point_load.jl
@@ -6,7 +6,7 @@ using JuliaFEM.Preprocess
 using JuliaFEM.Testing
 
 @testset "test Pyr5 elasticity with point load" begin
-    fn = Pkg.dir("JuliaFEM") * "/geometry/3d_pyr/Pyr5.med"
+    fn = dirname(@__DIR__) * "/geometry/3d_pyr/Pyr5.med"
     mesh = aster_read_mesh(fn)
     element_sets = join(keys(mesh.element_sets), ", ")
     info("element sets: $element_sets")
diff --git a/test/test_elasticity_plane_strain.jl b/test/test_elasticity_plane_strain.jl
index 33ed666..2dfdda7 100644
--- a/test/test_elasticity_plane_strain.jl
+++ b/test/test_elasticity_plane_strain.jl
@@ -8,7 +8,7 @@ using JuliaFEM.Testing
 
 @testset "test 2d linear elasticity with surface + volume load" begin
     meshfile = "/geometry/2d_block/BLOCK_1elem.med"
-    mesh = aster_read_mesh(Pkg.dir("JuliaFEM")*meshfile)
+    mesh = aster_read_mesh(dirname(@__DIR__)*meshfile)
 
     # field problem
     block = Problem(Elasticity, "BLOCK", 2)
diff --git a/test/test_elasticplastic_3d_linear_with_surface_load.jl b/test/test_elasticplastic_3d_linear_with_surface_load.jl
index 66f77ea..7bbe24b 100644
--- a/test/test_elasticplastic_3d_linear_with_surface_load.jl
+++ b/test/test_elasticplastic_3d_linear_with_surface_load.jl
@@ -62,7 +62,7 @@ using JuliaFEM.Testing
 #end
 
 # function solve_rod_model_elasticity(eltype)
-#     fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med"
+#     fn = @__DIR__() * "/testdata/rod_short.med"
 #     mesh = aster_read_mesh(fn, eltype)
 #     element_sets = join(keys(mesh.element_sets), ", ")
 #     info("element sets: $element_sets")
diff --git a/test/test_heat.jl b/test/test_heat.jl
index 06cb091..d80e7e8 100644
--- a/test/test_heat.jl
+++ b/test/test_heat.jl
@@ -8,7 +8,7 @@ using JuliaFEM.Postprocess
 
 @testset "Tet10 + convection" begin
     # For some reason Tet10 fails, maybe because of convection.
-    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh_file = @__DIR__() * "/testdata/primitives.med"
     mesh = aster_read_mesh(mesh_file, "TETRA_TET10_1")
     prob = Problem(Heat, "tet", 1)
     face = Problem(Heat, "face 4", 1)
@@ -146,7 +146,7 @@ end
 =#
 
 @testset "compare simple 3d heat problem to code aster solution" begin
-    fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med"
+    fn = @__DIR__() * "/testdata/rod_short.med"
     mesh = aster_read_mesh(fn, "Hex8")
     element_sets = join(keys(mesh.element_sets), ", ")
     info("element sets: $element_sets")
@@ -221,7 +221,7 @@ end
 
 @testset "compare simple 3d heat problem to analytical solution" begin
     function calc_3d_heat_model(mesh_name)
-        fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med"
+        fn = @__DIR__() * "/testdata/rod_short.med"
         mesh = aster_read_mesh(fn, mesh_name)
         p1 = Problem(Heat, "rod", 1)
         p2 = Problem(Dirichlet, "left support T=100", 1, "temperature")
@@ -248,7 +248,7 @@ end
 @testset "compare simple 3d heat problem to code aster solution" begin
 
     function calc_3d_heat_model(mesh_name)
-        fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med"
+        fn = @__DIR__() * "/testdata/rod_short.med"
         mesh = aster_read_mesh(fn, mesh_name)
         element_sets = join(keys(mesh.element_sets), ", ")
         info("element sets: $element_sets")
diff --git a/test/test_heat_2.jl b/test/test_heat_2.jl
index a57c8fa..8f2cdc5 100644
--- a/test/test_heat_2.jl
+++ b/test/test_heat_2.jl
@@ -7,7 +7,7 @@ using JuliaFEM.Postprocess
 using JuliaFEM.Testing
 
 @testset "3d rod" begin
-    mesh = aster_read_mesh(Pkg.dir("JuliaFEM")*"/test/testdata/primitives.med", "CYLINDER_20_TET4")
+    mesh = aster_read_mesh(@__DIR__()*"/testdata/primitives.med", "CYLINDER_20_TET4")
     problem = Problem(Heat, "rod of length 20", 1)
     problem.elements = create_elements(mesh, "CYLINDER")
     update!(problem, "temperature thermal conductivity", 200.0)
diff --git a/test/test_heat_3.jl b/test/test_heat_3.jl
index 82ab891..bd41186 100644
--- a/test/test_heat_3.jl
+++ b/test/test_heat_3.jl
@@ -10,7 +10,7 @@ using JuliaFEM.Testing
     # from FENiCS tutorial, u(x,y) = 1 + x² + 2y² on [0x1]×[0,1]
     # and u₀(x,y) = 1 + x² + 2y², f(x,y) = -6
 
-    mesh_file = Pkg.dir("JuliaFEM")*"/test/testdata/primitives.med"
+    mesh_file = @__DIR__()*"/testdata/primitives.med"
     mesh = aster_read_mesh(mesh_file, "UNITSQUARE_6X4")
 
     field = Problem(Heat, "unit square, 6x4 triangular mesh", 1)
diff --git a/test/test_heat_3d_two_rings.jl b/test/test_heat_3d_two_rings.jl
index 0d74776..d6efd8d 100644
--- a/test/test_heat_3d_two_rings.jl
+++ b/test/test_heat_3d_two_rings.jl
@@ -13,7 +13,7 @@ inner ring and T=2.0 for outer ring, measure temperature from middle of ring.
 Results are calculated using Code Aster for comparison.
 =#
 @testset "test 3d heat, two rings, and compare to CA solution" begin
-    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    meshfile = @__DIR__() * "/testdata/primitives.med"
     mesh = aster_read_mesh(meshfile, "RINGS_UNION")
 
     rings = Problem(Heat, "RINGS", 1)
@@ -33,7 +33,7 @@ Results are calculated using Code Aster for comparison.
 
     temp_jf = rings("temperature", 0.0)
 
-    fn = Pkg.dir("JuliaFEM") * "/test/testdata/rings.rmed"
+    fn = @__DIR__() * "/testdata/rings.rmed"
     results = RMEDFile(fn)
     nodes = aster_read_nodes(results)
     temp_ca = aster_read_data(results, "TEMP")
diff --git a/test/test_modal_analysis_elasticity.jl b/test/test_modal_analysis_elasticity.jl
index c0ec77c..dbf304e 100644
--- a/test/test_modal_analysis_elasticity.jl
+++ b/test/test_modal_analysis_elasticity.jl
@@ -8,7 +8,7 @@ using JuliaFEM.Testing
 
 #= this has nothing to do here
 @testset "calculate cross-sectional properties" begin
-    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh_file = @__DIR__() * "/testdata/primitives.med"
     mesh = aster_read_mesh(mesh_file, "CYLINDER_20_TET4")
     # calculate cross-sectional properties A and Iₓ
     fixed1 = Problem(Dirichlet, "left support", 3, "displacement")
@@ -65,7 +65,7 @@ numéro    fréquence (HZ)     norme d'erreur
 =#
 
 @testset "long rod natural frequencies" begin
-    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh_file = @__DIR__() * "/testdata/primitives.med"
     mesh = aster_read_mesh(mesh_file, "CYLINDER_20_TET10")
 #   for (id, coords) in mesh.nodes
 #       mesh.nodes[id][1] *= 5.0
@@ -179,7 +179,7 @@ numéro    fréquence (HZ)     norme d'erreur
 end
 
 @testset "eigenvalues of cube (tet4)" begin
-    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    meshfile = @__DIR__() * "/testdata/primitives.med"
     mesh = aster_read_mesh(meshfile, "CUBE_TET4")
     cube = Problem(mesh, Elasticity, "CUBE", 3)
     update!(cube.elements, "youngs modulus", 10000.0)
diff --git a/test/test_modal_analysis_elasticity_2.jl b/test/test_modal_analysis_elasticity_2.jl
index cfded39..bde3308 100644
--- a/test/test_modal_analysis_elasticity_2.jl
+++ b/test/test_modal_analysis_elasticity_2.jl
@@ -7,7 +7,7 @@ using JuliaFEM.Postprocess
 using JuliaFEM.Testing
 
 @testset "eigenvalues of CYLINDER1" begin
-    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    meshfile = @__DIR__() * "/testdata/primitives.med"
     mesh = aster_read_mesh(meshfile, "CYLINDER_1_TET4")
     cylinder = Problem(mesh, Elasticity, "CYLINDER", 3)
     update!(cylinder.elements, "youngs modulus", 10000.0)
@@ -35,7 +35,7 @@ using JuliaFEM.Testing
 end
 
 @testset "eigenvalues of CYLINDER20" begin
-    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    meshfile = @__DIR__() * "/testdata/primitives.med"
     mesh = aster_read_mesh(meshfile, "CYLINDER_20_TET4")
     cylinder = Problem(mesh, Elasticity, "CYLINDER", 3)
     #update!(cylinder.elements, "youngs modulus", 10.0e6)
diff --git a/test/test_modal_analysis_zero_eigenmodes.jl b/test/test_modal_analysis_zero_eigenmodes.jl
index d987b3a..543b47d 100644
--- a/test/test_modal_analysis_zero_eigenmodes.jl
+++ b/test/test_modal_analysis_zero_eigenmodes.jl
@@ -7,7 +7,7 @@ using JuliaFEM.Postprocess
 using JuliaFEM.Testing
 
 @testset "zero eigenmode model" begin
-    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    meshfile = @__DIR__() * "/testdata/primitives.med"
     mesh = aster_read_mesh(meshfile, "TETRA_TET10_1")
     model = Problem(mesh, Elasticity, "TET", 3)
     update!(model.elements, "youngs modulus", 10.0)
diff --git a/test/test_mortar_2d_mesh_tie.jl b/test/test_mortar_2d_mesh_tie.jl
index 54884bd..de17751 100644
--- a/test/test_mortar_2d_mesh_tie.jl
+++ b/test/test_mortar_2d_mesh_tie.jl
@@ -57,7 +57,7 @@ end
 end
 
 @testset "test that interface transfers constant field without error" begin
-    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d.med"
+    meshfile = @__DIR__() * "/testdata/block_2d.med"
     mesh = aster_read_mesh(meshfile)
 
     upper = Problem(Heat, "upper", 1)
@@ -132,7 +132,7 @@ end
 =#
 
 @testset "test mesh tie with splitted block and plane stress elasticity" begin
-    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d.med"
+    meshfile = @__DIR__() * "/testdata/block_2d.med"
     mesh = aster_read_mesh(meshfile)
 
     upper = Problem(Elasticity, "upper", 2)
diff --git a/test/test_mortar_2d_mesh_tie_forwarddiff.jl b/test/test_mortar_2d_mesh_tie_forwarddiff.jl
index 49d48cf..577df07 100644
--- a/test/test_mortar_2d_mesh_tie_forwarddiff.jl
+++ b/test/test_mortar_2d_mesh_tie_forwarddiff.jl
@@ -11,7 +11,7 @@ function get_model(::Type{Val{Symbol("mesh tie with curved 2d block")}};
             dual_basis=false, use_forwarddiff=true, finite_strain=false,
             geometric_stiffness=false)
 
-    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d_curved.med"
+    meshfile = @__DIR__() * "/testdata/block_2d_curved.med"
     mesh = aster_read_mesh(meshfile)
 
     upper = Problem(Elasticity, "upper", 2)
diff --git a/test/test_mortar_3d_mesh_tie_modal.jl b/test/test_mortar_3d_mesh_tie_modal.jl
index 8e331eb..bd5b821 100644
--- a/test/test_mortar_3d_mesh_tie_modal.jl
+++ b/test/test_mortar_3d_mesh_tie_modal.jl
@@ -122,7 +122,7 @@ FIN()
 
 @testset "splitted rod with tie contact" begin
     # CYLINDER_20_1_FACE1 -- CYLINDER_20_1_FACE2 -- CYLINDER_20_2_FACE_1 -- CYLINDER_20_2_FACE_2
-    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh_file = @__DIR__() * "/testdata/primitives.med"
     mesh = aster_read_mesh(mesh_file, "CYLINDER_20_SPLITTED")
     body1 = Problem(mesh, Elasticity, "CYLINDER_20_1", 3)
     body2 = Problem(mesh, Elasticity, "CYLINDER_20_2", 3)
diff --git a/test/test_mortar_3d_mesh_tie_two_rings.jl b/test/test_mortar_3d_mesh_tie_two_rings.jl
index 06e6067..2bc81d7 100644
--- a/test/test_mortar_3d_mesh_tie_two_rings.jl
+++ b/test/test_mortar_3d_mesh_tie_two_rings.jl
@@ -16,7 +16,7 @@ surface of outer ring. We should expect constant temperature in contact surface.
 This is conforming mesh so result should match to the conforming situation.
 =#
 @testset "test that curved interface transfers constant field without error, two rings problem" begin
-    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    meshfile = @__DIR__() * "/testdata/primitives.med"
     mesh = aster_read_mesh(meshfile, "RINGS")
 
     ring1 = Problem(Heat, "RING1", 1)
@@ -44,7 +44,7 @@ This is conforming mesh so result should match to the conforming situation.
     solver = LinearSolver(ring1, ring2, bc_inner, bc_outer, interface)
     solver()
 
-    fn = Pkg.dir("JuliaFEM") * "/test/testdata/rings.rmed"
+    fn = @__DIR__() * "/testdata/rings.rmed"
     results = RMEDFile(fn)
     nodes = aster_read_nodes(results)
     temp_ca = aster_read_data(results, "TEMP")