Rigid body connector

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;

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

src/Dofs/ConstraintHandler.jl

@@ -82,7 +82,7 @@
 `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 @@
     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

src/Export/VTK.jl

@@ -65,6 +65,7 @@
     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
 

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}

src/Grid/grid.jl

@@ -193,6 +193,11 @@
     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 @@
 ######################################################
 
 # 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 @@
     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}()

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,

src/interpolations.jl

@@ -1308,6 +1308,19 @@
     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 #
 ###################