From 852ee67f79bb67ad4b4f36ba65e45de3cb3a8500 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 9 May 2026 16:33:16 +0300 Subject: [PATCH] refactor(basis): remove legacy plate_elements DKT basis module Delete AbstractPlateBasis/DKT definitions from basis/; plate handling moves elsewhere. --- src/basis/plate_elements.jl | 416 ------------------------------------ 1 file changed, 416 deletions(-) delete mode 100644 src/basis/plate_elements.jl diff --git a/src/basis/plate_elements.jl b/src/basis/plate_elements.jl deleted file mode 100644 index ad28ba9..0000000 --- a/src/basis/plate_elements.jl +++ /dev/null @@ -1,416 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE - -# ============================================================================ -# Plate Element Basis Functions -# ============================================================================ -# -# This file defines basis functions for plate bending elements, which are -# typically NON-CONFORMING (C0 continuity for deflection w, discontinuous -# rotations θx, θy). -# -# These elements are special because: -# 1. Multiple DOF types per node: {w, θx, θy} -# 2. Mixed continuity: w is C0, rotations are discontinuous -# 3. Enforce Kirchhoff constraints (zero transverse shear) via basis construction -# 4. Often use higher-order basis internally (e.g., DKT uses Tri6 internally) -# -# Theory: -# Plate bending theory (Kirchhoff-Love or Mindlin-Reissner) -# Non-conforming finite elements -# Discrete Kirchhoff technique -# -# References: -# - Batoz, J.-L., et al. (1980). "A study of three-node triangular plate -# bending elements." Int. J. Numer. Methods Eng., 15(12), 1771-1812. -# - Zienkiewicz & Taylor, "The Finite Element Method" Vol. 2, Chapter 10 -# -# See also: docs/book/plate_element_basis.md (when created) -# ============================================================================ - -using Tensors -using LinearAlgebra - -""" - AbstractPlateBasis <: AbstractBasis - -Abstract base type for plate bending element basis functions. - -Plate elements have special characteristics: -- **Multiple DOF types per node:** Typically {w, θx, θy} at each node -- **Mixed continuity:** Deflection w is C0, rotations are discontinuous -- **Non-conforming:** Violate strict C1 requirement but converge correctly -- **Constraint-enforcing:** Kirchhoff constraint built into shape functions - -# Interface Requirements - -In addition to standard `AbstractBasis` interface, plate basis types should implement: -- `ndofs(::AbstractPlateBasis)` - Total number of DOFs (often > nnodes) -- `dof_types(::AbstractPlateBasis)` - Tuple of DOF symbols, e.g., (:w, :θx, :θy) - -# Concrete Types - -- [`DKT`](@ref) - Discrete Kirchhoff Triangle (3 nodes, 9 DOFs) -- DST - Discrete Shear Triangle (future) -- Morley - Morley triangle (future) -- DKQ - Discrete Kirchhoff Quadrilateral (future) - -See also: [`AbstractBasis`](@ref), [`DKT`](@ref) -""" -abstract type AbstractPlateBasis <: AbstractBasis end - -# ============================================================================ -# DKT (Discrete Kirchhoff Triangle) -# ============================================================================ - -""" - DKT <: AbstractPlateBasis - -Discrete Kirchhoff Triangle - Non-conforming plate bending element. - -# Mathematical Properties - -**Element Type:** 3-node triangular plate bending element - -**DOFs:** 9 total (3 per node) -- Node i: {wᵢ, θxᵢ, θyᵢ} - - wᵢ: Deflection (transverse displacement) - - θxᵢ: Rotation about x-axis (∂w/∂y) - - θyᵢ: Rotation about y-axis (-∂w/∂x) - -**Continuity:** -- C0 for deflection w (continuous across element boundaries) -- Discontinuous for rotations θx, θy (non-conforming!) - -**Kirchhoff Constraint:** -Zero transverse shear enforced at specific points on each edge: -- γxz = 0 (shear strain in xz plane) -- γyz = 0 (shear strain in yz plane) - -**Internal Basis:** -Uses Lagrange{Triangle, 2} (Tri6 quadratic) basis functions internally -to construct the 9 DOF shape functions. - -**Convergence:** -Despite being non-conforming (discontinuous rotations), DKT passes the -patch test and converges correctly to the Kirchhoff plate solution. - -# Topology Requirements - -- **Topology:** Triangle only -- **Nodes:** 3 (vertices of triangle) -- **Integration:** Typically Gauss{2} or Gauss{3} - -# DOF Ordering - -DOFs are ordered by node, then by type: -``` -DOF index: 1 2 3 4 5 6 7 8 9 -Node: 1 1 1 2 2 2 3 3 3 -Type: w θx θy w θx θy w θx θy -``` - -# Mathematical Background - -The DKT element constructs shape functions that satisfy: - -1. **Compatibility:** ∑ Nᵢ = 1 (partition of unity) -2. **Interpolation:** wᵢ(xⱼ) = δᵢⱼ for deflection -3. **Kirchhoff constraint:** γxz = γyz = 0 at 2 points per edge (6 total) - -The rotations are interpolated using modified shape functions: -- θx = ∑ Hxᵢ θxᵢ -- θy = ∑ Hyᵢ θyᵢ - -where Hx and Hy are constructed from Tri6 basis to enforce constraints. - -# References - -- Batoz, J.-L., Bathe, K.-J., & Ho, L.-W. (1980). - "A study of three-node triangular plate bending elements." - International Journal for Numerical Methods in Engineering, 15(12), 1771-1812. - DOI: 10.1002/nme.1620151206 - -- Batoz, J.-L., & Dhatt, G. (1990). - "Modélisation des structures par éléments finis, Vol. 2: Poutres et plaques." - Hermès, Paris. - -# Example - -```julia -using JuliaFEM - -# Create DKT basis -basis = DKT() - -# Query properties -nnodes(basis) # → 3 (triangle vertices) -ndofs(basis) # → 9 (3 DOFs per node) -ndims(basis) # → 2 (2D element) -dof_types(basis) # → (:w, :θx, :θy) - -# Evaluate basis functions at parametric point -xi = Vec(0.25, 0.25) -N = get_basis_functions(Triangle(), basis, xi) -# Returns: NTuple{9, Float64} - all 9 DOF shape functions - -# Evaluate derivatives (for stiffness matrix) -dN = get_basis_derivatives(Triangle(), basis, xi) -# Returns: NTuple{9, Vec{2, Float64}} - gradients of all 9 DOF shape functions - -# Use in element assembly -topology = Triangle() -element = Element(topology, basis, Gauss{2}(), (1, 2, 3)) -``` - -See also: [`AbstractPlateBasis`](@ref), [`get_basis_functions`](@ref), [`get_basis_derivatives`](@ref) -""" -struct DKT <: AbstractPlateBasis end - -# ============================================================================ -# Interface Implementation for DKT -# ============================================================================ - -""" - nnodes(::DKT) -> Int - -Number of nodes in DKT element (always 3 for triangle vertices). -""" -nnodes(::DKT) = 3 -nnodes(::Type{DKT}) = 3 - -""" - ndofs(::DKT) -> Int - -Total number of degrees of freedom in DKT element (9 = 3 nodes × 3 DOFs/node). -""" -ndofs(::DKT) = 9 -ndofs(::Type{DKT}) = 9 - -""" - Base.ndims(::DKT) -> Int - -Spatial dimension of DKT element (always 2 for 2D plate element). -""" -Base.ndims(::DKT) = 2 -Base.ndims(::Type{DKT}) = 2 - -""" - dof_types(::DKT) -> NTuple{3, Symbol} - -Types of DOFs at each node: deflection w, rotation θx, rotation θy. -""" -dof_types(::DKT) = (:w, :θx, :θy) -dof_types(::Type{DKT}) = (:w, :θx, :θy) - -# ============================================================================ -# DKT Basis Function Evaluation -# ============================================================================ - -""" - get_basis_functions(::Triangle, ::DKT, xi::Vec{2, T}) where T -> NTuple{9, T} - -Evaluate all 9 DKT basis functions at parametric point ξ. - -Returns shape function values for 9 DOFs in order: -[N_w1, N_θx1, N_θy1, N_w2, N_θx2, N_θy2, N_w3, N_θx3, N_θy3] - -# Algorithm - -1. Evaluate Tri6 (quadratic triangle) basis functions internally -2. Construct DKT shape functions using geometric coefficients -3. Return 9 DOF shape functions - -# Arguments -- `::Triangle`: Triangle topology (required for dispatch) -- `::DKT`: DKT basis type -- `xi::Vec{2, T}`: Parametric coordinates (ξ, η) in area coordinates - -# Returns -- `NTuple{9, T}`: Shape function values for all 9 DOFs - -# Notes - -The DKT element has a special DOF structure: -- DOFs 1, 4, 7: Deflection w at nodes 1, 2, 3 -- DOFs 2, 5, 8: Rotation θx at nodes 1, 2, 3 -- DOFs 3, 6, 9: Rotation θy at nodes 1, 2, 3 - -For assembling plate problems, you typically need: -- N values for mass matrix -- dN values for stiffness matrix (use get_basis_derivatives) - -# Example - -```julia -topology = Triangle() -basis = DKT() -xi = Vec(1/3, 1/3) # Element centroid - -N = get_basis_functions(topology, basis, xi) -# N[1] = shape function for w at node 1 -# N[2] = shape function for θx at node 1 -# N[3] = shape function for θy at node 1 -# ... and so on -``` - -See also: [`get_basis_derivatives`](@ref), [`DKT`](@ref) -""" -function get_basis_functions(::Triangle, ::DKT, xi::Vec{2,T}) where {T} - # Extract area coordinates - ξ = xi[1] - η = xi[2] - - # Evaluate Tri6 (quadratic triangle) basis internally - # This uses the existing Lagrange{Triangle, 2} implementation - N_tri6 = get_basis_functions(Triangle(), Lagrange{Triangle,2}(), xi) - - # N_tri6 = (N1, N2, N3, N4, N5, N6) for 6-node quadratic triangle: - # N1, N2, N3: Corner nodes (vertices) - # N4: Midpoint of edge 2-3 - # N5: Midpoint of edge 3-1 - # N6: Midpoint of edge 1-2 - - N1, N2, N3, N4, N5, N6 = N_tri6 - - # For DKT, we need to construct shape functions for 9 DOFs: - # - Deflection w at 3 nodes (uses standard Tri3 linear basis) - # - Rotations θx, θy at 3 nodes (constructed from Tri6 basis) - - # Deflection w uses linear (P1) basis functions - # ζ = 1 - ξ - η (third area coordinate) - ζ = one(T) - ξ - η - - # Shape functions for deflection w (simple linear interpolation) - Nw1 = ζ # Node 1 - Nw2 = ξ # Node 2 - Nw3 = η # Node 3 - - # Shape functions for rotations θx and θy are more complex - # They are constructed from the Tri6 basis to enforce Kirchhoff constraints - # The exact construction depends on element geometry (edge lengths, angles) - # which is handled in the element-specific assembly code - - # For a generic evaluation (without element geometry), we use the - # Tri6 basis directly as a placeholder. The actual DKT shape functions - # are geometry-dependent and computed during element assembly. - - # NOTE: This is a simplified implementation for the interface. - # Real DKT assembly uses element-specific shape functions that depend - # on edge lengths, angles, etc. (see src/plates/dkt.jl) - - # Rotation shape functions (simplified - geometry-independent approximation) - # In practice, these should be constructed from DKTShapeFunctions coefficients - - # For now, return a valid NTuple{9, T} structure - # TODO: Implement full geometry-dependent shape functions - - # Placeholder: Use Tri6 basis for rotations (not physically correct!) - # This needs to be refined based on element geometry - Nθx1 = T(1.5) * (N6 - N5) # Simplified - Nθy1 = T(1.5) * (N6 - N5) - - Nθx2 = T(1.5) * (N4 - N6) - Nθy2 = T(1.5) * (N4 - N6) - - Nθx3 = T(1.5) * (N5 - N4) - Nθy3 = T(1.5) * (N5 - N4) - - # Return all 9 DOF shape functions - # Order: [w1, θx1, θy1, w2, θx2, θy2, w3, θx3, θy3] - return (Nw1, Nθx1, Nθy1, Nw2, Nθx2, Nθy2, Nw3, Nθx3, Nθy3) -end - -""" - get_basis_derivatives(::Triangle, ::DKT, xi::Vec{2, T}) where T -> NTuple{9, Vec{2, T}} - -Evaluate derivatives of all 9 DKT basis functions at parametric point ξ. - -Returns gradients (∂N/∂ξ, ∂N/∂η) for 9 DOFs in order: -[∇N_w1, ∇N_θx1, ∇N_θy1, ∇N_w2, ∇N_θx2, ∇N_θy2, ∇N_w3, ∇N_θx3, ∇N_θy3] - -# Algorithm - -1. Evaluate Tri6 basis function derivatives internally -2. Construct DKT shape function derivatives using geometric coefficients -3. Return 9 DOF gradient vectors - -# Arguments -- `::Triangle`: Triangle topology (required for dispatch) -- `::DKT`: DKT basis type -- `xi::Vec{2, T}`: Parametric coordinates (ξ, η) in area coordinates - -# Returns -- `NTuple{9, Vec{2, T}}`: Gradient vectors for all 9 DOFs - -# Notes - -These derivatives are critical for: -- Computing curvatures κx, κy, κxy (second derivatives of w) -- Assembling element stiffness matrix -- Computing bending moments and shear forces - -The derivatives are in parametric coordinates. Transform to physical coordinates using: -```julia -dN_dx = inv(J) * dN_dxi # J = Jacobian matrix -``` - -# Example - -```julia -topology = Triangle() -basis = DKT() -xi = Vec(1/3, 1/3) - -dN = get_basis_derivatives(topology, basis, xi) -# dN[1] = ∇N_w1 = (∂N_w1/∂ξ, ∂N_w1/∂η) -# dN[2] = ∇N_θx1 = (∂N_θx1/∂ξ, ∂N_θx1/∂η) -# ... and so on -``` - -See also: [`get_basis_functions`](@ref), [`DKT`](@ref) -""" -function get_basis_derivatives(::Triangle, ::DKT, xi::Vec{2,T}) where {T} - # Extract area coordinates - ξ = xi[1] - η = xi[2] - - # Evaluate Tri6 basis derivatives internally - dN_tri6 = get_basis_derivatives(Triangle(), Lagrange{Triangle,2}(), xi) - - # Derivatives of linear basis for deflection w - # ∂/∂ξ: [-1, 1, 0] - # ∂/∂η: [-1, 0, 1] - dNw1 = Vec{2,T}((-one(T), -one(T))) # ∂ζ/∂ξ = -1, ∂ζ/∂η = -1 - dNw2 = Vec{2,T}((one(T), zero(T))) # ∂ξ/∂ξ = 1, ∂ξ/∂η = 0 - dNw3 = Vec{2,T}((zero(T), one(T))) # ∂η/∂ξ = 0, ∂η/∂η = 1 - - # Derivatives of rotation shape functions - # These are geometry-dependent and should be computed from - # DKTShapeFunctions coefficients during element assembly - - # For now, use placeholder derivatives from Tri6 basis - # TODO: Implement full geometry-dependent derivatives - - dN4, dN5, dN6 = dN_tri6[4], dN_tri6[5], dN_tri6[6] - - # Simplified rotation derivatives (placeholder) - dNθx1 = T(1.5) * (dN6 - dN5) - dNθy1 = T(1.5) * (dN6 - dN5) - - dNθx2 = T(1.5) * (dN4 - dN6) - dNθy2 = T(1.5) * (dN4 - dN6) - - dNθx3 = T(1.5) * (dN5 - dN4) - dNθy3 = T(1.5) * (dN5 - dN4) - - # Return all 9 DOF gradients - # Order: [∇w1, ∇θx1, ∇θy1, ∇w2, ∇θx2, ∇θy2, ∇w3, ∇θx3, ∇θy3] - return (dNw1, dNθx1, dNθy1, dNw2, dNθx2, dNθy2, dNw3, dNθx3, dNθy3) -end - -# ============================================================================ -# Export -# ============================================================================ - -export AbstractPlateBasis, DKT, dof_types, ndofs