From 358a8e3c4ef4b8b8800f2e5a0ffa909bfc2c44d8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Elias=20B=C3=B6rjesson?= Date: Sun, 8 Dec 2024 10:35:09 +0100 Subject: [PATCH 1/3] Add Point cell --- src/Export/VTK.jl | 1 + src/Ferrite.jl | 1 + src/Grid/grid.jl | 9 +++++++++ src/exports.jl | 2 ++ src/interpolations.jl | 13 +++++++++++++ 5 files changed, 26 insertions(+) diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl index 840817fd4b..4ad5084cb9 100644 --- a/src/Export/VTK.jl +++ b/src/Export/VTK.jl @@ -65,6 +65,7 @@ function Base.setindex!(pvd::WriteVTK.CollectionFile, datfile::VTKGridFile, time return WriteVTK.collection_add_timestep(pvd, datfile, time) end +cell_to_vtkcell(::Type{Point}) = VTKCellTypes.VTK_VERTEX cell_to_vtkcell(::Type{Line}) = VTKCellTypes.VTK_LINE cell_to_vtkcell(::Type{QuadraticLine}) = VTKCellTypes.VTK_QUADRATIC_EDGE diff --git a/src/Ferrite.jl b/src/Ferrite.jl index 7b358463be..b9e28af746 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -45,6 +45,7 @@ abstract type AbstractRefShape{refdim} end # See src/docs.jl for detailed documentation struct RefHypercube{refdim} <: AbstractRefShape{refdim} end struct RefSimplex{refdim} <: AbstractRefShape{refdim} end +const RefPoint = RefHypercube{0} const RefLine = RefHypercube{1} const RefQuadrilateral = RefHypercube{2} const RefHexahedron = RefHypercube{3} diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index 3d27431bd1..7cadbee883 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -193,6 +193,11 @@ function faces(c::AbstractCell{RefShape}) where {RefShape} end end +# RefPoint (refdim = 0) +reference_vertices(::Type{RefPoint}) = (1,) +reference_edges(::Type{RefPoint}) = () # - +reference_faces(::Type{RefPoint}) = () # - + # RefLine (refdim = 1) reference_vertices(::Type{RefLine}) = (1, 2) reference_edges(::Type{RefLine}) = ((1, 2),) # e1 @@ -263,6 +268,9 @@ end ###################################################### # Lagrange interpolation based cells +struct Point <: AbstractCell{RefPoint} + nodes::NTuple{1, Int} +end struct Line <: AbstractCell{RefLine} nodes::NTuple{2, Int} end @@ -300,6 +308,7 @@ struct Pyramid <: AbstractCell{RefPyramid} nodes::NTuple{5, Int} end +geometric_interpolation(::Type{Point}) = Lagrange{RefPoint, 1}() geometric_interpolation(::Type{Line}) = Lagrange{RefLine, 1}() geometric_interpolation(::Type{QuadraticLine}) = Lagrange{RefLine, 2}() geometric_interpolation(::Type{Triangle}) = Lagrange{RefTriangle, 1}() diff --git a/src/exports.jl b/src/exports.jl index c710861c15..f6ec96548a 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -4,6 +4,7 @@ export VectorInterpolation, ScalarInterpolation, VectorizedInterpolation, + #RefPoint, RefLine, RefQuadrilateral, RefHexahedron, @@ -58,6 +59,7 @@ export # Grid Grid, Node, + #Point, Line, QuadraticLine, Triangle, diff --git a/src/interpolations.jl b/src/interpolations.jl index 5aaed4fdfb..3bc2528299 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -1308,6 +1308,19 @@ function reference_shape_value(ip::Lagrange{RefPyramid, 2}, ξ::Vec{3, T}, i::In throw(ArgumentError("no shape function $i for interpolation $ip")) end +##################################### +# Lagrange dim 0 RefPoint order 0 # +##################################### +getnbasefunctions(::Lagrange{RefPoint}) = 1 +vertexdof_indices(::Lagrange{RefPoint}) = ((1,),) +function reference_coordinates(::Lagrange{RefPoint}) + return [Tensor{1, 0, Float64, 0}(())] # zero dim Vec{0} +end +function reference_shape_value(ip::Lagrange{RefPoint}, ξ::Vec{0}, i::Int) + return i == 1 && 1.0 + throw(ArgumentError("no shape function $i for interpolation $ip")) +end + ################### # Bubble elements # ################### From fa0fc88d0c2dec53f840004e93764c73fa1b07b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Elias=20B=C3=B6rjesson?= Date: Sun, 8 Dec 2024 10:35:42 +0100 Subject: [PATCH 2/3] Work on Rigid Connector constraint. --- src/Dofs/ConstraintHandler.jl | 81 ++++++++++++++++++++++++++++++++++- 1 file changed, 80 insertions(+), 1 deletion(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index a48cb3940f..64a87ff24f 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -82,7 +82,7 @@ A collection of constraints associated with the dof handler `dh`. `T` is the numeric type for stored values. """ mutable struct ConstraintHandler{DH <: AbstractDofHandler, T} - const dbcs::Vector{Dirichlet} + const dbcs::Vector{Dirichlet} #TODO: Vector{Constraint} const prescribed_dofs::Vector{Int} const free_dofs::Vector{Int} const inhomogeneities::Vector{T} @@ -1762,3 +1762,82 @@ function _condense_local!( end return end + + +##### New constraint + +struct RigidConnector # <: Constraint + rigidbody_cellid::Int #Cell id + rigidbody_field::Symbol + facets::Vector{FacetIndex} + field_name::Symbol +end + +function RigidConnector(; + rigidbody_cellid, + rigidbody_field = :rb, + facets, + field_name = :u + ) + return RigidConnector(rigidbody_cellid, rigidbody_field, collect(facets), field_name) +end + +function add!(ch::ConstraintHandler, c::RigidConnector) + (; dh) = ch + grid = get_grid(dh) + + xdof, ydof, θdof = celldofs(dh, c.rigidbody_cellid) + + _check_same_celltype(grid, c.facets) + @assert getcells(grid, c.rigidbody_cellid) isa Point + r = getcoordinates(grid, c.rigidbody_cellid)[1] + + cellid = first(c.facets)[1] + sdh_index = dh.cell_to_subdofhandler[cellid] + sdh = dh.subdofhandlers[sdh_index] + + field_idx = find_field(sdh, :u) #Hardcoded + CT = getcelltype(sdh) + ip = getfieldinterpolation(sdh, field_idx) + ip_geo = geometric_interpolation(CT) + offset = field_offset(sdh, field_idx) + n_comp = n_dbc_components(ip) + + local_facet_dofs, local_facet_dofs_offset = + _local_facet_dofs_for_bc(ip.ip, n_comp, 1:n_comp, offset) + + dofsadded = Set{Int}() + fv = BCValues(ip.ip, ip_geo) + + cc = CellCache(dh, UpdateFlags(; nodes = false, coords = true, dofs = true)) + for (cellidx, entityidx) in c.facets + reinit!(cc, cellidx) + + # no need to reinit!, enough to update current_entity since we only need geometric shape functions M + fv.current_entity = entityidx + + # local dof-range for this facet + erange = local_facet_dofs_offset[entityidx]:(local_facet_dofs_offset[entityidx + 1] - 1) + dofs = cc.dofs + #dofs[1] in dofsadded && continue + + c = 0 + for ipoint in getnquadpoints(fv) + x = spatial_coordinate(fv, ipoint, cc.coords) + + # r + A*ū == x + ū = x - r + uperp = rotate(ū, pi / 2) + + for d in 1:2 + c += 1 + globaldof = cc.dofs[local_facet_dofs[erange[c]]] + a1 = AffineConstraint(globaldof, [xdof => d == 1 ? 1.0 : 0.0, ydof => d == 2 ? 1.0 : 0.0, θdof => uperp[d]], 0.0) + add!(ch, a1) + push!(dofsadded, globaldof) + end + end + end + + return +end From be445f8f1af383e4b527032c4adec3484cf891d9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Elias=20B=C3=B6rjesson?= Date: Sun, 8 Dec 2024 10:35:57 +0100 Subject: [PATCH 3/3] Example demonstrating rigid connector --- docs/src/literate-tutorials/plate_hole.geo | 44 ++++++ .../src/literate-tutorials/rigid_connector.jl | 136 ++++++++++++++++++ 2 files changed, 180 insertions(+) create mode 100644 docs/src/literate-tutorials/plate_hole.geo create mode 100644 docs/src/literate-tutorials/rigid_connector.jl diff --git a/docs/src/literate-tutorials/plate_hole.geo b/docs/src/literate-tutorials/plate_hole.geo new file mode 100644 index 0000000000..6bd7ee7b63 --- /dev/null +++ b/docs/src/literate-tutorials/plate_hole.geo @@ -0,0 +1,44 @@ +// Parameters +L = 10; // Length of the plate +H = 10; // Height of the plate +r = 3; // Radius of the hole +lc = 0.2; // Mesh size parameter + +// Define the plate geometry +Point(1) = {0, 0, 0, lc}; +Point(2) = {L, 0, 0, lc}; +Point(3) = {L, H, 0, lc}; +Point(4) = {0, H, 0, lc}; + +// Define the hole geometry +Point(5) = {L/2, H/2, 0, lc}; // Center of the hole +Point(6) = {L/2 + r, H/2, 0, lc}; // Start point on the circle +Point(7) = {L/2, H/2 + r, 0, lc}; // Top point on the circle +Point(8) = {L/2 - r, H/2, 0, lc}; // Left point on the circle +Point(9) = {L/2, H/2 - r, 0, lc}; // Bottom point on the circle + +// Plate lines +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 4}; +Line(4) = {4, 1}; + +// Define circular arcs for the hole +Circle(5) = {6, 5, 7}; +Circle(6) = {7, 5, 8}; +Circle(7) = {8, 5, 9}; +Circle(8) = {9, 5, 6}; + +// Create surface with the hole +Curve Loop(1) = {1, 2, 3, 4}; +Curve Loop(2) = {5, 6, 7, 8}; +Plane Surface(1) = {1, 2}; + +// Define physical groups for boundary conditions +Physical Curve("PlateRightLeft") = {2, 4}; +Physical Curve("PlateExterior") = {1, 2, 3, 4}; +Physical Curve("HoleInterior") = {5, 6, 7, 8}; +Physical Surface("PlateWithHole") = {1}; + +// Mesh the geometry +Mesh 2; diff --git a/docs/src/literate-tutorials/rigid_connector.jl b/docs/src/literate-tutorials/rigid_connector.jl new file mode 100644 index 0000000000..9df57ae5b0 --- /dev/null +++ b/docs/src/literate-tutorials/rigid_connector.jl @@ -0,0 +1,136 @@ +using Ferrite +using FerriteGmsh + + +grid = togrid("docs/src/literate-tutorials/plate_hole.geo") + +grid = Grid( + Ferrite.AbstractCell[grid.cells...], + grid.nodes, + facetsets = grid.facetsets, + cellsets = grid.cellsets +) +rigidbody_node = Node(Vec((5.0, 5.0))) +rigidbody_cellid = getncells(grid) + 1 +push!(grid.nodes, rigidbody_node) +push!(grid.cells, Ferrite.Point((getnnodes(grid),))) + +addcellset!(grid, "rigidbody", [getncells(grid)]) +addvertexset!(grid, "rigidvertex", x -> x ≈ rigidbody_node.x) + +ip_u = Lagrange{RefTriangle, 1}()^2 +ip_rb_u = Lagrange{Ferrite.RefPoint, 1}()^2 +ip_rb_θ = Lagrange{Ferrite.RefPoint, 1}() + +qr = QuadratureRule{RefTriangle}(2) +cellvalues = CellValues(qr, ip_u) + +dh = DofHandler(grid) +sdh = SubDofHandler(dh, getcellset(grid, "PlateWithHole")) +add!(sdh, :u, ip_u) + +sdh = SubDofHandler(dh, getcellset(grid, "rigidbody")) +add!(sdh, :u, ip_rb_u) +add!(sdh, :θ, ip_rb_θ) +close!(dh) + +ch = ConstraintHandler(dh) +rb = Ferrite.RigidConnector(; + rigidbody_cellid = rigidbody_cellid, + facets = getfacetset(grid, "HoleInterior"), +) +add!(ch, rb) +add!(ch, Dirichlet(:u, getfacetset(grid, "PlateRightLeft"), x -> (0.0, 0.0))) +add!(ch, Dirichlet(:u, getvertexset(grid, "rigidvertex"), x -> (0.0, 0.0))) +add!(ch, Dirichlet(:θ, getvertexset(grid, "rigidvertex"), x -> (0.1))) +close!(ch) + +Emod = 200.0e3 # Young's modulus [MPa] +ν = 0.3 # Poisson's ratio [-] +Gmod = Emod / (2(1 + ν)) # Shear modulus +Kmod = Emod / (3(1 - 2ν)) # Bulk modulus +C = gradient(ϵ -> 2 * Gmod * dev(ϵ) + 3 * Kmod * vol(ϵ), zero(SymmetricTensor{2, 2})); + +function assemble_cell!(ke, cellvalues, C) + for q_point in 1:getnquadpoints(cellvalues) + ## Get the integration weight for the quadrature point + dΩ = getdetJdV(cellvalues, q_point) + for i in 1:getnbasefunctions(cellvalues) + ## Gradient of the test function + ∇Nᵢ = shape_gradient(cellvalues, q_point, i) + for j in 1:getnbasefunctions(cellvalues) + ## Symmetric gradient of the trial function + ∇ˢʸᵐNⱼ = shape_symmetric_gradient(cellvalues, q_point, j) + ke[i, j] += (∇Nᵢ ⊡ C ⊡ ∇ˢʸᵐNⱼ) * dΩ + end + end + end + return ke +end + +function assemble_global!(K, dh, cellvalues, C, cellset) + ## Allocate the element stiffness matrix + n_basefuncs = getnbasefunctions(cellvalues) + ke = zeros(n_basefuncs, n_basefuncs) + ## Create an assembler + assembler = start_assemble(K) + ## Loop over all cells + for cell in CellIterator(dh, cellset) + ## Update the shape function gradients based on the cell coordinates + reinit!(cellvalues, cell) + ## Reset the element stiffness matrix + fill!(ke, 0.0) + ## Compute element contribution + assemble_cell!(ke, cellvalues, C) + ## Assemble ke into K + assemble!(assembler, celldofs(cell), ke) + end + return K +end + +K = allocate_matrix(dh, ch) +f_ext = zeros(Float64, ndofs(dh)) + +assemble_global!(K, dh, cellvalues, C, getcellset(grid, "PlateWithHole")); + +apply!(K, f_ext, ch) +a = K \ f_ext; + +apply!(a, ch) + + +function calculate_stresses(grid, dh, cv, u, C, cellset) + qp_stresses = [ + [zero(SymmetricTensor{2, 2}) for _ in 1:getnquadpoints(cv)] + for _ in 1:getncells(grid) + ] + avg_cell_stresses = tuple((zeros(length(cellset)) for _ in 1:3)...) + for cell in CellIterator(dh, cellset) + reinit!(cv, cell) + cell_stresses = qp_stresses[cellid(cell)] + for q_point in 1:getnquadpoints(cv) + ε = function_symmetric_gradient(cv, q_point, u, celldofs(cell)) + cell_stresses[q_point] = C ⊡ ε + end + σ_avg = sum(cell_stresses) / getnquadpoints(cv) + avg_cell_stresses[1][cellid(cell)] = σ_avg[1, 1] + avg_cell_stresses[2][cellid(cell)] = σ_avg[2, 2] + avg_cell_stresses[3][cellid(cell)] = σ_avg[1, 2] + end + return qp_stresses, avg_cell_stresses +end + +qp_stresses, avg_cell_stresses = calculate_stresses(grid, dh, cellvalues, a, C, getcellset(grid, "PlateWithHole")); + +# We now use the the L2Projector to project the stress-field onto the piecewise linear +# finite element space that we used to solve the problem. +proj = L2Projector(grid) +add!(proj, getcellset(grid, "PlateWithHole"), ip_u; qr_rhs = qr) +close!(proj) + +projected = project(proj, qp_stresses) + +VTKGridFile("rigid_con", grid) do vtk + write_solution(vtk, dh, a) + write_projection(vtk, proj, projected, "stress") +end