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] 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