diff --git a/src/particle_refinement/refinement.jl b/src/particle_refinement/refinement.jl
index 0fd687b51..5ac228b27 100644
--- a/src/particle_refinement/refinement.jl
+++ b/src/particle_refinement/refinement.jl
@@ -196,6 +196,8 @@ function refine_particles!(system_parent::FluidSystem,
     (; candidates, candidates_mass, system_child) = particle_refinement
 
     if !isempty(candidates)
+        nhs = get_neighborhood_search(system_parent, semi)
+
         # Old storage
         v_parent = _wrap_v(_v_cache, system_parent, semi, callback)
         u_parent = _wrap_u(_u_cache, system_parent, semi, callback)
@@ -211,8 +213,8 @@ function refine_particles!(system_parent::FluidSystem,
         mass_index = 1
         for particle_parent in candidates
             mass_parent = candidates_mass[mass_index]
-            bear_childs!(system_child, system_parent, particle_parent, mass_parent,
-                         particle_refinement, v_parent, u_parent, v_child, u_child, semi)
+            bear_children!(system_child, system_parent, particle_parent, mass_parent, nhs,
+                           particle_refinement, v_parent, u_parent, v_child, u_child)
 
             particle_refinement.available_children -= nchilds(system_parent,
                                                               particle_refinement)
@@ -226,11 +228,10 @@ end
 #
 # Reducing the dof by using a fixed regular refinement pattern
 # (given: position and number of child particles)
-function bear_childs!(system_child, system_parent, particle_parent, mass_parent,
-                      particle_refinement, v_parent, u_parent, v_child, u_child, semi)
+function bear_children!(system_child, system_parent, particle_parent, mass_parent, nhs,
+                        particle_refinement, v_parent, u_parent, v_child, u_child)
     (; rel_position_children, available_children, mass_ratio) = particle_refinement
 
-    nhs = get_neighborhood_search(system_parent, system_parent, semi)
     parent_coords = current_coords(u_parent, system_parent, particle_parent)
 
     # Loop over all child particles of parent particle
diff --git a/src/particle_refinement/refinement_criteria.jl b/src/particle_refinement/refinement_criteria.jl
index 54ea2311c..b777543c6 100644
--- a/src/particle_refinement/refinement_criteria.jl
+++ b/src/particle_refinement/refinement_criteria.jl
@@ -55,7 +55,6 @@ struct RefinementZone{NDIMS, ELTYPE, ZO} <: RefinementCriteria{NDIMS, ELTYPE}
         return new{NDIMS, ELTYPE,
                    typeof(zone_origin_function)}(zone_origin_function, spanning_set)
     end
-
 end
 
 @inline Base.ndims(::RefinementCriteria{NDIMS}) where {NDIMS} = NDIMS
diff --git a/src/particle_refinement/refinement_pattern.jl b/src/particle_refinement/refinement_pattern.jl
index f7024bc3d..2d94c2a0d 100644
--- a/src/particle_refinement/refinement_pattern.jl
+++ b/src/particle_refinement/refinement_pattern.jl
@@ -1,4 +1,6 @@
 function mass_distribution(system, refinement_pattern)
+    # TODO:
+    #=
     if refinement_pattern.center_particle
         # solve minimisation problem
 
@@ -12,10 +14,11 @@ function mass_distribution(system, refinement_pattern)
 
         error("no mass conservation")
     else
-        lambda = 1 / nchilds(system, refinement_pattern)
+        =#
+    lambda = 1 / nchilds(system, refinement_pattern)
 
-        return fill(lambda, SVector{nchilds(system, refinement_pattern), eltype(system)})
-    end
+    return fill(lambda, SVector{nchilds(system, refinement_pattern), eltype(system)})
+    #end
 end
 
 struct CubicSplitting{ELTYPE}
diff --git a/src/particle_refinement/resize.jl b/src/particle_refinement/resize.jl
index d7db857c7..ccdafa32c 100644
--- a/src/particle_refinement/resize.jl
+++ b/src/particle_refinement/resize.jl
@@ -54,7 +54,7 @@ end
         capacity_parent = nparticles(system) - candidates_refine
         capacity_child = nparticles(system_child) + n_new_child
 
-        capacity_parent <= 0 && error("`RefinementCriteria` affects all particles")
+        capacity_parent < 0 && error("`RefinementCriteria` affects more than all particles")
 
         # Resize child system (extending)
         resize_system!(system_child, capacity_child)
diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl
index 28c07fd57..d58bdd257 100644
--- a/src/schemes/fluid/weakly_compressible_sph/system.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/system.jl
@@ -170,7 +170,8 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst
         summary_line(io, "acceleration", system.acceleration)
         summary_line(io, "source terms", system.source_terms |> typeof |> nameof)
         if !isnothing(system.particle_refinement)
-        summary_line(io, "refinement level", refinement_level(system.particle_refinement))
+            summary_line(io, "refinement level",
+                         refinement_level(system.particle_refinement))
         end
         summary_footer(io)
     end
diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl
index ac0f34300..a0be9f1b1 100644
--- a/test/general/semidiscretization.jl
+++ b/test/general/semidiscretization.jl
@@ -23,8 +23,8 @@
         semi = Semidiscretization(system1, system2, neighborhood_search=nothing)
 
         # Verification
-        @test semi.ranges_u == (1:6, 7:18)
-        @test semi.ranges_v == (1:6, 7:12)
+        @test semi.ranges_u == ([1:6], [7:18])
+        @test semi.ranges_v == ([1:6], [7:12])
 
         nhs = ((TrixiParticles.TrivialNeighborhoodSearch{3}(0.2, Base.OneTo(2)),
                 TrixiParticles.TrivialNeighborhoodSearch{3}(0.2, Base.OneTo(3))),
diff --git a/test/particle_refinement/particle_refinement.jl b/test/particle_refinement/particle_refinement.jl
new file mode 100644
index 000000000..899e54e30
--- /dev/null
+++ b/test/particle_refinement/particle_refinement.jl
@@ -0,0 +1,2 @@
+include("refinement.jl")
+include("refinement_resize.jl")
diff --git a/test/particle_refinement/refinement.jl b/test/particle_refinement/refinement.jl
new file mode 100644
index 000000000..c42ae1478
--- /dev/null
+++ b/test/particle_refinement/refinement.jl
@@ -0,0 +1,223 @@
+@trixi_testset "Create System" begin
+    particle_spacing = 0.1
+    smoothing_length = 1.2 * particle_spacing
+
+    refinement_criterion_1 = RefinementZone(edge_length_x=1.0, edge_length_y=0.5,
+                                            zone_origin=(0.0, 0.0))
+    refinement_criterion_2 = RefinementZone(edge_length_x=0.5, edge_length_y=0.75,
+                                            zone_origin=(0.5, 1.0))
+
+    ic = InitialCondition(; particle_spacing, coordinates=ones(2, 2), density=[1.0, 1.0])
+
+    @testset "Single Criteria" begin
+        particle_refinement = ParticleRefinement(refinement_criterion_1)
+
+        system_parent = WeaklyCompressibleSPHSystem(ic, ContinuityDensity(),
+                                                    nothing,
+                                                    SchoenbergCubicSplineKernel{2}(),
+                                                    smoothing_length;
+                                                    particle_refinement)
+
+        systems = TrixiParticles.create_child_systems((system_parent,))
+
+        system_child = systems[2]
+
+        particle_spacing_child = system_child.initial_condition.particle_spacing
+        factor = particle_refinement.refinement_pattern.separation_parameter
+
+        @test length(systems) == 2
+        @test nparticles(system_parent) == 2
+        @test nparticles(system_child) == 0
+        @test system_parent.particle_refinement.system_child == system_child
+        @test TrixiParticles.refinement_level(system_parent.particle_refinement) == 0
+        @test system_child.particle_refinement isa Nothing
+        @test particle_spacing_child == factor * smoothing_length
+        # TODO: Test mass distribution
+    end
+
+    @testset "Multilevel Refinement" begin
+        particle_refinement = ParticleRefinement(refinement_criterion_1,
+                                                 criteria_next_levels=[
+                                                     refinement_criterion_1,
+                                                 ])
+
+        system_parent = WeaklyCompressibleSPHSystem(ic, ContinuityDensity(),
+                                                    nothing,
+                                                    SchoenbergCubicSplineKernel{2}(),
+                                                    smoothing_length;
+                                                    particle_refinement)
+
+        systems = TrixiParticles.create_child_systems((system_parent,))
+
+        system_child_1 = systems[2]
+        system_child_2 = systems[3]
+
+        particle_spacing_child_1 = system_child_1.initial_condition.particle_spacing
+        particle_spacing_child_2 = system_child_2.initial_condition.particle_spacing
+        factor1 = particle_refinement.refinement_pattern.separation_parameter
+        factor2 = factor1^2
+
+        @test length(systems) == 3
+        @test nparticles(system_parent) == 2
+        @test nparticles(system_child_1) == 0
+        @test nparticles(system_child_2) == 0
+        @test system_parent.particle_refinement.system_child == system_child_1
+        @test system_child_1.particle_refinement.system_child == system_child_2
+        @test TrixiParticles.refinement_level(system_parent.particle_refinement) == 0
+        @test TrixiParticles.refinement_level(system_child_1.particle_refinement) == 1
+        @test system_child_2.particle_refinement isa Nothing
+        @test particle_spacing_child_1 == factor1 * smoothing_length
+        @test particle_spacing_child_2 == factor2 * smoothing_length
+        # TODO: Test mass distribution
+    end
+
+    @testset "Multiple Criteria" begin
+        particle_refinement = ParticleRefinement(refinement_criterion_1,
+                                                 refinement_criterion_2)
+
+        system_parent = WeaklyCompressibleSPHSystem(ic, ContinuityDensity(),
+                                                    nothing,
+                                                    SchoenbergCubicSplineKernel{2}(),
+                                                    smoothing_length;
+                                                    particle_refinement)
+
+        @test length(system_parent.particle_refinement.refinement_criteria) == 2
+
+        (; refinement_criteria) = system_parent.particle_refinement
+        @test refinement_criteria[1].zone_origin(0, 0, 0, 0, 0, 0, 0) == [0.0, 0.0]
+        @test refinement_criteria[2].zone_origin(0, 0, 0, 0, 0, 0, 0) == [0.5, 1.0]
+    end
+end
+
+@trixi_testset "Refinement Pattern" begin
+    struct MockSystem <: TrixiParticles.System{2}
+        smoothing_length::Any
+    end
+
+    Base.eltype(::MockSystem) = Float64
+
+    mock_system = MockSystem(1.0)
+
+    refinement_patterns = [
+        TriangularSplitting(),
+        CubicSplitting(),
+        HexagonalSplitting(),
+        TriangularSplitting(; center_particle=false),
+        CubicSplitting(; center_particle=false),
+        HexagonalSplitting(; center_particle=false),
+    ]
+
+    expected_positions = [
+        [[0.0, 0.5], [-0.4330127018922193, -0.25], [0.4330127018922193, -0.25], [0.0, 0.0]],
+        [
+            [0.35355339059327373, 0.35355339059327373],
+            [0.35355339059327373, -0.35355339059327373],
+            [-0.35355339059327373, -0.35355339059327373],
+            [-0.35355339059327373, 0.35355339059327373],
+            [0.0, 0.0],
+        ],
+        [
+            [0.5, 0.0],
+            [-0.5, 0.0],
+            [0.25, 0.4330127018922193],
+            [0.25, -0.4330127018922193],
+            [-0.25, 0.4330127018922193],
+            [-0.25, -0.4330127018922193],
+            [0.0, 0.0],
+        ],
+        [[0.0, 0.5], [-0.4330127018922193, -0.25], [0.4330127018922193, -0.25]],
+        [
+            [0.35355339059327373, 0.35355339059327373],
+            [0.35355339059327373, -0.35355339059327373],
+            [-0.35355339059327373, -0.35355339059327373],
+            [-0.35355339059327373, 0.35355339059327373],
+        ],
+        [
+            [0.5, 0.0],
+            [-0.5, 0.0],
+            [0.25, 0.4330127018922193],
+            [0.25, -0.4330127018922193],
+            [-0.25, 0.4330127018922193],
+            [-0.25, -0.4330127018922193],
+        ],
+    ]
+
+    @testset "$(refinement_patterns[i])" for i in eachindex(refinement_patterns)
+        positions = TrixiParticles.relative_position_children(mock_system,
+                                                              refinement_patterns[i])
+
+        @test expected_positions[i] ≈ positions
+    end
+end
+
+@trixi_testset "Refine Particle" begin
+    nx = 5
+    n_particles = nx^2
+    particle_parent = 12
+    particle_spacing = 0.1
+    smoothing_length = 1.2 * particle_spacing
+
+    nhs = TrixiParticles.TrivialNeighborhoodSearch{2}(2smoothing_length, 1:n_particles)
+
+    ic = RectangularShape(particle_spacing, (nx, nx), (0.0, 0.0), velocity=ones(2),
+                          mass=1.0, density=2.0)
+
+    system_child = WeaklyCompressibleSPHSystem(ic, ContinuityDensity(),
+                                               nothing,
+                                               SchoenbergCubicSplineKernel{2}(),
+                                               smoothing_length)
+
+    refinement_criterion = RefinementZone(edge_length_x=Inf, edge_length_y=Inf,
+                                          zone_origin=(0.0, 0.0))
+
+    refinement_patterns = [
+        TriangularSplitting(),
+        CubicSplitting(),
+        HexagonalSplitting(),
+        TriangularSplitting(; center_particle=false),
+        CubicSplitting(; center_particle=false),
+        HexagonalSplitting(; center_particle=false),
+    ]
+
+    @testset "$refinement_pattern" for refinement_pattern in refinement_patterns
+        particle_refinement = ParticleRefinement(refinement_criterion; refinement_pattern)
+
+        system_parent = WeaklyCompressibleSPHSystem(ic, ContinuityDensity(),
+                                                    nothing,
+                                                    SchoenbergCubicSplineKernel{2}(),
+                                                    smoothing_length;
+                                                    particle_refinement)
+
+        n_children = TrixiParticles.nchilds(system_parent, particle_refinement)
+
+        resize!(system_child.mass, n_children)
+
+        relative_positions = TrixiParticles.relative_position_children(system_parent,
+                                                                       refinement_pattern)
+        mass_ratios = TrixiParticles.mass_distribution(system_parent,
+                                                       refinement_pattern)
+
+        particle_refinement.rel_position_children = relative_positions
+        particle_refinement.available_children = n_children
+        particle_refinement.mass_ratio = mass_ratios
+
+        v_parent = vcat(ic.velocity, ic.density')
+        u_parent = copy(ic.coordinates)
+
+        v_child = Array{Float64, 2}(undef, 3, n_children)
+        u_child = Array{Float64, 2}(undef, 2, n_children)
+
+        mass_parent = system_parent.mass[particle_parent]
+
+        TrixiParticles.bear_children!(system_child, system_parent, particle_parent,
+                                      mass_parent, nhs, particle_refinement,
+                                      v_parent, u_parent, v_child, u_child)
+
+        parent_pos = u_parent[:, particle_parent]
+        for child in 1:n_children
+            @test u_child[:, child] ≈ parent_pos .+ relative_positions[child]
+            @test v_child[:, child] == [1.0, 1.0, 2.0]
+            @test system_child.mass[child] == mass_parent * mass_ratios[child]
+        end
+    end
+end
diff --git a/test/particle_refinement/refinement_resize.jl b/test/particle_refinement/refinement_resize.jl
new file mode 100644
index 000000000..0cfa01d1e
--- /dev/null
+++ b/test/particle_refinement/refinement_resize.jl
@@ -0,0 +1,155 @@
+@trixi_testset "Resize and Copy" begin
+    n_particles = 10
+    dict_ = Dict(12 => "Filled", 0 => "Empty")
+
+    @testset "$(dict_[n_particles_child]) Child System" for n_particles_child in [0, 12]
+        refinement_patterns = [
+            TriangularSplitting(),
+            CubicSplitting(),
+            HexagonalSplitting(),
+        ]
+
+        candidates = [[1, 10], [1, 3, 5, 7, 9], [2, 4, 6, 8], collect(1:n_particles)]
+
+        old_ranges_v = (UnitRange{Int64}[1:(3n_particles)],
+                        UnitRange{Int64}[(3n_particles + 1):(3n_particles + 3n_particles_child)])
+        old_ranges_u = (UnitRange{Int64}[1:(2n_particles)],
+                        UnitRange{Int64}[(2n_particles + 1):(2n_particles + 2n_particles_child)])
+
+        old_nparticles = (n_particles, n_particles_child)
+
+        ic = InitialCondition(; particle_spacing=0.1, coordinates=ones(2, n_particles),
+                              density=ones(n_particles))
+
+        refinement_criterion = RefinementZone(edge_length_x=Inf, edge_length_y=Inf,
+                                              zone_origin=(0.0, 0.0))
+
+        refinement_callback = ParticleRefinementCallback(interval=1)
+        callback = refinement_callback.affect!
+
+        @testset "Refinement Pattern $refinement_pattern" for refinement_pattern in refinement_patterns
+            @testset "Refinement Candidates $j" for j in eachindex(candidates)
+                particle_refinement = ParticleRefinement(refinement_criterion;
+                                                         refinement_pattern=refinement_pattern)
+
+                system_parent = WeaklyCompressibleSPHSystem(ic, ContinuityDensity(),
+                                                            nothing,
+                                                            SchoenbergCubicSplineKernel{2}(),
+                                                            1.0; particle_refinement)
+
+                semi = Semidiscretization(system_parent)
+
+                # Resize emtpy child system
+                resize!(semi.systems[2].mass, n_particles_child)
+
+                # Add candidates
+                system_parent.particle_refinement.candidates = candidates[j]
+
+                # Create vectors filled with the corresponding particle index
+                v_ode_parent = reshape(stack([i * ones(Int, 3) for i in 1:n_particles]),
+                                       (3 * n_particles,))
+                u_ode_parent = reshape(stack([i * ones(Int, 2) for i in 1:n_particles]),
+                                       (2 * n_particles,))
+
+                if n_particles_child > 0
+                    v_ode_child = reshape(stack([i * ones(Int, 3)
+                                                 for i in 1:n_particles_child]),
+                                          (3 * n_particles_child,))
+                    u_ode_child = reshape(stack([i * ones(Int, 2)
+                                                 for i in 1:n_particles_child]),
+                                          (2 * n_particles_child,))
+
+                    v_ode = vcat(v_ode_parent, v_ode_child)
+                    u_ode = vcat(u_ode_parent, u_ode_child)
+                else
+                    v_ode = v_ode_parent
+                    u_ode = u_ode_parent
+                end
+
+                _v_cache = copy(v_ode)
+                _u_cache = copy(u_ode)
+
+                TrixiParticles.resize_and_copy!(callback, semi, v_ode, u_ode,
+                                                _v_cache, _u_cache)
+
+                # Iterator without candidates
+                eachparticle_excluding_candidates = (setdiff(1:n_particles, candidates[j]),
+                                                     Base.OneTo(n_particles_child))
+
+                n_candidates = length(candidates[j])
+                n_children = TrixiParticles.nchilds(system_parent, particle_refinement)
+                n_total_particles = n_particles + n_particles_child - n_candidates +
+                                    n_candidates * n_children
+
+                # Ranges after resizing
+                new_ranges_v = (UnitRange{Int64}[1:(3nparticles(system_parent))],
+                                UnitRange{Int64}[(3nparticles(system_parent) + 1):(3n_total_particles)])
+                new_ranges_u = (UnitRange{Int64}[1:(2nparticles(system_parent))],
+                                UnitRange{Int64}[(2nparticles(system_parent) + 1):(2n_total_particles)])
+
+                @testset "Iterators And Ranges" begin
+                    @test callback.nparticles_cache == old_nparticles
+                    @test callback.ranges_v_cache == old_ranges_v
+                    @test callback.ranges_u_cache == old_ranges_u
+                    @test callback.eachparticle_cache == eachparticle_excluding_candidates
+
+                    @test nparticles(system_parent) == n_particles - n_candidates ==
+                          length(eachparticle_excluding_candidates[1])
+                    @test nparticles(semi.systems[2]) ==
+                          n_candidates * n_children + n_particles_child
+
+                    @test semi.ranges_v == new_ranges_v
+                    @test semi.ranges_u == new_ranges_u
+                end
+
+                @testset "Resized Integrator-Array" begin
+                    # Parent system
+                    v_parent = TrixiParticles.wrap_v(v_ode, system_parent, semi)
+
+                    @test size(v_parent, 2) == length(eachparticle_excluding_candidates[1])
+
+                    # Test `copy_values_v!`
+                    for particle in 1:nparticles(system_parent)
+                        for dim in 1:3
+                            @test v_parent[dim, particle] ==
+                                  eachparticle_excluding_candidates[1][particle]
+                        end
+                    end
+
+                    u_parent = TrixiParticles.wrap_u(u_ode, system_parent, semi)
+
+                    @test size(u_parent, 2) == length(eachparticle_excluding_candidates[1])
+
+                    # Test `copy_values_u!`
+                    for particle in 1:nparticles(system_parent)
+                        for dim in 1:ndims(system_parent)
+                            @test u_parent[dim, particle] ==
+                                  eachparticle_excluding_candidates[1][particle]
+                        end
+                    end
+
+                    # Child system
+                    v_child = TrixiParticles.wrap_v(v_ode, semi.systems[2], semi)
+                    @test size(v_child, 2) == n_candidates * n_children + n_particles_child
+
+                    # Test `copy_values_v!`
+                    for particle in 1:n_particles_child
+                        for dim in 1:3
+                            @test v_child[dim, particle] == particle
+                        end
+                    end
+
+                    u_child = TrixiParticles.wrap_u(u_ode, semi.systems[2], semi)
+                    @test size(u_child, 2) == n_candidates * n_children + n_particles_child
+
+                    # Test `copy_values_u!`
+                    for particle in 1:n_particles_child
+                        for dim in 1:2
+                            @test u_child[dim, particle] == particle
+                        end
+                    end
+                end
+            end
+        end
+    end
+end
diff --git a/test/systems/wcsph_system.jl b/test/systems/wcsph_system.jl
index 0749261bd..83bffd45b 100644
--- a/test/systems/wcsph_system.jl
+++ b/test/systems/wcsph_system.jl
@@ -193,7 +193,7 @@
                                              smoothing_length,
                                              density_diffusion=density_diffusion)
 
-        show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), [0.0, 0.0], nothing) with 2 particles"
+        show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), [0.0, 0.0], nothing, nothing) with 2 particles"
         @test repr(system) == show_compact
         show_box = """
         ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
diff --git a/test/unittest.jl b/test/unittest.jl
index 01308bd44..06878d714 100644
--- a/test/unittest.jl
+++ b/test/unittest.jl
@@ -7,4 +7,5 @@
     include("setups/setups.jl")
     include("systems/systems.jl")
     include("schemes/schemes.jl")
+    include("particle_refinement/particle_refinement.jl")
 end;