diff --git a/examples/fluid/vortex_street.jl b/examples/fluid/vortex_street.jl new file mode 100644 index 0000000000..fc33e9b994 --- /dev/null +++ b/examples/fluid/vortex_street.jl @@ -0,0 +1,193 @@ +using TrixiParticles +using OrdinaryDiffEqLowStorageRK + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.005 + +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 4 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +# Boundary geometry and initial fluid particle positions +tank_size = (1.0, 0.5) +# tank_size = (0.6, 0.6) +initial_fluid_size = tank_size + +Re = 10000 +initial_velocity = (0.0, 0.0) +nu = 1 * 1 / Re + +tspan = (0.0, 2.0) + +fluid_density = 1000.0 +sound_speed = 50.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1, background_pressure=0.0) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(false, false, true, true), velocity=initial_velocity) + +center = tank_size ./ 2 +hollow_sphere = SphereShape(fluid_particle_spacing, 0.1, center, fluid_density, + n_layers=4, sphere_type=RoundSphere()) + +filled_sphere = SphereShape(fluid_particle_spacing, 0.1, center, fluid_density, + sphere_type=RoundSphere()) + +# hollow_sphere = RectangularShape(fluid_particle_spacing, round.(Int, (0.1, 0.3) ./ fluid_particle_spacing), center .- 0.15, +# density=fluid_density) +# filled_sphere = hollow_sphere + +# using PythonCall +# using CondaPkg + +# CondaPkg.add("svgpathtools") +# CondaPkg.add("ezdxf") +# svgpath = pyimport("svgpathtools") + +# svg_path = "M509.299 100.016C507.966 91.5871 505.915 85.7111 503.145 82.3879C498.991 77.4031 479.883 60.7871 475.521 58.2947C471.16 55.8023 455.167 49.1559 448.936 47.702C442.705 46.2481 355.471 30.463 339.893 28.8014C329.508 27.6937 287.761 23.5397 214.651 16.3394C199.727 15.7768 185.488 15.1656 171.934 14.5056C151.602 13.5157 106.318 5.19544 82.4982 0.830064C58.6781 -3.53532 3.26262 9.37214 0.279574 38.2664C-2.54326 65.6089 16.5085 89.4186 34.9345 99.2843C49.9801 107.34 58.7166 107.403 71.5338 114.275C101.912 130.564 100.849 169.8 123.353 175.939C145.856 182.078 146.637 180.9 155.986 179.279C162.22 178.199 199.267 168.32 267.129 149.644L359.355 117.398L405.801 102.863C446.849 101.008 472.412 100.061 482.489 100.022C497.605 99.9635 507.524 100.062 509.299 100.016Z" + +# path = svgpath.parse_path(svg_path) +# t_range = range(0, 1, length=50 * length(path)) +# points = [(pyconvert(Float64, p.real), -pyconvert(Float64, p.imag)) +# for p in (path.point(t) for t in t_range)] + +# # ezdxf = pyimport("ezdxf") +# # doc = ezdxf.new(dxfversion="R2010") +# # msp = doc.modelspace() +# # msp.add_polyline2d(points) +# # doc.saveas("output.dxf") + +# center = tank_size ./ 2 +# points_matrix = reinterpret(reshape, Float64, points) +# scaling = 0.3 / maximum(points_matrix, dims=2)[1] +# points_matrix .*= scaling +# points_matrix .+= (-0.15, -points_matrix[2, 1]) + +# geometry = TrixiParticles.Polygon(points_matrix) + +# # trixi2vtk(geometry) + +# point_in_geometry_algorithm = WindingNumberJacobson(; geometry, +# winding_number_factor=0.4, +# hierarchical_winding=true) + +# # Returns `InitialCondition` +# shape_sampled = ComplexShape(geometry; particle_spacing=fluid_particle_spacing, density=fluid_density, +# store_winding_number=true, +# point_in_geometry_algorithm) + +# # angle = pi / 4 +# # using StaticArrays +# # rotation_matrix = @SMatrix [cos(angle) -sin(angle) +# # sin(angle) cos(angle)] +# # shape_sampled.initial_condition.coordinates .= rotation_matrix * shape_sampled.initial_condition.coordinates +# shape_sampled.initial_condition.coordinates .+= center + +# hollow_sphere = shape_sampled.initial_condition +# filled_sphere = hollow_sphere + +# n_particles = round(Int, 0.12 / fluid_particle_spacing) +# cylinder = RectangularShape(fluid_particle_spacing, (n_particles, n_particles), (0.2 - 1.0, 0.24 - 1.0), density=fluid_density) + +fluid = setdiff(tank.fluid, filled_sphere) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 2 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +viscosity_fluid = ViscosityAdami(; nu) +# viscosity_fluid = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +viscosity_wall = ViscosityAdami(; nu) +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +fluid_density_calculator = ContinuityDensity() +fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + density_diffusion=density_diffusion, + pressure_acceleration=TrixiParticles.tensile_instability_control, + particle_shifting=TrixiParticles.ParticleShiftingSun2019(5.0), + smoothing_length, viscosity=viscosity_fluid) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation(pressure_offset=10000.0) +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length) + +# Movement function +# https://en.wikipedia.org/wiki/Triangle_wave#Harmonics +# triangle_series(t, N) = 8 / pi^2 * sum(i -> (-1)^i / (2i + 1)^2 * sin(2pi * (2i + 1) * t), 0:(N-1)) +# movement_function(x, t) = x + SVector(0.4 * triangle_series(2 * t, 5), 0.0) +# is_moving(t) = true +# boundary_movement = BoundaryMovement(movement_function, is_moving) + +boundary_movement = TrixiParticles.oscillating_movement(1.0, + SVector(0.4, 0.0), + 0.0, center; + ramp_up=0.5) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +boundary_model_cylinder = BoundaryModelDummyParticles(hollow_sphere.density, + hollow_sphere.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity_wall) + +boundary_system_cylinder = BoundarySPHSystem(hollow_sphere, boundary_model_cylinder, + movement=boundary_movement) + +# boundary_system_cylinder = TotalLagrangianSPHSystem(hollow_sphere, smoothing_kernel, smoothing_length, +# 1e5, 0.0; +# n_fixed_particles=nparticles(hollow_sphere), movement=boundary_movement, +# boundary_model=boundary_model) + +# ========================================================================================== +# ==== Simulation +min_corner = minimum(fluid.coordinates, dims=2) .- fluid_particle_spacing / 2 +max_corner = maximum(fluid.coordinates, dims=2) .+ fluid_particle_spacing / 2 +periodic_box = PeriodicBox(; min_corner, max_corner) +cell_list = FullGridCellList(; min_corner, max_corner) +neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list) + +semi = Semidiscretization(fluid_system, boundary_system_cylinder; + neighborhood_search) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=10) +saving_callback = SolutionSavingCallback(dt=0.01, prefix="") +shifting_callback = ParticleShiftingCallback() + +stepsize_callback = StepsizeCallback(cfl=1.5) + +callbacks = CallbackSet(info_callback, saving_callback, shifting_callback, + stepsize_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so because forces become extremely large when +# fluid particles are very close to boundary particles, and the time integration method +# interprets this as an instability. +# fluid_dt = 1e-3 +# sol = solve(ode, RDPK3SpFSAL49(), +# # abstol=1.0e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) +# # reltol=1.0e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) +# adaptive=false, dt=fluid_dt, +# save_everystep=false, callback=callbacks); + +time_step = 1.0 +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=time_step, # This is overwritten by the stepsize callback + save_everystep=false, callback=callbacks); + +# plane = interpolate_plane([0.0, -0.25], [1.0, 0.75], 0.0025, semi, fluid_system, sol) diff --git a/examples/fluid/vortex_street2.jl b/examples/fluid/vortex_street2.jl new file mode 100644 index 0000000000..8f89583693 --- /dev/null +++ b/examples/fluid/vortex_street2.jl @@ -0,0 +1,233 @@ +using TrixiParticles +using OrdinaryDiffEqLowStorageRK + +# ========================================================================================== +# ==== Resolution +fluid_particle_spacing = 0.002 + +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 4 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +# Boundary geometry and initial fluid particle positions +tank_size = (1.0, 0.5) +# tank_size = (0.6, 0.6) +initial_fluid_size = tank_size + +initial_velocity = (1.0, 0.0) +nu = 1e-4 + +tspan = (0.0, 2.0) + +fluid_density = 1000.0 +sound_speed = 20.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1, background_pressure=0.0) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(false, false, true, true), velocity=initial_velocity) + +center = tank_size ./ 2 +# hollow_sphere = SphereShape(fluid_particle_spacing, 0.1, center, fluid_density, +# n_layers=4, sphere_type=RoundSphere()) + +# filled_sphere = SphereShape(fluid_particle_spacing, 0.1, center, fluid_density, +# sphere_type=RoundSphere()) + +# fluid = setdiff(tank.fluid, filled_sphere) + + +# ========================================================================================== +# ==== Sample and pack foot pocket +using PythonCall +using CondaPkg + +CondaPkg.add("svgpathtools") +CondaPkg.add("ezdxf") +svgpath = pyimport("svgpathtools") + +svg_path = "M509.299 100.016C507.966 91.5871 505.915 85.7111 503.145 82.3879C498.991 77.4031 479.883 60.7871 475.521 58.2947C471.16 55.8023 455.167 49.1559 448.936 47.702C442.705 46.2481 355.471 30.463 339.893 28.8014C329.508 27.6937 287.761 23.5397 214.651 16.3394C199.727 15.7768 185.488 15.1656 171.934 14.5056C151.602 13.5157 106.318 5.19544 82.4982 0.830064C58.6781 -3.53532 3.26262 9.37214 0.279574 38.2664C-2.54326 65.6089 16.5085 89.4186 34.9345 99.2843C49.9801 107.34 58.7166 107.403 71.5338 114.275C101.912 130.564 100.849 169.8 123.353 175.939C145.856 182.078 146.637 180.9 155.986 179.279C162.22 178.199 199.267 168.32 267.129 149.644L359.355 117.398L405.801 102.863C446.849 101.008 472.412 100.061 482.489 100.022C497.605 99.9635 507.524 100.062 509.299 100.016Z" + +path = svgpath.parse_path(svg_path) +t_range = range(0, 1, length=50 * length(path)) +points = [(pyconvert(Float64, p.real), -pyconvert(Float64, p.imag)) + for p in (path.point(t) for t in t_range)] + +points_matrix = reinterpret(reshape, Float64, points) +scaling = 0.3 / maximum(points_matrix, dims=2)[1] +points_matrix .*= scaling +points_matrix .+= (-0.15, -points_matrix[2, 1]) .+ center .- (0.0, 1e-4) +# # Clamp the blade in one layer of particles by moving the foot down by a particle spacing +# points_matrix .-= (0.0, particle_spacing) +geometry = TrixiParticles.Polygon(points_matrix) + +point_in_geometry_algorithm = WindingNumberJacobson(; geometry, + winding_number_factor=0.4, + hierarchical_winding=true) + +# Returns `InitialCondition` +shape_sampled = ComplexShape(geometry; particle_spacing=fluid_particle_spacing, density=fluid_density, + point_in_geometry_algorithm) + +foot_sdf = SignedDistanceField(geometry, fluid_particle_spacing; + max_signed_distance=4 * fluid_particle_spacing, + use_for_boundary_packing=true) + +boundary_packing = sample_boundary(foot_sdf; boundary_density=fluid_density, + boundary_thickness=4 * fluid_particle_spacing) + +background_pressure = 1.0 +smoothing_length_packing = 0.8 * fluid_particle_spacing +foot_packing_system = ParticlePackingSystem(shape_sampled; smoothing_length=smoothing_length_packing, + signed_distance_field=foot_sdf, background_pressure) + +fluid_packing_system = ParticlePackingSystem(boundary_packing; smoothing_length=smoothing_length_packing, + signed_distance_field=foot_sdf, is_boundary=true, background_pressure, + boundary_compress_factor=0.8) + +min_corner = minimum(tank.fluid.coordinates, dims=2) .- fluid_particle_spacing / 2 +max_corner = maximum(tank.fluid.coordinates, dims=2) .+ fluid_particle_spacing / 2 +periodic_box = PeriodicBox(; min_corner, max_corner) +cell_list = FullGridCellList(; min_corner, max_corner) +neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list, update_strategy=ParallelUpdate()) + +semi_packing = Semidiscretization(foot_packing_system, fluid_packing_system; + neighborhood_search) + +ode_packing = semidiscretize(semi_packing, (0.0, 10.0)) + +sol_packing = solve(ode_packing, RDPK3SpFSAL35(); + save_everystep=false, + callback=CallbackSet(InfoCallback(interval=50), + # SolutionSavingCallback(interval=50, prefix="packing"), + UpdateCallback()), + dtmax=1e-2) + +hollow_sphere = InitialCondition(sol_packing, foot_packing_system, semi_packing) +fluid = setdiff(tank.fluid, hollow_sphere) + +fluid_packing_system = ParticlePackingSystem(fluid; smoothing_length=smoothing_length_packing, + signed_distance_field=nothing, background_pressure) + +boundary_packing_system = ParticlePackingSystem(hollow_sphere; smoothing_length=smoothing_length_packing, + fixed_system=true, signed_distance_field=nothing, background_pressure) + +semi_packing = Semidiscretization(fluid_packing_system, boundary_packing_system; + neighborhood_search) + +ode_packing = semidiscretize(semi_packing, (0.0, 2.0)) + +sol_packing = solve(ode_packing, RDPK3SpFSAL35(); + save_everystep=false, + callback=CallbackSet(InfoCallback(interval=50), + # SolutionSavingCallback(interval=50, prefix="packing"), + UpdateCallback()), + dtmax=1e-2) + +fluid = InitialCondition(sol_packing, fluid_packing_system, semi_packing) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 2 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +viscosity_fluid = ViscosityAdami(; nu) +# viscosity_fluid = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) +viscosity_wall = ViscosityAdami(; nu) +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +fluid_density_calculator = ContinuityDensity() +fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + density_diffusion=density_diffusion, + pressure_acceleration=TrixiParticles.tensile_instability_control, + particle_shifting=TrixiParticles.ParticleShiftingSun2019(2.0), + smoothing_length, viscosity=viscosity_fluid) + +# ========================================================================================== +# ==== Boundary +# Movement function +frequency = 1.3 # Hz +amplitude = 0.18 # m +rotation_deg = 25 # degrees +rotation_phase_offset = 0.12 # periods +translation_vector = SVector(0.0, amplitude) +rotation_angle = rotation_deg * pi / 180 + +boundary_movement = TrixiParticles.oscillating_movement(frequency, + SVector(0.0, amplitude), + rotation_angle, center; + rotation_phase_offset, ramp_up=0.5) + +boundary_density_calculator = BernoulliPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +boundary_model_cylinder = BoundaryModelDummyParticles(hollow_sphere.density, + hollow_sphere.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity_wall) + +# boundary_system_cylinder = BoundarySPHSystem(hollow_sphere, boundary_model_cylinder, +# movement=boundary_movement) + +boundary_system_cylinder = TotalLagrangianSPHSystem(hollow_sphere, smoothing_kernel, smoothing_length, + 1e9, 0.3; + n_fixed_particles=nparticles(hollow_sphere), movement=boundary_movement, + boundary_model=boundary_model_cylinder, + penalty_force=PenaltyForceGanzenmueller(alpha=0.1)) + +# boundary_system_cylinder = TotalLagrangianSPHSystem(hollow_sphere, smoothing_kernel, smoothing_length, +# 1e5, 0.0; +# n_fixed_particles=nparticles(hollow_sphere), movement=boundary_movement, +# boundary_model=boundary_model) + +# ========================================================================================== +# ==== Simulation +min_corner = minimum(tank.fluid.coordinates, dims=2) .- fluid_particle_spacing / 2 +max_corner = maximum(tank.fluid.coordinates, dims=2) .+ fluid_particle_spacing / 2 +periodic_box = PeriodicBox(; min_corner, max_corner) +cell_list = FullGridCellList(; min_corner, max_corner) +neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list) + +semi = Semidiscretization(fluid_system, boundary_system_cylinder; + neighborhood_search) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=10) +saving_callback = SolutionSavingCallback(dt=0.01, prefix="") +shifting_callback = ParticleShiftingCallback() + +stepsize_callback = StepsizeCallback(cfl=1.5) + +callbacks = CallbackSet(info_callback, saving_callback, shifting_callback, + stepsize_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so because forces become extremely large when +# fluid particles are very close to boundary particles, and the time integration method +# interprets this as an instability. +# fluid_dt = 1e-3 +# sol = solve(ode, RDPK3SpFSAL49(), +# # abstol=1.0e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) +# # reltol=1.0e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) +# adaptive=false, dt=fluid_dt, +# save_everystep=false, callback=callbacks); + +time_step = 1.0 +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=time_step, # This is overwritten by the stepsize callback + save_everystep=false, callback=callbacks); + +# plane = interpolate_plane([0.0, -0.25], [1.0, 0.75], 0.0025, semi, fluid_system, sol) diff --git a/examples/fluid/vortex_street_postprocess.jl b/examples/fluid/vortex_street_postprocess.jl new file mode 100644 index 0000000000..155d592907 --- /dev/null +++ b/examples/fluid/vortex_street_postprocess.jl @@ -0,0 +1,168 @@ +using TrixiParticles +using OrdinaryDiffEqLowStorageRK + +# ========================================================================================== +# ==== Resolution +const fluid_particle_spacing = 0.2 + +# Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model +boundary_layers = 4 +spacing_ratio = 1 + +# ========================================================================================== +# ==== Experiment Setup +# Boundary geometry and initial fluid particle positions +tank_size = (65.0, 20.0) +initial_fluid_size = tank_size + +Re = 200 +initial_velocity = (1.0, 0.0) +nu = 1 * 1 / Re + +strouhal_number = 0.198 * (1 - 19.7 / Re) +frequency = strouhal_number * initial_velocity[1] / 1 + +tspan = (0.0, 50.0) + +fluid_density = 1000.0 +sound_speed = 10initial_velocity[1] +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7) + +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, spacing_ratio=spacing_ratio, + faces=(false, false, true, true), velocity=initial_velocity) + +const radius = 2.0 +const center = SVector(5.0, 10.0) +hollow_sphere = SphereShape(fluid_particle_spacing, radius, Tuple(center), fluid_density, + n_layers=4, sphere_type=RoundSphere()) + +filled_sphere = SphereShape(fluid_particle_spacing, radius, Tuple(center), fluid_density, + sphere_type=RoundSphere()) + +# n_particles = round(Int, 0.12 / fluid_particle_spacing) +# cylinder = RectangularShape(fluid_particle_spacing, (n_particles, n_particles), (0.2 - 1.0, 0.24 - 1.0), density=fluid_density) + +fluid = setdiff(tank.fluid, filled_sphere) + +# ========================================================================================== +# ==== Fluid +smoothing_length = 2 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +viscosity = ViscosityAdami(; nu) +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +fluid_density_calculator = ContinuityDensity() +fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + density_diffusion=density_diffusion, + pressure_acceleration=TrixiParticles.tensile_instability_control, + # transport_velocity=TransportVelocityAdami(50_000.0), + smoothing_length, viscosity=viscosity) + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length) + +boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) + +boundary_model_cylinder = BoundaryModelDummyParticles(hollow_sphere.density, + hollow_sphere.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length, + viscosity=viscosity) + +boundary_system_cylinder = BoundarySPHSystem(hollow_sphere, boundary_model_cylinder) + +# ========================================================================================== +# ==== Simulation +periodic_box = PeriodicBox(min_corner=[0.0, -1.0], max_corner=[65.0, 21.0]) +cell_list = FullGridCellList(min_corner=[0.0, -1.0], max_corner=[65.0, 21.0]) +neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list) + +semi = Semidiscretization(fluid_system, boundary_system, boundary_system_cylinder; + neighborhood_search) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=10) +saving_callback = SolutionSavingCallback(dt=1.0, prefix="tvf_old") +shifting_callback = ParticleShiftingCallback() + +# ========================================================================================== +# ==== Postprocessing +circle = SphereShape(fluid_particle_spacing, (2 * radius + fluid_particle_spacing) / 2, + Tuple(center), fluid_density, n_layers=1, + sphere_type=RoundSphere()) + +# Points for pressure interpolation, located at the wall interface +const data_points = copy(circle.coordinates) + +calculate_lift_force(system, v_ode, u_ode, semi, t) = nothing +function calculate_lift_force(system::TrixiParticles.FluidSystem, v_ode, u_ode, semi, t) + force = zero(SVector{ndims(system), eltype(system)}) + + values = interpolate_points(data_points, semi, system, v_ode, u_ode; cut_off_bnd=false, + clip_negative_pressure=false) + pressure = Array(values.pressure) + + for i in axes(data_points, 2) + point = TrixiParticles.current_coords(data_points, system, i) + + # F = ∑ -p_i * A_i * n_i + force -= pressure[i] * fluid_particle_spacing .* + TrixiParticles.normalize(point - center) + end + + return 2 * force[2] / (fluid_density * 1^2 * 2 * radius) +end + +calculate_drag_force(system, v_ode, u_ode, semi, t) = nothing +function calculate_drag_force(system::TrixiParticles.FluidSystem, v_ode, u_ode, semi, t) + force = zero(SVector{ndims(system), eltype(system)}) + + values = interpolate_points(data_points, semi, system, v_ode, u_ode; cut_off_bnd=false, + clip_negative_pressure=false) + pressure = Array(values.pressure) + + for i in axes(data_points, 2) + point = TrixiParticles.current_coords(data_points, system, i) + + # F = ∑ -p_i * A_i * n_i + force -= pressure[i] * fluid_particle_spacing .* + TrixiParticles.normalize(point - center) + end + + return 2 * force[1] / (fluid_density * 1^2 * 2 * radius) +end + +pp_callback = PostprocessCallback(; dt=0.5, + f_l=calculate_lift_force, f_d=calculate_drag_force, + filename="resulting_force_pst", + write_csv=true, write_file_interval=10) + +callbacks = CallbackSet(info_callback, saving_callback, shifting_callback, + UpdateCallback(), pp_callback) + +# Use a Runge-Kutta method with automatic (error based) time step size control. +# Limiting of the maximum stepsize is necessary to prevent crashing. +# When particles are approaching a wall in a uniform way, they can be advanced +# with large time steps. Close to the wall, the stepsize has to be reduced drastically. +# Sometimes, the method fails to do so because forces become extremely large when +# fluid particles are very close to boundary particles, and the time integration method +# interprets this as an instability. +sol = solve(ode, RDPK3SpFSAL49(), + abstol=1.0e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1.0e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + save_everystep=false, callback=callbacks); + +# sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), +# dt=1.0, # This is overwritten by the stepsize callback +# save_everystep=false, callback=callbacks); + +# plane = interpolate_plane([0.0, -0.25], [1.0, 0.75], 0.0025, semi, fluid_system, sol) diff --git a/examples/fsi/fin_2d.jl b/examples/fsi/fin_2d.jl new file mode 100644 index 0000000000..06ff920ed1 --- /dev/null +++ b/examples/fsi/fin_2d.jl @@ -0,0 +1,420 @@ +using TrixiParticles +using OrdinaryDiffEqLowStorageRK +using OrdinaryDiffEqSymplecticRK + +function convert_ic(ic, T) + return InitialCondition{ndims(ic)}(ic.coordinates, ic.velocity, ic.mass, ic.density, + ic.pressure, T(ic.particle_spacing)) +end + +# ========================================================================================== +# ==== Resolution +n_particles_y = 4 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 2.0) + +fin_length = 0.6 +fin_thickness = 30e-3 +real_thickness = 1e-3 +real_modulus = 125e9 +poisson_ratio = 0.3 +flexural_rigidity = real_modulus * real_thickness^3 / (1 - poisson_ratio^2) / 12 +modulus = 12 * (1 - poisson_ratio^2) * flexural_rigidity / fin_thickness^3 + +fiber_volume_fraction = 0.6 +fiber_density = 1800.0 +epoxy_density = 1250.0 +density = fiber_volume_fraction * fiber_density + + (1 - fiber_volume_fraction) * epoxy_density + +clamp_radius = 0.05 + +tank_size = (2.0, 1.0) +center = (tank_size[2] / 2, tank_size[2] / 2) +initial_fluid_size = tank_size +initial_velocity = (1.0, 0.0) + +# The structure starts at the position of the first particle and ends +# at the position of the last particle. +particle_spacing = fin_thickness / (n_particles_y - 1) +fluid_particle_spacing = particle_spacing + +smoothing_length_structure = sqrt(2) * particle_spacing +smoothing_length_fluid = 1.5 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +file = joinpath(examples_dir(), "preprocessing", "data", "fin.dxf") +geometry = load_geometry(file) + +# trixi2vtk(geometry) + +point_in_geometry_algorithm = WindingNumberJacobson(; geometry, + winding_number_factor=0.4, + hierarchical_winding=true) + +# Returns `InitialCondition` +shape_sampled = ComplexShape(geometry; particle_spacing, density=density, + grid_offset=center, point_in_geometry_algorithm) +shape_sampled = TrixiParticles.@set shape_sampled.coordinates = Float64.(shape_sampled.coordinates) + +# Beam and clamped particles +length_clamp = round(Int, 0.15 / particle_spacing) * particle_spacing # m +n_particles_per_dimension = (round(Int, (fin_length + length_clamp) / particle_spacing) + 2,# + n_particles_clamp_x, + n_particles_y) + +# Note that the `RectangularShape` puts the first particle half a particle spacing away +# from the boundary, which is correct for fluids, but not for structures. +# We therefore need to pass `place_on_shell=true`. +beam = RectangularShape(particle_spacing, n_particles_per_dimension, + (-length_clamp, 0.0), density=density, place_on_shell=true) + +fixed_particles = setdiff(shape_sampled, beam) + +# structure = union(beam, fixed_particles) + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# Make sure that the kernel support of fluid particles at an open boundary is always +# fully sampled. +# Note: Due to the dynamics at the inlets and outlets of open boundaries, +# it is recommended to use `open_boundary_layers > boundary_layers` +open_boundary_layers = 10 + +fluid_density = 1000.0 +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, + faces=(false, false, true, true), velocity=initial_velocity) +# fluid = setdiff(tank.fluid, structure) + +open_boundary_size = (fluid_particle_spacing * open_boundary_layers, tank_size[2]) + +min_coords_inlet = (-open_boundary_layers * fluid_particle_spacing, 0.0) +inlet = RectangularTank(fluid_particle_spacing, open_boundary_size, open_boundary_size, + fluid_density, n_layers=boundary_layers, + min_coordinates=min_coords_inlet, + velocity=initial_velocity, + faces=(false, false, true, true)) + +min_coords_outlet = (tank.fluid_size[1], 0.0) +outlet = RectangularTank(fluid_particle_spacing, open_boundary_size, open_boundary_size, + fluid_density, n_layers=boundary_layers, + min_coordinates=min_coords_outlet, + velocity=initial_velocity, + faces=(false, false, true, true)) + + +NDIMS = ndims(tank.fluid) +n_buffer_particles = 20 * tank.n_particles_per_dimension[2]^(NDIMS - 1) + +# ========================================================================================== +# ==== Packing +packing = false +if packing + foot_sdf = SignedDistanceField(geometry, particle_spacing; + max_signed_distance=4 * particle_spacing, + use_for_boundary_packing=true) + + boundary_packing = sample_boundary(foot_sdf; boundary_density=density, + boundary_thickness=4 * particle_spacing) + boundary_packing = TrixiParticles.@set boundary_packing.coordinates = Float64.(boundary_packing.coordinates) + boundary_packing = setdiff(boundary_packing, beam) + + background_pressure = 1.0 + smoothing_length_packing = 0.8 * particle_spacing + foot_packing_system = ParticlePackingSystem(fixed_particles; smoothing_length=smoothing_length_packing, + signed_distance_field=foot_sdf, background_pressure) + + fluid_packing_system = ParticlePackingSystem(boundary_packing; smoothing_length=smoothing_length_packing, + signed_distance_field=foot_sdf, is_boundary=true, background_pressure, + boundary_compress_factor=0.8) + + blade_packing_system = ParticlePackingSystem(beam; smoothing_length=smoothing_length_packing, + fixed_system=true, signed_distance_field=nothing, background_pressure) + + min_corner = minimum(tank.boundary.coordinates, dims=2) .- fluid_particle_spacing / 2 + max_corner = maximum(tank.boundary.coordinates, dims=2) .+ fluid_particle_spacing / 2 + min_corner .-= center + max_corner .-= center + cell_list = FullGridCellList(; min_corner, max_corner) + neighborhood_search = GridNeighborhoodSearch{2}(; cell_list, update_strategy=ParallelUpdate()) + + semi_packing = Semidiscretization(foot_packing_system, fluid_packing_system, + blade_packing_system; neighborhood_search) + + ode_packing = semidiscretize(semi_packing, (0.0, 10.0)) + + sol_packing = solve(ode_packing, RDPK3SpFSAL35(); + abstol=1e-8, + save_everystep=false, + callback=CallbackSet(InfoCallback(interval=50), + # SolutionSavingCallback(interval=50, prefix="packing_foot"), + UpdateCallback()), + dtmax=1e-2) + + packed_foot = InitialCondition(sol_packing, foot_packing_system, semi_packing) + + # Move the fin to the center of the tank + packed_foot.coordinates .+= center + beam.coordinates .+= center + + # `union(packed_foot, beam)`, but when particles are too close together, keep the ones + # from `beam` instead of `packed_foot` to ensure that the blade doesn't have holes. + structure = union(setdiff(packed_foot, beam), beam) + fluid = setdiff(tank.fluid, structure) + + # Pack the fluid against the fin and the tank boundary + pack_window = TrixiParticles.Polygon(stack([ + [0.15, 0.42], + [0.3, 0.42], + [0.44, 0.48], + [1.12, 0.48], + [1.12, 0.52], + [0.55, 0.52], + [0.5, 0.56], + [0.24, 0.6], + [0.15, 0.6], + [0.15, 0.42] + ])) + + # Then, we extract the particles that fall inside this window + pack_fluid = intersect(fluid, pack_window) + # and those outside the window + fixed_fluid = setdiff(fluid, pack_fluid) + fixed_union = union(fixed_fluid, structure) + + fluid_packing_system = ParticlePackingSystem(pack_fluid; smoothing_length=smoothing_length_packing, + signed_distance_field=nothing, background_pressure) + + fixed_packing_system = ParticlePackingSystem(fixed_union; smoothing_length=smoothing_length_packing, + fixed_system=true, signed_distance_field=nothing, background_pressure) + + min_corner = minimum(tank.boundary.coordinates, dims=2) .- fluid_particle_spacing / 2 + max_corner = maximum(tank.boundary.coordinates, dims=2) .+ fluid_particle_spacing / 2 + cell_list = FullGridCellList(; min_corner, max_corner) + neighborhood_search = GridNeighborhoodSearch{2}(; cell_list, update_strategy=ParallelUpdate()) + + semi_packing = Semidiscretization(fluid_packing_system, fixed_packing_system; + neighborhood_search) + + ode_packing = semidiscretize(semi_packing, (0.0, 2.0)) + + sol_packing = solve(ode_packing, RDPK3SpFSAL35(); + save_everystep=false, + callback=CallbackSet(InfoCallback(interval=50), + # SolutionSavingCallback(interval=50, prefix="packing"), + UpdateCallback()), + abstol=1e-8, + dtmax=1e-2) + + fluid = InitialCondition(sol_packing, fluid_packing_system, semi_packing) + fluid = union(fluid, fixed_fluid) +else + structure = union(fixed_particles, beam) + # Move the fin to the center of the tank + structure.coordinates .+= center + + fluid = setdiff(tank.fluid, structure) +end + +n_clamped_particles = nparticles(structure) - nparticles(beam) + +# Movement function +frequency = 1.3 # Hz +amplitude = 0.18 # m +rotation_deg = 25 # degrees +rotation_phase_offset = 0.12 # periods +translation_vector = SVector(0.0, amplitude) +rotation_angle = rotation_deg * pi / 180 + +boundary_motion = OscillatingMotion2D(; frequency, + translation_vector=SVector(0.0, amplitude), + rotation_angle, rotation_center=center, + rotation_phase_offset, ramp_up_tspan=(0.0, 0.5)) + +sound_speed = 40.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1, background_pressure=0.0) + +# ========================================================================================== +# ==== Structure +boundary_density_calculator = AdamiPressureExtrapolation() +viscosity_fluid = ViscosityAdami(nu=1e-4) +viscosity_fin = ViscosityAdami(nu=1e-4) + +# For the FSI we need the hydrodynamic masses and densities in the structure boundary model +hydrodynamic_densites = fluid_density * ones(size(structure.density)) +hydrodynamic_masses = hydrodynamic_densites * particle_spacing^2 + +boundary_model_structure = BoundaryModelDummyParticles(hydrodynamic_densites, + hydrodynamic_masses, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length_fluid, + viscosity=viscosity_fin) + +# k_structure = 1.0 +# beta_structure = fluid_particle_spacing / particle_spacing +# boundary_model_structure = BoundaryModelMonaghanKajtar(k_structure, beta_structure, +# particle_spacing, +# hydrodynamic_masses) + +viscosity_structure = ArtificialViscosityMonaghan(alpha=0.2) +structure_system = TotalLagrangianSPHSystem(structure; smoothing_kernel, smoothing_length=smoothing_length_structure, + young_modulus=modulus, poisson_ratio, + clamped_particles=1:n_clamped_particles, + clamped_particles_motion=boundary_motion, + boundary_model=boundary_model_structure, + velocity_averaging=TrixiParticles.VelocityAveraging(time_constant=5e-4), + viscosity=viscosity_structure, + penalty_force=PenaltyForceGanzenmueller(alpha=0.1)) + +# ========================================================================================== +# ==== Fluid +fluid_density_calculator = ContinuityDensity() +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) + +fluid_system = WeaklyCompressibleSPHSystem(fluid; density_calculator=fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length=smoothing_length_fluid, + viscosity=viscosity_fluid, + density_diffusion, + shifting_technique=ParticleShiftingTechnique(sound_speed_factor=0.2, v_max_factor=0.0), + pressure_acceleration=tensile_instability_control, + buffer_size=n_buffer_particles) +# fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, +# sound_speed, viscosity=ViscosityAdami(; nu), +# transport_velocity=TransportVelocityAdami(10 * sound_speed^2 * fluid_density)) + +# ========================================================================================== +# ==== Open Boundaries +periodic = false +if periodic + min_corner = minimum(tank.boundary.coordinates, dims=2) .- fluid_particle_spacing / 2 + max_corner = maximum(tank.boundary.coordinates, dims=2) .+ fluid_particle_spacing / 2 + min_corner = convert.(typeof(fluid_particle_spacing), min_corner) + max_corner = convert.(typeof(fluid_particle_spacing), max_corner) + periodic_box = PeriodicBox(; min_corner, max_corner) + open_boundary_system = nothing + wall = tank.boundary +else + periodic_box = nothing + + open_boundary_model = BoundaryModelDynamicalPressureZhang() + # open_boundary_model = BoundaryModelMirroringTafuni(; mirror_method=ZerothOrderMirroring()) + reference_velocity_in = SVector(1.0, 0.0) + reference_pressure_in = 0.0 + reference_density_in = nothing + boundary_type_in = InFlow() + face_in = ([0.0, 0.0], [0.0, tank_size[2]]) + flow_direction = [1.0, 0.0] + inflow = BoundaryZone(; boundary_face=face_in, face_normal=flow_direction, + open_boundary_layers, density=fluid_density, particle_spacing, + reference_density=reference_density_in, + reference_pressure=reference_pressure_in, + reference_velocity=reference_velocity_in, + initial_condition=inlet.fluid, boundary_type=boundary_type_in) + + reference_velocity_out = SVector(1.0, 0.0) + reference_pressure_out = nothing + reference_density_out = nothing + boundary_type_out = OutFlow() + face_out = ([min_coords_outlet[1], 0.0], [min_coords_outlet[1], tank_size[2]]) + outflow = BoundaryZone(; boundary_face=face_out, face_normal=(-flow_direction), + open_boundary_layers, density=fluid_density, particle_spacing, + reference_density=reference_density_out, + reference_pressure=reference_pressure_out, + reference_velocity=reference_velocity_out, + initial_condition=outlet.fluid, boundary_type=boundary_type_out) + + open_boundary_system = OpenBoundarySystem(inflow, outflow; fluid_system, + boundary_model=open_boundary_model, + buffer_size=n_buffer_particles) + + wall = union(tank.boundary, inlet.boundary, outlet.boundary) + min_corner = minimum(wall.coordinates, dims=2) .- 5 * fluid_particle_spacing + max_corner = maximum(wall.coordinates, dims=2) .+ 5 * fluid_particle_spacing +end + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(wall.density, wall.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length_fluid) + +boundary_system = WallBoundarySystem(wall, boundary_model) + +# ========================================================================================== +# ==== Simulation +cell_list = FullGridCellList(; min_corner, max_corner) +neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list, + update_strategy=ParallelUpdate()) + +semi = Semidiscretization(fluid_system, boundary_system, open_boundary_system, structure_system; neighborhood_search, + parallelization_backend=PolyesterBackend()) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +prefix = "" +saving_callback = SolutionSavingCallback(dt=0.01; prefix) + +split_cfl = 1.5 +# SSPRK104 CFL = 2.5, 15k RHS evaluations +# CarpenterKennedy2N54 CFL = 1.6, 11k RHS evaluations +# RK4 CFL = 1.2, 12k RHS evaluations +# VerletLeapfrog CFL = 0.5, 6.75k RHS evaluations +# VelocityVerlet CFL = 0.5, 6.75k RHS evaluations +# DPRKN4 CFL = 1.7, 9k RHS evaluations + +# function tip_velocity(system::TotalLagrangianSPHSystem, data, t) +# return data.velocity[2254] +# end +# pp_tip = PostprocessCallback(; tip_velocity, interval=1, +# filename="$(prefix)_tip_velocity", write_file_interval=10_000) +split_integration = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=false), adaptive=false, + stage_coupling=true, + dt=1e-5, # This is overwritten by the stepsize callback + callback=StepsizeCallback(cfl=split_cfl), + maxiters=10^8) + +fluid_cfl = 1.2 +stepsize_callback = StepsizeCallback(cfl=fluid_cfl) + +function total_volume(system::WeaklyCompressibleSPHSystem, data, t) + return sum(data.mass ./ data.density) +end +function total_volume(system, data, t) + return nothing +end +pp_cb = PostprocessCallback(; total_volume, interval=100, + filename="$(prefix)_total_volume", write_file_interval=50) + +function plane_vtk(system, dv_ode, du_ode, v_ode, u_ode, semi, t) + return nothing +end +function plane_vtk(system::WeaklyCompressibleSPHSystem, dv_ode, du_ode, v_ode, u_ode, semi, t) + resolution = fluid_particle_spacing / 6 + pvd = TrixiParticles.paraview_collection("out/$(prefix)_plane"; append=t > 0) + interpolate_plane_2d_vtk(min_corner, max_corner, resolution, + semi, semi.systems[1], v_ode, u_ode, include_wall_velocity=true, + filename="$(prefix)_plane_$(round(Int, t * 1000))", pvd=pvd, t=t) + TrixiParticles.vtk_save(pvd) + return nothing +end +interpolate_cb = PostprocessCallback(; plane_vtk, dt=0.01, filename="plane") + +callbacks = CallbackSet(info_callback, saving_callback, + stepsize_callback, split_integration, pp_cb, interpolate_cb, + UpdateCallback(), SortingCallback(interval=10_000)) + +dt_fluid = 1.25e-4 +sol = solve(ode, + # RDPK3SpFSAL35(), + CarpenterKennedy2N54(williamson_condition=false), + dt=dt_fluid, # This is overwritten by the stepsize callback + # reltol=1e-5, abstol=1e-7, + save_everystep=false, callback=callbacks, maxiters=10^8); diff --git a/examples/fsi/fin_2d_simplified.jl b/examples/fsi/fin_2d_simplified.jl new file mode 100644 index 0000000000..2234fe1df8 --- /dev/null +++ b/examples/fsi/fin_2d_simplified.jl @@ -0,0 +1,253 @@ +using TrixiParticles +using OrdinaryDiffEqLowStorageRK +using OrdinaryDiffEqSymplecticRK + +# ========================================================================================== +# ==== Resolution +n_particles_y = 4 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 2.0) + +fin_length = 0.6 +fin_thickness = 2e-3 +real_thickness = 1e-3 +real_modulus = 125e9 +poisson_ratio = 0.3 +flexural_rigidity = real_modulus * real_thickness^3 / (1 - poisson_ratio^2) / 12 +modulus = 12 * (1 - poisson_ratio^2) * flexural_rigidity / fin_thickness^3 + +fiber_volume_fraction = 0.6 +fiber_density = 1800.0 +epoxy_density = 1250.0 +density = fiber_volume_fraction * fiber_density + + (1 - fiber_volume_fraction) * epoxy_density + +clamp_radius = 0.05 + +tank_size = (0.8, 0.2) +center = (0.05, 0.1) +initial_fluid_size = tank_size +initial_velocity = (1.0, 0.0) + +# The structure starts at the position of the first particle and ends +# at the position of the last particle. +particle_spacing = fin_thickness / (n_particles_y - 1) +fluid_particle_spacing = particle_spacing + +smoothing_length_structure = sqrt(2) * particle_spacing +smoothing_length_fluid = 1.5 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +# Beam and clamped particles +length_clamp = round(Int, 0.01 / particle_spacing) * particle_spacing # m +n_particles_per_dimension = (round(Int, (length_clamp) / particle_spacing) + 2,# + n_particles_clamp_x, + n_particles_y) +shape_sampled = RectangularShape(particle_spacing, n_particles_per_dimension, + (-length_clamp, 0.0), density=density, place_on_shell=true) +length_clamp = 0.0 +n_particles_per_dimension = (round(Int, (fin_length + length_clamp) / particle_spacing) + 2,# + n_particles_clamp_x, + n_particles_y) + +# Note that the `RectangularShape` puts the first particle half a particle spacing away +# from the boundary, which is correct for fluids, but not for structures. +# We therefore need to pass `place_on_shell=true`. +beam = RectangularShape(particle_spacing, n_particles_per_dimension, + (-length_clamp, 0.0), density=density, place_on_shell=true) + +fixed_particles = setdiff(shape_sampled, beam) + +# structure = union(beam, fixed_particles) + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# Make sure that the kernel support of fluid particles at an open boundary is always +# fully sampled. +# Note: Due to the dynamics at the inlets and outlets of open boundaries, +# it is recommended to use `open_boundary_layers > boundary_layers` +open_boundary_layers = 6 + +fluid_density = 1000.0 +tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density, + n_layers=boundary_layers, + faces=(false, false, true, true), velocity=initial_velocity) +# fluid = setdiff(tank.fluid, structure) + +open_boundary_size = (fluid_particle_spacing * open_boundary_layers, tank_size[2]) + +min_coords_inlet = (-open_boundary_layers * fluid_particle_spacing, 0.0) +inlet = RectangularTank(fluid_particle_spacing, open_boundary_size, open_boundary_size, + fluid_density, n_layers=boundary_layers, + min_coordinates=min_coords_inlet, + faces=(false, false, true, true)) + +min_coords_outlet = (tank.fluid_size[1], 0.0) +outlet = RectangularTank(fluid_particle_spacing, open_boundary_size, open_boundary_size, + fluid_density, n_layers=boundary_layers, + min_coordinates=min_coords_outlet, + faces=(false, false, true, true)) + + +NDIMS = ndims(tank.fluid) +n_buffer_particles = 10 * tank.n_particles_per_dimension[2]^(NDIMS - 1) + + +structure = union(beam, fixed_particles) +# Move the fin to the center of the tank +structure.coordinates .+= center .+ (fluid_particle_spacing / 2, fluid_particle_spacing / 2) + +fluid = setdiff(tank.fluid, structure) + +n_clamped_particles = nparticles(structure) - nparticles(beam) + +# Movement function +frequency = 1.3 # Hz +amplitude = 0.18 # m +rotation_deg = 25 # degrees +rotation_phase_offset = 0.12 # periods +translation_vector = SVector(0.0, amplitude) +rotation_angle = rotation_deg * pi / 180 + +boundary_motion = OscillatingMotion2D(; frequency, + translation_vector=SVector(0.0, amplitude), + rotation_angle, rotation_center=center, + rotation_phase_offset, ramp_up_tspan=(0.0, 0.5)) + +sound_speed = 40.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=1, background_pressure=0.0) + +# ========================================================================================== +# ==== Structure +boundary_density_calculator = AdamiPressureExtrapolation() +viscosity_fluid = ViscosityAdami(nu=1e-4) +viscosity_fin = ViscosityAdami(nu=1e-4) + +# For the FSI we need the hydrodynamic masses and densities in the structure boundary model +hydrodynamic_densites = fluid_density * ones(size(structure.density)) +hydrodynamic_masses = hydrodynamic_densites * particle_spacing^2 + +boundary_model_structure = BoundaryModelDummyParticles(hydrodynamic_densites, + hydrodynamic_masses, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length_fluid, + viscosity=viscosity_fin) + +# k_structure = 1.0 +# beta_structure = fluid_particle_spacing / particle_spacing +# boundary_model_structure = BoundaryModelMonaghanKajtar(k_structure, beta_structure, +# particle_spacing, +# hydrodynamic_masses) + +structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothing_length_structure, + modulus, poisson_ratio; + n_clamped_particles, #clamped_particles_motion=boundary_motion, + velocity_averaging=nothing, + boundary_model=boundary_model_structure, + viscosity=ArtificialViscosityMonaghan(alpha=0.1), + penalty_force=PenaltyForceGanzenmueller(alpha=0.1)) + +# ========================================================================================== +# ==== Fluid +fluid_density_calculator = ContinuityDensity() +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) +# density_diffusion = DensityDiffusionAntuono(fluid, delta=0.1) + +fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length_fluid, viscosity=viscosity_fluid, + density_diffusion=density_diffusion, + shifting_technique=ParticleShiftingTechnique(sound_speed_factor=0.2, v_max_factor=0.0), + pressure_acceleration=tensile_instability_control, + buffer_size=n_buffer_particles) +# fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, +# sound_speed, viscosity=ViscosityAdami(; nu), +# transport_velocity=TransportVelocityAdami(10 * sound_speed^2 * fluid_density)) + +min_corner = minimum(tank.boundary.coordinates, dims=2) .- fluid_particle_spacing / 2 +max_corner = maximum(tank.boundary.coordinates, dims=2) .+ fluid_particle_spacing / 2 +min_corner = convert.(typeof(fluid_particle_spacing), min_corner) +max_corner = convert.(typeof(fluid_particle_spacing), max_corner) +periodic_box = PeriodicBox(; min_corner, max_corner) +open_boundary_system = nothing +wall = tank.boundary + +# ========================================================================================== +# ==== Boundary +boundary_density_calculator = AdamiPressureExtrapolation() +boundary_model = BoundaryModelDummyParticles(wall.density, wall.mass, + state_equation=state_equation, + boundary_density_calculator, + smoothing_kernel, smoothing_length_fluid) + +boundary_system = WallBoundarySystem(wall, boundary_model) + +# ========================================================================================== +# ==== Simulation +min_corner = minimum(wall.coordinates, dims=2) .- fluid_particle_spacing / 2 +max_corner = maximum(wall.coordinates, dims=2) .+ fluid_particle_spacing / 2 +cell_list = FullGridCellList(; min_corner, max_corner) +neighborhood_search = GridNeighborhoodSearch{2}(; periodic_box, cell_list, + update_strategy=ParallelUpdate()) + +semi = Semidiscretization(fluid_system, boundary_system, open_boundary_system, structure_system; neighborhood_search, + parallelization_backend=PolyesterBackend()) +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +prefix = "simplified" +saving_callback = SolutionSavingCallback(dt=0.01, prefix=prefix) + +split_cfl = 1.6 +# SSPRK104 CFL = 2.5, 15k RHS evaluations +# CarpenterKennedy2N54 CFL = 1.6, 11k RHS evaluations +# RK4 CFL = 1.2, 12k RHS evaluations +# VerletLeapfrog CFL = 0.5, 6.75k RHS evaluations +# VelocityVerlet CFL = 0.5, 6.75k RHS evaluations +# DPRKN4 CFL = 1.7, 9k RHS evaluations + +split_integration = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=false), adaptive=false, + stage_coupling=true, + dt=1e-5, # This is overwritten by the stepsize callback + callback=StepsizeCallback(cfl=split_cfl), + maxiters=10^8) + +fluid_cfl = 1.2 +stepsize_callback = StepsizeCallback(cfl=fluid_cfl) + +function total_volume(system::WeaklyCompressibleSPHSystem, data, t) + return sum(data.mass ./ data.density) +end +function total_volume(system, data, t) + return nothing +end +pp_cb = PostprocessCallback(; total_volume, interval=100, + filename="$(prefix)_total_volume", write_file_interval=50) + +function plane_vtk(system, dv_ode, du_ode, v_ode, u_ode, semi, t) + return nothing +end +function plane_vtk(system::WeaklyCompressibleSPHSystem, dv_ode, du_ode, v_ode, u_ode, semi_, t) + resolution = fluid_particle_spacing / 6 + pvd = TrixiParticles.paraview_collection("out/$(prefix)_plane"; append=t > 0) + interpolate_plane_2d_vtk(min_corner, max_corner, resolution, + semi_, semi_.systems[1], v_ode, u_ode, include_wall_velocity=true, + filename="$(prefix)_plane_$(round(Int, t * 1000))", pvd=pvd, t=t) + TrixiParticles.vtk_save(pvd) + return nothing +end +interpolate_cb = PostprocessCallback(; plane_vtk, dt=0.01, filename="plane") + +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), + stepsize_callback, split_integration, pp_cb, interpolate_cb) + +dt_fluid = 1.25e-4 +sol = solve(ode, + # RDPK3SpFSAL35(), + CarpenterKennedy2N54(williamson_condition=false), + dt=dt_fluid, # This is overwritten by the stepsize callback + # reltol=1e-5, abstol=1e-7, + save_everystep=false, callback=callbacks, maxiters=10^8); diff --git a/examples/preprocessing/data/fin.dxf b/examples/preprocessing/data/fin.dxf new file mode 100644 index 0000000000..3f56143d5a --- /dev/null +++ b/examples/preprocessing/data/fin.dxf @@ -0,0 +1,6758 @@ + 0 +SECTION + 2 +HEADER + 9 +$ACADVER + 1 +AC1024 + 9 +$ACADMAINTVER + 70 +6 + 9 +$DWGCODEPAGE + 3 +ANSI_1252 + 9 +$LASTSAVEDBY + 1 +ezdxf + 9 +$INSBASE + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$EXTMIN + 10 +-100 + 20 +-100 + 30 +-100 + 9 +$EXTMAX + 10 +100 + 20 +100 + 30 +100 + 9 +$LIMMIN + 10 +0.0 + 20 +0.0 + 9 +$LIMMAX + 10 +420.0 + 20 +297.0 + 9 +$ORTHOMODE + 70 +0 + 9 +$REGENMODE + 70 +1 + 9 +$FILLMODE + 70 +1 + 9 +$QTEXTMODE + 70 +0 + 9 +$MIRRTEXT + 70 +1 + 9 +$LTSCALE + 40 +1.0 + 9 +$ATTMODE + 70 +1 + 9 +$TEXTSIZE + 40 +2.5 + 9 +$TRACEWID + 40 +1.0 + 9 +$TEXTSTYLE + 7 +Standard + 9 +$CLAYER + 8 +0 + 9 +$CELTYPE + 6 +ByLayer + 9 +$CECOLOR + 62 +256 + 9 +$CELTSCALE + 40 +1.0 + 9 +$DISPSILH + 70 +0 + 9 +$DIMSCALE + 40 +1.0 + 9 +$DIMASZ + 40 +2.5 + 9 +$DIMEXO + 40 +0.625 + 9 +$DIMDLI + 40 +3.75 + 9 +$DIMRND + 40 +0.0 + 9 +$DIMDLE + 40 +0.0 + 9 +$DIMEXE + 40 +1.25 + 9 +$DIMTP + 40 +0.0 + 9 +$DIMTM + 40 +0.0 + 9 +$DIMTXT + 40 +2.5 + 9 +$DIMCEN + 40 +2.5 + 9 +$DIMTSZ + 40 +0.0 + 9 +$DIMTOL + 70 +0 + 9 +$DIMLIM + 70 +0 + 9 +$DIMTIH + 70 +0 + 9 +$DIMTOH + 70 +0 + 9 +$DIMSE1 + 70 +0 + 9 +$DIMSE2 + 70 +0 + 9 +$DIMTAD + 70 +1 + 9 +$DIMZIN + 70 +8 + 9 +$DIMBLK + 1 + + 9 +$DIMASO + 70 +1 + 9 +$DIMSHO + 70 +1 + 9 +$DIMPOST + 1 + + 9 +$DIMAPOST + 1 + + 9 +$DIMALT + 70 +0 + 9 +$DIMALTD + 70 +3 + 9 +$DIMALTF + 40 +0.03937007874 + 9 +$DIMLFAC + 40 +1.0 + 9 +$DIMTOFL + 70 +1 + 9 +$DIMTVP + 40 +0.0 + 9 +$DIMTIX + 70 +0 + 9 +$DIMSOXD + 70 +0 + 9 +$DIMSAH + 70 +0 + 9 +$DIMBLK1 + 1 + + 9 +$DIMBLK2 + 1 + + 9 +$DIMSTYLE + 2 +ISO-25 + 9 +$DIMCLRD + 70 +0 + 9 +$DIMCLRE + 70 +0 + 9 +$DIMCLRT + 70 +0 + 9 +$DIMTFAC + 40 +1.0 + 9 +$DIMGAP + 40 +0.625 + 9 +$DIMJUST + 70 +0 + 9 +$DIMSD1 + 70 +0 + 9 +$DIMSD2 + 70 +0 + 9 +$DIMTOLJ + 70 +0 + 9 +$DIMTZIN + 70 +8 + 9 +$DIMALTZ + 70 +0 + 9 +$DIMALTTZ + 70 +0 + 9 +$DIMUPT + 70 +0 + 9 +$DIMDEC + 70 +2 + 9 +$DIMTDEC + 70 +2 + 9 +$DIMALTU + 70 +2 + 9 +$DIMALTTD + 70 +3 + 9 +$DIMTXSTY + 7 +Standard + 9 +$DIMAUNIT + 70 +0 + 9 +$DIMADEC + 70 +0 + 9 +$DIMALTRND + 40 +0.0 + 9 +$DIMAZIN + 70 +0 + 9 +$DIMDSEP + 70 +44 + 9 +$DIMATFIT + 70 +3 + 9 +$DIMFRAC + 70 +0 + 9 +$DIMLDRBLK + 1 + + 9 +$DIMLUNIT + 70 +2 + 9 +$DIMLWD + 70 +-2 + 9 +$DIMLWE + 70 +-2 + 9 +$DIMTMOVE + 70 +0 + 9 +$DIMFXL + 40 +1.0 + 9 +$DIMFXLON + 70 +0 + 9 +$DIMJOGANG + 40 +0.785398163397 + 9 +$DIMTFILL + 70 +0 + 9 +$DIMTFILLCLR + 70 +0 + 9 +$DIMARCSYM + 70 +0 + 9 +$DIMLTYPE + 6 + + 9 +$DIMLTEX1 + 6 + + 9 +$DIMLTEX2 + 6 + + 9 +$DIMTXTDIRECTION + 70 +0 + 9 +$LUNITS + 70 +2 + 9 +$LUPREC + 70 +4 + 9 +$SKETCHINC + 40 +1.0 + 9 +$FILLETRAD + 40 +10.0 + 9 +$AUNITS + 70 +0 + 9 +$AUPREC + 70 +2 + 9 +$MENU + 1 +. + 9 +$ELEVATION + 40 +0.0 + 9 +$PELEVATION + 40 +0.0 + 9 +$THICKNESS + 40 +0.0 + 9 +$LIMCHECK + 70 +0 + 9 +$CHAMFERA + 40 +0.0 + 9 +$CHAMFERB + 40 +0.0 + 9 +$CHAMFERC + 40 +0.0 + 9 +$CHAMFERD + 40 +0.0 + 9 +$SKPOLY + 70 +0 + 9 +$TDCREATE + 40 +2461046.8383217594 + 9 +$TDUCREATE + 40 +2458532.153996898 + 9 +$TDUPDATE + 40 +2461046.8383217594 + 9 +$TDUUPDATE + 40 +2458532.1544311 + 9 +$TDINDWG + 40 +0.0 + 9 +$TDUSRTIMER + 40 +0.0 + 9 +$USRTIMER + 70 +1 + 9 +$ANGBASE + 50 +0.0 + 9 +$ANGDIR + 70 +0 + 9 +$PDMODE + 70 +0 + 9 +$PDSIZE + 40 +0.0 + 9 +$PLINEWID + 40 +0.0 + 9 +$SPLFRAME + 70 +0 + 9 +$SPLINETYPE + 70 +6 + 9 +$SPLINESEGS + 70 +8 + 9 +$HANDSEED + 5 +DA + 9 +$SURFTAB1 + 70 +6 + 9 +$SURFTAB2 + 70 +6 + 9 +$SURFTYPE + 70 +6 + 9 +$SURFU + 70 +6 + 9 +$SURFV + 70 +6 + 9 +$UCSBASE + 2 + + 9 +$UCSNAME + 2 + + 9 +$UCSORG + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSXDIR + 10 +1.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSYDIR + 10 +0.0 + 20 +1.0 + 30 +0.0 + 9 +$UCSORTHOREF + 2 + + 9 +$UCSORTHOVIEW + 70 +0 + 9 +$UCSORGTOP + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSORGBOTTOM + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSORGLEFT + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSORGRIGHT + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSORGFRONT + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$UCSORGBACK + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSBASE + 2 + + 9 +$PUCSNAME + 2 + + 9 +$PUCSORG + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSXDIR + 10 +1.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSYDIR + 10 +0.0 + 20 +1.0 + 30 +0.0 + 9 +$PUCSORTHOREF + 2 + + 9 +$PUCSORTHOVIEW + 70 +0 + 9 +$PUCSORGTOP + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSORGBOTTOM + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSORGLEFT + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSORGRIGHT + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSORGFRONT + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PUCSORGBACK + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$USERI1 + 70 +0 + 9 +$USERI2 + 70 +0 + 9 +$USERI3 + 70 +0 + 9 +$USERI4 + 70 +0 + 9 +$USERI5 + 70 +0 + 9 +$USERR1 + 40 +0.0 + 9 +$USERR2 + 40 +0.0 + 9 +$USERR3 + 40 +0.0 + 9 +$USERR4 + 40 +0.0 + 9 +$USERR5 + 40 +0.0 + 9 +$WORLDVIEW + 70 +1 + 9 +$SHADEDGE + 70 +3 + 9 +$SHADEDIF + 70 +70 + 9 +$TILEMODE + 70 +1 + 9 +$MAXACTVP + 70 +64 + 9 +$PINSBASE + 10 +0.0 + 20 +0.0 + 30 +0.0 + 9 +$PLIMCHECK + 70 +0 + 9 +$PEXTMIN + 10 +1e+20 + 20 +1e+20 + 30 +1e+20 + 9 +$PEXTMAX + 10 +-1e+20 + 20 +-1e+20 + 30 +-1e+20 + 9 +$PLIMMIN + 10 +0.0 + 20 +0.0 + 9 +$PLIMMAX + 10 +420.0 + 20 +297.0 + 9 +$UNITMODE + 70 +0 + 9 +$VISRETAIN + 70 +1 + 9 +$PLINEGEN + 70 +0 + 9 +$PSLTSCALE + 70 +1 + 9 +$TREEDEPTH + 70 +3020 + 9 +$CMLSTYLE + 2 +Standard + 9 +$CMLJUST + 70 +0 + 9 +$CMLSCALE + 40 +20.0 + 9 +$PROXYGRAPHICS + 70 +1 + 9 +$MEASUREMENT + 70 +1 + 9 +$CELWEIGHT +370 +-1 + 9 +$ENDCAPS +280 +0 + 9 +$JOINSTYLE +280 +0 + 9 +$LWDISPLAY +290 +0 + 9 +$INSUNITS + 70 +6 + 9 +$HYPERLINKBASE + 1 + + 9 +$STYLESHEET + 1 + + 9 +$XEDIT +290 +1 + 9 +$CEPSNTYPE +380 +0 + 9 +$PSTYLEMODE +290 +1 + 9 +$FINGERPRINTGUID + 2 +BD1C87EE-EA69-11F0-BD04-A088C2115B2C + 9 +$VERSIONGUID + 2 +BD1D174A-EA69-11F0-BD04-A088C2115B2C + 9 +$EXTNAMES +290 +1 + 9 +$PSVPSCALE + 40 +0.0 + 9 +$OLESTARTUP +290 +0 + 9 +$SORTENTS +280 +127 + 9 +$INDEXCTL +280 +0 + 9 +$HIDETEXT +280 +1 + 9 +$XCLIPFRAME +280 +2 + 9 +$HALOGAP +280 +0 + 9 +$OBSCOLOR + 70 +257 + 9 +$OBSLTYPE +280 +0 + 9 +$INTERSECTIONDISPLAY +280 +0 + 9 +$INTERSECTIONCOLOR + 70 +257 + 9 +$DIMASSOC +280 +2 + 9 +$PROJECTNAME + 1 + + 9 +$CAMERADISPLAY +290 +0 + 9 +$LENSLENGTH + 40 +50.0 + 9 +$CAMERAHEIGHT + 40 +0.0 + 9 +$STEPSPERSEC + 40 +24.0 + 9 +$STEPSIZE + 40 +100.0 + 9 +$3DDWFPREC + 40 +2.0 + 9 +$PSOLWIDTH + 40 +0.005 + 9 +$PSOLHEIGHT + 40 +0.08 + 9 +$LOFTANG1 + 40 +1.570796326795 + 9 +$LOFTANG2 + 40 +1.570796326795 + 9 +$LOFTMAG1 + 40 +0.0 + 9 +$LOFTMAG2 + 40 +0.0 + 9 +$LOFTPARAM + 70 +7 + 9 +$LOFTNORMALS +280 +1 + 9 +$LATITUDE + 40 +37.795 + 9 +$LONGITUDE + 40 +-122.394 + 9 +$NORTHDIRECTION + 40 +0.0 + 9 +$TIMEZONE + 70 +-8000 + 9 +$LIGHTGLYPHDISPLAY +280 +1 + 9 +$TILEMODELIGHTSYNCH +280 +1 + 9 +$CMATERIAL +347 +20 + 9 +$SOLIDHIST +280 +0 + 9 +$SHOWHIST +280 +1 + 9 +$DWFFRAME +280 +2 + 9 +$DGNFRAME +280 +2 + 9 +$REALWORLDSCALE +290 +1 + 9 +$INTERFERECOLOR + 62 +256 + 9 +$CSHADOW +280 +0 + 9 +$SHADOWPLANELOCATION + 40 +0.0 + 0 +ENDSEC + 0 +SECTION + 2 +CLASSES + 0 +CLASS + 1 +ACDBDICTIONARYWDFLT + 2 +AcDbDictionaryWithDefault + 3 +ObjectDBX Classes + 90 +0 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +SUN + 2 +AcDbSun + 3 +SCENEOE + 90 +1153 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +VISUALSTYLE + 2 +AcDbVisualStyle + 3 +ObjectDBX Classes + 90 +4095 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +MATERIAL + 2 +AcDbMaterial + 3 +ObjectDBX Classes + 90 +1153 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +SCALE + 2 +AcDbScale + 3 +ObjectDBX Classes + 90 +1153 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +TABLESTYLE + 2 +AcDbTableStyle + 3 +ObjectDBX Classes + 90 +4095 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +MLEADERSTYLE + 2 +AcDbMLeaderStyle + 3 +ACDB_MLEADERSTYLE_CLASS + 90 +4095 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +DICTIONARYVAR + 2 +AcDbDictionaryVar + 3 +ObjectDBX Classes + 90 +0 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +CELLSTYLEMAP + 2 +AcDbCellStyleMap + 3 +ObjectDBX Classes + 90 +1152 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +MENTALRAYRENDERSETTINGS + 2 +AcDbMentalRayRenderSettings + 3 +SCENEOE + 90 +1024 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +ACDBDETAILVIEWSTYLE + 2 +AcDbDetailViewStyle + 3 +ObjectDBX Classes + 90 +1025 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +ACDBSECTIONVIEWSTYLE + 2 +AcDbSectionViewStyle + 3 +ObjectDBX Classes + 90 +1025 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +RASTERVARIABLES + 2 +AcDbRasterVariables + 3 +ISM + 90 +0 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +LAYOUT + 2 +AcDbLayout + 3 +ObjectDBX Classes + 90 +0 + 91 +0 +280 +0 +281 +0 + 0 +CLASS + 1 +ACDBPLACEHOLDER + 2 +AcDbPlaceHolder + 3 +ObjectDBX Classes + 90 +0 + 91 +0 +280 +0 +281 +0 + 0 +ENDSEC + 0 +SECTION + 2 +TABLES + 0 +TABLE + 2 +VPORT + 5 +8 +330 +0 +100 +AcDbSymbolTable + 70 +1 + 0 +VPORT + 5 +23 +330 +8 +100 +AcDbSymbolTableRecord +100 +AcDbViewportTableRecord + 2 +*Active + 70 +0 + 10 +0.0 + 20 +0.0 + 11 +1.0 + 21 +1.0 + 12 +70.0 + 22 +50.0 + 13 +0.0 + 23 +0.0 + 14 +0.5 + 24 +0.5 + 15 +0.5 + 25 +0.5 + 16 +0.0 + 26 +0.0 + 36 +1.0 + 17 +0.0 + 27 +0.0 + 37 +0.0 + 40 +1.0 + 41 +1.34 + 42 +50.0 + 43 +0.0 + 44 +0.0 + 50 +0.0 + 51 +0.0 + 71 +0 + 72 +1000 + 73 +1 + 74 +3 + 75 +0 + 76 +0 + 77 +0 + 78 +0 +281 +0 + 65 +0 +146 +0.0 + 0 +ENDTAB + 0 +TABLE + 2 +LTYPE + 5 +2 +330 +0 +100 +AcDbSymbolTable + 70 +3 + 0 +LTYPE + 5 +24 +330 +2 +100 +AcDbSymbolTableRecord +100 +AcDbLinetypeTableRecord + 2 +ByBlock + 70 +0 + 3 + + 72 +65 + 73 +0 + 40 +0.0 + 0 +LTYPE + 5 +25 +330 +2 +100 +AcDbSymbolTableRecord +100 +AcDbLinetypeTableRecord + 2 +ByLayer + 70 +0 + 3 + + 72 +65 + 73 +0 + 40 +0.0 + 0 +LTYPE + 5 +26 +330 +2 +100 +AcDbSymbolTableRecord +100 +AcDbLinetypeTableRecord + 2 +Continuous + 70 +0 + 3 + + 72 +65 + 73 +0 + 40 +0.0 + 0 +ENDTAB + 0 +TABLE + 2 +LAYER + 5 +1 +330 +0 +100 +AcDbSymbolTable + 70 +2 + 0 +LAYER + 5 +27 +330 +1 +100 +AcDbSymbolTableRecord +100 +AcDbLayerTableRecord + 2 +0 + 70 +0 + 62 +7 + 6 +Continuous +370 +-3 +390 +13 +347 +21 + 0 +LAYER + 5 +28 +330 +1 +100 +AcDbSymbolTableRecord +100 +AcDbLayerTableRecord + 2 +Defpoints + 70 +0 + 62 +7 + 6 +Continuous +290 +0 +370 +-3 +390 +13 +347 +21 + 0 +ENDTAB + 0 +TABLE + 2 +STYLE + 5 +5 +330 +0 +100 +AcDbSymbolTable + 70 +1 + 0 +STYLE + 5 +29 +330 +5 +100 +AcDbSymbolTableRecord +100 +AcDbTextStyleTableRecord + 2 +Standard + 70 +0 + 40 +0.0 + 41 +1.0 + 50 +0.0 + 71 +0 + 42 +2.5 + 3 +txt + 4 + + 0 +ENDTAB + 0 +TABLE + 2 +VIEW + 5 +7 +330 +0 +100 +AcDbSymbolTable + 70 +0 + 0 +ENDTAB + 0 +TABLE + 2 +UCS + 5 +6 +330 +0 +100 +AcDbSymbolTable + 70 +0 + 0 +ENDTAB + 0 +TABLE + 2 +APPID + 5 +3 +330 +0 +100 +AcDbSymbolTable + 70 +2 + 0 +APPID + 5 +2A +330 +3 +100 +AcDbSymbolTableRecord +100 +AcDbRegAppTableRecord + 2 +ACAD + 70 +0 + 0 +APPID + 5 +D9 +330 +3 +100 +AcDbSymbolTableRecord +100 +AcDbRegAppTableRecord + 2 +HATCHBACKGROUNDCOLOR + 70 +0 + 0 +ENDTAB + 0 +TABLE + 2 +DIMSTYLE + 5 +4 +330 +0 +100 +AcDbSymbolTable + 70 +1 +100 +AcDbDimStyleTable + 0 +DIMSTYLE +105 +2B +330 +4 +100 +AcDbSymbolTableRecord +100 +AcDbDimStyleTableRecord + 2 +Standard + 70 +0 + 40 +1.0 + 41 +2.5 + 42 +0.625 + 43 +3.75 + 44 +1.25 + 45 +0.0 + 46 +0.0 + 47 +0.0 + 48 +0.0 + 49 +2.5 +140 +2.5 +141 +2.5 +142 +0.0 +143 +0.03937007874 +144 +1.0 +145 +0.0 +146 +1.0 +147 +0.625 +148 +0.0 + 69 +0 + 70 +0 + 71 +0 + 72 +0 + 73 +0 + 74 +0 + 75 +0 + 76 +0 + 77 +1 + 78 +8 + 79 +0 +170 +0 +171 +3 +172 +1 +173 +0 +174 +0 +175 +0 +176 +0 +177 +0 +178 +0 +179 +0 +271 +0 +272 +2 +273 +2 +274 +3 +275 +0 +276 +0 +277 +2 +278 +44 +279 +0 +280 +0 +281 +0 +282 +0 +283 +0 +284 +8 +285 +0 +286 +0 +288 +0 +289 +3 +290 +0 +371 +-2 +372 +-2 + 0 +ENDTAB + 0 +TABLE + 2 +BLOCK_RECORD + 5 +9 +330 +0 +100 +AcDbSymbolTable + 70 +2 + 0 +BLOCK_RECORD + 5 +17 +330 +9 +100 +AcDbSymbolTableRecord +100 +AcDbBlockTableRecord + 2 +*Model_Space +340 +1A + 70 +0 +280 +1 +281 +0 + 0 +BLOCK_RECORD + 5 +1B +330 +9 +100 +AcDbSymbolTableRecord +100 +AcDbBlockTableRecord + 2 +*Paper_Space +340 +1E + 70 +0 +280 +1 +281 +0 + 0 +ENDTAB + 0 +ENDSEC + 0 +SECTION + 2 +BLOCKS + 0 +BLOCK + 5 +18 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbBlockBegin + 2 +*Model_Space + 70 +0 + 10 +0.0 + 20 +0.0 + 30 +0.0 + 3 +*Model_Space + 1 + + 0 +ENDBLK + 5 +19 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbBlockEnd + 0 +BLOCK + 5 +1C +330 +1B +100 +AcDbEntity + 8 +0 +100 +AcDbBlockBegin + 2 +*Paper_Space + 70 +0 + 10 +0.0 + 20 +0.0 + 30 +0.0 + 3 +*Paper_Space + 1 + + 0 +ENDBLK + 5 +1D +330 +1B +100 +AcDbEntity + 8 +0 +100 +AcDbBlockEnd + 0 +ENDSEC + 0 +SECTION + 2 +ENTITIES + 0 +POLYLINE + 5 +2D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDb2dPolyline + 66 +1 + 10 +0.0 + 20 +0.0 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +2F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +0.0 + 20 +5.551115123125783e-17 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +30 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0010647730582502013 + 20 +-2.8770736736127844e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +31 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0034971934090743395 + 20 +-4.953574980059994e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +32 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.007027977637432459 + 20 +-6.326524879191053e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +33 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.011387842328284381 + 20 +-7.092944330866491e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +34 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.01630750406658993 + 20 +-7.34985429493018e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +35 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.02151767943730898 + 20 +-7.19427573124265e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +36 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.026749085025401576 + 20 +-6.723229599658875e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +37 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.031732437415827486 + 20 +-6.0337368600282826e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +38 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.03619845319354664 + 20 +-5.222818472205848e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +39 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.03987784894351887 + 20 +-4.387495396040997e-05 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.04378446456114665 + 20 +-0.0008783586002438781 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.048103522430103185 + 20 +-0.0025049262890193824 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.051718721680916624 + 20 +-0.003948651160244676 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0549385925607315 + 20 +-0.0052322536218777915 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.05807166531669217 + 20 +-0.006378454081876428 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +3F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.061532460806013545 + 20 +-0.007673859629111657 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +40 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.06516986230184571 + 20 +-0.009171395348674838 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +41 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.06880726379767793 + 20 +-0.010668931068237963 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +42 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.07244466529351018 + 20 +-0.0121664667878012 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +43 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.07608206678934243 + 20 +-0.013664002507364326 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +44 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0797194682851746 + 20 +-0.015161538226927451 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +45 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.08335686978100681 + 20 +-0.016659073946490688 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +46 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.08697507316482567 + 20 +-0.01820113484491409 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +47 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.09056097837690974 + 20 +-0.019818103206756144 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +48 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.09414688358899384 + 20 +-0.02143507156859814 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +49 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.09773278880107797 + 20 +-0.023052039930440138 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1013186940131621 + 20 +-0.02466900829228219 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10490459922524614 + 20 +-0.026285976654124188 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10849050443733027 + 20 +-0.027902945015966185 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1120764096494144 + 20 +-0.029519913377808238 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.11566231486149847 + 20 +-0.031136881739650235 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +4F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.11924822007358254 + 20 +-0.03275385010149223 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +50 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.12283412528566659 + 20 +-0.03437081846333434 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +51 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.12642003049775075 + 20 +-0.03598778682517634 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +52 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1300059357098348 + 20 +-0.037604755187018335 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +53 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1335918409219189 + 20 +-0.03922172354886039 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +54 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.13717774613400305 + 20 +-0.040838691910702385 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +55 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.14372856522610838 + 20 +-0.04327496838266798 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +56 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.14995755911305256 + 20 +-0.04558405899092027 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +57 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.15582399389417462 + 20 +-0.04775445959307795 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +58 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.16132787111901833 + 20 +-0.04978616888891935 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +59 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.16646919233712754 + 20 +-0.05167918557822265 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.17124795909804597 + 20 +-0.053433508360766435 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1756641729513177 + 20 +-0.05504913593632882 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.17971783544648637 + 20 +-0.05652606700468815 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18340894813309583 + 20 +-0.05786430026562289 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18673751256069004 + 20 +-0.05906383441891133 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +5F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18970353027881282 + 20 +-0.06012466816433182 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +60 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1923070028370079 + 20 +-0.06104680020166264 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +61 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.19454793178481922 + 20 +-0.0618302292306821 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +62 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.19642631867179056 + 20 +-0.06247495395116859 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +63 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1979421650474658 + 20 +-0.06298097306290046 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +64 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1990954724613888 + 20 +-0.063348285265656 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +65 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2000151253740754 + 20 +-0.06361126741164846 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +66 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.20276874829217467 + 20 +-0.06435484737989106 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +67 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.205246815645514 + 20 +-0.06479549956536129 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +68 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.20836029970179235 + 20 +-0.06480866119589457 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +69 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.21302017272870882 + 20 +-0.06426976949932639 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.21934980512311386 + 20 +-0.06318850973792922 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2223415227029405 + 20 +-0.062085340964722235 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.22495973231109506 + 20 +-0.060292609672437714 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2273012549555563 + 20 +-0.057927912877090926 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.22946291164430282 + 20 +-0.05510884759469742 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +6F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2315415233853133 + 20 +-0.05195301084127257 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +70 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2336339111865664 + 20 +-0.04857799963283188 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +71 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.23583689605604086 + 20 +-0.04510141098539078 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +72 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.23824729900171532 + 20 +-0.04164084191496459 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +73 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.24096194103156843 + 20 +-0.038313889437569026 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +74 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.24407764315357894 + 20 +-0.03523815056921925 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +75 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.24769122637572547 + 20 +-0.032531222325930864 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +76 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.25189527154487523 + 20 +-0.03031250748991371 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +77 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2556710587823801 + 20 +-0.028897726461885542 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +78 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.259176350329518 + 20 +-0.027895338568229233 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +79 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.26261507934256245 + 20 +-0.027062875912520212 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2661911789777869 + 20 +-0.026157870598334243 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2701085823914646 + 20 +-0.02493785472924681 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2742767322269781 + 20 +-0.023287611479191006 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2773537617614356 + 20 +-0.02171288824969564 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.28041119952074683 + 20 +-0.019761313469147057 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +7F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.283386801947182 + 20 +-0.01745533373208813 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +80 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.28621832548301135 + 20 +-0.014817395633061614 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +81 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.28884352657050505 + 20 +-0.011869945766610546 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +82 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.29120016165193335 + 20 +-0.008635430727277682 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +83 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.29322598716956644 + 20 +-0.005136297109605947 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +84 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2948587595656746 + 20 +-0.0013949915081380992 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +85 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.296036235282528 + 20 +0.0025660394825827715 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +86 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2966961707623968 + 20 +0.006724349268014074 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +87 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.29677131237060905 + 20 +0.010823895979494569 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +88 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.29606911246228473 + 20 +0.014181500952964832 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +89 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2945903225789705 + 20 +0.01733619705562711 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.29243224340788126 + 20 +0.020277968983556938 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.28969217563623184 + 20 +0.022996801432829728 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.28646741995123726 + 20 +0.025482679099520955 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2828552770401123 + 20 +0.02772558667970615 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.278953047590072 + 20 +0.029715508869460727 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +8F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.27485803228833117 + 20 +0.03144243036486016 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +90 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.27066753182210485 + 20 +0.032896335861979986 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +91 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2664788468786078 + 20 +0.034067210056895614 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +92 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.26238927814505497 + 20 +0.034945037645682575 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +93 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2584961263086613 + 20 +0.03551980332441629 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +94 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2548966920566418 + 20 +0.035781491789172226 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +95 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2516882760762112 + 20 +0.03572008773602592 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +96 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.24842500237321694 + 20 +0.03541667366752854 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +97 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.24477200340806995 + 20 +0.035075636354958084 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +98 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.24082262719541409 + 20 +0.03471373073340123 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +99 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.23664618716419447 + 20 +0.034341520250685376 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +9A +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.23231199674335595 + 20 +0.03396956835463805 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +9B +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.22788936936184365 + 20 +0.033608438493086756 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +9C +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.22344761844860245 + 20 +0.03326869411385891 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +9D +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.21905605743257744 + 20 +0.03296089866478208 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +9E +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.21478399974271362 + 20 +0.03269561559368367 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +9F +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2107007588079559 + 20 +0.03248340834839114 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A0 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.20687564805724937 + 20 +0.03233484037673212 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A1 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.2033779809195389 + 20 +0.032260475126533905 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A2 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.20027707082376967 + 20 +0.03227087604562412 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A3 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.19698654858711234 + 20 +0.03238592413156666 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A4 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.19318285904709753 + 20 +0.0325402622832367 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A5 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18931970887041638 + 20 +0.03270409348936071 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A6 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18539709781429586 + 20 +0.032877414975028385 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A7 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.18141502563596268 + 20 +0.033060223965329594 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A8 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.17737349209264375 + 20 +0.03325251768535398 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +A9 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1726381225467503 + 20 +0.03338911516149823 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +AA +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1659840849214246 + 20 +0.033316250989470886 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +AB +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.15964888467818186 + 20 +0.03324643511774267 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +AC +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1536325212730734 + 20 +0.033179667593903006 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +AD +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.14793499416215058 + 20 +0.03311594846554122 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +AE +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.14255630280146458 + 20 +0.03305527778024664 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +AF +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.13749644664706695 + 20 +0.032997655585608576 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B0 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.13275542515500882 + 20 +0.03294308192921647 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B1 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.12833323778134154 + 20 +0.032891556858659576 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B2 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.12422988398211648 + 20 +0.032843080421527227 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B3 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.12044536321338492 + 20 +0.03279765266540885 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B4 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.11697967493119824 + 20 +0.03275527363789377 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B5 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1138328185916076 + 20 +0.032715943386571245 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B6 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.11100479365066446 + 20 +0.03267966195903066 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B7 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10849559956442015 + 20 +0.03264642940286139 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B8 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.1063052357889259 + 20 +0.03261624576565281 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +B9 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10443370178023303 + 20 +0.03258911109499418 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +BA +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10288099699439296 + 20 +0.03256502543847484 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +BB +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.10164712088745684 + 20 +0.03254398884368426 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +BC +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.09994337892822275 + 20 +0.03249160766593884 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +BD +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.09722703434661115 + 20 +0.03234615182366568 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +BE +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.093771592897088 + 20 +0.032118852135510856 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +BF +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.08970115316225025 + 20 +0.03182178565127325 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C0 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0851398137246947 + 20 +0.03146702942075169 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C1 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.08021167316701835 + 20 +0.031066660493744958 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C2 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.07504083007181808 + 20 +0.03063275592005199 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C3 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.06975138302169073 + 20 +0.03017739274947162 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C4 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.06446743059923335 + 20 +0.029712648031802624 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C5 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0593130713870427 + 20 +0.029250598816843887 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C6 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.05441240396771574 + 20 +0.028803322154394295 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C7 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.04988952692384935 + 20 +0.028382895094252625 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C8 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.045868538838040485 + 20 +0.028001394686217707 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +C9 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.04247353829288597 + 20 +0.027670897980088482 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +CA +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.03982862387098274 + 20 +0.027403482025663728 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +CB +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.03805789415492772 + 20 +0.02721122387274233 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +CC +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.036264793274075446 + 20 +0.0269216897066466 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +CD +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.032322371975447184 + 20 +0.025941611713292345 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +CE +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.027830560327799014 + 20 +0.0246016350498387 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +CF +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.023812359863461197 + 20 +0.02324602396599168 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D0 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.021294698693667236 + 20 +0.022218126153409845 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D1 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.01872618634040374 + 20 +0.02064227908050753 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D2 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.015093256261121235 + 20 +0.018146755728090103 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D3 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.011096281688223297 + 20 +0.01523911360159308 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D4 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.007435635854113387 + 20 +0.012426910206452146 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D5 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.004811691991195022 + 20 +0.010217703048102877 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D6 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.003023760328472236 + 20 +0.008011242975078714 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D7 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +-0.0014008764041681387 + 20 +0.004592179313590128 + 30 +0.0 + 70 +0 + 0 +VERTEX + 5 +D8 +330 +17 +100 +AcDbEntity + 8 +0 +100 +AcDbVertex +100 +AcDb2dVertex + 10 +0.0 + 20 +0.0 + 30 +0.0 + 70 +0 + 0 +SEQEND + 5 +2E +330 +17 +100 +AcDbEntity + 8 +0 + 0 +ENDSEC + 0 +SECTION + 2 +OBJECTS + 0 +DICTIONARY + 5 +A +330 +0 +100 +AcDbDictionary +281 +1 + 3 +ACAD_COLOR +350 +B + 3 +ACAD_GROUP +350 +C + 3 +ACAD_LAYOUT +350 +D + 3 +ACAD_MATERIAL +350 +E + 3 +ACAD_MLEADERSTYLE +350 +F + 3 +ACAD_MLINESTYLE +350 +10 + 3 +ACAD_PLOTSETTINGS +350 +11 + 3 +ACAD_PLOTSTYLENAME +350 +12 + 3 +ACAD_SCALELIST +350 +14 + 3 +ACAD_TABLESTYLE +350 +15 + 3 +ACAD_VISUALSTYLE +350 +16 + 0 +DICTIONARY + 5 +B +330 +A +100 +AcDbDictionary +281 +1 + 0 +DICTIONARY + 5 +C +330 +A +100 +AcDbDictionary +281 +1 + 0 +DICTIONARY + 5 +D +330 +A +100 +AcDbDictionary +281 +1 + 3 +Model +350 +1A + 3 +Layout1 +350 +1E + 0 +DICTIONARY + 5 +E +330 +A +100 +AcDbDictionary +281 +1 + 3 +ByBlock +350 +1F + 3 +ByLayer +350 +20 + 3 +Global +350 +21 + 0 +DICTIONARY + 5 +F +330 +A +100 +AcDbDictionary +281 +1 + 3 +Standard +350 +2C + 0 +DICTIONARY + 5 +10 +330 +A +100 +AcDbDictionary +281 +1 + 3 +Standard +350 +22 + 0 +DICTIONARY + 5 +11 +330 +A +100 +AcDbDictionary +281 +1 + 0 +ACDBDICTIONARYWDFLT + 5 +12 +330 +A +100 +AcDbDictionary +281 +1 + 3 +Normal +350 +13 +100 +AcDbDictionaryWithDefault +340 +13 + 0 +ACDBPLACEHOLDER + 5 +13 +330 +12 + 0 +DICTIONARY + 5 +14 +330 +A +100 +AcDbDictionary +281 +1 + 0 +DICTIONARY + 5 +15 +330 +A +100 +AcDbDictionary +281 +1 + 0 +DICTIONARY + 5 +16 +330 +A +100 +AcDbDictionary +281 +1 + 0 +LAYOUT + 5 +1A +330 +D +100 +AcDbPlotSettings + 1 + + 2 +Adobe PDF + 4 +A4 + 6 + + 40 +3.175 + 41 +3.175 + 42 +3.175 + 43 +3.175 + 44 +209.91 + 45 +297.03 + 46 +0.0 + 47 +0.0 + 48 +0.0 + 49 +0.0 +140 +0.0 +141 +0.0 +142 +1.0 +143 +1.0 + 70 +1024 + 72 +0 + 73 +1 + 74 +5 + 7 + + 75 +16 + 76 +0 + 77 +2 + 78 +300 +147 +1.0 +148 +0.0 +149 +0.0 +100 +AcDbLayout + 1 +Model + 70 +1 + 71 +0 + 10 +-3.175 + 20 +-3.175 + 11 +293.857 + 21 +206.735 + 12 +0.0 + 22 +0.0 + 32 +0.0 + 14 +29.068 + 24 +20.356 + 34 +0.0 + 15 +261.614 + 25 +183.204 + 35 +0.0 +146 +0.0 + 13 +0.0 + 23 +0.0 + 33 +0.0 + 16 +1.0 + 26 +0.0 + 36 +0.0 + 17 +0.0 + 27 +1.0 + 37 +0.0 + 76 +1 +330 +17 + 0 +LAYOUT + 5 +1E +330 +D +100 +AcDbPlotSettings + 1 + + 2 +Adobe PDF + 4 +A4 + 6 + + 40 +3.175 + 41 +3.175 + 42 +3.175 + 43 +3.175 + 44 +209.91 + 45 +297.03 + 46 +0.0 + 47 +0.0 + 48 +0.0 + 49 +0.0 +140 +0.0 +141 +0.0 +142 +1.0 +143 +1.0 + 70 +0 + 72 +0 + 73 +1 + 74 +5 + 7 + + 75 +16 + 76 +0 + 77 +2 + 78 +300 +147 +1.0 +148 +0.0 +149 +0.0 +100 +AcDbLayout + 1 +Layout1 + 70 +1 + 71 +1 + 10 +-3.175 + 20 +-3.175 + 11 +293.857 + 21 +206.735 + 12 +0.0 + 22 +0.0 + 32 +0.0 + 14 +29.068 + 24 +20.356 + 34 +0.0 + 15 +261.614 + 25 +183.204 + 35 +0.0 +146 +0.0 + 13 +0.0 + 23 +0.0 + 33 +0.0 + 16 +1.0 + 26 +0.0 + 36 +0.0 + 17 +0.0 + 27 +1.0 + 37 +0.0 + 76 +1 +330 +1B + 0 +MATERIAL + 5 +1F +102 +{ACAD_REACTORS +330 +E +102 +} +330 +E +100 +AcDbMaterial + 1 +ByBlock + 2 + + 70 +0 + 40 +1.0 + 71 +1 + 41 +1.0 + 91 +-1023410177 + 42 +1.0 + 72 +1 + 3 + + 73 +1 + 74 +1 + 75 +1 + 44 +0.5 + 73 +0 + 45 +1.0 + 46 +1.0 + 77 +1 + 4 + + 78 +1 + 79 +1 +170 +1 + 48 +1.0 +171 +1 + 6 + +172 +1 +173 +1 +174 +1 +140 +1.0 +141 +1.0 +175 +1 + 7 + +176 +1 +177 +1 +178 +1 +143 +1.0 +179 +1 + 8 + +270 +1 +271 +1 +272 +1 +145 +1.0 +146 +1.0 +273 +1 + 9 + +274 +1 +275 +1 +276 +1 + 42 +1.0 + 72 +1 + 3 + + 73 +1 + 74 +1 + 75 +1 + 94 +63 + 0 +MATERIAL + 5 +20 +102 +{ACAD_REACTORS +330 +E +102 +} +330 +E +100 +AcDbMaterial + 1 +ByLayer + 2 + + 70 +0 + 40 +1.0 + 71 +1 + 41 +1.0 + 91 +-1023410177 + 42 +1.0 + 72 +1 + 3 + + 73 +1 + 74 +1 + 75 +1 + 44 +0.5 + 73 +0 + 45 +1.0 + 46 +1.0 + 77 +1 + 4 + + 78 +1 + 79 +1 +170 +1 + 48 +1.0 +171 +1 + 6 + +172 +1 +173 +1 +174 +1 +140 +1.0 +141 +1.0 +175 +1 + 7 + +176 +1 +177 +1 +178 +1 +143 +1.0 +179 +1 + 8 + +270 +1 +271 +1 +272 +1 +145 +1.0 +146 +1.0 +273 +1 + 9 + +274 +1 +275 +1 +276 +1 + 42 +1.0 + 72 +1 + 3 + + 73 +1 + 74 +1 + 75 +1 + 94 +63 + 0 +MATERIAL + 5 +21 +102 +{ACAD_REACTORS +330 +E +102 +} +330 +E +100 +AcDbMaterial + 1 +Global + 2 + + 70 +0 + 40 +1.0 + 71 +1 + 41 +1.0 + 91 +-1023410177 + 42 +1.0 + 72 +1 + 3 + + 73 +1 + 74 +1 + 75 +1 + 44 +0.5 + 73 +0 + 45 +1.0 + 46 +1.0 + 77 +1 + 4 + + 78 +1 + 79 +1 +170 +1 + 48 +1.0 +171 +1 + 6 + +172 +1 +173 +1 +174 +1 +140 +1.0 +141 +1.0 +175 +1 + 7 + +176 +1 +177 +1 +178 +1 +143 +1.0 +179 +1 + 8 + +270 +1 +271 +1 +272 +1 +145 +1.0 +146 +1.0 +273 +1 + 9 + +274 +1 +275 +1 +276 +1 + 42 +1.0 + 72 +1 + 3 + + 73 +1 + 74 +1 + 75 +1 + 94 +63 + 0 +MLINESTYLE + 5 +22 +102 +{ACAD_REACTORS +330 +10 +102 +} +330 +10 +100 +AcDbMlineStyle + 2 +Standard + 70 +0 + 3 + + 62 +256 + 51 +90.0 + 52 +90.0 + 71 +2 + 49 +0.5 + 62 +256 + 6 +BYLAYER + 49 +-0.5 + 62 +256 + 6 +BYLAYER + 0 +MLEADERSTYLE + 5 +2C +102 +{ACAD_REACTORS +330 +F +102 +} +330 +F +100 +AcDbMLeaderStyle +179 +2 +170 +2 +171 +1 +172 +0 + 90 +2 + 40 +0.0 + 41 +0.0 +173 +1 + 91 +-1056964608 + 92 +-2 +290 +1 + 42 +2.0 +291 +1 + 43 +8.0 + 3 +Standard + 44 +4.0 +300 + +342 +29 +174 +1 +175 +1 +176 +0 +178 +1 + 93 +-1056964608 + 45 +4.0 +292 +0 +297 +0 + 46 +4.0 + 94 +-1056964608 + 47 +1.0 + 49 +0.0 +140 +1.0 +294 +1 +141 +0.0 +177 +0 +142 +1.0 +295 +0 +296 +0 +143 +3.75 +271 +0 +272 +9 +272 +9 + 0 +ENDSEC + 0 +EOF diff --git a/src/callbacks/info.jl b/src/callbacks/info.jl index 69d2d38567..cb494e12ea 100644 --- a/src/callbacks/info.jl +++ b/src/callbacks/info.jl @@ -1,5 +1,7 @@ mutable struct InfoCallback start_time::Float64 + last_performance_print_time::Float64 + ncalls_at_last_print::Int interval::Int end @@ -33,7 +35,7 @@ beginning of a simulation and then resets the timer. When the returned callback directly, the current timer values are shown. """ function InfoCallback(; interval=0, reset_threads=true) - info_callback = InfoCallback(0.0, interval) + info_callback = InfoCallback(0.0, 0.0, 0, interval) function initialize(cb, u, t, integrator) initialize_info_callback(cb, u, t, integrator; @@ -65,6 +67,17 @@ function (info_callback::InfoCallback)(integrator) integrator.stats.naccept, integrator.dt, t, sim_time_percentage), 71) * @sprintf("│ run time: %.4e s", runtime_absolute)) + + if condition_integrator_interval(integrator, 1 * info_callback.interval) + runtime_since_last_print = (time_ns() - info_callback.last_performance_print_time) + info_callback.last_performance_print_time = time_ns() + n_calls = integrator.stats.nf - info_callback.ncalls_at_last_print + info_callback.ncalls_at_last_print = integrator.stats.nf + pid = runtime_since_last_print / n_calls / sum(nparticles.(integrator.p.semi.systems)) + println("─"^100) + println(@sprintf("time/particle/rhs: %.4e ns", pid)) + println("─"^100) + end end # Tell OrdinaryDiffEq that u has not been modified @@ -143,6 +156,8 @@ function initialize_info_callback(discrete_callback, u, t, integrator; # Save current time as start_time info_callback.start_time = time_ns() + info_callback.last_performance_print_time = time_ns() + info_callback.ncalls_at_last_print = 0 return nothing end @@ -256,9 +271,16 @@ end function print_summary(integrator) println("─"^100) - println("Trixi simulation finished. Final time: ", integrator.t, - " Time steps: ", integrator.stats.naccept, " (accepted), ", - integrator.iter, " (total)") + # Print the final time as string to let Julia decide on the formatting + @printf(" %-31s %10s\n", "Final time:", string(integrator.t)) + @printf(" %-31s %10d (accepted) %10d (total)\n", + "Time steps:", integrator.stats.naccept, integrator.iter) + if !isnothing(integrator.p.split_integration_data) + @printf(" %-31s %10d (accepted) %10d (total)\n", + "Split integration time steps:", + integrator.p.split_integration_data.integrator.stats.naccept, + integrator.p.split_integration_data.integrator.iter) + end println("─"^100) println() diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index b78ae8013f..d49553bf47 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -95,7 +95,7 @@ initial_condition = InitialCondition(; coordinates, velocity=x -> 2x, mass=1.0, └──────────────────────────────────────────────────────────────────────────────────────────────────┘ ``` """ -struct InitialCondition{ELTYPE, C, MATRIX, VECTOR} +struct InitialCondition{ELTYPE, C, MATRIX <: AbstractArray{ELTYPE}, VECTOR <: AbstractVector{ELTYPE}} particle_spacing :: ELTYPE coordinates :: C # Array{coordinates_eltype, 2} velocity :: MATRIX # Array{ELTYPE, 2} diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl index a15bf88f07..d3cc223b88 100644 --- a/src/general/interpolation.jl +++ b/src/general/interpolation.jl @@ -155,7 +155,7 @@ function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_ v_ode, u_ode; include_wall_velocity=false, smoothing_length=initial_smoothing_length(ref_system), cut_off_bnd=true, clip_negative_pressure=false, - output_directory="out", filename="plane") + output_directory="out", filename="plane", pvd=nothing, t=-1) # Don't filter out particles without neighbors to keep 2D grid structure filter_no_neighbors = false @trixi_timeit timer() "interpolate plane" begin @@ -177,6 +177,10 @@ function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_ vtk["density"] = density vtk["velocity"] = velocity vtk["pressure"] = pressure + + if !isnothing(pvd) && t >= 0 + pvd[t] = vtk + end end end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index f2f651fa91..1551fb09a9 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -135,6 +135,10 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; calculate_flow_rate, cache) end +function Base.:(==)(system1::OpenBoundarySystem, system2::OpenBoundarySystem) + return system1.mass === system2.mass +end + function initialize!(system::OpenBoundarySystem, semi) (; boundary_zones) = system diff --git a/src/schemes/boundary/prescribed_motion.jl b/src/schemes/boundary/prescribed_motion.jl index 6e5fff7272..4ea5b19ef2 100644 --- a/src/schemes/boundary/prescribed_motion.jl +++ b/src/schemes/boundary/prescribed_motion.jl @@ -83,9 +83,10 @@ function initialize_prescribed_motion!(::Nothing, initial_condition, end function (prescribed_motion::PrescribedMotion)(coordinates, velocity, acceleration, - ismoving, system, semi, t) + ismoving, system, semi, t_) (; movement_function, is_moving, moving_particles) = prescribed_motion + t = convert(Float64, t_) ismoving[] = is_moving(t) is_moving(t) || return nothing @@ -160,19 +161,22 @@ PrescribedMotion{...} function OscillatingMotion2D(; frequency, translation_vector, rotation_angle, rotation_center, rotation_phase_offset=0, tspan=(-Inf, Inf), ramp_up_tspan=(0.0, 0.0), moving_particles=nothing) - translation_vector_ = SVector{2}(translation_vector) - rotation_center_ = SVector{2}(rotation_center) + translation_vector_ = SVector{2, Float64}(translation_vector) + rotation_center_ = SVector{2, Float64}(rotation_center) + frequency_ = convert(Float64, frequency) + rotation_angle_ = convert(Float64, rotation_angle) + rotation_phase_offset_ = convert(Float64, rotation_phase_offset) @inline function movement_function(x, t) if isfinite(tspan[1]) t = t - tspan[1] end - sin_scaled = sin(frequency * 2pi * t) + sin_scaled = sin(frequency_ * 2pi * t) translation = sin_scaled * translation_vector_ x_centered = x .- rotation_center_ - sin_scaled_offset = sin(2pi * (frequency * t - rotation_phase_offset)) - angle = rotation_angle * sin_scaled_offset + sin_scaled_offset = sin(2pi * (frequency_ * t - rotation_phase_offset_)) + angle = rotation_angle_ * sin_scaled_offset rotated = SVector(x_centered[1] * cos(angle) - x_centered[2] * sin(angle), x_centered[1] * sin(angle) + x_centered[2] * cos(angle)) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index eb6338f1cc..a81807f5a4 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -176,7 +176,7 @@ end # Artificial density diffusion should only be applied to systems representing a fluid # with the same physical properties i.e. density and viscosity. # TODO: shouldn't be applied to particles on the interface (depends on PR #539) - if particle_system === neighbor_system + if particle_system == neighbor_system density_diffusion!(drho_particle, density_diffusion(particle_system), particle_system, particle, neighbor, pos_diff, distance, m_b, rho_a, rho_b, grad_kernel) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 2e5474a33c..00eb47409a 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -173,6 +173,11 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, cache) end +@inline function Base.:(==)(system1::WeaklyCompressibleSPHSystem, + system2::WeaklyCompressibleSPHSystem) + return system1.mass === system2.mass +end + function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) @nospecialize system # reduce precompilation time diff --git a/src/schemes/structure/total_lagrangian_sph/rhs.jl b/src/schemes/structure/total_lagrangian_sph/rhs.jl index d3ee370237..232f011838 100644 --- a/src/schemes/structure/total_lagrangian_sph/rhs.jl +++ b/src/schemes/structure/total_lagrangian_sph/rhs.jl @@ -108,6 +108,8 @@ function interact!(dv, v_particle_system, u_particle_system, neighbor_system::AbstractFluidSystem, semi; integrate_tlsph=semi.integrate_tlsph[], eachparticle=each_integrated_particle(particle_system)) + (; boundary_model) = particle_system + # Skip interaction if TLSPH systems are integrated separately integrate_tlsph || return dv diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 82f02dcc40..eb19636745 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -454,6 +454,8 @@ function update_quantities!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode # Precompute PK1 stress tensor @trixi_timeit timer() "stress tensor" compute_pk1_corrected!(system, semi) + # @trixi_timeit timer() "KV damping" apply_kelvin_voigt_damping!(system.pk1_rho2, v, system, semi) + return system end @@ -552,6 +554,54 @@ end return deformation_grad end +@inline function apply_kelvin_voigt_damping!(pk1_rho2, v, system, semi) + (; mass, material_density) = system + + # Compute bulk modulus from Young's modulus and Poisson's ratio. + # See the table at the end of https://en.wikipedia.org/wiki/Lam%C3%A9_parameters + # TODO Should we compute the sound speed per particle and then use the maximum? + E = maximum(system.young_modulus) + K = E / (ndims(system) * (1 - 2 * maximum(system.poisson_ratio))) + + # Newton–Laplace equation + sound_speed = sqrt(K / minimum(system.material_density)) + + # Loop over all pairs of particles and neighbors within the kernel cutoff + initial_coords = initial_coordinates(system) + foreach_point_neighbor(system, system, initial_coords, initial_coords, + semi) do particle, neighbor, initial_pos_diff, initial_distance + # Only consider particles with a distance > 0. + # See `src/general/smoothing_kernels.jl` for more details. + initial_distance^2 < eps(initial_smoothing_length(system)^2) && return + + volume = @inbounds mass[neighbor] / material_density[neighbor] + v_diff = @inbounds current_velocity(v, system, particle) - + current_velocity(v, system, neighbor) + + grad_kernel = smoothing_kernel_grad(system, initial_pos_diff, + initial_distance, particle) + + # Multiply by L_{0a} + L = @inbounds correction_matrix(system, particle) + + dF = -volume * v_diff * grad_kernel' * L' + + F = deformation_gradient(system, particle) + + alpha = 1.0 + rho = material_density[particle] + + h = smoothing_length(system, particle) + P_rho2 = alpha / rho * sound_speed * h / 2 * F * (dF' * F + F' * dF) + + for j in 1:ndims(system), i in 1:ndims(system) + @inbounds pk1_rho2[i, j, particle] += P_rho2[i, j] + end + end + + return pk1_rho2 +end + # First Piola-Kirchhoff stress tensor @propagate_inbounds function pk1_stress_tensor(system, particle) (; lame_lambda, lame_mu) = system diff --git a/src/setups/complex_shape.jl b/src/setups/complex_shape.jl index 6a78b412e2..fd7e821ad1 100644 --- a/src/setups/complex_shape.jl +++ b/src/setups/complex_shape.jl @@ -46,16 +46,12 @@ function ComplexShape(geometry; particle_spacing, density, point_in_geometry_algorithm=WindingNumberJacobson(; geometry, hierarchical_winding=true, winding_number_factor=sqrt(eps())), - store_winding_number=false, grid_offset::Real=0.0, + store_winding_number=false, grid_offset=0.0, max_nparticles=10^7, pad_initial_particle_grid=2particle_spacing) if ndims(geometry) == 3 && point_in_geometry_algorithm isa WindingNumberHormann throw(ArgumentError("`WindingNumberHormann` only supports 2D geometries")) end - if grid_offset < 0.0 - throw(ArgumentError("only a positive `grid_offset` is supported")) - end - grid = particle_grid(geometry, particle_spacing; padding=pad_initial_particle_grid, grid_offset, max_nparticles) @@ -145,7 +141,12 @@ function particle_grid(geometry, particle_spacing; padding=2particle_spacing, grid_offset=0.0, max_nparticles=10^7) (; max_corner) = geometry + # First subtract the grid offset, then add it again after rounding. + # This is making sure that `min_corner` is `n * particle_spacing` + # away from `grid_offset`, so that a particle is placed at `grid_offset`. min_corner = geometry.min_corner .- grid_offset .- padding + min_corner = floor.(Int, min_corner ./ particle_spacing) .* particle_spacing + min_corner = min_corner .+ grid_offset n_particles_per_dimension = Tuple(ceil.(Int, (max_corner .- min_corner .+ 2padding) ./