diff --git a/docs/src/general/initial_condition.md b/docs/src/general/initial_condition.md index 4f6ba6f723..28b1ae15af 100644 --- a/docs/src/general/initial_condition.md +++ b/docs/src/general/initial_condition.md @@ -11,3 +11,8 @@ Pages = [joinpath("general", "initial_condition.jl")] Modules = [TrixiParticles] Pages = map(file -> joinpath("setups", file), readdir(joinpath("..", "src", "setups"))) ``` + +```@autodocs +Modules = [TrixiParticles] +Filter = t -> t === TrixiParticles.OrientedBoundingBox +``` diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 3a34b10914..145e3fcaaf 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -330,6 +330,15 @@ @Article{Ganzenmueller2015 publisher = {Elsevier {BV}}, } +@Book{Gasser2021, + author = {Gasser, T. Christian}, + title = {Vascular Biomechanics: Concepts, Models, and Applications}, + year = {2021}, + isbn = {9783030709662}, + doi = {10.1007/978-3-030-70966-2}, + publisher = {Springer International Publishing}, +} + @Article{Giles1990, author = {Giles, Michael B.}, title = {Nonreflecting boundary conditions for {Euler} equation calculations}, @@ -828,6 +837,20 @@ @Article{Wendland1995 publisher = {Springer Science and Business Media LLC}, } +@Article{Westerhof2008, + author = {Westerhof, Nico and Lankhaar, Jan-Willem and Westerhof, Berend E.}, + journal = {Medical \& Biological Engineering \& Computing}, + title = {The arterial Windkessel}, + year = {2008}, + issn = {1741-0444}, + month = jun, + number = {2}, + pages = {131--141}, + volume = {47}, + doi = {10.1007/s11517-008-0359-2}, + publisher = {Springer Science and Business Media LLC}, +} + @Article{Zhang2025, author = {Zhang, Shuoguo and Fan, Yu and Wu, Dong and Zhang, Chi and Hu, Xiangyu}, journal = {Physics of Fluids}, diff --git a/docs/src/systems/boundary.md b/docs/src/systems/boundary.md index 45c3e602be..95a6f2bc63 100644 --- a/docs/src/systems/boundary.md +++ b/docs/src/systems/boundary.md @@ -342,3 +342,9 @@ without the need to specifically identify those near the free surface. To further handle incomplete kernel support, for example in the viscous term of the momentum equation, the updated velocity of particles within the [`BoundaryZone`](@ref) is projected onto the face normal, so that only the component in flow direction is kept. + +# Pressure Models +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "boundary", "open_boundary", "pressure_model.jl")] +``` diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index c7b88a8982..dc0c19c6ea 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -90,13 +90,14 @@ export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureEx export FirstOrderMirroring, ZerothOrderMirroring, SimpleMirroring export HertzContactModel, LinearContactModel export PrescribedMotion, OscillatingMotion2D +export RCRWindkesselModel export examples_dir, validation_dir export trixi2vtk, vtk2trixi export RectangularTank, RectangularShape, SphereShape, ComplexShape export ParticlePackingSystem, SignedDistanceField export WindingNumberHormann, WindingNumberJacobson export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry, - sample_boundary, planar_geometry_to_face + sample_boundary, planar_geometry_to_face, OrientedBoundingBox, is_in_oriented_box export SourceTermDamping export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection, GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 1c136ebc69..44af5fdb6e 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -1149,5 +1149,6 @@ function set_system_links(system::OpenBoundarySystem, semi) system.buffer, system.pressure_acceleration_formulation, system.shifting_technique, + system.pressure_model_values, system.cache) end diff --git a/src/preprocessing/geometries/geometries.jl b/src/preprocessing/geometries/geometries.jl index 01fbc69eed..dab118fe59 100644 --- a/src/preprocessing/geometries/geometries.jl +++ b/src/preprocessing/geometries/geometries.jl @@ -92,7 +92,7 @@ face, face_normal = planar_geometry_to_face(planar_geometry) function planar_geometry_to_face(planar_geometry::TriangleMesh) face_normal = normalize(sum(planar_geometry.face_normals) / nfaces(planar_geometry)) - face_vertices = oriented_bounding_box(stack(planar_geometry.vertices)) + face_vertices, _, _ = oriented_bounding_box(stack(planar_geometry.vertices)) # Vectors spanning the face edge1 = face_vertices[:, 2] - face_vertices[:, 1] @@ -122,11 +122,247 @@ function oriented_bounding_box(point_cloud) aligned_coords = eigen_vectors' * centered_data - min_corner = minimum(aligned_coords, dims=2) - max_corner = maximum(aligned_coords, dims=2) + min_corner = minimum(aligned_coords, dims=2) .- eps() + max_corner = maximum(aligned_coords, dims=2) .+ eps() - face_vertices = hcat(min_corner, max_corner, - [min_corner[1], max_corner[2], min_corner[3]]) + if length(min_corner) == 2 + rect_coords = hcat(min_corner, max_corner, [min_corner[1], max_corner[2]]) + else + rect_coords = hcat(min_corner, max_corner, + [min_corner[1], max_corner[2], min_corner[3]]) + end + + rotated_rect_coords = eigen_vectors * rect_coords .+ means + + return rotated_rect_coords, eigen_vectors, (min_corner, max_corner) +end + +""" + OrientedBoundingBox(; box_origin, orientation_matrix, edge_lengths, local_axis_scale::Tuple)) + OrientedBoundingBox(point_cloud::AbstractMatrix; local_axis_scale::Tuple) + OrientedBoundingBox(geometry::Union{Polygon, TriangleMesh}; local_axis_scale::Tuple) + +Creates an oriented bounding box (rectangle in 2D or cuboid in 3D) that can be +rotated and positioned arbitrarily in space. + +The box is defined either by explicit parameters +or by automatically fitting it around an existing geometry or a point cloud with optional scaling. + +# Arguments +- `point_cloud`: An array where the ``i``-th column holds the coordinates of a point ``i``. +- `geometry`: Geometry returned by [`load_geometry`](@ref). + +# Keywords +- `box_origin`: The corner point from which the box is constructed. +- `orientation_matrix`: A matrix defining the orientation of the box in space. + Each column of the matrix represents one of the local axes of the box (e.g., x, y, z in 3D). + The matrix must be orthogonal, ensuring that the local axes are perpendicular to each other + and normalized. +- `edge_lengths`: The lengths of the edges of the box: + - In 2D: `(width, height)` + - In 3D: `(width, height, depth)` +- `local_axis_scale`: Allows for anisotropic scaling along the oriented axes of the `OrientedBoundingBox` + (the eigenvectors of the geometry's covariance matrix). Default is no scaling. + The tuple components correspond to: + - first element: scaling along the first eigenvector (local x-axis), + - second element: scaling along the second eigenvector (local y-axis), + - third element (only in 3D): scaling along the third eigenvector (local z-axis). + + **Note:** Scaling is always applied in the local `OrientedBoundingBox` + coordinate system, i.e. along its oriented axes. + Scaling along arbitrary world directions is not supported, + as this would break the orthogonality of the spanning vectors. + +# Examples +```jldoctest; output=false, setup = :(using LinearAlgebra: I) +# 2D axis aligned +OrientedBoundingBox(box_origin=[0.0, 0.0], orientation_matrix=I(2), edge_lengths=[2.0, 1.0]) + +# 2D rotated +orientation_matrix = [cos(π/4) -sin(π/4); + sin(π/4) cos(π/4)] +OrientedBoundingBox(; box_origin=[0.0, 0.0], orientation_matrix, edge_lengths=[2.0, 1.0]) + +# 2D point cloud +# Create a point cloud from a spherical shape (quarter circle sector) +shape = SphereShape(0.1, 1.5, (0.2, 0.4), 1.0, n_layers=4, + sphere_type=RoundSphere(; start_angle=0, end_angle=π/4)) +OrientedBoundingBox(shape.coordinates) + +# 3D rotated +orientation_matrix = [cos(π/4) -sin(π/4) 0.0; # x-axis after rotation + sin(π/4) cos(π/4) 0.0; # y-axis after rotation + 0.0 0.0 1.0] # z-axis remains unchanged +OrientedBoundingBox(; box_origin=[0.5, -0.2, 0.0], orientation_matrix, + edge_lengths=[1.0, 2.0, 3.0]) + +# output +┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ +│ OrientedBoundingBox (3D) │ +│ ════════════════════════ │ +│ box origin: ………………………………………………… [0.5, -0.2, 0.0] │ +│ edge lengths: …………………………………………… (1.0, 2.0, 3.0) │ +└──────────────────────────────────────────────────────────────────────────────────────────────────┘ +``` +""" +struct OrientedBoundingBox{NDIMS, ELTYPE <: Real, SV} + box_origin :: SVector{NDIMS, ELTYPE} + spanning_vectors :: SV +end + +function OrientedBoundingBox(; box_origin, orientation_matrix, edge_lengths, + local_axis_scale::Tuple=ntuple(_ -> 0, length(box_origin))) + box_origin_ = SVector(box_origin...) + NDIMS = length(box_origin_) + + @assert size(orientation_matrix) == (NDIMS, NDIMS) + @assert length(edge_lengths) == NDIMS + @assert length(local_axis_scale) == NDIMS + + spanning_vectors = ntuple(i -> SVector{NDIMS}(orientation_matrix[:, i] * + edge_lengths[i]), NDIMS) + + prod(local_axis_scale) > 0 || return OrientedBoundingBox(box_origin_, spanning_vectors) + + return scaled_bbox(SVector(local_axis_scale), box_origin_, spanning_vectors) +end + +function OrientedBoundingBox(point_cloud::AbstractMatrix; + local_axis_scale::Tuple=ntuple(_ -> 0, size(point_cloud, 1))) + NDIMS = size(point_cloud, 1) + vertices, eigen_vectors, (min_corner, max_corner) = oriented_bounding_box(point_cloud) + + # Use the first vertex as box origin (bottom-left-front corner) + box_origin = SVector{NDIMS}(vertices[:, 1]) + + # Calculate edge lengths from min/max corners + edge_lengths = max_corner - min_corner + + return OrientedBoundingBox(; box_origin, orientation_matrix=eigen_vectors, edge_lengths, + local_axis_scale) +end + +function OrientedBoundingBox(geometry::Union{Polygon, TriangleMesh}; + local_axis_scale::Tuple=ntuple(_ -> 0, ndims(geometry))) + point_cloud = stack(geometry.vertices) + + return OrientedBoundingBox(point_cloud; local_axis_scale) +end + +@inline Base.ndims(::OrientedBoundingBox{NDIMS}) where {NDIMS} = NDIMS + +function Base.show(io::IO, ::MIME"text/plain", box::OrientedBoundingBox) + @nospecialize box # reduce precompilation time + + if get(io, :compact, false) + show(io, box) + else + summary_header(io, "OrientedBoundingBox ($(ndims(box))D)") + summary_line(io, "box origin", box.box_origin) + summary_line(io, "edge lengths", norm.(box.spanning_vectors)) + summary_footer(io) + end +end + +function scaled_bbox(local_axis_scale::SVector{2}, box_origin, spanning_vectors) + # Uniform scaling about the center, center remains unchanged + v1, v2 = spanning_vectors + center = box_origin + (v1 + v2) / 2 + + # Scaling factor per oriented axis + s1, s2 = local_axis_scale + + # New spanning vectors + v1p, v2p = s1 * v1, s2 * v2 + + new_origin = center - (v1p + v2p) / 2 + + return OrientedBoundingBox(new_origin, (v1p, v2p)) +end + +function scaled_bbox(local_axis_scale::SVector{3}, box_origin, spanning_vectors) + # Uniform scaling about the center, center remains unchanged + v1, v2, v3 = spanning_vectors + center = box_origin + (v1 + v2 + v3) / 2 + + # Scaling factor per oriented axis + s1, s2, s3 = local_axis_scale + + # New spanning vectors + v1p, v2p, v3p = s1 * v1, s2 * v2, s3 * v3 + + new_origin = center - (v1p + v2p + v3p) / 2 + + return OrientedBoundingBox(new_origin, (v1p, v2p, v3p)) +end + +function is_in_oriented_box(coordinates::AbstractArray, box) + is_in_box = fill(false, size(coordinates, 2)) + @threaded default_backend(coordinates) for particle in axes(coordinates, 2) + particle_coords = current_coords(coordinates, box, particle) + is_in_box[particle] = is_in_oriented_box(particle_coords, box) + end + + return findall(is_in_box) +end + +@inline function is_in_oriented_box(particle_coords::SVector{NDIMS}, box) where {NDIMS} + (; spanning_vectors, box_origin) = box + relative_position = particle_coords - box_origin + + for dim in 1:NDIMS + span_dim = spanning_vectors[dim] + # Checks whether the projection of the particle position + # falls within the range of the zone + if !(0 <= dot(relative_position, span_dim) <= dot(span_dim, span_dim)) + + # Particle is not in box + return false + end + end + + # Particle is in box + return true +end + +@inline function Base.intersect(initial_condition::InitialCondition, + boxes::OrientedBoundingBox...) + (; coordinates, density, mass, velocity, pressure, particle_spacing) = initial_condition + box = first(boxes) + + if ndims(box) != ndims(initial_condition) + throw(ArgumentError("all passed `OrientedBoundingBox`s must have the same dimensionality as the initial condition")) + end + + keep_indices = is_in_oriented_box(coordinates, box) + + result = InitialCondition(; coordinates=coordinates[:, keep_indices], + density=density[keep_indices], mass=mass[keep_indices], + velocity=velocity[:, keep_indices], + pressure=pressure[keep_indices], + particle_spacing) + + return intersect(result, Base.tail(boxes)...) +end + +@inline function Base.setdiff(initial_condition::InitialCondition, + boxes::OrientedBoundingBox...) + (; coordinates, density, mass, velocity, pressure, particle_spacing) = initial_condition + box = first(boxes) + + if ndims(box) != ndims(initial_condition) + throw(ArgumentError("all passed `OrientedBoundingBox`s must have the same dimensionality as the initial condition")) + end + + delete_indices = is_in_oriented_box(coordinates, box) + keep_indices = setdiff(eachparticle(initial_condition), delete_indices) + + result = InitialCondition(; coordinates=coordinates[:, keep_indices], + density=density[keep_indices], + mass=mass[keep_indices], + velocity=velocity[:, keep_indices], + pressure=pressure[keep_indices], + particle_spacing) - return eigen_vectors * face_vertices .+ means + return setdiff(result, Base.tail(boxes)...) end diff --git a/src/preprocessing/geometries/io.jl b/src/preprocessing/geometries/io.jl index d910e4ed6e..72a1d929b8 100644 --- a/src/preprocessing/geometries/io.jl +++ b/src/preprocessing/geometries/io.jl @@ -132,7 +132,8 @@ end # FileIO.jl docs: # https://juliaio.github.io/FileIO.jl/stable/implementing/#All-at-once-I/O:-implementing-load-and-save -function load(fn::FileIO.File{FileIO.format"STL_BINARY"}; element_types...) +function load(fn::Union{FileIO.File{FileIO.format"STL_BINARY"}, + FileIO.File{FileIO.format"STL_ASCII"}}; element_types...) open(fn) do s FileIO.skipmagic(s) # skip over the magic bytes load(s; element_types...) @@ -155,6 +156,81 @@ function load(fs::FileIO.Stream{FileIO.format"STL_BINARY"}; ELTYPE=Float64) return TriangleMesh(face_vertices, normals, vertices) end +function load(fs::FileIO.Stream{FileIO.format"STL_ASCII"}; ELTYPE=Float64) + # ASCII STL (solid ... endsolid) + io = FileIO.stream(fs) + + # Unlike the binary STL format, the ASCII STL does not include a face count up front; + # faces must be discovered by parsing until the end of the file. + face_vertices = Tuple{SVector{3, ELTYPE}, SVector{3, ELTYPE}, SVector{3, ELTYPE}}[] + vertices = SVector{3, ELTYPE}[] + normals = SVector{3, ELTYPE}[] + # ASCII STL files may contain multiple "solid ... endsolid" patches. + # Collect each patch as a separate `TriangleMesh` and return a Vector. + # If the file contains only one patch we return the single `TriangleMesh`. + geometries = TriangleMesh[] + + load_data_ascii!(geometries, face_vertices, vertices, normals, io) + + length(geometries) == 1 && return first(geometries) + + return geometries +end + +function load_data_ascii!(geometries, + face_vertices::Vector{Tuple{SVector{3, T}, SVector{3, T}, + SVector{3, T}}}, + vertices, normals, io) where {T} + while !eof(io) + line = strip(readline(io)) + low = lowercase(line) + + if startswith(low, "endsolid") + # First patch finished + push!(geometries, + TriangleMesh(copy(face_vertices), copy(normals), copy(vertices))) + + empty!(face_vertices) + empty!(normals) + empty!(vertices) + end + + # look for facet normal ... + if startswith(low, "facet") + parts = split(line) + # Expect: `["facet", "normal", nx, ny, nz]` + @assert length(parts)>=5 "Unexpected facet normal line: $line" + n = SVector{3, T}(parse(T, parts[3]), parse(T, parts[4]), parse(T, parts[5])) + push!(normals, n) + + # Consume "outer loop" + while !eof(io) + lowercase(strip(readline(io))) == "outer loop" && break + end + + # Read three vertex lines + for _ in 1:3 + l = strip(readline(io)) + partv = split(l) + @assert lowercase(partv[1])=="vertex" "Unexpected vertex line: $l" + vertex = SVector{3, T}(parse(T, partv[2]), parse(T, partv[3]), + parse(T, partv[4])) + + push!(vertices, vertex) + end + + push!(face_vertices, (vertices[end - 2], vertices[end - 1], vertices[end])) + + # Consume "endloop" and "endfacet" + while !eof(io) + lowercase(strip(readline(io))) == "endfacet" && break + end + end + end + + return geometries +end + function load_data!(face_vertices::Vector{Tuple{SVector{3, T}, SVector{3, T}, SVector{3, T}}}, vertices, normals, io) where {T} diff --git a/src/preprocessing/geometries/triangle_mesh.jl b/src/preprocessing/geometries/triangle_mesh.jl index b11cc95ef4..a72492fd81 100644 --- a/src/preprocessing/geometries/triangle_mesh.jl +++ b/src/preprocessing/geometries/triangle_mesh.jl @@ -254,3 +254,155 @@ function volume(mesh::TriangleMesh) return volume end + +function extrude_geometry(geometry_bottom::TriangleMesh{NDIMS, ELTYPE}, + extrude_length::Real; omit_top_face=false, + omit_bottom_face=false) where {NDIMS, ELTYPE} + face_normal = normalize(sum(geometry_bottom.face_normals) / nfaces(geometry_bottom)) + + shift = face_normal * extrude_length + + # Bottom (original) geometry data + vertices_bottom = copy(geometry_bottom.vertices) + # Flip the winding of the bottom faces so that after extrusion the bottom faces' + # normals point outward from the resulting solid (consistent with the top faces). + # We reverse the triangle vertex order (v2 <-> v3) and negate the face normals to + # reflect that orientation change. + flipped_faces = copy(geometry_bottom.face_vertices) + face_vertices_bottom = [(v1, v3, v2) for (v1, v2, v3) in flipped_faces] + normals_bottom = .-copy(geometry_bottom.face_normals) + + # Top (shifted) geometry data. + # Shift every bottom vertex by the same vector to create the top vertices. + vertices_top = [v .+ shift for v in vertices_bottom] + normals_top = copy(geometry_bottom.face_normals) + # Shift each face-tuple component-wise to create top face tuples + face_vertices_top = [(v1 .+ shift, v2 .+ shift, v3 .+ shift) + for (v1, v2, v3) in flipped_faces] + + # Build a temporary `TriangleMesh` for the top layer to reuse its indexing helpers + geometry_top = TriangleMesh(face_vertices_top, normals_top, vertices_top) + + # Find boundary edges on bottom to generate side faces only for exterior edges. + directed_edges = zeros(Int, length(geometry_bottom.edge_normals)) + for face in eachindex(geometry_bottom.face_vertices) + v1 = geometry_bottom.face_vertices_ids[face][1] + v2 = geometry_bottom.face_vertices_ids[face][2] + v3 = geometry_bottom.face_vertices_ids[face][3] + + edge1 = geometry_bottom.face_edges_ids[face][1] + edge2 = geometry_bottom.face_edges_ids[face][2] + edge3 = geometry_bottom.face_edges_ids[face][3] + + if geometry_bottom.edge_vertices_ids[edge1] == (v1, v2) + directed_edges[edge1] += 1 + else + directed_edges[edge1] -= 1 + end + if geometry_bottom.edge_vertices_ids[edge2] == (v2, v3) + directed_edges[edge2] += 1 + else + directed_edges[edge2] -= 1 + end + if geometry_bottom.edge_vertices_ids[edge3] == (v3, v1) + directed_edges[edge3] += 1 + else + directed_edges[edge3] -= 1 + end + end + + boundary_edges = findall(!iszero, directed_edges) + + @inline function triangle_normal(v1, v2, v3) + n = cross(v2 - v1, v3 - v1) + norm(n) > 0 && return normalize(n) + + # Return zero for degenerate triangles + return zero(n) + end + + # Compute approximate geometric center (used to ensure face winding produces outward normals) + center = (reduce(+, vertices_bottom) / length(vertices_bottom)) .+ (shift / 2) + + faces_vertices_closing = Tuple{SVector{3, ELTYPE}, SVector{3, ELTYPE}, + SVector{3, ELTYPE}}[] + normals_closing = SVector{3, ELTYPE}[] + for edge in boundary_edges + # Bottom edge endpoints (in the bottom mesh vertex list) + v1b = geometry_bottom.vertices[geometry_bottom.edge_vertices_ids[edge][1]] + v2b = geometry_bottom.vertices[geometry_bottom.edge_vertices_ids[edge][2]] + + # Corresponding top edge endpoints (in the top mesh vertex list) + v1t = geometry_top.vertices[geometry_top.edge_vertices_ids[edge][1]] + v2t = geometry_top.vertices[geometry_top.edge_vertices_ids[edge][2]] + + # Split the quad (v1b, v2b, v2t, v1t) into two triangles + face_1 = (v1b, v2b, v2t) + face_2 = (v1b, v2t, v1t) + + n1 = triangle_normal(face_1[1], face_1[2], face_1[3]) + n2 = triangle_normal(face_2[1], face_2[2], face_2[3]) + + # Compute correct winding so normals point away from the body center + c1 = (face_1[1] + face_1[2] + face_1[3]) / 3 + if dot(n1, c1 - center) < 0 + face_1 = (face_1[1], face_1[3], face_1[2]) + n1 = -n1 + end + + c2 = (face_2[1] + face_2[2] + face_2[3]) / 3 + if dot(n2, c2 - center) < 0 + face_2 = (face_2[1], face_2[3], face_2[2]) + n2 = -n2 + end + + push!(faces_vertices_closing, face_1) + push!(normals_closing, n1) + push!(faces_vertices_closing, face_2) + push!(normals_closing, n2) + end + + vertices_new = vcat(vertices_bottom, vertices_top) + if omit_bottom_face && omit_top_face + # Only sides + face_vertices_new = vcat(faces_vertices_closing) + normals_new = vcat(normals_closing) + elseif omit_bottom_face + # top + sides + face_vertices_new = vcat(face_vertices_top, faces_vertices_closing) + normals_new = vcat(normals_top, normals_closing) + elseif omit_top_face + # bottom + sides + face_vertices_new = vcat(face_vertices_bottom, faces_vertices_closing) + normals_new = vcat(normals_bottom, normals_closing) + else + # bottom + top + sides (default) + face_vertices_new = vcat(face_vertices_bottom, face_vertices_top, + faces_vertices_closing) + normals_new = vcat(normals_bottom, normals_top, normals_closing) + end + + return TriangleMesh(face_vertices_new, normals_new, vertices_new) +end + +function Base.union(geometry::TriangleMesh, geometries::TriangleMesh...) + other_geometry = first(geometries) + + vertices_other = copy(other_geometry.vertices) + face_vertices_other = copy(other_geometry.face_vertices) + normals_other = copy(other_geometry.face_normals) + + vertices = copy(geometry.vertices) + face_vertices = copy(geometry.face_vertices) + normals = copy(geometry.face_normals) + + vertices_new = vcat(vertices, vertices_other) + face_vertices_new = vcat(face_vertices, face_vertices_other) + normals_new = vcat(normals, normals_other) + + result = TriangleMesh(face_vertices_new, normals_new, vertices_new) + + return union(result, Base.tail(geometries)...) +end + +Base.union(geometry::TriangleMesh) = geometry diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index a4fead937e..ea8598023a 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -147,7 +147,7 @@ bidirectional_flow = BoundaryZone(; boundary_face=face_vertices, face_normal, !!! warning "Experimental Implementation" This is an experimental feature and may change in any future releases. """ -struct BoundaryZone{IC, S, ZO, ZW, FD, FN, ELTYPE, R} +struct BoundaryZone{IC, S, ZO, ZW, FD, FN, ELTYPE, R, PM} initial_condition :: IC spanning_set :: S zone_origin :: ZO @@ -156,6 +156,7 @@ struct BoundaryZone{IC, S, ZO, ZW, FD, FN, ELTYPE, R} face_normal :: FN rest_pressure :: ELTYPE # Only required for `BoundaryModelDynamicalPressureZhang` reference_values :: R + pressure_model :: PM # Note that the following can't be static type parameters, as all boundary zones in a system # must have the same type, so that we can loop over them in a type-stable way. average_inflow_velocity :: Bool @@ -168,6 +169,7 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing, initial_condition=nothing, extrude_geometry=nothing, open_boundary_layers::Integer, average_inflow_velocity=true, boundary_type=BidirectionalFlow(), + pressure_model=nothing, rest_pressure=zero(eltype(density)), reference_density=nothing, reference_pressure=nothing, reference_velocity=nothing) @@ -246,7 +248,11 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing, coordinates_svector = reinterpret(reshape, SVector{NDIMS, ELTYPE}, ic.coordinates) if prescribed_pressure - ic.pressure .= pressure_ref.(coordinates_svector, 0) + if isnothing(pressure_model) + ic.pressure .= pressure_ref.(coordinates_svector, 0) + else + throw(ArgumentError("Setting prescribed pressure together with a pressure model is not supported.")) + end end if prescribed_density ic.density .= density_ref.(coordinates_svector, 0) @@ -258,8 +264,8 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing, return BoundaryZone(ic, spanning_set_, zone_origin, zone_width, flow_direction, face_normal_, rest_pressure, reference_values, - average_inflow_velocity, prescribed_density, prescribed_pressure, - prescribed_velocity) + pressure_model, average_inflow_velocity, prescribed_density, + prescribed_pressure, prescribed_velocity) end function boundary_type_name(boundary_zone::BoundaryZone) @@ -291,6 +297,10 @@ function Base.show(io::IO, ::MIME"text/plain", boundary_zone::BoundaryZone) summary_line(io, "boundary type", boundary_type_name(boundary_zone)) summary_line(io, "#particles", nparticles(boundary_zone.initial_condition)) summary_line(io, "width", round(boundary_zone.zone_width, digits=6)) + if !isnothing(boundary_zone.pressure_model) && + boundary_zone.pressure_model.is_prescribed + summary_line(io, "pressure model", type2string(boundary_zone.pressure_model)) + end summary_footer(io) end end @@ -494,7 +504,11 @@ function reference_pressure(boundary_zone, v, system, particle, pos, t) # `pressure_reference_values[zone_id](pos, t)`, but in a type-stable way return apply_ith_function(pressure_reference_values, zone_id, pos, t) else - return current_pressure(v, system, particle) + pressure = current_pressure(v, system, particle) + + # Return either the current pressure or a pressure computed by a pressure model + return imposed_pressure(system, system.pressure_model_values, boundary_zone, + pressure, particle) end end diff --git a/src/schemes/boundary/open_boundary/dynamical_pressure.jl b/src/schemes/boundary/open_boundary/dynamical_pressure.jl index e7cb837d90..70acf3873a 100644 --- a/src/schemes/boundary/open_boundary/dynamical_pressure.jl +++ b/src/schemes/boundary/open_boundary/dynamical_pressure.jl @@ -120,7 +120,9 @@ function reference_pressure(boundary_zone, v, # `pressure_reference_values[zone_id](pos, t)`, but in a type-stable way return apply_ith_function(pressure_reference_values, zone_id, pos, t) else - return rest_pressure + # Return either the rest pressure or a pressure computed by a pressure model + return imposed_pressure(system, system.pressure_model_values, boundary_zone, + rest_pressure, particle) end end diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 32209965e5..abcf3450cf 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -95,14 +95,14 @@ function update_boundary_quantities!(system, boundary_model::BoundaryModelMirror @threaded semi for particle in each_integrated_particle(system) boundary_zone = current_boundary_zone(system, particle) - (; prescribed_density, prescribed_pressure, prescribed_velocity) = boundary_zone + (; prescribed_density, prescribed_velocity) = boundary_zone particle_coords = current_coords(u, system, particle) - if prescribed_pressure - pressure[particle] = reference_pressure(boundary_zone, v, system, particle, - particle_coords, t) - end + # The pressure can be prescribed, left unprescribed (in which case the current + # pressure is returned), or determined by a boundary model. + pressure[particle] = reference_pressure(boundary_zone, v, system, particle, + particle_coords, t) if prescribed_density density[particle] = reference_density(boundary_zone, v, system, particle, diff --git a/src/schemes/boundary/open_boundary/open_boundary.jl b/src/schemes/boundary/open_boundary/open_boundary.jl index 8f5c770c94..5f83dfa735 100644 --- a/src/schemes/boundary/open_boundary/open_boundary.jl +++ b/src/schemes/boundary/open_boundary/open_boundary.jl @@ -4,3 +4,4 @@ include("method_of_characteristics.jl") include("system.jl") include("dynamical_pressure.jl") include("rhs.jl") +include("pressure_model.jl") diff --git a/src/schemes/boundary/open_boundary/pressure_model.jl b/src/schemes/boundary/open_boundary/pressure_model.jl new file mode 100644 index 0000000000..40b539eebd --- /dev/null +++ b/src/schemes/boundary/open_boundary/pressure_model.jl @@ -0,0 +1,138 @@ +""" + RCRWindkesselModel(; characteristic_resistance, peripheral_resistance, compliance) + +The `RCRWindkessel` model is a biomechanical lumped-parameter representation +that captures the relationship between pressure and flow in pulsatile systems, e.g., vascular systems. +It is derived from an electrical circuit analogy and consists of three elements: + +- characteristic resistance (``R_1``): Represents the proximal resistance at the vessel entrance. + It models the immediate pressure drop that arises at the entrance of a vessel segment, + due either to a geometric narrowing or to a mismatch in characteristic impedance between adjacent segments. + A larger ``R_1`` produces a sharper initial pressure rise at the onset of flow. +- peripheral resistance (``R_2``): Represents the distal resistance, + which controls the sustained outflow into the peripheral circulation and thereby determines the level of the mean pressure. + A high ``R_2`` maintains a higher pressure (reduced outflow), whereas a low ``R_2`` allows a faster pressure decay. +- compliance (``C``): Connected in parallel with ``R_2`` and represents the capacity of elastic walls + to store and release volume; in other words, it models the "stretchiness" of walls. + Analogous to a capacitor in an electrical circuit, it absorbs blood when pressure rises and releases it during diastole. + The presence of ``C`` smooths pulsatile flow and produces a more uniform outflow profile. + + Lumped-parameter models for the vascular system are well described in the literature (e.g. [Westerhof2008](@cite)). + A practical step-by-step procedure for identifying the corresponding model parameters is provided by [Gasser2021](@cite). + +# Keywords +- `characteristic_resistance`: characteristic resistance (``R_1``) +- `peripheral_resistance`: peripheral resistance (``R_2``) +- `compliance`: compliance (``C``) +""" +struct RCRWindkesselModel{ELTYPE <: Real} + characteristic_resistance :: ELTYPE + peripheral_resistance :: ELTYPE + compliance :: ELTYPE + is_prescribed :: Bool +end + +function RCRWindkesselModel(; characteristic_resistance, peripheral_resistance, compliance) + return RCRWindkesselModel(characteristic_resistance, peripheral_resistance, compliance, + true) +end + +function Base.show(io::IO, ::MIME"text/plain", pressure_model::RCRWindkesselModel) + @nospecialize pressure_model # reduce precompilation time + + if get(io, :compact, false) + show(io, pressure_model) + else + summary_header(io, "RCRWindkesselModel") + summary_line(io, "characteristic_resistance", + pressure_model.characteristic_resistance) + summary_line(io, "peripheral_resistance", + pressure_model.peripheral_resistance) + summary_line(io, "compliance", pressure_model.compliance) + summary_footer(io) + end +end + +function update_pressure_model!(system::OpenBoundarySystem, v, u, semi, dt) + isnothing(system.pressure_model_values) && return system + + calculate_flow_rate_and_pressure!(system, v, u, dt) + + return system +end + +function calculate_flow_rate_and_pressure!(system, v, u, dt) + for (zone_id, boundary_zone) in enumerate(system.boundary_zones) + if boundary_zone.pressure_model.is_prescribed + calculate_flow_rate_and_pressure!(boundary_zone.pressure_model, system, + boundary_zone, zone_id, v, u, dt) + end + end + + return system +end + +function calculate_flow_rate_and_pressure!(pressure_model, system, boundary_zone, + zone_id, v, u, dt) + dt < sqrt(eps()) && return pressure_model + (; particle_spacing) = system.initial_condition + (; characteristic_resistance, peripheral_resistance, compliance) = pressure_model + (; flow_rate, pressure) = system.pressure_model_values[zone_id] + (; face_normal, zone_origin) = boundary_zone + + # Use a thin slice for the flow rate calculation + slice = particle_spacing * 125 / 100 + + # Find particles within a thin slice near the boundary face for flow rate computation + candidates = findall(x -> dot(x - zone_origin, -face_normal) <= slice, + reinterpret(reshape, SVector{ndims(system), eltype(u)}, u)) + + particles_in_zone = findall(particle -> boundary_zone == + current_boundary_zone(system, particle), + each_integrated_particle(system)) + + intersect!(candidates, particles_in_zone) + + cross_sectional_area = length(candidates) * particle_spacing^(ndims(system) - 1) + + # Division inside the `sum` closure to maintain GPU compatibility + velocity_avg = sum(candidates) do particle + return dot(current_velocity(v, system, particle), -face_normal) / length(candidates) + end + + # Compute volumetric flow rate: Q = v * A + current_flow_rate = velocity_avg * cross_sectional_area + + previous_pressure = pressure[] + previous_flow_rate = flow_rate[] + flow_rate[] = current_flow_rate + + # Calculate new pressure according to eq. 22 in Zhang et al. (2025) + R1 = characteristic_resistance + R2 = peripheral_resistance + C = compliance + + term_1 = (1 + R1 / R2) * flow_rate[] + term_2 = C * R1 * (flow_rate[] - previous_flow_rate) / dt + term_3 = C * previous_pressure / dt + divisor = C / dt + 1 / R2 + + pressure_new = (term_1 + term_2 + term_3) / divisor + + pressure[] = pressure_new + + return system +end + +function imposed_pressure(system, pressure_model_values, boundary_zone, + pressure, particle) + boundary_zone.pressure_model.is_prescribed || return pressure + + zone_id = system.boundary_zone_indices[particle] + return pressure_model_values[zone_id].pressure[] +end + +function imposed_pressure(system, pressure_model_values::Nothing, boundary_zone, + pressure, particle) + return pressure +end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 2e8222faed..8bd43ada1a 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -23,7 +23,7 @@ Open boundary system for in- and outflow particles. This is an experimental feature and may change in any future releases. """ struct OpenBoundarySystem{BM, ELTYPE, NDIMS, IC, FS, FSI, K, ARRAY1D, BC, FC, BZI, BZ, - B, PF, ST, C} <: AbstractSystem{NDIMS} + B, PF, ST, PMV, C} <: AbstractSystem{NDIMS} boundary_model :: BM initial_condition :: IC fluid_system :: FS @@ -39,6 +39,7 @@ struct OpenBoundarySystem{BM, ELTYPE, NDIMS, IC, FS, FSI, K, ARRAY1D, BC, FC, BZ buffer :: B pressure_acceleration_formulation :: PF shifting_technique :: ST + pressure_model_values :: PMV cache :: C end @@ -46,19 +47,21 @@ function OpenBoundarySystem(boundary_model, initial_condition, fluid_system, fluid_system_index, smoothing_kernel, smoothing_length, mass, volume, boundary_candidates, fluid_candidates, boundary_zone_indices, boundary_zone, buffer, - pressure_acceleration, shifting_technique, cache) + pressure_acceleration, shifting_technique, + pressure_model_values, cache) OpenBoundarySystem{typeof(boundary_model), eltype(mass), ndims(initial_condition), typeof(initial_condition), typeof(fluid_system), typeof(fluid_system_index), typeof(smoothing_kernel), typeof(mass), typeof(boundary_candidates), typeof(fluid_candidates), typeof(boundary_zone_indices), typeof(boundary_zone), typeof(buffer), typeof(pressure_acceleration), typeof(shifting_technique), + typeof(pressure_model_values), typeof(cache)}(boundary_model, initial_condition, fluid_system, fluid_system_index, smoothing_kernel, smoothing_length, mass, volume, boundary_candidates, fluid_candidates, boundary_zone_indices, boundary_zone, buffer, pressure_acceleration, - shifting_technique, cache) + shifting_technique, pressure_model_values, cache) end function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; @@ -72,9 +75,8 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; BoundaryModelDynamicalPressureZhang ? shifting_technique(fluid_system) : nothing) boundary_zones_ = filter(bz -> !isnothing(bz), boundary_zones) - reference_values_ = map(bz -> bz.reference_values, boundary_zones_) - initial_conditions = union((bz.initial_condition for bz in boundary_zones)...) + initial_conditions = union((bz.initial_condition for bz in boundary_zones_)...) buffer = SystemBuffer(nparticles(initial_conditions), buffer_size) @@ -86,7 +88,7 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; cache = (; create_cache_shifting(initial_conditions, shifting_technique)..., create_cache_open_boundary(boundary_model, fluid_system, initial_conditions, - reference_values_)...) + boundary_zones_)...) fluid_system_index = Ref(0) @@ -98,28 +100,32 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; boundary_zone_indices = zeros(Int, nparticles(initial_conditions)) + pressure_models, pressure_model_values = extract_pressure_models(boundary_zones_) + # Create new `BoundaryZone`s with `reference_values` set to `nothing` for type stability. # `reference_values` are only used as API feature to temporarily store the reference values # in the `BoundaryZone`, but they are not used in the actual simulation. - boundary_zones_new = map(zone -> BoundaryZone(zone.initial_condition, - zone.spanning_set, - zone.zone_origin, - zone.zone_width, - zone.flow_direction, - zone.face_normal, - zone.rest_pressure, - nothing, - zone.average_inflow_velocity, - zone.prescribed_density, - zone.prescribed_pressure, - zone.prescribed_velocity), - boundary_zones) + boundary_zones_new = (map((zone, + pressure_model) -> BoundaryZone(zone.initial_condition, + zone.spanning_set, + zone.zone_origin, + zone.zone_width, + zone.flow_direction, + zone.face_normal, + zone.rest_pressure, + nothing, pressure_model, + zone.average_inflow_velocity, + zone.prescribed_density, + zone.prescribed_pressure, + zone.prescribed_velocity), + boundary_zones_, pressure_models)...,) return OpenBoundarySystem(boundary_model, initial_conditions, fluid_system, fluid_system_index, smoothing_kernel, smoothing_length, mass, volume, boundary_candidates, fluid_candidates, boundary_zone_indices, boundary_zones_new, buffer, - pressure_acceleration, shifting_technique, cache) + pressure_acceleration, shifting_technique, + pressure_model_values, cache) end function initialize!(system::OpenBoundarySystem, semi) @@ -130,8 +136,9 @@ function initialize!(system::OpenBoundarySystem, semi) return system end -function create_cache_open_boundary(boundary_model, fluid_system, - initial_condition, reference_values) +function create_cache_open_boundary(boundary_model, fluid_system, initial_condition, + boundary_zones) + reference_values = map(bz -> bz.reference_values, boundary_zones) ELTYPE = eltype(initial_condition) # Separate `reference_values` into pressure, density and velocity reference values @@ -191,6 +198,34 @@ function create_cache_open_boundary(boundary_model, fluid_system, end end +function extract_pressure_models(boundary_zones) + zone_indices = findall(zone -> !isnothing(zone.pressure_model), boundary_zones) + if isempty(zone_indices) + pressure_models = fill(nothing, length(boundary_zones)) + pressure_model_values = nothing + + else + # Build a vector of pressure models, using a dummy instance for boundary zones + # that lack a pressure model to maintain type stability. + # This vector is then passed to the new boundary zone instances. + zero_ = zero(eltype(first(boundary_zones).initial_condition)) + dummy_model = RCRWindkesselModel(zero_, zero_, zero_, false) + pressure_models = fill(dummy_model, length(boundary_zones)) + pressure_models[zone_indices] .= map(zone -> RCRWindkesselModel(zone.pressure_model.characteristic_resistance, + zone.pressure_model.peripheral_resistance, + zone.pressure_model.compliance, + true), + boundary_zones[zone_indices]) + + # Store values that are updated during the simulation. + # These must be kept separately to ensure type stability of the boundary zone. + pressure_model_values = map(bz -> (pressure=Ref(bz.rest_pressure), + flow_rate=Ref(zero_)), boundary_zones) + end + + return pressure_models, pressure_model_values +end + timer_name(::OpenBoundarySystem) = "open_boundary" vtkname(system::OpenBoundarySystem) = "open_boundary" @@ -217,6 +252,12 @@ function Base.show(io::IO, ::MIME"text/plain", system::OpenBoundarySystem) summary_line(io, "density diffusion", density_diffusion(system)) summary_line(io, "shifting technique", shifting_technique(system)) end + if !isnothing(system.pressure_model_values) + summary_line(io, "pressure model", + type2string(first(system.boundary_zones).pressure_model) * + " (in boundary zones: " * + "$(findall(zone -> zone.pressure_model.is_prescribed, system.boundary_zones)))") + end summary_footer(io) end end @@ -300,6 +341,8 @@ function update_open_boundary_eachstep!(system::OpenBoundarySystem, v_ode, u_ode @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) + update_pressure_model!(system, v, u, semi, integrator.dt) + # Update density, pressure and velocity based on the specific boundary model @trixi_timeit timer() "update boundary quantities" begin update_boundary_quantities!(system, boundary_model, v, u, v_ode, u_ode, semi, t) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index c9ecb1ca82..15efb91fa0 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -15,6 +15,16 @@ all relevant quantities and operators are measured with respect to the initial configuration (O’Connor & Rogers 2021, Belytschko et al. 2000). See [Total Lagrangian SPH](@ref tlsph) for more details on the method. +There are two approaches to specify clamped particles in the system: + +1. **Manual ordering**: Provide the number of clamped particles via `n_clamped_particles`, + ensuring that these particles are the last entries in the `InitialCondition` + (see the info box below for details). + +2. **Automatic ordering**: Pass the indices of the particles to be clamped via `clamped_particles`. + In this case, the specified particles will be internally reordered so that they become + the last particles in the `InitialCondition`. + # Arguments - `initial_condition`: Initial condition representing the system's particles. - `young_modulus`: Young's modulus. diff --git a/src/visualization/recipes_plots.jl b/src/visualization/recipes_plots.jl index ad70124b44..d562fbdf5e 100644 --- a/src/visualization/recipes_plots.jl +++ b/src/visualization/recipes_plots.jl @@ -170,3 +170,41 @@ RecipesBase.@recipe function f(geometry::Polygon) return (x, y) end end + +RecipesBase.@recipe function f(box::OrientedBoundingBox{2}) + # Get the box origin and spanning vectors + origin = box.box_origin + v1, v2 = box.spanning_vectors + + # Calculate all 4 corners of the rectangle + # Corner 1: origin + # Corner 2: origin + v1 + # Corner 3: origin + v1 + v2 (diagonal corner) + # Corner 4: origin + v2 + corners = hcat(origin, + origin + v1, + origin + v1 + v2, + origin + v2, + origin) # Close the rectangle by returning to start + + x = corners[1, :] + y = corners[2, :] + + aspect_ratio --> :equal + grid --> false + + # First plot the vertices as a scatter plot + @series begin + seriestype --> :scatter + return (x, y) + end + + # Now plot the edges as a line plot + @series begin + # Ignore series in legend and color cycling. Note that `:=` forces the attribute, + # whereas `-->` would only set it if it is not already set. + primary := false + + return (x, y) + end +end diff --git a/test/preprocessing/data/cuboid.stl b/test/preprocessing/data/cuboid.stl new file mode 100644 index 0000000000..dc5390fb07 --- /dev/null +++ b/test/preprocessing/data/cuboid.stl @@ -0,0 +1,101 @@ +solid front + facet normal 0 0 1 + outer loop + vertex 0 0 1 + vertex 1 0 1 + vertex 1 1 1 + endloop + endfacet + facet normal 0 0 1 + outer loop + vertex 0 0 1 + vertex 1 1 1 + vertex 0 1 1 + endloop + endfacet +endsolid front + +solid back + facet normal 0 0 -1 + outer loop + vertex 0 0 0 + vertex 1 1 0 + vertex 1 0 0 + endloop + endfacet + facet normal 0 0 -1 + outer loop + vertex 0 0 0 + vertex 0 1 0 + vertex 1 1 0 + endloop + endfacet +endsolid back + +solid left + facet normal -1 0 0 + outer loop + vertex 0 0 0 + vertex 0 0 1 + vertex 0 1 1 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 0 0 + vertex 0 1 1 + vertex 0 1 0 + endloop + endfacet +endsolid left + +solid right + facet normal 1 0 0 + outer loop + vertex 1 0 0 + vertex 1 1 1 + vertex 1 0 1 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 1 0 0 + vertex 1 1 0 + vertex 1 1 1 + endloop + endfacet +endsolid right + +solid top + facet normal 0 1 0 + outer loop + vertex 0 1 0 + vertex 0 1 1 + vertex 1 1 1 + endloop + endfacet + facet normal 0 1 0 + outer loop + vertex 0 1 0 + vertex 1 1 1 + vertex 1 1 0 + endloop + endfacet +endsolid top + +solid bottom + facet normal 0 -1 0 + outer loop + vertex 0 0 0 + vertex 1 0 1 + vertex 0 0 1 + endloop + endfacet + facet normal 0 -1 0 + outer loop + vertex 0 0 0 + vertex 1 0 0 + vertex 1 0 1 + endloop + endfacet +endsolid bottom diff --git a/test/preprocessing/geometries/geometries.jl b/test/preprocessing/geometries/geometries.jl index 4096f42027..08a46c343d 100644 --- a/test/preprocessing/geometries/geometries.jl +++ b/test/preprocessing/geometries/geometries.jl @@ -129,6 +129,123 @@ end end + @testset verbose=true "Multiple Patches" begin + file = pkgdir(TrixiParticles, "test", "preprocessing", "data") + geometries = load_geometry(joinpath(file, "cuboid.stl")) + + @test isa(geometries, Vector{TrixiParticles.TriangleMesh}) + @test length(geometries) == 6 # six solids in the file + total_faces = sum(TrixiParticles.nfaces(m) for m in geometries) + @test total_faces == 12 + + expected_normals = [ + SVector(0.0, 0.0, 1.0), # front + SVector(0.0, 0.0, -1.0), # back + SVector(-1.0, 0.0, 0.0), # left + SVector(1.0, 0.0, 0.0), # right + SVector(0.0, 1.0, 0.0), # top + SVector(0.0, -1.0, 0.0) # bottom + ] + + # Each patch should contain exactly two triangles and their normals should + # match the facet normals. + @testset "patch $i" for i in eachindex(geometries) + mesh = geometries[i] + @test TrixiParticles.nfaces(mesh) == 2 + @test all(isapprox.(mesh.face_normals, Ref(expected_normals[i]); + atol=1e-12)) + end + end + + @testset verbose=true "Union" begin + # Build a single geometry by uniting multiple STL patches (cuboid.stl contains separate solids). + # The union should produce a closed volume. + file = pkgdir(TrixiParticles, "test", "preprocessing", "data") + geometries = load_geometry(joinpath(file, "cuboid.stl")) + geometry_union = union(geometries...) + + # Employing `winding_number_factor = sqrt(eps())` serves as + # a quantitative criterion for geometric watertightness. + # Thus, we can check the resulting particle count with validated reference values. + winding_number_factor = sqrt(eps()) + ic = ComplexShape(geometry_union; particle_spacing=0.05, density=1.0, + point_in_geometry_algorithm=WindingNumberJacobson(; + geometry=geometry_union, + winding_number_factor)) + + @test nparticles(ic) == 9261 + end + + @testset verbose=true "Extrude Geometry" begin + file = pkgdir(TrixiParticles, "test", "preprocessing", "data") + geometry = load_geometry(joinpath(file, "inflow_geometry.stl")) + + @testset verbose=true "Watertightness" begin + geometry_extruded = extrude_geometry(geometry, 0.8) + + # Employing `winding_number_factor = sqrt(eps())` serves as + # a quantitative criterion for geometric watertightness. + # Thus, we can check the resulting particle count with validated reference values. + winding_number_factor = sqrt(eps()) + ic = ComplexShape(geometry_extruded; particle_spacing=0.03, density=1.0, + point_in_geometry_algorithm=WindingNumberJacobson(; + geometry=geometry_extruded, + winding_number_factor)) + + @test nparticles(ic) == 2692 + end + + @testset verbose=true "Omitting Top/Bottom Face" begin + geometry_extruded_1 = extrude_geometry(geometry, 0.8; omit_top_face=true) + geometry_extruded_2 = extrude_geometry(geometry, 0.8; omit_bottom_face=true) + geometry_extruded_3 = extrude_geometry(geometry, 0.8; omit_top_face=true, + omit_bottom_face=true) + + winding_number_factor = 0.2 + @testset verbose=true "Omit Top Face" begin + expected_min_corner = [-0.036399998962879196; 0.24624998748302457; -0.5233639197487431;;] + expected_max_corner = [0.38360000103712083; 1.1462499874830245; -0.07336391974874301;;] + + ic_1 = ComplexShape(geometry_extruded_1; particle_spacing=0.03, density=1.0, + point_in_geometry_algorithm=WindingNumberJacobson(; + geometry=geometry_extruded_1, + winding_number_factor)) + + @test nparticles(ic_1) == 2994 + @test isapprox(maximum(ic_1.coordinates, dims=2), expected_max_corner) + @test isapprox(minimum(ic_1.coordinates, dims=2), expected_min_corner) + end + + @testset verbose=true "Omit Bottom Face" begin + expected_min_corner = [-0.0663999989628792; 0.1562499874830246; -0.49336391974874305;;] + expected_max_corner = [0.38360000103712083; 1.0562499874830245; -0.07336391974874301;;] + + ic_2 = ComplexShape(geometry_extruded_2; particle_spacing=0.03, density=1.0, + point_in_geometry_algorithm=WindingNumberJacobson(; + geometry=geometry_extruded_2, + winding_number_factor)) + + @test nparticles(ic_2) == 2988 + @test isapprox(maximum(ic_2.coordinates, dims=2), expected_max_corner) + @test isapprox(minimum(ic_2.coordinates, dims=2), expected_min_corner) + end + + @testset verbose=true "Omit Both" begin + expected_min_corner = [-0.0663999989628792; 0.1562499874830246; -0.5233639197487431;;] + expected_max_corner = [0.38360000103712083; 1.1462499874830245; -0.07336391974874301;;] + + ic_3 = ComplexShape(geometry_extruded_3; particle_spacing=0.03, density=1.0, + point_in_geometry_algorithm=WindingNumberJacobson(; + geometry=geometry_extruded_3, + winding_number_factor)) + + @test nparticles(ic_3) == 3258 + @test isapprox(maximum(ic_3.coordinates, dims=2), expected_max_corner) + @test isapprox(minimum(ic_3.coordinates, dims=2), expected_min_corner) + end + end + end + @testset verbose=true "Boundary Face" begin file = pkgdir(TrixiParticles, "test", "preprocessing", "data") planar_geometry = load_geometry(joinpath(file, "inflow_geometry.stl")) @@ -190,4 +307,250 @@ @test TrixiParticles.unique_sorted(copy(x)) == sort(unique(x)) end + + @testset verbose=true "OrientedBoundingBox" begin + @testset verbose=true "2D" begin + @testset verbose=true "Manual Construction" begin + # Create a 2D oriented bounding box and test its spanning vectors + orientation_matrix = [cos(pi / 4) -sin(pi / 4); + sin(pi / 4) cos(pi / 4)] + box_2d = OrientedBoundingBox(; box_origin=[0.1, -2.0], + orientation_matrix, + edge_lengths=(2.0, 1.0)) + expected = [0.1+sqrt(2) 0.1-sqrt(2) / 2; + -2.0+sqrt(2) -2.0+sqrt(2) / 2] + @test isapprox(stack(box_2d.spanning_vectors) .+ box_2d.box_origin, + expected) + end + + @testset verbose=true "From Geometry" begin + # Load 2D geometry from file and test oriented bounding box calculation + data_dir = pkgdir(TrixiParticles, "examples", "preprocessing", "data") + geometry = load_geometry(joinpath(data_dir, "potato.asc")) + + box_1 = OrientedBoundingBox(geometry) + expected = [1.925952836831006 0.07595980222547905; + -1.026422999553135 0.14290580311134243] + @test isapprox(stack(box_1.spanning_vectors) .+ box_1.box_origin, expected) + + # Test oriented bounding box with custom local axis scale + box_2 = OrientedBoundingBox(geometry, local_axis_scale=(1.5, 2)) + expected = [2.8819385332081096 -0.8800258941516246; + -1.3869998124508005 0.503482616009008] + @test isapprox(stack(box_2.spanning_vectors) .+ box_2.box_origin, expected) + + box_3 = OrientedBoundingBox(geometry, local_axis_scale=(3, 0.5)) + expected = [1.3085086828079229 0.6934039562485623; + -1.8545287410598816 0.9710115446180893] + @test isapprox(stack(box_3.spanning_vectors) .+ box_3.box_origin, expected) + end + + @testset verbose=true "From Point-Cloud" begin + # Create a point cloud from a spherical shape (quarter circle sector) + shape = SphereShape(0.1, 1.5, (0.2, 0.4), 1.0, n_layers=4, + sphere_type=RoundSphere(; start_angle=0, + end_angle=pi / 4)) + + box = OrientedBoundingBox(shape.coordinates) + expected = [1.3939339828220176 1.3264241636171004; + 0.29393398282201755 1.4671898309959612] + @test isapprox(stack(box.spanning_vectors) .+ box.box_origin, expected) + end + end + + @testset verbose=true "3D" begin + @testset verbose=true "Manual Construction" begin + # Create a 3D oriented bounding box and test its spanning vectors + orientation_matrix = [1/sqrt(2) -1/sqrt(2) 0.0; + 1/sqrt(2) 1/sqrt(2) 0.0; + 0.0 0.0 1.0] + box_3d = OrientedBoundingBox(; box_origin=[0.5, -0.2, 0.0], + orientation_matrix, + edge_lengths=(1.0, 2.0, 3.0)) + expected = [0.5+1 / sqrt(2) 0.5-2 / sqrt(2) 0.5; + -0.2+1 / sqrt(2) -0.2+2 / sqrt(2) -0.2; + 0.0 0.0 3.0] + @test isapprox(stack(box_3d.spanning_vectors) .+ box_3d.box_origin, + expected) + end + + @testset verbose=true "From Geometry" begin + # Load 3D geometry from file and test oriented bounding box calculation + file = pkgdir(TrixiParticles, "test", "preprocessing", "data") + geometry = load_geometry(joinpath(file, "inflow.stl")) + box_1 = OrientedBoundingBox(geometry) + expected = [-0.1095304728105212 -0.03438646676380468 0.2379628221247777; + 0.16873336880727502 0.2999997960771134 0.20335870225239414; + -0.36117732916596945 -0.04334770439536223 -0.42937034485453873] + @test isapprox(stack(box_1.spanning_vectors) .+ box_1.box_origin, expected) + + # Test oriented bounding box with custom local axis scale + box_2 = OrientedBoundingBox(geometry, local_axis_scale=(1.5, 1.5, 0.5)) + expected = [-0.044939194583412626 0.06777681448666216 0.1427916288643849; + 0.12137652133720928 0.3182761622419668 0.23147548411531255; + -0.4543422598660782 0.022402177289832625 -0.5018016853691041] + @test isapprox(stack(box_2.spanning_vectors) .+ box_2.box_origin, expected) + end + end + + @testset verbose=true "`is_in_oriented_box`" begin + @testset verbose=true "Manual Construction 2D" begin + # Simple axis-aligned box + orientation_matrix = [1.0 0.0; + 0.0 1.0] + box_2d = OrientedBoundingBox(; box_origin=[0.0, 0.0], + orientation_matrix, + edge_lengths=(2.0, 1.0)) + # Test points + query_points = Dict( + "Inside center" => ([1.0, 0.5], [1]), + "On corner" => ([0.0, 0.0], [1]), + "On edge" => ([1.0, 0.0], [1]), + "Just outside left" => ([-eps(), 0.5], []), + "Just outside right" => ([2.0 + eps(2.0), 0.5], []), + "Just outside bottom" => ([1.0, -eps()], []), + "Just outside top" => ([1.0, 1.0 + eps()], []), + "Far outside" => ([-5.0, -5.0], []) + ) + @testset "$k" for k in keys(query_points) + (point, expected) = query_points[k] + @test expected == TrixiParticles.is_in_oriented_box(point, box_2d) + end + end + + @testset verbose=true "Rotated Box 2D" begin + # 45° rotated box + orientation_matrix = [1/sqrt(2) -1/sqrt(2); + 1/sqrt(2) 1/sqrt(2)] + box_rotated = OrientedBoundingBox(; box_origin=[0.0, 0.0], + orientation_matrix, + edge_lengths=(sqrt(2), 1.0)) + v1_normalized = [1 / sqrt(2), 1 / sqrt(2)] # First edge direction + v2_normalized = [-1 / sqrt(2), 1 / sqrt(2)] # Second edge direction (perpendicular) + + # Test points + query_points = Dict( + "Inside center" => ([0.5, 0.5], [1]), + "Origin" => ([0.0, 0.0], [1]), + "Just outside along v1" => ([0.0, 0.0] + + (sqrt(2) + eps()) * v1_normalized, []), + "Just outside along v2 positive" => ([0.0, 0.0] + + (1.0 + eps()) * v2_normalized, []), + "Just outside along v2 negative" => ([0.0, 0.0] - + eps() * v2_normalized, []), + "Just inside along v1" => ([0.0, 0.0] + sqrt(2) * v1_normalized, [1]), + "Just inside along v1 negative" => ([0.0, 0.0] + + eps() * v1_normalized, [1]), + "Just inside along v2 positive" => ([0.0, 0.0] + v2_normalized, [1]), + "Just inside along v2 negative" => ([0.0, 0.0] + + eps() * v2_normalized, [1]), + "Far outside rotated" => ([2.0, 2.0], []) + ) + + @testset "$k" for k in keys(query_points) + (point, expected) = query_points[k] + @test expected == TrixiParticles.is_in_oriented_box(point, box_rotated) + end + end + + @testset verbose=true "Manual Construction 3D" begin + # Simple axis-aligned box + box_3d = OrientedBoundingBox(; box_origin=[0.0, 0.0, 0.0], + orientation_matrix=I(3), + edge_lengths=(2.0, 1.0, 0.5)) + + # Test points + query_points = Dict( + "Inside center" => ([1.0, 0.5, 0.25], [1]), + "On corner origin" => ([0.0, 0.0, 0.0], [1]), + "On corner opposite" => ([2.0, 1.0, 0.5], [1]), + "On face" => ([1.0, 0.0, 0.25], [1]), + "Just outside x-negative" => ([-eps(), 0.5, 0.25], []), + "Just outside x-positive" => ([2.0 + eps(2.0), 0.5, 0.25], []), + "Just outside y-negative" => ([1.0, -eps(), 0.25], []), + "Just outside y-positive" => ([1.0, 1.0 + eps(), 0.25], []), + "Just outside z-negative" => ([1.0, 0.5, -eps()], []), + "Just outside z-positive" => ([1.0, 0.5, 0.5 + eps()], []), + "Far outside" => ([-5.0, -5.0, -5.0], []) + ) + + @testset "$k" for k in keys(query_points) + (point, expected) = query_points[k] + @test expected == TrixiParticles.is_in_oriented_box(point, box_3d) + end + end + + @testset verbose=true "Rotated Box 3D" begin + # Box oriented along space diagonal + orientation_matrix = [1/sqrt(3) -1/sqrt(2) -1/sqrt(6); + 1/sqrt(3) 1/sqrt(2) -1/sqrt(6); + 1/sqrt(3) 0.0 2/sqrt(6)] + box_rotated = OrientedBoundingBox(; box_origin=[0.0, 0.0, 0.0], + orientation_matrix, + edge_lengths=(2.0, 1.0, 0.5)) + + # Test points + query_points = Dict( + "Origin" => ([0.0, 0.0, 0.0], [1]), + "Inside diagonal center" => ([0.0, 0.0, 0.0] + + orientation_matrix * [1.0, 0.5, 0.25], + [1]), + "Just outside along diagonal" => ([0.0, 0.0, 0.0] + + orientation_matrix * + ([2.0, 1.0, 0.5] .+ eps()), []), + "Just inside along diagonal" => ([0.0, 0.0, 0.0] + + orientation_matrix * [2.0, 1.0, 0.5], + [1]), + "Far outside diagonal" => ([0.0, 0.0, 0.0] + + orientation_matrix * [3.0, 2.0, 2.0], []) + ) + + @testset "$k" for k in keys(query_points) + (point, expected) = query_points[k] + @test expected == TrixiParticles.is_in_oriented_box(point, box_rotated) + end + end + + @testset verbose=true "Set Operations" begin + shape = RectangularShape(0.1, (10, 10), (0.0, 0.0), density=1.0) + box_1 = OrientedBoundingBox(box_origin=[0.0, 0.0], + orientation_matrix=I(2), + edge_lengths=(1.0, 1.0)) + box_2 = OrientedBoundingBox(box_origin=[0.0, 0.0], + orientation_matrix=I(2), + edge_lengths=(1.0, 0.5)) + + @test nparticles(intersect(shape, box_1)) == 100 + @test nparticles(setdiff(shape, box_1)) == 0 + @test nparticles(intersect(shape, box_2)) == 50 + @test nparticles(setdiff(shape, box_2)) == 50 + end + end + + @testset verbose=true "show" begin + box_2d = OrientedBoundingBox(box_origin=[0.1, -2.0], + orientation_matrix=I(2), + edge_lengths=(2.0, 1.0)) + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ OrientedBoundingBox (2D) │ + │ ════════════════════════ │ + │ box origin: ………………………………………………… [0.1, -2.0] │ + │ edge lengths: …………………………………………… (2.0, 1.0) │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", box_2d) == show_box + + box_3d = OrientedBoundingBox(box_origin=[0.5, -0.2, 0.0], + orientation_matrix=I(3), + edge_lengths=(1.0, 2.0, 3.0)) + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ OrientedBoundingBox (3D) │ + │ ════════════════════════ │ + │ box origin: ………………………………………………… [0.5, -0.2, 0.0] │ + │ edge lengths: …………………………………………… (1.0, 2.0, 3.0) │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", box_3d) == show_box + end + end end diff --git a/test/schemes/boundary/open_boundary/open_boundary.jl b/test/schemes/boundary/open_boundary/open_boundary.jl index 981facc427..78c73345e4 100644 --- a/test/schemes/boundary/open_boundary/open_boundary.jl +++ b/test/schemes/boundary/open_boundary/open_boundary.jl @@ -1,3 +1,4 @@ include("characteristic_variables.jl") include("mirroring.jl") include("boundary_zone.jl") +include("pressure_model.jl") diff --git a/test/schemes/boundary/open_boundary/pressure_model.jl b/test/schemes/boundary/open_boundary/pressure_model.jl new file mode 100644 index 0000000000..90bbac5825 --- /dev/null +++ b/test/schemes/boundary/open_boundary/pressure_model.jl @@ -0,0 +1,158 @@ +@testset verbose=true "`RCRWindkesselModel`" begin + @testset verbose=true "Show" begin + pressure_model = RCRWindkesselModel(; peripheral_resistance=1.2, + compliance=0.5, + characteristic_resistance=2.3) + + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ RCRWindkesselModel │ + │ ══════════════════ │ + │ characteristic_resistance: ………… 2.3 │ + │ peripheral_resistance: …………………… 1.2 │ + │ compliance: ………………………………………………… 0.5 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + + @test repr("text/plain", pressure_model) == show_box + end + + @testset verbose=true "Validation" begin + # The following test example is adapted from a case study presented + # in Gasser (2021, https://link.springer.com/book/10.1007/978-3-030-70966-2), + # which uses data from Burattini et al. (2002, https://doi.org/10.1152/ajpheart.2002.282.1.h244) + # to simulate a ferret vascular system. + function pulsatile_flow(t) + c_q = [ + 0.0136667 - 0.00144338im, # k = -10 + 0.00722562 - 0.0347752im, # k = -9 + -0.0593984 + 0.0217432im, # k = -8 + -0.0233298 + 0.0691505im, # k = -7 + 0.0250477 + 0.0231058im, # k = -6 + 0.0369504 + 0.0725im, # k = -5 + 0.216635 - 0.0557698im, # k = -4 + -0.0896779 - 0.408525im, # k = -3 + -0.555816 + 0.112961im, # k = -2 + 0.0276938 + 0.646413im, # k = -1 + 0.786333, # k = 0 + 0.0276938 - 0.646413im, # k = 1 + -0.555816 - 0.112961im, # k = 2 + -0.0896779 + 0.408525im, # k = 3 + 0.216635 + 0.0557698im, # k = 4 + 0.0369504 - 0.0725im, # k = 5 + 0.0250477 - 0.0231058im, # k = 6 + -0.0233298 - 0.0691505im, # k = 7 + -0.0593984 - 0.0217432im, # k = 8 + 0.00722562 + 0.0347752im, # k = 9 + 0.0136667 + 0.00144338im # k = 10 + ] + + signal = 0.0 + for (k, c) in zip(-10:10, c_q) + signal += real(c * exp(1im * k * 2 * pi / T * t)) + end + + return signal + end + + # Mock fluid system + struct FluidSystemMockRCR <: TrixiParticles.AbstractFluidSystem{2} + pressure_acceleration_formulation::Nothing + density_diffusion::Nothing + end + TrixiParticles.initial_smoothing_length(system::FluidSystemMockRCR) = 1.0 + TrixiParticles.nparticles(system::FluidSystemMockRCR) = 1 + TrixiParticles.system_smoothing_kernel(system::FluidSystemMockRCR) = nothing + + # Encapsulated simulation updating only `calculate_flow_rate_and_pressure!` + function simulate_rcr(R1, R2, C, tspan, dt, p_0, func) + pressure_model = RCRWindkesselModel(; peripheral_resistance=R2, + compliance=C, + characteristic_resistance=R1) + # Define a boundary zone with height=1.0 to ensure a unit volume, + # so velocity directly corresponds to flow rate. + reference_velocity(pos, t) = SVector(func(t), 0.0) + boundary_zone = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]), + particle_spacing=0.1, face_normal=(-1.0, 0.0), + density=1000.0, reference_velocity, pressure_model, + open_boundary_layers=1, rest_pressure=p_0) + + system = OpenBoundarySystem(boundary_zone; buffer_size=0, + boundary_model=nothing, + fluid_system=FluidSystemMockRCR(nothing, nothing)) + system.boundary_zone_indices .= 1 + + u = system.initial_condition.coordinates + v = system.initial_condition.velocity + + times = collect(tspan[1]:dt:tspan[2]) + p_calculated = empty(times) + for t in times + v[1, :] .= func(t) + TrixiParticles.calculate_flow_rate_and_pressure!(system, v, u, dt) + + # Store only values after the seventh cycle + if t >= 7T + push!(p_calculated, system.pressure_model_values[1].pressure[]) + end + end + + return p_calculated + end + + T = 0.375 + dt = T / 500 + tspan = (0.0, 8T) + p0 = 78.0 + + R1 = 1.7714 + R2 = 106.66 + C = 1.1808e-2 + + # The reference pressure values are computed using an ODE that describes + # the behavior of the RCR Windkessel model. The governing equation is: + # + # dp/dt + p / (R_1 * C) = R_2 * dq/dt + (R_1 + R_2) / (R_1 * C) * q + # + # where + # - p: pressure + # - q: time-dependent flow, provided by `pulsatile_flow` + # - R_1: characteristic resistance + # - R_2: peripheral resistance + # - C: compliance + # + # The function `pressure_RCR_ode!` implements this equation for numerical solution: + # function pressure_RCR_ode!(dp, p, params, t) + # (; R_1, R_2, C, q_func, dt) = params + # dq_dt = (q_func(t) - q_func(t - dt)) / dt # numerical derivative of inflow + # dp[1] = -p[1] / (R_2 * C) + R_1 * dq_dt + q_func(t) * (R_2 + R_1) / (R_2 * C) + # end + # + # The solution is obtained using simple Euler integration over the time interval `tspan`. + # Reference pressure values are evaluated at discrete time points (`validation_times`): + # params_RCR = (R_2=R2, R_1=R1, C=C, q_func=pulsatile_flow, dt=dt) + # sol_RCR = solve(ODEProblem(pressure_RCR_ode!, [p0], tspan, params_RCR), Euler(), + # dt=dt, adaptive=false) + # validation_times = collect(range(7T, 8T, step=50 * dt)) + # pressures_ref = vec(stack(sol_RCR(validation_times))) + # + # The values listed below are taken from this ODE simulation and serve as a baseline + # for validating the implementation of `RCRWindkesselModel` in the test. + pressures_ref = [ + 78.23098908459171, + 76.0357395466389, + 86.63098598789874, + 96.04303886358863, + 91.35701317358787, + 88.89163496483422, + 87.1572215816713, + 84.95271784314382, + 82.79365887019439, + 80.25102392189406, + 78.23943618409275 + ] + + pressures = simulate_rcr(R1, R2, C, tspan, dt, p0, pulsatile_flow) + + @test isapprox(pressures_ref, pressures[1:50:end], rtol=5e-3) + end +end