bugfix for add_sphere!, more tests

Author: boriskaus

src/Setup_geometry.jl

@@ -420,13 +420,15 @@ function add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid;    # required input
     # Set phase number & thermal structure in the full domain
     ind = findall(((X .- cen[1]).^2 + (Y .- cen[2]).^2 + (Z .- cen[3]).^2).^0.5 .< radius)
 
+    ind_flat = flatten_index_dimensions(Phase, ind)
+
     # Compute thermal structure accordingly. See routines below for different options
     if T != nothing
-        Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T)
+        Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
     end
 
     # Set the phase. Different routines are available for that - see below.
-    Phase[ind] = compute_phase(Phase[ind], Temp[ind], X[ind], Y[ind], Z[ind], phase)
+    Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase)
 
     return nothing
 end

test/test_Chmy.jl

@@ -68,4 +68,9 @@ grid   = UniformGrid(arch;
 Temp2D_C   =  Field(backend, grid, Center())
 Phases2D_C =  Field(backend, grid, Center(), Int32)
 
+# check
 add_box!(Phases2D_C,Temp2D_C,grid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(2))
+@test extrema(Phases2D_C) == (0,2)
+
+add_sphere!(Phases2D_C,Temp2D_C,grid,  cen=(0,0,-1), radius=2.5, phase=ConstantPhase(3), T=ConstantTemp(800))
+@test extrema(Phases2D_C) == (0,3)