Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-17 01:02:13 +00:00
Code Issues Packages Projects Releases Wiki Activity

chore: remove obsolete and legacy files

Browse Source 
Major cleanup: remove obsolete files that have been replaced or
moved to new locations. Files removed include legacy assembly
implementations, old API files, and deprecated test files.

Removed files:

Assemblers and assembly:
- src/assemblers/nodal_based.jl
- src/assemblers/nodal_cache.jl
- src/assemblers/node_based_coo.jl
- src/assembly/assembly.jl
- src/assembly/element_structures.jl
- src/assembly/framework.jl
- src/assembly/nodal_structures.jl
- src/assembly/problems.jl
- src/element_assembly_structures.jl
- src/nodal_assembly_structures.jl

Legacy API and structure files:
- src/beams/api.jl
- src/formulations/api.jl
- src/gpu_elasticity.jl
- src/io.jl
- src/materials_plasticity.jl
- src/postprocess_utils.jl
- src/preprocess.jl
- src/quadrature.jl
- src/readers.jl
- src/shells/api.jl
- src/trusses/api.jl

Elements and domains:
- src/domains/continuum/assemble_v2.jl
- src/elements/integrate.jl

Quadrature legacy files:
- src/quadrature/gauss_points.jl
- src/quadrature/integration.jl

Test files:
- test/runtests_new.jl
- test/runtests.jl.old
- test/test_problems_elasticity_assemble_3d_seg3.jl
This commit is contained in:
Jukka Aho
2025-12-12 22:45:49 +02:00
parent cba8d19c1c
commit ff3c26d2c3
28 changed files with 0 additions and 6419 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/assemblers/nodal_based.jl
-177
View File
@@ -1,177 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
Nodal-based assembly (future implementation).
Node-by-node assembly using inverse connectivity (node-to-elements map).
Natural for GPU parallelization (one thread per node).
**Performance**: Expected 2-10x speedup on GPU for large problems (> 100k nodes).
**Best for**: GPU acceleration, very large problems.
**Status**: Planned for future implementation.
"""
using SparseArrays
"""
assemble!(
cache::NodalCache,
assembler::NodalAssembler,
kernel::AbstractKernel,
mesh::AbstractMesh
) -> Nothing
Assemble global system using nodal traversal **in-place, zero allocations**.
**Status**: Not yet implemented. Placeholder for future GPU-based assembly.
# Algorithm (Planned)
1. Reset cache
2. Loop over nodes (parallelizable on GPU):
a. Get all elements touching this node
b. For each touching element:
- Compute full element stiffness (or use cached value)
- Extract only rows/columns for this node
- Accumulate to global system (atomic add on GPU)
# GPU Parallelization
Each node processed by one GPU thread:
```cuda
__global__ void assemble_nodal(nodes, elements, K, f) {
int node_id = blockIdx.x * blockDim.x + threadIdx.x;
if (node_id >= nnodes) return;
// Get touching elements
for (elem in touching_elements[node_id]) {
// Compute element contribution for this node
// Atomic add to K, f
}
}
```
# Arguments
- `cache`: Pre-allocated nodal cache
- `assembler`: Nodal assembler
- `kernel`: Domain kernel
- `mesh`: Finite element mesh
# Zero-Allocation Guarantee
All arrays pre-allocated. Atomic operations on GPU ensure thread-safety.
"""
function assemble!(
cache::NodalCache,
assembler::NodalAssembler,
kernel::AbstractKernel,
mesh::AbstractMesh
)
error("NodalAssembler not yet implemented. Use COOAssembler or CSCAssembler.")
# Planned implementation:
# 1. Reset cache
# 2. Loop over nodes
# 3. For each node, get touching elements from cache.node_to_elements
# 4. Accumulate contributions from all touching elements
# 5. Write to global K, f (atomic on GPU)
return nothing
end
# ============================================================================
# HELPER FUNCTIONS (for future implementation)
# ============================================================================
"""
create_cache(assembler::NodalAssembler, mesh::AbstractMesh, kernel::AbstractKernel) -> NodalCache
Create pre-allocated cache for nodal assembly.
Builds node-to-elements map (inverse connectivity) for efficient nodal traversal.
# Arguments
- `assembler`: Nodal assembler
- `mesh`: Finite element mesh
- `kernel`: Domain kernel
# Returns
- Pre-allocated nodal cache with inverse connectivity
"""
function create_cache(assembler::NodalAssembler, mesh::AbstractMesh, kernel::AbstractKernel)
return NodalCache(mesh, kernel)
end
"""
compute_node_contributions!(
node_cache::NodeCache,
element_cache::ElementCache,
node_id::Int,
kernel::AbstractKernel,
mesh::AbstractMesh,
node_to_elements::NodeToElementsMap
) -> Nothing
Compute contributions to this node from all touching elements **in-place**.
# Algorithm
1. Get all elements touching this node
2. For each element:
a. Compute full element stiffness (using `compute_element_stiffness!`)
b. Find local node index within element
c. Extract rows/columns corresponding to this node's DOFs
d. Accumulate to node contributions
# Arguments
- `node_cache`: Node workspace (output)
- `element_cache`: Element workspace (for kernel calls)
- `node_id`: Node index
- `kernel`: Domain kernel
- `mesh`: Finite element mesh
- `node_to_elements`: Inverse connectivity
# Returns
Nothing - writes to `node_cache` in-place.
"""
function compute_node_contributions!(
node_cache::NodeCache,
element_cache::ElementCache,
node_id::Int,
kernel::AbstractKernel,
mesh::AbstractMesh,
node_to_elements::NodeToElementsMap
)
error("compute_node_contributions! not yet implemented")
# Planned implementation:
# 1. Get touching elements: get_node_spider(node_to_elements, node_id)
# 2. For each element:
# - Compute element stiffness: compute_element_stiffness!(element_cache, ...)
# - Find local node index in element
# - Extract node DOF rows/columns from Ke, fe
# - Accumulate to node_cache
# 3. Return node_cache (contains all contributions for this node)
return nothing
end
# ============================================================================
# GPU KERNEL STUBS (for future CUDA implementation)
# ============================================================================
# """
# assemble_nodal_gpu!(K, f, mesh, kernel, node_to_elements)
#
# GPU kernel for nodal assembly.
#
# Launches one thread per node. Each thread:
# 1. Gets touching elements for its node
# 2. Computes contributions from all elements
# 3. Atomically adds to global K, f
#
# Requires CUDA.jl or similar GPU framework.
# """
# function assemble_nodal_gpu! end
src/assemblers/nodal_cache.jl
-125
View File
@@ -1,125 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
Nodal cache for node-based assembly.
Nodal assembly iterates over nodes rather than elements. Each node assembles
contributions from all touching elements. This approach has advantages for:
- Contact mechanics (contact is inherently nodal)
- Domain decomposition (clear node ownership)
- Matrix-free operations (natural node-based matvec)
- Adaptive refinement (local node operations)
# Structure
Pre-builds node-to-elements map (inverse connectivity) once. During assembly,
each node visits all touching elements and accumulates contributions.
# Performance
Similar to CSC for standard problems. Better locality for nodal operations
like contact and matrix-free solvers.
# Use case
Best for problems with nodal phenomena (contact, nodal plasticity) or
matrix-free iterative solvers.
"""
using SparseArrays
"""
NodalCache <: AbstractAssemblerCache
Cache for nodal-based assembly.
Pre-allocates sparse matrix, force vector, and node-to-elements map.
Each node assembles contributions from all touching elements.
# Fields
- `K::SparseMatrixCSC{Float64,Int}`: Sparse matrix
- `f::Vector{Float64}`: Global force vector
- `node_cache::NodeCache`: Per-node workspace
- `element_cache::ElementCache`: Per-element workspace (for kernel calls)
- `node_to_elements::NodeToElementsMap`: Inverse connectivity
# Zero-Allocation Usage
```julia
cache = NodalCache(mesh, kernel)
fill!(cache)
assemble!(cache, assembler, kernel, mesh) # No allocations
K, f = extract_system(cache)
```
"""
mutable struct NodalCache <: AbstractAssemblerCache
K::SparseMatrixCSC{Float64,Int} # Sparse matrix
f::Vector{Float64} # Force vector
node_cache::NodeCache # Node workspace
element_cache::ElementCache # Element workspace
node_to_elements::NodeToElementsMap # Inverse connectivity
end
"""
NodalCache(mesh::AbstractMesh, kernel::AbstractKernel) -> NodalCache
Create pre-allocated nodal cache.
Builds node-to-elements map (inverse connectivity) for efficient nodal traversal.
# Arguments
- `mesh`: Finite element mesh
- `kernel`: Domain kernel defining DOF structure
# Returns
- `NodalCache` with pre-allocated workspace and inverse connectivity
"""
function NodalCache(mesh::AbstractMesh, kernel::AbstractKernel)
ndofs_per_node = dofs_per_node(kernel)
nnodes_mesh = nnodes_total(mesh)
ndofs = nnodes_mesh * ndofs_per_node
# Build sparsity pattern (same as CSC)
K = build_sparsity_pattern(mesh, kernel)
f = zeros(Float64, ndofs)
# Create caches
node_cache = create_node_cache(mesh, kernel)
element_cache = create_element_cache(mesh, kernel)
# Build inverse connectivity
node_to_elements = NodeToElementsMap(mesh)
return NodalCache(K, f, node_cache, element_cache, node_to_elements)
end
"""
reset!(cache::NodalCache)
Reset nodal cache for new assembly.
Zeros out matrix values and force vector.
**Zero allocations** - reuses existing arrays.
"""
function reset!(cache::NodalCache)
fill!(cache.K.nzval, 0.0)
fill!(cache.f, 0.0)
return nothing
end
"""
extract_system(cache::NodalCache) -> (K, f)
Extract global system from nodal cache.
Returns references to sparse matrix and force vector.
**Zero allocations** - no copying.
# Arguments
- `cache`: Nodal cache after assembly
# Returns
- `K`: Sparse matrix (reference, no copy)
- `f`: Force vector (reference, no copy)
"""
function extract_system(cache::NodalCache)
return cache.K, cache.f
end
src/assemblers/node_based_coo.jl
-484
View File
@@ -1,484 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
Node-based COO assembly using block integration.
**NODAL ASSEMBLY PARADIGM**: Loop over nodes, not elements!
Each node:
1. Finds all elements touching it (via inverse connectivity)
2. For each touching element:
- Prepares element geometry once (PreparedElement)
- Computes only needed 3×3 blocks (compute_block!)
3. Scatters blocks to COO triplets
# Key Differences from Element-Based Assembly
**Element-Based (traditional):**
```julia
for element in elements
K_e = compute_element_stiffness(element) # Full N×N matrix of 3×3 blocks
scatter(K_e) # Scatter all entries
end
```
**Node-Based (this file):**
```julia
for node_i in nodes
for element in elements_touching(node_i)
prepared = prepare_element(element) # Geometry preprocessing
for node_j in element.nodes
K_ij = compute_block!(prepared, i, j) # Single 3×3 block
scatter(K_ij, i, j) # Scatter one block
end
end
end
```
# Advantages
1. **GPU-friendly**: One thread per node, no race conditions
2. **Contact-ready**: Contact is naturally node-based
3. **Matrix-free ready**: Can compute K*v without forming K
4. **Cache-friendly**: Reuses PreparedElement for multiple blocks
5. **Adaptive-ready**: Easy to refine/coarsen at node level
# Performance Expectations
- **CPU Single-thread**: ~1.5-2x slower than element-based (more kernel calls)
- **CPU Multi-thread**: ~1.5-2x faster (better parallelization)
- **GPU**: ~10-50x faster (massive parallelization, no atomics needed)
# References
- Golden standard: `docs/src/book/multigpu_nodal_assembly.md`
- PreparedElement: `src/domains/continuum/integration.jl`
- Block kernel: `src/domains/continuum/kernel.jl`
# Example
```julia
# Setup
mesh = create_cantilever_mesh(50, 10, 10)
material = LinearElastic(E=210e9, ν=0.3)
kernel = ContinuumKernel(
ContinuumFormulation{FullThreeD}(),
material,
Displacement{3}()
)
# Create node-based assembler and cache
assembler = NodeBasedCOOAssembler()
cache = create_cache(assembler, mesh, kernel)
# Assemble (zero allocations after warmup!)
assemble!(cache, assembler, kernel, mesh)
# Extract system
K, f = extract_system(cache)
# Solve
apply_dirichlet_bcs!(K, f, kernel, mesh, bc_dirichlet)
u = K \\ f
```
"""
using SparseArrays
using Tensors
"""
NodeBasedCOOCache
Pre-allocated cache for node-based COO assembly.
Similar to COOCache but includes inverse connectivity mapping.
# Fields
- `I::Vector{Int}`: Row indices (COO format)
- `J::Vector{Int}`: Column indices (COO format)
- `V::Vector{Float64}`: Values (COO format)
- `f::Vector{Float64}`: Global force vector
- `counter::Ref{Int}`: Current triplet count
- `capacity::Int`: Maximum triplet capacity
- `node_to_elements::NodeToElementsMap`: Inverse connectivity
- `element_cache::ElementCache`: Cache for element operations
- `ndofs::Int`: Total DOFs in system
"""
struct NodeBasedCOOCache{T<:AbstractTopology,B<:AbstractBasis,IPS}
I::Vector{Int}
J::Vector{Int}
V::Vector{Float64}
f::Vector{Float64}
counter::Ref{Int}
capacity::Int
node_to_elements::NodeToElementsMap
element_cache::ElementCache{T,B,IPS}
ndofs::Int
end
"""
NodeBasedCOOCache(mesh::AbstractMesh, kernel::ContinuumKernel)
Create cache for node-based assembly.
Builds inverse connectivity and allocates buffers.
# Arguments
- `mesh`: Finite element mesh
- `kernel`: Continuum kernel
# Returns
- Pre-allocated node-based COO cache
"""
function NodeBasedCOOCache(mesh::AbstractMesh, kernel::ContinuumKernel)
# Build inverse connectivity
node_to_elements = NodeToElementsMap(mesh.connectivity)
# Estimate triplet count (same as element-based)
ndofs_per_node = dofs_per_node(kernel)
nnodes = length(mesh.nodes)
ndofs = ndofs_per_node * nnodes
# Estimate: For each node, sum over touching elements
# Each element contributes N blocks (N = nodes per element)
# Each block = 3×3 = 9 triplets
avg_elements_per_node = node_to_elements.nelements / nnodes
N = length(first(mesh.connectivity)) # Nodes per element
estimated_triplets = Int(ceil(1.2 * nnodes * avg_elements_per_node * N * 9))
# Allocate triplet arrays
I = Vector{Int}(undef, estimated_triplets)
J = Vector{Int}(undef, estimated_triplets)
V = Vector{Float64}(undef, estimated_triplets)
f = zeros(Float64, ndofs)
counter = Ref(0)
# Create element cache (for prepare_element! and compute_block!)
element_cache = ElementCache(mesh, kernel)
return NodeBasedCOOCache(I, J, V, f, counter, estimated_triplets,
node_to_elements, element_cache, ndofs)
end
"""
reset!(cache::NodeBasedCOOCache)
Reset cache for new assembly (zero force vector, reset counter).
Does NOT clear inverse connectivity (that's permanent structure).
"""
function reset!(cache::NodeBasedCOOCache)
fill!(cache.f, 0.0)
cache.counter[] = 0
return nothing
end
"""
assemble!(
cache::NodeBasedCOOCache,
assembler::NodeBasedCOOAssembler,
kernel::ContinuumKernel,
mesh::AbstractMesh
) -> Nothing
Assemble global system using **node-based traversal**.
# Algorithm
```julia
for node_i in 1:nnodes
# Get all elements touching this node
for elem_info in node_to_elements[node_i]
element_id = elem_info.element_id
local_i = elem_info.local_node_idx
# Prepare element geometry ONCE
prepared = prepare_element!(cache.element_cache, kernel, element_id, mesh)
# Compute blocks for all nodes in this element
for local_j in 1:N
global_j = connectivity[element_id][local_j]
# Compute single 3×3 block
K_ij = compute_block!(prepared, kernel.material, local_i, local_j)
# Scatter to triplets
scatter_block_to_triplets!(cache, K_ij, node_i, global_j)
end
end
end
```
# Key Operations
1. **prepare_element!** - Precompute geometry (Jacobian, gradients) once per element
2. **compute_block!** - Compute single 3×3 stiffness block using prepared geometry
3. **scatter_block_to_triplets!** - Add 9 triplets (i,j,value) for 3×3 block
# Zero-Allocation (After Warmup)
All arrays pre-allocated. Element preparation reuses cache buffers.
# Arguments
- `cache`: Pre-allocated node-based COO cache
- `assembler`: Node-based COO assembler
- `kernel`: Continuum kernel
- `mesh`: Finite element mesh
"""
function assemble!(
cache::NodeBasedCOOCache,
assembler::NodeBasedCOOAssembler,
kernel::ContinuumKernel,
mesh::AbstractMesh
)
# Reset cache
reset!(cache)
nnodes = length(mesh.nodes)
ndofs_per_node = dofs_per_node(kernel)
# NODAL LOOP: One iteration per node (GPU: one thread per node!)
for node_i in 1:nnodes
# Get all elements touching this node
touching_elements = cache.node_to_elements.node_to_elements[node_i]
# Loop over touching elements
for elem_info in touching_elements
element_id = elem_info.element_id
local_i = elem_info.local_node_idx # Position of node_i in element
# Prepare element geometry ONCE (reuses cache.element_cache)
prepared = prepare_element!(cache.element_cache, kernel, element_id, mesh)
# Get element connectivity
conn = mesh.connectivity[element_id]
N = length(conn) # Nodes per element
# Compute blocks for all nodes j in this element
for local_j in 1:N
global_j = conn[local_j]
# Compute single 3×3 block K[i,j]
# This is THE KEY OPERATION: block-based integration
K_ij = compute_block!(
prepared,
kernel.material,
local_i,
local_j
)
# Scatter 3×3 block to triplets (adds 9 entries)
scatter_block_to_triplets!(
cache,
K_ij,
node_i,
global_j,
ndofs_per_node
)
end
end
end
return nothing
end
"""
scatter_block_to_triplets!(
cache::NodeBasedCOOCache,
K_block::Tensor{2,3},
node_i::Int,
node_j::Int,
ndofs_per_node::Int
)
Scatter single 3×3 block to COO triplets **in-place**.
Maps block[α,β] → triplet at DOF indices:
- Row: 3*(node_i-1) + α
- Col: 3*(node_j-1) + β
- Val: K_block[α,β]
# Arguments
- `cache`: Node-based COO cache
- `K_block`: 3×3 stiffness block (Tensor{2,3})
- `node_i`: Global row node index
- `node_j`: Global column node index
- `ndofs_per_node`: DOFs per node (typically 3)
# Zero-Allocation
Writes to pre-allocated triplet arrays, updates counter.
"""
function scatter_block_to_triplets!(
cache::NodeBasedCOOCache,
K_block::Tensor{2,3,Float64},
node_i::Int,
node_j::Int,
ndofs_per_node::Int
)
counter = cache.counter[]
# Check capacity
new_triplets = ndofs_per_node * ndofs_per_node # 3×3 = 9
if counter + new_triplets > cache.capacity
error("Node-based COO cache overflow: need $(counter + new_triplets) triplets, " *
"capacity is $(cache.capacity). Increase cache size.")
end
# DOF offsets for nodes i and j
row_offset = ndofs_per_node * (node_i - 1)
col_offset = ndofs_per_node * (node_j - 1)
# Scatter 3×3 block to triplets
for β in 1:ndofs_per_node # Column (node j DOF)
j_global = col_offset + β
for α in 1:ndofs_per_node # Row (node i DOF)
i_global = row_offset + α
counter += 1
cache.I[counter] = i_global
cache.J[counter] = j_global
cache.V[counter] = K_block[α, β]
end
end
cache.counter[] = counter
return nothing
end
"""
extract_system(cache::NodeBasedCOOCache) -> (K, f)
Build sparse matrix from triplets and return system.
Calls `sparse(I, J, V)` to build CSC matrix. Duplicates are summed automatically.
# Arguments
- `cache`: Assembled node-based COO cache
# Returns
- `K::SparseMatrixCSC`: Global stiffness matrix
- `f::Vector`: Global force vector
# Allocation
Allocates sparse matrix structure (CSC format). This is the only allocation
outside cache construction.
"""
function extract_system(cache::NodeBasedCOOCache)
ntriplets = cache.counter[]
# Build sparse matrix (duplicates are summed automatically)
I_used = @view cache.I[1:ntriplets]
J_used = @view cache.J[1:ntriplets]
V_used = @view cache.V[1:ntriplets]
K = sparse(I_used, J_used, V_used, cache.ndofs, cache.ndofs)
return K, cache.f
end
# ============================================================================
# HELPER FUNCTIONS
# ============================================================================
"""
create_cache(
assembler::NodeBasedCOOAssembler,
mesh::AbstractMesh,
kernel::ContinuumKernel
) -> NodeBasedCOOCache
Create pre-allocated cache for node-based COO assembly.
Convenience function that wraps `NodeBasedCOOCache(mesh, kernel)`.
# Example
```julia
assembler = NodeBasedCOOAssembler()
cache = create_cache(assembler, mesh, kernel)
assemble!(cache, assembler, kernel, mesh)
K, f = extract_system(cache)
```
"""
function create_cache(
assembler::NodeBasedCOOAssembler,
mesh::AbstractMesh,
kernel::ContinuumKernel
)
return NodeBasedCOOCache(mesh, kernel)
end
"""
dofs_per_node(kernel::ContinuumKernel) -> Int
Return DOFs per node for continuum kernel (always 3 for displacement).
Dispatches on kernel field dimension.
"""
function dofs_per_node(kernel::ContinuumKernel{Theory,Mat}) where {Theory,Mat}
field = kernel.field
return field.dim # Displacement{3} → 3
end
# ============================================================================
# PERFORMANCE NOTES
# ============================================================================
#=
# CPU Performance Comparison (Estimated)
**Element-Based Assembly:**
- Elements: 1000 Tet4
- Operations: 1000 elements × 4×4 blocks × 3×3 entries = 48,000 block computations
- Time: ~5ms (baseline)
**Node-Based Assembly:**
- Nodes: 500 nodes
- Operations: 500 nodes × 8 elements/node × 4 blocks/element = 16,000 block computations
- But: 3× more kernel calls due to overlaps
- Time: ~7-10ms (1.5-2× slower single-threaded)
**Why slower on CPU?**
- Each block computed once in element assembly
- Each block computed 2× on average in nodal assembly (shared between 2 elements)
- More function call overhead
**Why faster on GPU?**
- Element assembly: Sequential (can't parallelize over elements efficiently)
- Nodal assembly: Massive parallelism (one thread per node)
- GPU speedup: ~10-50× depending on problem size
**Multi-threaded CPU (Threads.@threads):**
- Can parallelize outer node loop
- Expected speedup: 1.5-2× over element-based
- No race conditions (each node writes different triplets)
# Memory Comparison
**Element-Based:**
- Triplet storage: ~50 KB per 1000 elements
- Element cache: ~2 KB per thread
**Node-Based:**
- Triplet storage: Same (~50 KB)
- Element cache: ~2 KB per thread
- Inverse connectivity: ~10-20 KB (one-time)
→ Nearly identical memory usage!
# When to Use Node-Based Assembly
**Use when:**
- ✅ GPU acceleration needed
- ✅ Contact mechanics (naturally nodal)
- ✅ Matrix-free methods (K*v without forming K)
- ✅ Adaptive refinement (local node operations)
- ✅ Multi-threading on CPU
**Don't use when:**
- ❌ Single-threaded CPU only
- ❌ Simple problems (< 1000 nodes)
- ❌ Prototyping/debugging (element-based is clearer)
=#
src/assembly/assembly.jl
-201
View File
@@ -1,201 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
function isapprox(a1::Assembly, a2::Assembly)
T = isapprox(a1.K, a2.K)
T &= isapprox(a1.C1, a2.C1)
T &= isapprox(a1.C2, a2.C2)
T &= isapprox(a1.D, a2.D)
T &= isapprox(a1.f, a2.f)
T &= isapprox(a1.g, a2.g)
return T
end
function assemble_prehook!(::Problem, ::T) where T<:Number end
function assemble_posthook!(::Problem, ::T) where T<:Number end
"""
assemble_elements!(problem, assembly, elements, time)
Assemble elements for problem.
This should be overridden with own `assemble_elements!`-implementation.
"""
function assemble_elements!(problem::Problem, assembly::Assembly,
elements::Vector{T}, time) where T<:AbstractElement{E} where E
elements2 = convert(Vector{Element}, elements)
assemble!(assembly, problem, elements2, time)
end
function assemble!(problem::Problem, time)
assemble_prehook!(problem, time)
elements = get_elements(problem)
assembly = get_assembly(problem)
if !isempty(assembly)
@warn("Problem assembly is not empty before assembling. This is probably " *
"causing unexpected results. To remove old assembly, use " *
"`empty!(problem.assembly)`", typeof(problem), problem.name)
assemble_posthook!(problem, time)
return nothing
end
if isempty(elements)
@warn("There is no elements defined in problem. Before assembling a " *
"problem, elements must be added using " *
"`add_elements!(problem, elements)`.", typeof(problem), problem.name)
assemble_posthook!(problem, time)
return nothing
end
first_element = first(elements)
unknown_field_name = get_unknown_field_name(problem)
if !haskey(first_element, unknown_field_name)
#=
warn("Assembling elements for problem $(problem.name): seems that ",
"problem is uninitialized. To initialize problem, use ",
"`initialize!(problem, time)`.")
info("Initializing problem $(problem.name) at time $time automatically.")
=#
initialize!(problem, time)
end
for (element_type, elements) in group_by_element_type(elements)
assemble_elements!(problem, assembly, elements, time)
end
assemble_posthook!(problem, time)
return nothing
end
function assemble!(problem::Problem)
@warn("assemble!(problem) will be deprecated. Use assemble!(problem, time)")
assemble!(problem, 0.0)
end
function assemble_mass_matrix!(problem::Problem, time::Float64)
if !isempty(problem.assembly.M)
@info("Mass matrix for is already assembled, not assembling.",
problem.name)
return nothing
end
elements = get_elements(problem)
for (element_type, elements) in group_by_element_type(get_elements(problem))
assemble_mass_matrix!(problem::Problem, elements, time)
end
return nothing
end
function assemble_mass_matrix!(problem::Problem, elements::Vector{E}, time) where E<:AbstractElement{M_,B} where {M_,B}
nnodes = length(first(elements))
dim = get_unknown_field_dimension(problem)
M = zeros(nnodes, nnodes)
N = zeros(1, nnodes)
NtN = zeros(nnodes, nnodes)
ldofs = zeros(Int, nnodes)
for element in elements
fill!(M, 0.0)
for ip in get_integration_points(element, 2)
detJ = element(ip, time, Val{:detJ})
rho = element("density", ip, time)
w = ip.weight * rho * detJ
eval_basis!(B, N, ip)
N = element(ip, time)
mul!(NtN, transpose(N), N)
rmul!(NtN, w)
for i = 1:nnodes^2
M[i] += NtN[i]
end
end
for (i, j) in enumerate(get_connectivity(element))
@inbounds ldofs[i] = (j - 1) * dim
end
for i = 1:dim
add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, M)
end
end
return
end
# TODO (Phase 1B): Re-enable after resolving Tet10 topology vs Tet10Basis name conflict
# This specialized method assumes Tet10 <: AbstractBasis, but Tet10 is now a topology type.
# Need to refactor to use Tet10Basis or accept AbstractElement{M, T, Tet10Basis}
#=
"""
assemble_mass_matrix!(problem, elements::Vector{Element{Tet10}}, time)
Assemble Tet10 mass matrices using special method. If Tet10 has constant metric
if can be integrated analytically to gain performance.
"""
function assemble_mass_matrix!(problem::Problem, elements::Vector{E}, time) where E<:AbstractElement{M_, Tet10} where M_
nnodes = length(Tet10)
dim = get_unknown_field_dimension(problem)
M = zeros(nnodes, nnodes)
N = zeros(1, nnodes)
NtN = zeros(nnodes, nnodes)
ldofs = zeros(Int, nnodes)
M_CM = 1.0/2520.0 * [
6 1 1 1 -4 -6 -4 -4 -6 -6
1 6 1 1 -4 -4 -6 -6 -4 -6
1 1 6 1 -6 -4 -4 -6 -6 -4
1 1 1 6 -6 -6 -6 -4 -4 -4
-4 -4 -6 -6 32 16 16 16 16 8
-6 -4 -4 -6 16 32 16 8 16 16
-4 -6 -4 -6 16 16 32 16 8 16
-4 -6 -6 -4 16 8 16 32 16 16
-6 -4 -6 -4 16 16 8 16 32 16
-6 -6 -4 -4 8 16 16 16 16 32]
function is_CM(::AbstractElement{M, Tet10}, X; rtol=1.0e-6) where M
isapprox(X[5], 1/2*(X[1]+X[2]); rtol=rtol) || return false
isapprox(X[6], 1/2*(X[2]+X[3]); rtol=rtol) || return false
isapprox(X[7], 1/2*(X[3]+X[1]); rtol=rtol) || return false
isapprox(X[8], 1/2*(X[1]+X[4]); rtol=rtol) || return false
isapprox(X[9], 1/2*(X[2]+X[4]); rtol=rtol) || return false
isapprox(X[10], 1/2*(X[3]+X[4]); rtol=rtol) || return false
return true
end
n_CM = 0
for element in elements
for (i, j) in enumerate(get_connectivity(element))
@inbounds ldofs[i] = (j-1)*dim
end
X = element("geometry", time)
rho = element("density", time)
if is_CM(element, X) && length(rho) == 1
ip = (1.0/3.0, 1.0/3.0, 1.0/3.0)
detJ = element(ip, time, Val{:detJ})
rho = element("density", ip, time)
CM_s = detJ*rho
n_CM += 1
for i=1:dim
add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, CM_s * M_CM)
end
else
fill!(M, 0.0)
for ip in get_integration_points(element, 2)
detJ = element(ip, time, Val{:detJ})
rho = element("density", ip, time)
w = ip.weight*rho*detJ
eval_basis!(Tet10, N, ip)
N = element(ip, time)
mul!(NtN, transpose(N), N)
rmul!(NtN, w)
for i=1:nnodes^2
M[i] += NtN[i]
end
end
for i=1:dim
add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, M)
end
end
end
@info("$n_CM of $(length(elements)) was constant metric.")
return
end
=#
src/assembly/element_structures.jl
-341
View File
@@ -1,341 +0,0 @@
# Traditional Element Assembly
#
# This module provides the standard element-by-element assembly approach
# for comparison with nodal assembly. Builds global tangent stiffness matrix
# and residual force vector using sparse matrix formats.
using Tensors
using SparseArrays
using LinearAlgebra
"""
ElementAssemblyData{T}
Storage for element assembly using traditional (element-by-element) approach.
# Fields
- `K_global::SparseMatrixCSC{T}`: Global tangent stiffness matrix
- `r_global::Vector{T}`: Global residual force vector (r = f_int - f_ext)
- `f_int_global::Vector{T}`: Global internal force vector
- `f_ext_global::Vector{T}`: Global external force vector
- `ndof::Int`: Total number of degrees of freedom
# Notes
- Assembly uses COO (coordinate) format, then converts to CSC
- Multiple elements can write to same global DOF (summed automatically)
"""
mutable struct ElementAssemblyData{T}
K_global::SparseMatrixCSC{T,Int}
r_global::Vector{T}
f_int_global::Vector{T}
f_ext_global::Vector{T}
ndof::Int
end
"""
ElementAssemblyData(ndof::Int, ::Type{T}=Float64)
Allocate storage for traditional element assembly.
# Arguments
- `ndof`: Total degrees of freedom (nnodes × 3 for 3D)
- `T`: Floating point type (default Float64)
# Example
```julia
nnodes = 100
assembly = ElementAssemblyData(3 * nnodes, Float64)
```
"""
function ElementAssemblyData(ndof::Int, ::Type{T}=Float64) where T
# Pre-allocate empty sparse matrix (will fill during assembly)
K_global = spzeros(T, ndof, ndof)
r_global = zeros(T, ndof)
f_int_global = zeros(T, ndof)
f_ext_global = zeros(T, ndof)
return ElementAssemblyData{T}(K_global, r_global, f_int_global, f_ext_global, ndof)
end
"""
reset!(assembly::ElementAssemblyData)
Reset assembly data to zero (for incremental/iterative solvers).
"""
function reset!(assembly::ElementAssemblyData{T}) where T
assembly.K_global = spzeros(T, assembly.ndof, assembly.ndof)
fill!(assembly.r_global, 0.0)
fill!(assembly.f_int_global, 0.0)
fill!(assembly.f_ext_global, 0.0)
end
"""
ElementContribution{T}
Local element contribution before scattering to global.
# Fields
- `element_id::Int`: Element ID
- `gdofs::Vector{Int}`: Global DOF indices (e.g., [1,2,3,4,5,6,...] for nodes)
- `K_local::Matrix{T}`: Local stiffness matrix (ndofs_local × ndofs_local)
- `f_int_local::Vector{T}`: Local internal force vector
- `f_ext_local::Vector{T}`: Local external force vector
# Notes
- For Tet4: ndofs_local = 12 (4 nodes × 3 DOF)
- For Tet10: ndofs_local = 30 (10 nodes × 3 DOF)
"""
struct ElementContribution{T}
element_id::Int
gdofs::Vector{Int}
K_local::Matrix{T}
f_int_local::Vector{T}
f_ext_local::Vector{T}
end
"""
ElementContribution(element_id::Int, gdofs::Vector{Int}, ::Type{T}=Float64)
Allocate storage for element contribution.
# Arguments
- `element_id`: Element ID
- `gdofs`: Global DOF indices
- `T`: Floating point type
# Example
```julia
# Tet4 element connecting nodes [5, 7, 12, 15]
gdofs = [13,14,15, 19,20,21, 34,35,36, 43,44,45] # 3 DOF per node
contrib = ElementContribution(1, gdofs, Float64)
```
"""
function ElementContribution(element_id::Int, gdofs::Vector{Int}, ::Type{T}=Float64) where T
ndofs = length(gdofs)
K_local = zeros(T, ndofs, ndofs)
f_int_local = zeros(T, ndofs)
f_ext_local = zeros(T, ndofs)
return ElementContribution{T}(element_id, gdofs, K_local, f_int_local, f_ext_local)
end
"""
scatter_to_global!(assembly::ElementAssemblyData, contrib::ElementContribution)
Scatter element contribution to global matrices/vectors (traditional assembly).
This is the key operation in element assembly: add local element quantities
to global system. Uses COO format (accumulates into lists).
# Arguments
- `assembly`: Global assembly data
- `contrib`: Element contribution
# Notes
- Multiple elements can contribute to same global DOF (summed)
- For GPU: Would require atomic operations (slow!)
- For CPU: Direct scatter-add works fine
"""
function scatter_to_global!(assembly::ElementAssemblyData{T},
contrib::ElementContribution{T}) where T
# Scatter forces (simple vector addition)
for (local_i, global_i) in enumerate(contrib.gdofs)
assembly.f_int_global[global_i] += contrib.f_int_local[local_i]
assembly.f_ext_global[global_i] += contrib.f_ext_local[local_i]
end
# Scatter stiffness (matrix addition)
# Build list of (I, J, V) triplets for sparse matrix
I_rows = Int[]
J_cols = Int[]
values = T[]
ndofs_local = length(contrib.gdofs)
for i in 1:ndofs_local, j in 1:ndofs_local
if abs(contrib.K_local[i, j]) > 1e-14 # Skip near-zeros
push!(I_rows, contrib.gdofs[i])
push!(J_cols, contrib.gdofs[j])
push!(values, contrib.K_local[i, j])
end
end
# Add to existing sparse matrix
K_elem = sparse(I_rows, J_cols, values, assembly.ndof, assembly.ndof)
assembly.K_global += K_elem
end
"""
compute_residual!(assembly::ElementAssemblyData)
Compute residual force vector: r = f_int - f_ext
Should be called after all elements have been assembled.
"""
function compute_residual!(assembly::ElementAssemblyData{T}) where T
assembly.r_global .= assembly.f_int_global .- assembly.f_ext_global
end
"""
assemble_elements!(assembly::ElementAssemblyData,
contributions::Vector{ElementContribution})
Assemble all element contributions to global system.
# Arguments
- `assembly`: Global assembly data (modified in-place)
- `contributions`: Vector of element contributions
# Example
```julia
assembly = ElementAssemblyData(ndof)
contributions = compute_all_element_contributions(elements, u, time)
assemble_elements!(assembly, contributions)
compute_residual!(assembly)
# Now solve: K_global * Δu = -r_global
```
"""
function assemble_elements!(assembly::ElementAssemblyData{T},
contributions::Vector{ElementContribution{T}}) where T
reset!(assembly)
# Loop over elements and scatter (element assembly)
for contrib in contributions
scatter_to_global!(assembly, contrib)
end
# Compute residual
compute_residual!(assembly)
end
"""
apply_dirichlet_bc!(assembly::ElementAssemblyData,
fixed_dofs::Vector{Int},
prescribed_values::Vector{T}=zeros(length(fixed_dofs)))
Apply Dirichlet (essential) boundary conditions by penalty method.
# Arguments
- `assembly`: Global assembly data (modified in-place)
- `fixed_dofs`: DOF indices to fix
- `prescribed_values`: Prescribed displacement values (default: zeros)
# Method
Uses penalty method: adds large stiffness to diagonal and corresponding RHS.
For DOF i with prescribed value u_prescribed:
- K[i,i] += penalty (e.g., 1e10 * max_K)
- r[i] = penalty * (u_current - u_prescribed)
# Example
```julia
# Fix nodes 1 and 2 in all directions (zero displacement)
fixed_dofs = [1,2,3, 4,5,6] # Nodes 1,2 × 3 DOF
apply_dirichlet_bc!(assembly, fixed_dofs)
```
"""
function apply_dirichlet_bc!(assembly::ElementAssemblyData{T},
fixed_dofs::Vector{Int},
prescribed_values::Vector{T}=zeros(T, length(fixed_dofs))) where T
# Penalty parameter (large relative to stiffness)
max_K = maximum(abs, assembly.K_global)
penalty = 1e10 * max_K
for (idx, dof) in enumerate(fixed_dofs)
# Add penalty stiffness to diagonal
assembly.K_global[dof, dof] += penalty
# Modify residual (assuming current displacement is zero for now)
# In full Newton: r[i] += penalty * (u_current[i] - u_prescribed[i])
assembly.r_global[dof] = penalty * prescribed_values[idx]
end
end
"""
get_dof_indices(connectivity::NTuple{N,Int}, dim::Int=3) -> Vector{Int}
Get global DOF indices for an element given node connectivity.
# Arguments
- `connectivity`: Element node IDs (e.g., (5, 7, 12, 15) for Tet4)
- `dim`: Dimension (3 for 3D elasticity)
# Returns
- `gdofs::Vector{Int}`: Global DOF indices
# Example
```julia
# Element with nodes [5, 7, 12, 15]
gdofs = get_dof_indices((5, 7, 12, 15), 3)
# Returns: [13,14,15, 19,20,21, 34,35,36, 43,44,45]
```
"""
function get_dof_indices(connectivity::NTuple{N,Int}, dim::Int=3) where N
nnodes = length(connectivity)
gdofs = zeros(Int, dim * nnodes)
for (local_i, global_node) in enumerate(connectivity)
for d in 1:dim
gdofs[dim*(local_i-1)+d] = dim * (global_node - 1) + d
end
end
return gdofs
end
"""
matrix_vector_product(assembly::ElementAssemblyData, v::Vector{T}) -> Vector{T}
Compute matrix-vector product: w = K * v using assembled sparse matrix.
# Arguments
- `assembly`: Assembly data (contains K_global)
- `v`: Input vector (ndof)
# Returns
- `w`: Output vector w = K * v
# Example
```julia
# GMRES matrix-free operator
function matvec(v)
return matrix_vector_product(assembly, v)
end
Δu = gmres(matvec, -r, tol=1e-6)
```
"""
function matrix_vector_product(assembly::ElementAssemblyData{T}, v::Vector{T}) where T
return assembly.K_global * v
end
"""
print_assembly_stats(assembly::ElementAssemblyData)
Print statistics about assembled system (for debugging).
"""
function print_assembly_stats(assembly::ElementAssemblyData)
nnz_K = nnz(assembly.K_global)
ndof = assembly.ndof
fill_ratio = nnz_K / (ndof * ndof)
println("="^60)
println("Traditional Element Assembly Statistics")
println("="^60)
println(" Total DOF: ", ndof)
println(" K matrix size: ", size(assembly.K_global))
println(" K non-zeros: ", nnz_K)
println(" K fill ratio: ", round(fill_ratio, sigdigits=2))
println(" K memory (MB): ", round(nnz_K * 16 / 1024^2, digits=2))
println(" ||f_int||: ", round(norm(assembly.f_int_global), sigdigits=2))
println(" ||f_ext||: ", round(norm(assembly.f_ext_global), sigdigits=2))
println(" ||residual||: ", round(norm(assembly.r_global), sigdigits=2))
println(" K symmetric: ", issymmetric(assembly.K_global))
println("="^60)
end
# Export main types and functions
export ElementAssemblyData, ElementContribution
export reset!, scatter_to_global!, compute_residual!
export assemble_elements!, apply_dirichlet_bc!
export get_dof_indices, matrix_vector_product
export print_assembly_stats
src/assembly/framework.jl
-201
View File
@@ -1,201 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
function isapprox(a1::Assembly, a2::Assembly)
T = isapprox(a1.K, a2.K)
T &= isapprox(a1.C1, a2.C1)
T &= isapprox(a1.C2, a2.C2)
T &= isapprox(a1.D, a2.D)
T &= isapprox(a1.f, a2.f)
T &= isapprox(a1.g, a2.g)
return T
end
function assemble_prehook!(::Problem, ::T) where T<:Number end
function assemble_posthook!(::Problem, ::T) where T<:Number end
"""
assemble_elements!(problem, assembly, elements, time)
Assemble elements for problem.
This should be overridden with own `assemble_elements!`-implementation.
"""
function assemble_elements!(problem::Problem, assembly::Assembly,
elements::Vector{T}, time) where T<:AbstractElement{E} where E
elements2 = convert(Vector{Element}, elements)
assemble!(assembly, problem, elements2, time)
end
function assemble!(problem::Problem, time)
assemble_prehook!(problem, time)
elements = get_elements(problem)
assembly = get_assembly(problem)
if !isempty(assembly)
@warn("Problem assembly is not empty before assembling. This is probably " *
"causing unexpected results. To remove old assembly, use " *
"`empty!(problem.assembly)`", typeof(problem), problem.name)
assemble_posthook!(problem, time)
return nothing
end
if isempty(elements)
@warn("There is no elements defined in problem. Before assembling a " *
"problem, elements must be added using " *
"`add_elements!(problem, elements)`.", typeof(problem), problem.name)
assemble_posthook!(problem, time)
return nothing
end
first_element = first(elements)
unknown_field_name = get_unknown_field_name(problem)
if !haskey(first_element, unknown_field_name)
#=
warn("Assembling elements for problem $(problem.name): seems that ",
"problem is uninitialized. To initialize problem, use ",
"`initialize!(problem, time)`.")
info("Initializing problem $(problem.name) at time $time automatically.")
=#
initialize!(problem, time)
end
for (element_type, elements) in group_by_element_type(elements)
assemble_elements!(problem, assembly, elements, time)
end
assemble_posthook!(problem, time)
return nothing
end
function assemble!(problem::Problem)
@warn("assemble!(problem) will be deprecated. Use assemble!(problem, time)")
assemble!(problem, 0.0)
end
function assemble_mass_matrix!(problem::Problem, time::Float64)
if !isempty(problem.assembly.M)
@info("Mass matrix for is already assembled, not assembling.",
problem.name)
return nothing
end
elements = get_elements(problem)
for (element_type, elements) in group_by_element_type(get_elements(problem))
assemble_mass_matrix!(problem::Problem, elements, time)
end
return nothing
end
function assemble_mass_matrix!(problem::Problem, elements::Vector{E}, time) where E<:AbstractElement{M_,B} where {M_,B}
nnodes = length(first(elements))
dim = get_unknown_field_dimension(problem)
M = zeros(nnodes, nnodes)
N = zeros(1, nnodes)
NtN = zeros(nnodes, nnodes)
ldofs = zeros(Int, nnodes)
for element in elements
fill!(M, 0.0)
for ip in get_integration_points(element, 2)
detJ = element(ip, time, Val{:detJ})
rho = element("density", ip, time)
w = ip.weight * rho * detJ
eval_basis!(B, N, ip)
N = element(ip, time)
mul!(NtN, transpose(N), N)
rmul!(NtN, w)
for i = 1:nnodes^2
M[i] += NtN[i]
end
end
for (i, j) in enumerate(get_connectivity(element))
@inbounds ldofs[i] = (j - 1) * dim
end
for i = 1:dim
add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, M)
end
end
return
end
# TODO (Phase 1B): Re-enable after resolving Tet10 topology vs Tet10Basis name conflict
# This specialized method assumes Tet10 <: AbstractBasis, but Tet10 is now a topology type.
# Need to refactor to use Tet10Basis or accept AbstractElement{M, T, Tet10Basis}
#=
"""
assemble_mass_matrix!(problem, elements::Vector{Element{Tet10}}, time)
Assemble Tet10 mass matrices using special method. If Tet10 has constant metric
if can be integrated analytically to gain performance.
"""
function assemble_mass_matrix!(problem::Problem, elements::Vector{E}, time) where E<:AbstractElement{M_, Tet10} where M_
nnodes = length(Tet10)
dim = get_unknown_field_dimension(problem)
M = zeros(nnodes, nnodes)
N = zeros(1, nnodes)
NtN = zeros(nnodes, nnodes)
ldofs = zeros(Int, nnodes)
M_CM = 1.0/2520.0 * [
6 1 1 1 -4 -6 -4 -4 -6 -6
1 6 1 1 -4 -4 -6 -6 -4 -6
1 1 6 1 -6 -4 -4 -6 -6 -4
1 1 1 6 -6 -6 -6 -4 -4 -4
-4 -4 -6 -6 32 16 16 16 16 8
-6 -4 -4 -6 16 32 16 8 16 16
-4 -6 -4 -6 16 16 32 16 8 16
-4 -6 -6 -4 16 8 16 32 16 16
-6 -4 -6 -4 16 16 8 16 32 16
-6 -6 -4 -4 8 16 16 16 16 32]
function is_CM(::AbstractElement{M, Tet10}, X; rtol=1.0e-6) where M
isapprox(X[5], 1/2*(X[1]+X[2]); rtol=rtol) || return false
isapprox(X[6], 1/2*(X[2]+X[3]); rtol=rtol) || return false
isapprox(X[7], 1/2*(X[3]+X[1]); rtol=rtol) || return false
isapprox(X[8], 1/2*(X[1]+X[4]); rtol=rtol) || return false
isapprox(X[9], 1/2*(X[2]+X[4]); rtol=rtol) || return false
isapprox(X[10], 1/2*(X[3]+X[4]); rtol=rtol) || return false
return true
end
n_CM = 0
for element in elements
for (i, j) in enumerate(get_connectivity(element))
@inbounds ldofs[i] = (j-1)*dim
end
X = element("geometry", time)
rho = element("density", time)
if is_CM(element, X) && length(rho) == 1
ip = (1.0/3.0, 1.0/3.0, 1.0/3.0)
detJ = element(ip, time, Val{:detJ})
rho = element("density", ip, time)
CM_s = detJ*rho
n_CM += 1
for i=1:dim
add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, CM_s * M_CM)
end
else
fill!(M, 0.0)
for ip in get_integration_points(element, 2)
detJ = element(ip, time, Val{:detJ})
rho = element("density", ip, time)
w = ip.weight*rho*detJ
eval_basis!(Tet10, N, ip)
N = element(ip, time)
mul!(NtN, transpose(N), N)
rmul!(NtN, w)
for i=1:nnodes^2
M[i] += NtN[i]
end
end
for i=1:dim
add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, M)
end
end
end
@info("$n_CM of $(length(elements)) was constant metric.")
return
end
=#
src/assembly/nodal_structures.jl
-234
View File
@@ -1,234 +0,0 @@
# Nodal Assembly Data Structures
#
# This module provides the inverse mapping needed for efficient nodal assembly:
# Given a node, find all elements touching it and the local node index within each element.
using Tensors
"""
ElementNodeInfo
Information about how a node appears in an element.
# Fields
- `element_id::Int`: Global element ID
- `local_node_idx::Int`: Local node index within the element (1-based)
"""
struct ElementNodeInfo
element_id::Int
local_node_idx::Int
end
"""
NodeToElementsMap
Inverse connectivity mapping: for each node, lists all elements touching it.
# Fields
- `node_to_elements::Vector{Vector{ElementNodeInfo}}`: For node j, gives all elements touching it
- `nnodes::Int`: Total number of nodes in mesh
- `nelements::Int`: Total number of elements in mesh
# Example
```julia
map = NodeToElementsMap(connectivity)
# Get all elements touching node 5
elements_touching_5 = map.node_to_elements[5]
for info in elements_touching_5
println("Node 5 is local node ", info.local_node_idx, " in element ", info.element_id)
end
```
"""
struct NodeToElementsMap
node_to_elements::Vector{Vector{ElementNodeInfo}}
nnodes::Int
nelements::Int
end
"""
NodeToElementsMap(connectivity::Vector{NTuple{N,Int}}) where N
Build inverse mapping from element connectivity.
# Arguments
- `connectivity`: Vector of element connectivity tuples, e.g., [(1,2,3,4), (2,3,5,6), ...]
# Returns
- `NodeToElementsMap`: Inverse mapping structure
# Example
```julia
# Tet4 mesh with 2 elements
connectivity = [(1,2,3,4), (2,3,4,5)]
map = NodeToElementsMap(connectivity)
# Node 2 appears in both elements
@assert length(map.node_to_elements[2]) == 2
```
"""
function NodeToElementsMap(connectivity::Vector{NTuple{N,Int}}) where N
nelements = length(connectivity)
# Find maximum node ID to determine array size
nnodes = maximum(maximum(conn) for conn in connectivity)
# Pre-allocate vectors for each node
node_to_elements = [Vector{ElementNodeInfo}() for _ in 1:nnodes]
# Build inverse mapping
for (elem_id, conn) in enumerate(connectivity)
for (local_idx, global_node_id) in enumerate(conn)
push!(node_to_elements[global_node_id],
ElementNodeInfo(elem_id, local_idx))
end
end
return NodeToElementsMap(node_to_elements, nnodes, nelements)
end
"""
get_node_spider(map::NodeToElementsMap, node_id::Int) -> Vector{Int}
Get the "spider" of a node - all nodes that couple with it (including itself).
This is the union of all nodes in elements touching `node_id`. These are exactly
the nodes for which we need to compute 3×3 stiffness blocks.
# Arguments
- `map`: Node-to-elements mapping
- `node_id`: Node for which to find the spider
# Returns
- `spider_nodes::Vector{Int}`: Sorted unique list of node IDs in the spider
# Example
```julia
# For node j, find all nodes it couples with
spider = get_node_spider(map, j)
# Now compute K_blocks[k] for each k in spider
```
"""
function get_node_spider(map::NodeToElementsMap, node_id::Int,
connectivity::Vector{NTuple{N,Int}}) where N
spider = Set{Int}()
# For each element touching this node
for elem_info in map.node_to_elements[node_id]
# Add all nodes in that element
for node in connectivity[elem_info.element_id]
push!(spider, node)
end
end
return sort(collect(spider))
end
"""
NodalStiffnessContribution{T}
Storage for nodal assembly contribution at a single node.
# Fields
- `node_id::Int`: Global node ID
- `spider_nodes::Vector{Int}`: Node IDs that couple with this node
- `K_blocks::Vector{Tensor{2,3,T}}`: 3×3 stiffness blocks for each spider node
- `f_int::Vec{3,T}`: Internal force at this node
- `f_ext::Vec{3,T}`: External force at this node
# Notes
- `K_blocks[k]` corresponds to `spider_nodes[k]`
- Diagonal block (self-coupling) is included in spider
- All quantities use Tensors.jl types (zero-allocation)
"""
struct NodalStiffnessContribution{T}
node_id::Int
spider_nodes::Vector{Int}
K_blocks::Vector{Tensor{2,3,T,9}}
f_int::Vec{3,T}
f_ext::Vec{3,T}
end
"""
NodalStiffnessContribution(node_id::Int, spider_nodes::Vector{Int}, ::Type{T}=Float64)
Allocate storage for nodal assembly contribution.
# Example
```julia
spider = get_node_spider(map, 5, connectivity)
contrib = NodalStiffnessContribution(5, spider, Float64)
# Now fill in K_blocks, f_int, f_ext during assembly
```
"""
function NodalStiffnessContribution(node_id::Int, spider_nodes::Vector{Int},
::Type{T}=Float64) where T
nspider = length(spider_nodes)
K_blocks = [zero(Tensor{2,3,T}) for _ in 1:nspider]
f_int = zero(Vec{3,T})
f_ext = zero(Vec{3,T})
return NodalStiffnessContribution{T}(node_id, spider_nodes, K_blocks, f_int, f_ext)
end
"""
matrix_vector_product_nodal(contrib::NodalStiffnessContribution,
u::Vector{Vec{3,T}}) -> Vec{3,T}
Compute the matrix-vector product for one node using nodal assembly.
This computes: w_i = sum_j K_ij * u_j for node i
# Arguments
- `contrib`: Nodal stiffness contribution (contains K_blocks for all j in spider)
- `u`: Displacement field at all nodes (Vec{3} per node)
# Returns
- `w_i::Vec{3}`: Result of K_i * u at this node
# Example
```julia
# Assemble contribution for node i
contrib = assemble_nodal_contribution(element_set, node_i, u, time)
# Matrix-free matvec: w_i = K_i * u
w_i = matrix_vector_product_nodal(contrib, u)
```
"""
function matrix_vector_product_nodal(contrib::NodalStiffnessContribution{T},
u::Vector{Vec{3,T}}) where T
w = zero(Vec{3,T})
# Loop over spider nodes (only non-zero columns)
for (k, node_j) in enumerate(contrib.spider_nodes)
K_ij = contrib.K_blocks[k] # 3×3 block
u_j = u[node_j] # 3×1 displacement
# Block matrix-vector product: K_ij is Tensor{2,3}, u_j is Vec{3}
# Use regular matrix-vector multiplication (single contraction)
w += K_ij ⋅ u_j # Tensor{2,3} ⋅ Vec{3} → Vec{3}
end
return w
end
"""
print_spider_info(map::NodeToElementsMap, node_id::Int,
connectivity::Vector{NTuple{N,Int}}) where N
Print diagnostic information about a node's spider for debugging.
"""
function print_spider_info(map::NodeToElementsMap, node_id::Int,
connectivity::Vector{NTuple{N,Int}}) where N
println("Node $node_id Spider Analysis:")
println(" Touches $(length(map.node_to_elements[node_id])) elements")
for elem_info in map.node_to_elements[node_id]
println(" Element $(elem_info.element_id): local node $(elem_info.local_node_idx)")
println(" Connectivity: $(connectivity[elem_info.element_id])")
end
spider = get_node_spider(map, node_id, connectivity)
println(" Spider has $(length(spider)) nodes: $spider")
println(" → Need to compute $(length(spider)) 3×3 blocks")
println(" → Diagonal block at node $node_id")
end
src/assembly/problems.jl
-478
View File
@@ -1,478 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
abstract type AbstractProblem end
abstract type FieldProblem <: AbstractProblem end
abstract type BoundaryProblem <: AbstractProblem end
abstract type MixedProblem <: AbstractProblem end
"""
General linearized problem to solve
(K₁+K₂)Δu + C1'*Δλ = f₁+f₂
C2Δu + D*Δλ = g
"""
mutable struct Assembly
M::SparseMatrixCOO # mass matrix
# for field assembly
K::SparseMatrixCOO # stiffness matrix
Kg::SparseMatrixCOO # geometric stiffness matrix
f::SparseMatrixCOO # force vector
fg::SparseMatrixCOO #
# for boundary assembly
C1::SparseMatrixCOO
C2::SparseMatrixCOO
D::SparseMatrixCOO
g::SparseMatrixCOO
c::SparseMatrixCOO
u::Vector{Float64} # solution vector u
u_prev::Vector{Float64} # previous solution vector u
u_norm_change::Real # change of norm in u
la::Vector{Float64} # solution vector la
la_prev::Vector{Float64} # previous solution vector u
la_norm_change::Real # change of norm in la
removed_dofs::Vector{Int} # manually remove dofs from assembly
end
function Assembly()
return Assembly(
SparseMatrixCOO(),
SparseMatrixCOO(),
SparseMatrixCOO(),
SparseMatrixCOO(),
SparseMatrixCOO(),
SparseMatrixCOO(),
SparseMatrixCOO(),
SparseMatrixCOO(),
SparseMatrixCOO(),
SparseMatrixCOO(),
[], [], Inf,
[], [], Inf,
[])
end
function empty!(assembly::Assembly)
empty!(assembly.M)
empty!(assembly.K)
empty!(assembly.Kg)
empty!(assembly.f)
empty!(assembly.fg)
empty!(assembly.C1)
empty!(assembly.C2)
empty!(assembly.D)
empty!(assembly.g)
empty!(assembly.c)
end
function Base.isempty(assembly::Assembly)
T = Base.isempty(assembly.M)
T &= Base.isempty(assembly.K)
T &= Base.isempty(assembly.Kg)
T &= Base.isempty(assembly.f)
T &= Base.isempty(assembly.fg)
T &= Base.isempty(assembly.C1)
T &= Base.isempty(assembly.C2)
T &= Base.isempty(assembly.D)
T &= isempty(assembly.g)
T &= isempty(assembly.c)
return T
end
"""
Problem{P<:AbstractProblem}
Defines a new problem of type `P`, where `P` characterizes the physics of the
problem. `P` can be for example `Elasticity`, if the physics of the system is
described by Cauchy's stress equation ∇⋅σ + b = ̈ρu, or `Heat`, if the physics
of the problem is described by heat equation -∇⋅(k∇u) = f.
"""
mutable struct Problem{P<:AbstractProblem}
name::AbstractString # descriptive name for the problem
dimension::Int # degrees of freedom per node
parent_field_name::AbstractString # (optional) name of the parent field e.g. "displacement"
elements::Vector{Element}
dofmap::Dict{Element,Vector{Int}} # connects the element local dofs to the global dofs
assembly::Assembly
fields::Dict{String,AbstractField}
postprocess_fields::Vector{String}
properties::P
end
"""
Problem(problem_type, problem_name, problem_dimension)
Construct a new field problem.
`problem_type` must be a subtype of `FieldProblem` (`Elasticity`, `Heat`, etc..).
`problem_dimensions` is the number of degrees of freedom each node is containing.
# Examples
To create vector-valued elasticity problem, having 3 dofs / node:
```julia
problem1 = Problem(Elasticity, "test problem", 3)
```
To create scalar-valued Poisson problem:
```julia
problem2 = Problem(Heat, "test problem 2", 1)
```
"""
function Problem(::Type{P}, name::AbstractString, dimension::Int) where P<:FieldProblem
parent_field_name = "none"
elements = []
dofmap = Dict()
assembly = Assembly()
fields = Dict()
postprocess_fields = Vector()
properties = P()
problem = Problem{P}(name, dimension, parent_field_name, elements, dofmap,
assembly, fields, postprocess_fields, properties)
@info("Creating a new problem of type $P, having name `$name` and " *
"dimension $dimension dofs/node.")
return problem
end
"""
Problem(problem_type, problem_name, problem_dimension, parent_field_name)
Construct a new boundary problem.
`problem_type` must be a subtype of `BoundaryProblem` (`Dirichlet`, `Contact`,
etc..). `problem_dimensions` is the number of degrees of freedom each node is
containing. `parent_field_name` is describing the field, where the boundary
problem is affecting.
# Examples
To create a Dirichlet boundary condition for a vector-valued elasticity problem,
having 3 dofs / node:
```julia
bc1 = Problem(Dirichlet, "fix displacement on support", 3, "displacement")
```
To create a Dirichlet boundary condition for scalar-valued Poisson problem:
```julia
bc2 = Problem(Dirichlet, "fix surface temperature", 1, "temperature")
```
"""
function Problem(::Type{P}, name, dimension, parent_field_name) where P<:BoundaryProblem
elements = []
dofmap = Dict()
assembly = Assembly()
fields = Dict()
postprocess_fields = Vector()
properties = P()
problem = Problem{P}(name, dimension, parent_field_name, elements, dofmap,
assembly, fields, postprocess_fields, properties)
@info("Creating a new boundary problem of type $P, having name `$name` and " *
"dimension $dimension dofs/node. This boundary problems fixes field " *
"`$parent_field_name`.")
return problem
end
function get_formulation_type(::Problem)
return :incremental
end
"""
get_unknown_field_dimension(problem)
Return the dimension of the unknown field of this problem.
"""
function get_unknown_field_dimension(problem::Problem)
return problem.dimension
end
"""
get_unknown_field_name(problem)
Default function if unknown field name is not defined for some problem.
"""
function get_unknown_field_name(::P) where P<:AbstractProblem
@warn("The name of unknown field (e.g. displacement, temperature, ...) of the " *
"problem type must be given by defining a function " *
"`get_unknown_field_name(::$P)`")
return "N/A"
end
""" Return the name of the unknown field of this problem. """
function get_unknown_field_name(problem::Problem{P}) where P
return get_unknown_field_name(problem.properties)
end
""" Return the name of the parent field of this (boundary) problem. """
function get_parent_field_name(problem::Problem{P}) where P<:BoundaryProblem
return problem.parent_field_name
end
function get_unknown_field_name(::P) where P<:BoundaryProblem
return "lambda"
end
is_field_problem(::Problem) = false
is_field_problem(::Problem{P}) where {P<:FieldProblem} = true
is_boundary_problem(::Problem) = false
is_boundary_problem(::Problem{P}) where {P<:BoundaryProblem} = true
function get_elements(problem::Problem)
return problem.elements
end
function update!(problem::P, attr::Pair{String,String}...) where P<:AbstractProblem
for (name, value) in attr
setfield!(problem, Meta.parse(name), Meta.parse(value))
end
end
"""
function initialize!(problem_type, element_name, time)
Initialize the element ready for calculation, where `problem_type` is the type
of the problem (Elasticity, Dirichlet, etc.), `element_name` is the name of a
constructed element (see Element(element_type, connectivity_vector)) and `time`
is the starting time of the initializing process.
"""
function initialize!(problem::Problem, element::AbstractElement, time::Float64)
field_name = get_unknown_field_name(problem)
field_dim = get_unknown_field_dimension(problem)
nnodes = length(element)
if field_dim == 1 # scalar field
empty_field = tuple(zeros(nnodes)...)
else # vector field
# FIXME: the most effective way to do
# ([0.0,0.0], [0.0,0.0], ..., [0.0,0.0]) ?
empty_field = tuple(map((x) -> zeros(field_dim) * x, 1:nnodes)...)
end
# initialize primary field
if !haskey(element, field_name)
update!(element, field_name, time => empty_field)
end
# if a boundary problem, initialize also a field for the main problem
is_boundary_problem(problem) || return
field_name = get_parent_field_name(problem)
if !haskey(element, field_name)
update!(element, field_name, time => empty_field)
end
end
function initialize!(problem::Problem, time::Float64=0.0)
for element in get_elements(problem)
initialize!(problem, element, time)
end
end
function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector)
# resize & fill with zeros vectors if length mismatch with current solution
if length(u) != length(assembly.u)
resize!(assembly.u, length(u))
fill!(assembly.u, 0.0)
end
if length(la) != length(assembly.la)
resize!(assembly.la, length(la))
fill!(assembly.la, 0.0)
end
# copy current solutions to previous ones and add/replace new solution
# TODO: here we have couple of options and they need to be clarified
# for total formulation we are solving total quantity Ku = f while in
# incremental formulation we solve KΔu = f and u = u + Δu
assembly.u_prev = copy(assembly.u)
assembly.la_prev = copy(assembly.la)
if get_formulation_type(problem) == :total
assembly.u = u
assembly.la = la
elseif get_formulation_type(problem) == :incremental
assembly.u += u
assembly.la = la
elseif get_formulation_type(problem) == :forwarddiff
assembly.u += u
assembly.la += la
else
@info("$(problem.name): unknown formulation type, don't know what to do with results")
error("serious failure with problem formulation: $(get_formulation_type(problem))")
end
# calculate change of norm
assembly.u_norm_change = norm(assembly.u - assembly.u_prev)
assembly.la_norm_change = norm(assembly.la - assembly.la_prev)
return assembly.u, assembly.la
end
"""
get_global_solution(problem, assembly)
Return a global solution (u, la) for a problem.
Notes
-----
If the length of solution vector != number of nodes, i.e. the field dimension is
something else than 1, reshape vectors so that their length matches to the
number of nodes. This helps to get nodal results easily.
"""
function get_global_solution(problem::Problem, assembly::Assembly)
u = assembly.u
la = assembly.la
field_dim = get_unknown_field_dimension(problem)
if field_dim == 1
return u, la
else
nnodes = round(Int, length(u) / field_dim)
u = reshape(u, field_dim, nnodes)
u = Vector{Float64}[u[:, i] for i in 1:nnodes]
la = reshape(la, field_dim, nnodes)
la = Vector{Float64}[la[:, i] for i in 1:nnodes]
return u, la
end
end
function update!(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64) where P<:FieldProblem
u, la = get_global_solution(problem, assembly)
field_name = get_unknown_field_name(problem)
# update solution u for elements
for element in elements
connectivity = get_connectivity(element)
update!(element, field_name, time => tuple(u[connectivity]...))
end
end
function update!(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64) where P<:BoundaryProblem
u, la = get_global_solution(problem, assembly)
parent_field_name = get_parent_field_name(problem) # displacement
field_name = get_unknown_field_name(problem) # lambda
# update solution and lagrange multipliers for boundary elements
for element in elements
connectivity = get_connectivity(element)
update!(element, parent_field_name, time => tuple(u[connectivity]...))
update!(element, field_name, time => tuple(la[connectivity]...))
end
end
"""
add_element!(problem, element1, element2, ...)
Add element(s) to the problem.
"""
function add_element!(problem, elements...)
for element in elements
push!(problem.elements, element)
end
return nothing
end
"""
add_elements!(problem, element_set_1, element_set_2, ...)
Add vectors/tuples of element(s) to the problem.
"""
function add_elements!(problem, element_sets::Union{Vector,Tuple}...)
for elements in element_sets
nelements = length(elements)
@info("Adding $nelements elements to problem `$(problem.name)`")
add_element!(problem, elements...)
end
return nothing
end
add_elements!(problem, elements::Element...) = add_element!(problem, elements...)
function add_elements!(problem, elements_or_lists_of_elements...)
for item in elements_or_lists_of_elements
add_elements!(problem, item)
end
end
get_assembly(problem::Problem) = problem.assembly
Base.length(problem::Problem) = length(problem.elements)
function update!(problem::Problem, field_name::AbstractString, data)
#if haskey(problem.fields, field_name)
# update!(problem.fields[field_name], field_name::AbstractString, data)
#else
# problem.fields[field_name] = Field(data)
#end
update!(problem.elements, field_name::AbstractString, data)
end
function haskey(problem::Problem, field_name::AbstractString)
return haskey(problem.fields, field_name)
end
function getindex(problem::Problem, field_name::String)
return problem.fields[field_name]
end
#""" Return field calculated to nodal points for elements in problem p. """
function (problem::Problem)(field_name::String, time::Float64)
#if haskey(problem, field_name)
# return problem[field_name](time)
#end
f = Dict{Int,Any}()
for element in get_elements(problem)
haskey(element, field_name) || continue
for (c, v) in zip(get_connectivity(element), element(field_name, time))
if haskey(f, c)
if !isapprox(f[c], v)
@info("several values for single node when returning field $field_name")
@info("already have: $(f[c]), and trying to set $v")
end
else
f[c] = v
end
end
end
#f == nothing && return f
#update!(problem, field_name, time => f)
return f
end
# REMOVED: push!(problem, elements) - Use add_elements!(problem, elements) instead
# This violated Julia semantics (push! should be for collections, not domain logic)
# The modern Physics API uses add_elements!, add_dirichlet!, add_neumann!
"""
set_gdofs!(problem, element)
Set element global degrees of freedom.
"""
function set_gdofs!(problem, element, dofs)
problem.dofmap[element] = dofs
end
"""
get_gdofs(problem, element)
Return the global degrees of freedom for element.
First make lookup from problem dofmap. If not defined there, make implicit
assumption that dofs follow formula `gdofs = [dim*(nid-1)+j for j=1:dim]`,
where `nid` is node id and `dim` is the dimension of problem. This formula
arranges dofs so that first comes all dofs of node 1, then node 2 and so on:
(u11, u12, u13, u21, u22, u23, ..., un1, un2, un3) for 3 dofs/node setting.
"""
function get_gdofs(problem::Problem, element::AbstractElement)
if haskey(problem.dofmap, element)
return problem.dofmap[element]
end
conn = get_connectivity(element)
if length(conn) == 0
error("element connectivity not defined, cannot determine global ",
"degrees of freedom for element #: $(element.id)")
end
dim = get_unknown_field_dimension(problem)
gdofs = [dim * (i - 1) + j for i in conn for j = 1:dim]
return gdofs
end
src/beams/api.jl
-98
View File
@@ -1,98 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
Beam formulation API definitions.
This file defines beam-specific abstract types and formulation theories.
Must be included after core api.jl.
"""
# ============================================================================
# BEAM FORMULATION THEORIES
# ============================================================================
"""
AbstractBeamTheory
Abstract type for beam theory variants.
Beam theories differ in how they model shear deformation and cross-section kinematics.
# Concrete Theories
- `EulerBernoulli`: Classical beam theory (no shear deformation)
- `Timoshenko`: Includes shear deformation (thick beams)
# See Also
- [`BeamFormulation`](@ref)
"""
abstract type AbstractBeamTheory end
"""
EulerBernoulli <: AbstractBeamTheory
Euler-Bernoulli beam theory (classical, no shear deformation).
Assumptions:
- Plane sections remain plane and perpendicular to neutral axis
- No transverse shear deformation
- Valid for slender beams (L/h > 10)
# Usage
```julia
formulation = BeamFormulation{EulerBernoulli}()
physics = Physics(
formulation=formulation,
field=DisplacementRotation{3}(),
mesh=beam_mesh,
material=steel
)
```
"""
struct EulerBernoulli <: AbstractBeamTheory end
"""
Timoshenko <: AbstractBeamTheory
Timoshenko beam theory (includes shear deformation).
Assumptions:
- Plane sections remain plane but NOT perpendicular to neutral axis
- Transverse shear deformation included
- Valid for thick beams and higher frequencies
# Usage
```julia
formulation = BeamFormulation{Timoshenko}()
physics = Physics(
formulation=formulation,
field=DisplacementRotation{3}(),
mesh=beam_mesh,
material=steel
)
```
"""
struct Timoshenko <: AbstractBeamTheory end
"""
BeamFormulation{Theory<:AbstractBeamTheory} <: AbstractFormulation
Beam element formulation with theory variant.
# Type Parameter
- `Theory`: Beam theory type (EulerBernoulli or Timoshenko)
# Examples
```julia
# Slender beam (classical theory)
BeamFormulation{EulerBernoulli}()
# Thick beam (includes shear)
BeamFormulation{Timoshenko}()
```
# Fields per Node
Typically 6 DOFs in 3D: (ux, uy, uz, θx, θy, θz)
Use with `DisplacementRotation{3}` field type.
"""
struct BeamFormulation{Theory<:AbstractBeamTheory} <: AbstractFormulation end
src/domains/continuum/assemble_v2.jl
-522
View File
@@ -1,522 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
Assembly for 3D Continuum Elasticity - Ferrite-Style Two-Pointer Merge (V2)
This is a CLEAN implementation using the Ferrite two-pointer merge algorithm
for 4.1x faster assembly compared to COO triplets.
Key differences from continuum_3d.jl:
- Uses pre-built CSC structure (K_csc) instead of COO triplets
- Sorts element DOFs for linear-time merge
- Zero allocations in assembly loop
- 4.1x faster assembly (2.36ms vs 9.71ms for 2500 Tet4 elements)
- 16.6x less memory (506 KB vs 8.4 MB)
References:
- Design: experiments/FERRITE_DATA_STRUCTURES_DESIGN.md
- Proof of concept: experiments/ferrite_style_assembly.jl
- Explanation: experiments/ferrite_sorteddofs_explained.jl
"""
using Tensors
using SparseArrays
using LinearAlgebra
# Import basevec for creating unit vectors
using Tensors: basevec
# ============================================================================
# Data Structures
# ============================================================================
"""
AssemblyCacheFerrite{N,T,Mat}
Pre-allocated buffers for zero-allocation Ferrite-style assembly.
Replaces COO triplets (I_rows, J_cols, K_values) with:
- Pre-built CSC structure (K_csc) - reused every assembly
- Sorted DOF buffers (sorteddofs, permutation) - for two-pointer merge
# Type Parameters
- `N`: Number of nodes per element (from topology)
- `T`: Topology type (Hexahedron{8}, Tet10, etc.)
- `Mat`: Material type
# Fields (Ferrite-specific)
- `K_csc`: Pre-built sparse matrix structure (zeros, reused)
- `sorteddofs`: Sorted element DOF buffer [3N]
- `permutation`: Sortperm buffer for DOF remapping [3N]
# Fields (Standard)
- `f`: Global force vector
- `K_blocks`: Element blocked stiffness matrix [N×N of Tensor{2,3}]
- `K_e`: Element stiffness in Float64 matrix form [3N×3N]
- `X_buffer`: Element coordinate buffer [N of Vec{3}]
- `gdofs`: Global DOF indices buffer (unsorted) [3N]
- `elements`: Element IDs to assemble (Vector{UInt32})
- `C`: Elasticity tensor (pre-computed)
- `topology`: Topology instance
- `basis`: Basis function instance
- `ips`: Integration points (pre-computed)
# Performance
- Memory: ~506 KB for typical mesh (vs 8.4 MB COO triplets)
- Time: 2.36ms for 2500 Tet4 (vs 9.71ms COO)
- Allocations: Zero after warmup
"""
struct AssemblyCacheFerrite{N,T<:AbstractTopology{N},Mat<:AbstractMaterial}
# Ferrite-style CSC storage (NEW!)
K_csc::SparseMatrixCSC{Float64,Int}
sorteddofs::Vector{Int}
permutation::Vector{Int}
# Global force vector
f::Vector{Float64}
# Element assembly buffers (same as V1)
K_blocks::Matrix{Tensor{2,3,Float64,9}}
K_e::Matrix{Float64}
X_buffer::Vector{Vec{3,Float64}}
gdofs::Vector{Int} # Unsorted DOFs
# Element set to assemble
elements::Vector{UInt32}
# Pre-computed material/topology data
C::Tensor{4,3,Float64,81}
topology::T
basis::Lagrange{1} # Basis order only (topology passed separately)
ips::Any # Integration points tuple
end
# ============================================================================
# Sparsity Pattern Construction
# ============================================================================
"""
Build sparsity pattern from mesh connectivity.
Returns (I, J) vectors for sparse matrix construction.
"""
function build_sparsity_pattern(mesh::M) where {M<:AbstractMesh}
# Get N from mesh type parameters
N = typeof(mesh).parameters[1]
ndofs_total = 3 * length(mesh.nodes)
# Pre-count entries (3N × 3N per element)
ndofs_per_elem = 3 * N
capacity = length(mesh.connectivity) * ndofs_per_elem * ndofs_per_elem
I = Vector{Int}()
J = Vector{Int}()
sizehint!(I, capacity)
sizehint!(J, capacity)
# Loop over elements
for conn in mesh.connectivity
# Global DOFs for this element
elem_dofs = Int[]
sizehint!(elem_dofs, ndofs_per_elem)
for node_id in conn
for α in 1:3
push!(elem_dofs, 3 * (node_id - 1) + α)
end
end
# All pairs (i,j) in elem_dofs
for i in elem_dofs
for j in elem_dofs
push!(I, i)
push!(J, j)
end
end
end
return I, J
end
# ============================================================================
# Cache Construction
# ============================================================================
"""
AssemblyCacheFerrite(physics::Physics{ContinuumFormulation{FullThreeD},
Displacement{3}, M, Mat})
Construct Ferrite-style assembly cache with pre-built CSC structure.
This is where ALL allocations happen. After construction, assembly is zero-allocation.
# Key Steps
1. Build sparsity pattern from mesh connectivity
2. Create K_csc with sparse(I, J, ones(...))
3. Zero K_csc.nzval for reuse
4. Allocate sorteddofs and permutation buffers
5. Pre-compute material and topology data
# Performance
- One-time cost: ~10-20ms for typical mesh
- Pays off after ~1-2 assemblies vs COO method
"""
function AssemblyCacheFerrite(
physics::Physics{ContinuumFormulation{FullThreeD},
Displacement{3},
M,
Mat}) where {M<:AbstractMesh,Mat<:AbstractMaterial}
mesh = physics.mesh
material = physics.material
element_set = physics.element_set
# Get topology type and N from mesh type parameters
N_param = typeof(mesh).parameters[1] # N (8 for Hex8)
T = typeof(mesh).parameters[2] # Hexahedron{8}
# Global system dimensions
nnodes = length(mesh.nodes)
ndofs = 3 * nnodes
# Element set
elem_set = get_element_set(mesh, element_set)
elements = collect(elem_set)
# Build sparsity pattern ONCE
I, J = build_sparsity_pattern(mesh)
K_csc = sparse(I, J, ones(length(I)), ndofs, ndofs)
fill!(K_csc.nzval, 0.0) # Zero values for reuse
# Ferrite buffers for sorting DOFs
max_ndofs = 3 * N_param
sorteddofs = Vector{Int}(undef, max_ndofs)
permutation = Vector{Int}(undef, max_ndofs)
# Global force vector
f = zeros(ndofs)
# Element assembly buffers (same as V1)
max_nnodes = N_param
K_blocks = Matrix{Tensor{2,3,Float64,9}}(undef, max_nnodes, max_nnodes)
K_e = zeros(max_ndofs, max_ndofs)
X_buffer = Vector{Vec{3,Float64}}(undef, max_nnodes)
gdofs = Vector{Int}(undef, max_ndofs)
# Pre-compute material and topology data
C = elasticity_tensor(material)
topology = T()
basis = Lagrange{1}() # Basis order only (topology passed separately)
integration_scheme = default_integration(T)
ips = integration_points(integration_scheme, topology)
return AssemblyCacheFerrite{N_param,T,Mat}(
K_csc, # Pre-built structure
sorteddofs, # Sorted DOF buffer
permutation, # Sortperm buffer
f,
K_blocks, K_e, X_buffer, gdofs,
elements,
C, topology, basis, ips
)
end
# ============================================================================
# Ferrite Two-Pointer Merge Assembly
# ============================================================================
"""
assemble_elements_ferrite!(cache::AssemblyCacheFerrite, mesh::M)
Assemble elements using Ferrite two-pointer merge (ZERO allocations).
This is the core Ferrite algorithm:
1. Zero K_csc.nzval once at start
2. For each element:
a. Compute element stiffness K_e
b. Get global DOFs (unsorted)
c. Sort DOFs → sorteddofs, permutation
d. Two-pointer merge: Linear scan through sorted lists
e. Accumulate to K_csc using permutation for correct K_e indices
# Algorithm Detail
For each column i_global in element DOFs:
- Get CSC column range: K_csc.colptr[i_global]:(colptr[i_global+1]-1)
- K_csc.rowval[range] is SORTED
- sorteddofs is SORTED
- Two pointers: Ri (CSC), ri (element)
- Advance Ri until K_csc.rowval[Ri] == sorteddofs[ri]
- Accumulate: K_csc.nzval[Ri] += K_e[permutation[ri], permutation[i_local]]
# Performance
- Time: 2.36ms for 2500 Tet4 elements
- Allocations: Zero after warmup
- 4.1x faster than COO method
- Linear-time merge vs O(log n) binary search
"""
function assemble_elements_ferrite!(
cache::AssemblyCacheFerrite{N,T,Mat},
mesh::M) where {M<:AbstractMesh,N,T,Mat}
# Zero K_csc once at start (reuse structure!)
fill!(cache.K_csc.nzval, 0.0)
ndofs_elem = 3 * N
# Loop over elements (ZERO allocations target!)
@inbounds for i in eachindex(cache.elements)
elem_id = cache.elements[i]
@inbounds conn = mesh.connectivity[elem_id]
# 1. Fill coordinate buffer (in-place)
@inbounds for j in 1:N
cache.X_buffer[j] = mesh.nodes[conn[j]]
end
# 2. Compute element stiffness
fill!(cache.K_blocks, zero(Tensor{2,3}))
compute_element_stiffness!(cache.K_blocks, cache.X_buffer,
cache.C, cache.topology, cache.basis, cache.ips)
blocked_tensor_to_matrix!(cache.K_e, cache.K_blocks)
# 3. Global DOFs (UNSORTED, follows connectivity)
@inbounds for (local_idx, node_id) in enumerate(conn)
for α in 1:3
cache.gdofs[3*(local_idx-1)+α] = 3 * (node_id - 1) + α
end
end
# 4. Sort DOFs (sortperm! is in-place, zero allocation)
sortperm!(cache.permutation, cache.gdofs)
@inbounds for k in 1:ndofs_elem
cache.sorteddofs[k] = cache.gdofs[cache.permutation[k]]
end
# 5. TWO-POINTER MERGE (Ferrite method!)
for i_local in 1:ndofs_elem
i_global = cache.sorteddofs[i_local]
# Column range in K_csc for column i_global
col_start = cache.K_csc.colptr[i_global]
col_end = cache.K_csc.colptr[i_global+1] - 1
# Two pointers: Ri (CSC), ri (element)
Ri = col_start
for ri in 1:ndofs_elem
row_i_sorted = cache.sorteddofs[ri]
# Advance Ri until K_csc.rowval[Ri] >= row_i_sorted
while Ri <= col_end && cache.K_csc.rowval[Ri] < row_i_sorted
Ri += 1
end
# If found, accumulate
if Ri <= col_end && cache.K_csc.rowval[Ri] == row_i_sorted
# Use permutation to get correct K_e indices!
orig_row = cache.permutation[ri]
orig_col = cache.permutation[i_local]
cache.K_csc.nzval[Ri] += cache.K_e[orig_row, orig_col]
end
end
end
end
nothing
end
# ============================================================================
# Assembly Functions (User-Facing API)
# ============================================================================
"""
assemble_v2!(physics::Physics{ContinuumFormulation{FullThreeD},
Displacement{3}, M, Mat}) -> (K, f)
Assemble using Ferrite two-pointer merge method (V2).
This is the user-facing function that:
1. Creates AssemblyCacheFerrite (allocates all buffers, builds CSC structure)
2. Calls assemble_elements_ferrite! (zero allocations)
3. Applies boundary conditions
4. Returns (K, f)
# Performance
- 4.1x faster assembly than V1 (COO method)
- 16.6x less memory
- Zero allocations in hot loop
# Example
```julia
physics = Physics(...)
K, f = assemble_v2!(physics) # Ferrite method
u = K \\ f
```
"""
function assemble_v2!(
physics::Physics{ContinuumFormulation{FullThreeD},
Displacement{3},
M,
Mat}) where {M<:AbstractMesh,Mat<:AbstractMaterial}
# Create cache (ALL allocations here!)
cache = AssemblyCacheFerrite(physics)
# Assemble (ZERO allocations!)
return _assemble_ferrite!(physics, cache)
end
"""
_assemble_ferrite!(physics, cache::AssemblyCacheFerrite) -> (K, f)
Zero-allocation assembly using Ferrite cache.
This is the internal function that performs the actual assembly.
Use `assemble_v2!` for the public API.
"""
function _assemble_ferrite!(
physics::Physics{ContinuumFormulation{FullThreeD},
Displacement{3},
M,
Mat},
cache::AssemblyCacheFerrite{N,T,Mat2}) where {M<:AbstractMesh,Mat<:AbstractMaterial,N,T,Mat2}
mesh = physics.mesh
bc_dirichlet = physics.bc_dirichlet
bc_neumann = physics.bc_neumann
# Get dimensions
nnodes = length(mesh.nodes)
ndofs = 3 * nnodes
# Clear force vector (in-place, zero allocation)
fill!(cache.f, 0.0)
# Ferrite assembly (ZERO allocations!)
assemble_elements_ferrite!(cache, mesh)
# Apply Neumann BCs (add forces to f)
for (surf_id, force) in zip(bc_neumann.surface_ids, bc_neumann.values)
# For now, interpret surface_ids as node_ids (simplified)
node = surf_id
if node <= nnodes
for α in 1:3
cache.f[3*(node-1)+α] += force[α]
end
end
end
# Copy K_csc (structure already correct, just copy!)
K = copy(cache.K_csc)
# Apply Dirichlet BCs (modify K and f)
for i in 1:length(bc_dirichlet.node_ids)
node = bc_dirichlet.node_ids[i]
components = bc_dirichlet.components[i]
value = bc_dirichlet.values[i]
for comp in components
dof = 3 * (node - 1) + comp
if dof <= ndofs # Safety check
K[dof, :] .= 0.0
K[:, dof] .= 0.0
K[dof, dof] = 1.0
cache.f[dof] = value
end
end
end
return (K, cache.f)
end
# ============================================================================
# Element Stiffness Computation (Reuse from V1)
# ============================================================================
"""
compute_stiffness_block(grad_k, grad_l, C) -> Tensor{2,3}
Compute single 3×3 stiffness block (reused from continuum_3d.jl).
See continuum_3d.jl for detailed documentation.
"""
@inline function compute_stiffness_block(
grad_k::Vec{3,Float64},
grad_l::Vec{3,Float64},
C::Tensor{4,3,Float64,81}
)::Tensor{2,3,Float64,9}
K_kl = zero(Tensor{2,3})
for α in 1:3, β in 1:3
e_α = basevec(Vec{3}, α)
e_β = basevec(Vec{3}, β)
B_k_α = 0.5 * (grad_k ⊗ e_α + e_α ⊗ grad_k)
B_l_β = 0.5 * (grad_l ⊗ e_β + e_β ⊗ grad_l)
k_αβ = dcontract(B_k_α, dcontract(C, B_l_β))
K_kl += k_αβ * (e_α ⊗ e_β)
end
return K_kl
end
"""
blocked_tensor_to_matrix!(K_e, K_blocks)
Convert blocked tensor to Float64 matrix (reused from continuum_3d.jl).
"""
function blocked_tensor_to_matrix!(
K_e::AbstractMatrix{Float64},
K_blocks::AbstractMatrix{Tensor{2,3,Float64,9}})
nnodes = size(K_blocks, 1)
for i in 1:nnodes, j in 1:nnodes
for α in 1:3, β in 1:3
K_e[3*(i-1)+α, 3*(j-1)+β] = K_blocks[i, j][α, β]
end
end
nothing
end
"""
compute_element_stiffness!(K_blocks, X, C, topology, basis, ips)
Compute element stiffness (reused from continuum_3d.jl).
See continuum_3d.jl for detailed documentation.
"""
function compute_element_stiffness!(
K_blocks::AbstractMatrix{Tensor{2,3,Float64,9}},
X::Vector{Vec{3,Float64}},
C::Tensor{4,3,Float64,81},
topology::T,
basis::B,
ips) where {T<:AbstractTopology{N},B<:AbstractBasis} where N
for k in 1:N, l in 1:N
for ip in ips
ξ = ip.ξ
w = ip.weight
dN_dξ = get_basis_derivatives(topology, basis, ξ)
J = X[1] ⊗ dN_dξ[1]
for i in 2:N
J += X[i] ⊗ dN_dξ[i]
end
detJ = det(J)
J_inv = inv(J)
J_inv_T = transpose(J_inv)
grad_k = J_inv_T ⋅ dN_dξ[k]
grad_l = J_inv_T ⋅ dN_dξ[l]
K_kl = compute_stiffness_block(grad_k, grad_l, C)
K_blocks[k, l] += K_kl * detJ * w
end
end
nothing
end
src/element_assembly_structures.jl
-341
View File
@@ -1,341 +0,0 @@
# Traditional Element Assembly
#
# This module provides the standard element-by-element assembly approach
# for comparison with nodal assembly. Builds global tangent stiffness matrix
# and residual force vector using sparse matrix formats.
using Tensors
using SparseArrays
using LinearAlgebra
"""
ElementAssemblyData{T}
Storage for element assembly using traditional (element-by-element) approach.
# Fields
- `K_global::SparseMatrixCSC{T}`: Global tangent stiffness matrix
- `r_global::Vector{T}`: Global residual force vector (r = f_int - f_ext)
- `f_int_global::Vector{T}`: Global internal force vector
- `f_ext_global::Vector{T}`: Global external force vector
- `ndof::Int`: Total number of degrees of freedom
# Notes
- Assembly uses COO (coordinate) format, then converts to CSC
- Multiple elements can write to same global DOF (summed automatically)
"""
mutable struct ElementAssemblyData{T}
K_global::SparseMatrixCSC{T,Int}
r_global::Vector{T}
f_int_global::Vector{T}
f_ext_global::Vector{T}
ndof::Int
end
"""
ElementAssemblyData(ndof::Int, ::Type{T}=Float64)
Allocate storage for traditional element assembly.
# Arguments
- `ndof`: Total degrees of freedom (nnodes × 3 for 3D)
- `T`: Floating point type (default Float64)
# Example
```julia
nnodes = 100
assembly = ElementAssemblyData(3 * nnodes, Float64)
```
"""
function ElementAssemblyData(ndof::Int, ::Type{T}=Float64) where T
# Pre-allocate empty sparse matrix (will fill during assembly)
K_global = spzeros(T, ndof, ndof)
r_global = zeros(T, ndof)
f_int_global = zeros(T, ndof)
f_ext_global = zeros(T, ndof)
return ElementAssemblyData{T}(K_global, r_global, f_int_global, f_ext_global, ndof)
end
"""
reset!(assembly::ElementAssemblyData)
Reset assembly data to zero (for incremental/iterative solvers).
"""
function reset!(assembly::ElementAssemblyData{T}) where T
assembly.K_global = spzeros(T, assembly.ndof, assembly.ndof)
fill!(assembly.r_global, 0.0)
fill!(assembly.f_int_global, 0.0)
fill!(assembly.f_ext_global, 0.0)
end
"""
ElementContribution{T}
Local element contribution before scattering to global.
# Fields
- `element_id::Int`: Element ID
- `gdofs::Vector{Int}`: Global DOF indices (e.g., [1,2,3,4,5,6,...] for nodes)
- `K_local::Matrix{T}`: Local stiffness matrix (ndofs_local × ndofs_local)
- `f_int_local::Vector{T}`: Local internal force vector
- `f_ext_local::Vector{T}`: Local external force vector
# Notes
- For Tet4: ndofs_local = 12 (4 nodes × 3 DOF)
- For Tet10: ndofs_local = 30 (10 nodes × 3 DOF)
"""
struct ElementContribution{T}
element_id::Int
gdofs::Vector{Int}
K_local::Matrix{T}
f_int_local::Vector{T}
f_ext_local::Vector{T}
end
"""
ElementContribution(element_id::Int, gdofs::Vector{Int}, ::Type{T}=Float64)
Allocate storage for element contribution.
# Arguments
- `element_id`: Element ID
- `gdofs`: Global DOF indices
- `T`: Floating point type
# Example
```julia
# Tet4 element connecting nodes [5, 7, 12, 15]
gdofs = [13,14,15, 19,20,21, 34,35,36, 43,44,45] # 3 DOF per node
contrib = ElementContribution(1, gdofs, Float64)
```
"""
function ElementContribution(element_id::Int, gdofs::Vector{Int}, ::Type{T}=Float64) where T
ndofs = length(gdofs)
K_local = zeros(T, ndofs, ndofs)
f_int_local = zeros(T, ndofs)
f_ext_local = zeros(T, ndofs)
return ElementContribution{T}(element_id, gdofs, K_local, f_int_local, f_ext_local)
end
"""
scatter_to_global!(assembly::ElementAssemblyData, contrib::ElementContribution)
Scatter element contribution to global matrices/vectors (traditional assembly).
This is the key operation in element assembly: add local element quantities
to global system. Uses COO format (accumulates into lists).
# Arguments
- `assembly`: Global assembly data
- `contrib`: Element contribution
# Notes
- Multiple elements can contribute to same global DOF (summed)
- For GPU: Would require atomic operations (slow!)
- For CPU: Direct scatter-add works fine
"""
function scatter_to_global!(assembly::ElementAssemblyData{T},
contrib::ElementContribution{T}) where T
# Scatter forces (simple vector addition)
for (local_i, global_i) in enumerate(contrib.gdofs)
assembly.f_int_global[global_i] += contrib.f_int_local[local_i]
assembly.f_ext_global[global_i] += contrib.f_ext_local[local_i]
end
# Scatter stiffness (matrix addition)
# Build list of (I, J, V) triplets for sparse matrix
I_rows = Int[]
J_cols = Int[]
values = T[]
ndofs_local = length(contrib.gdofs)
for i in 1:ndofs_local, j in 1:ndofs_local
if abs(contrib.K_local[i, j]) > 1e-14 # Skip near-zeros
push!(I_rows, contrib.gdofs[i])
push!(J_cols, contrib.gdofs[j])
push!(values, contrib.K_local[i, j])
end
end
# Add to existing sparse matrix
K_elem = sparse(I_rows, J_cols, values, assembly.ndof, assembly.ndof)
assembly.K_global += K_elem
end
"""
compute_residual!(assembly::ElementAssemblyData)
Compute residual force vector: r = f_int - f_ext
Should be called after all elements have been assembled.
"""
function compute_residual!(assembly::ElementAssemblyData{T}) where T
assembly.r_global .= assembly.f_int_global .- assembly.f_ext_global
end
"""
assemble_elements!(assembly::ElementAssemblyData,
contributions::Vector{ElementContribution})
Assemble all element contributions to global system.
# Arguments
- `assembly`: Global assembly data (modified in-place)
- `contributions`: Vector of element contributions
# Example
```julia
assembly = ElementAssemblyData(ndof)
contributions = compute_all_element_contributions(elements, u, time)
assemble_elements!(assembly, contributions)
compute_residual!(assembly)
# Now solve: K_global * Δu = -r_global
```
"""
function assemble_elements!(assembly::ElementAssemblyData{T},
contributions::Vector{ElementContribution{T}}) where T
reset!(assembly)
# Loop over elements and scatter (element assembly)
for contrib in contributions
scatter_to_global!(assembly, contrib)
end
# Compute residual
compute_residual!(assembly)
end
"""
apply_dirichlet_bc!(assembly::ElementAssemblyData,
fixed_dofs::Vector{Int},
prescribed_values::Vector{T}=zeros(length(fixed_dofs)))
Apply Dirichlet (essential) boundary conditions by penalty method.
# Arguments
- `assembly`: Global assembly data (modified in-place)
- `fixed_dofs`: DOF indices to fix
- `prescribed_values`: Prescribed displacement values (default: zeros)
# Method
Uses penalty method: adds large stiffness to diagonal and corresponding RHS.
For DOF i with prescribed value u_prescribed:
- K[i,i] += penalty (e.g., 1e10 * max_K)
- r[i] = penalty * (u_current - u_prescribed)
# Example
```julia
# Fix nodes 1 and 2 in all directions (zero displacement)
fixed_dofs = [1,2,3, 4,5,6] # Nodes 1,2 × 3 DOF
apply_dirichlet_bc!(assembly, fixed_dofs)
```
"""
function apply_dirichlet_bc!(assembly::ElementAssemblyData{T},
fixed_dofs::Vector{Int},
prescribed_values::Vector{T}=zeros(T, length(fixed_dofs))) where T
# Penalty parameter (large relative to stiffness)
max_K = maximum(abs, assembly.K_global)
penalty = 1e10 * max_K
for (idx, dof) in enumerate(fixed_dofs)
# Add penalty stiffness to diagonal
assembly.K_global[dof, dof] += penalty
# Modify residual (assuming current displacement is zero for now)
# In full Newton: r[i] += penalty * (u_current[i] - u_prescribed[i])
assembly.r_global[dof] = penalty * prescribed_values[idx]
end
end
"""
get_dof_indices(connectivity::NTuple{N,Int}, dim::Int=3) -> Vector{Int}
Get global DOF indices for an element given node connectivity.
# Arguments
- `connectivity`: Element node IDs (e.g., (5, 7, 12, 15) for Tet4)
- `dim`: Dimension (3 for 3D elasticity)
# Returns
- `gdofs::Vector{Int}`: Global DOF indices
# Example
```julia
# Element with nodes [5, 7, 12, 15]
gdofs = get_dof_indices((5, 7, 12, 15), 3)
# Returns: [13,14,15, 19,20,21, 34,35,36, 43,44,45]
```
"""
function get_dof_indices(connectivity::NTuple{N,Int}, dim::Int=3) where N
nnodes = length(connectivity)
gdofs = zeros(Int, dim * nnodes)
for (local_i, global_node) in enumerate(connectivity)
for d in 1:dim
gdofs[dim*(local_i-1)+d] = dim * (global_node - 1) + d
end
end
return gdofs
end
"""
matrix_vector_product(assembly::ElementAssemblyData, v::Vector{T}) -> Vector{T}
Compute matrix-vector product: w = K * v using assembled sparse matrix.
# Arguments
- `assembly`: Assembly data (contains K_global)
- `v`: Input vector (ndof)
# Returns
- `w`: Output vector w = K * v
# Example
```julia
# GMRES matrix-free operator
function matvec(v)
return matrix_vector_product(assembly, v)
end
Δu = gmres(matvec, -r, tol=1e-6)
```
"""
function matrix_vector_product(assembly::ElementAssemblyData{T}, v::Vector{T}) where T
return assembly.K_global * v
end
"""
print_assembly_stats(assembly::ElementAssemblyData)
Print statistics about assembled system (for debugging).
"""
function print_assembly_stats(assembly::ElementAssemblyData)
nnz_K = nnz(assembly.K_global)
ndof = assembly.ndof
fill_ratio = nnz_K / (ndof * ndof)
println("="^60)
println("Traditional Element Assembly Statistics")
println("="^60)
println(" Total DOF: ", ndof)
println(" K matrix size: ", size(assembly.K_global))
println(" K non-zeros: ", nnz_K)
println(" K fill ratio: ", round(fill_ratio, sigdigits=2))
println(" K memory (MB): ", round(nnz_K * 16 / 1024^2, digits=2))
println(" ||f_int||: ", round(norm(assembly.f_int_global), sigdigits=2))
println(" ||f_ext||: ", round(norm(assembly.f_ext_global), sigdigits=2))
println(" ||residual||: ", round(norm(assembly.r_global), sigdigits=2))
println(" K symmetric: ", issymmetric(assembly.K_global))
println("="^60)
end
# Export main types and functions
export ElementAssemblyData, ElementContribution
export reset!, scatter_to_global!, compute_residual!
export assemble_elements!, apply_dirichlet_bc!
export get_dof_indices, matrix_vector_product
export print_assembly_stats
src/elements/integrate.jl
-57
View File
@@ -1,57 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
# Default number of integration points for each element. First rule is the
# default integration rule returned by `get_integration_points(element)`.
# Sometimes we want to increase integration order, e.g. when integrating mass
# matrix or boundary conditions. For that reason, additional rules are provied
# in list, so e.g. `get_integration_points(element, 1)` returns the second rule,
# `get_integration_points(element, 2)` third rule and so on. Rules should be
# ordered so that picking next one integrates more accurately.
integration_rule_mapping = (
:Seg2 => (:GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5),
:Seg3 => (:GLSEG3, :GLSEG4, :GLSEG5),
:NSeg => (:GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5),
:Quad4 => (:GLQUAD4, :GLQUAD9, :GLQUAD16, :GLQUAD25),
:Quad8 => (:GLQUAD9, :GLQUAD16, :GLQUAD25),
:Quad9 => (:GLQUAD9, :GLQUAD16, :GLQUAD25),
:NSurf => (:GLQUAD9, :GLQUAD16, :GLQUAD25),
:Hex8 => (:GLHEX8, :GLHEX27, :GLHEX64, :GLHEX125),
:Hex20 => (:GLHEX27, :GLHEX64, :GLHEX125),
:Hex27 => (:GLHEX27, :GLHEX64, :GLHEX125),
:NSolid => (:GLHEX27, :GLHEX64, :GLHEX125),
:Tri3 => (:GLTRI1, :GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12),
:Tri6 => (:GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12),
:Tri7 => (:GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12),
:Tet4 => (:GLTET1, :GLTET4, :GLTET5, :GLTET15),
:Tet10 => (:GLTET4, :GLTET5, :GLTET15),
:Pyr5 => (:GLPYR5,),
:Wedge6 => (:GLWED6, :GLWED21),
:Wedge15 => (:GLWED21,))
for (E, R) in integration_rule_mapping
for i in 1:length(R)
P = Val{R[i]}
order = Val{i - 1}
local code # Explicitly declare as local to avoid warning
if isequal(i, 1)
code = quote
function get_integration_points(element::$E)
return get_quadrature_points($P)
end
end
else
code = quote
function get_integration_points(element::$E, ::Type{$order})
return get_quadrature_points($P)
end
end
end
eval(code)
end
end
# All good codes needs a special case. Here we have it: Poi1
function get_integration_points(::Poi1)
[(1.0, (0.0,))]
end
src/formulations/api.jl
-307
View File
@@ -1,307 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
Formulation API definitions.
This file defines formulation abstractions - the mathematical discretization strategies
for different types of FEM problems.
Must be included after fields/api.jl (formulations work with fields).
"""
# ============================================================================
# FORMULATION INTERFACE
# ============================================================================
"""
AbstractFormulation
Abstract type for discretization formulations.
Formulation defines HOW we discretize the governing equations. Different formulations
exist for different physics domains:
- **Continuum formulations** (this file) - Standard FEM for solid/fluid mechanics
- **Beam formulations** (src/beams/api.jl) - 1D structural elements
- **Shell formulations** (src/shells/api.jl) - 2D structural elements
- **Truss formulations** (src/trusses/api.jl) - 1D axial elements
# Type Hierarchy
- `ContinuumFormulation{Theory}` - Standard continuum FEM (here)
- `BeamFormulation{Theory}` - Beam elements (src/beams/api.jl)
- `ShellFormulation{Theory}` - Shell elements (src/shells/api.jl)
- `TrussFormulation{Theory}` - Truss elements (src/trusses/api.jl)
# Design Philosophy
**Formulation + Field = Dispatch pattern**
The combination of formulation and field type determines:
- Assembly method dispatch
- Element stiffness computation
- Stress/strain tensor dimensions
- DOF coupling patterns
# Examples
```julia
# 3D solid mechanics
physics = Physics(
formulation = ContinuumFormulation{FullThreeD}(),
field = Displacement{3}(),
mesh = mesh,
material = steel
)
# 2D plane stress
physics_2d = Physics(
formulation = ContinuumFormulation{PlaneStress}(),
field = Displacement{2}(),
mesh = mesh_2d,
material = aluminum
)
# Beam structure
physics_beam = Physics(
formulation = BeamFormulation{Timoshenko}(),
field = DisplacementRotation{3}(),
mesh = beam_mesh,
material = steel
)
```
# Assembly Dispatch
Specialized assembly methods dispatch on formulation × field:
```julia
# 3D continuum mechanics
function assemble!(physics::Physics{ContinuumFormulation{FullThreeD}, Displacement{3}, M, Mat})
# Standard 3D displacement-based assembly
# Implementation in src/assembly/continuum_3d.jl
end
# 2D plane stress
function assemble!(physics::Physics{ContinuumFormulation{PlaneStress}, Displacement{2}, M, Mat})
# 2D assembly with plane stress assumptions
# Implementation in src/assembly/continuum_2d.jl
end
# Beam elements
function assemble!(physics::Physics{BeamFormulation{Timoshenko}, DisplacementRotation{3}, M, Mat})
# Beam-specific assembly (6 DOFs per node)
# Implementation in src/assembly/beams.jl
end
```
# See Also
- Field types: src/fields/api.jl (Displacement, Temperature, DisplacementRotation)
- Physics coupling: src/physics/api.jl (AbstractPhysics)
- Domain-specific formulations: src/beams/api.jl, src/shells/api.jl, src/trusses/api.jl
- Assembly implementations: src/assembly/continuum_3d.jl, src/assembly/beams.jl, etc.
"""
abstract type AbstractFormulation end
# ============================================================================
# CONTINUUM FORMULATION (Standard FEM)
# ============================================================================
"""
AbstractContinuumTheory
Theory variants for continuum formulation.
Controls dimensionality reduction and stress/strain assumptions for continuum
mechanics problems.
# Concrete Theories
- `FullThreeD` - Full 3D analysis (no simplifications)
- `PlaneStress` - 2D plane stress (σ_zz = 0, thin plates)
- `PlaneStrain` - 2D plane strain (ε_zz = 0, thick plates)
- `Axisymmetric` - Axisymmetric analysis (rotation around z-axis)
# Theory Selection Guidelines
**FullThreeD (σ_xx, σ_yy, σ_zz, σ_xy, σ_yz, σ_xz):**
- General 3D solid mechanics
- No simplifying assumptions
- Most accurate but most expensive
**PlaneStress (σ_xx, σ_yy, σ_xy, σ_zz = 0):**
- Thin plates and membranes (thickness << length/width)
- Out-of-plane stress σ_zz = 0
- Examples: Sheet metal, aircraft skin, thin-walled structures
**PlaneStrain (ε_xx, ε_yy, ε_xy, ε_zz = 0):**
- Thick sections with no variation in z-direction
- Out-of-plane strain ε_zz = 0
- Examples: Dams, tunnels, retaining walls, long cylinders
**Axisymmetric (σ_rr, σ_θθ, σ_zz, σ_rz):**
- Geometry and loading symmetric about z-axis
- No circumferential variations
- Examples: Pressure vessels, pipes, rotating disks
# Usage
```julia
# Full 3D solid mechanics
formulation = ContinuumFormulation{FullThreeD}()
# 2D plane stress (thin plate)
formulation = ContinuumFormulation{PlaneStress}()
# 2D plane strain (thick section)
formulation = ContinuumFormulation{PlaneStrain}()
# Axisymmetric (cylinder, sphere)
formulation = ContinuumFormulation{Axisymmetric}()
```
# Mathematical Details
**Plane Stress (thin plate):**
- Stress state: σ_zz = σ_xz = σ_yz = 0
- Strain: ε_zz ≠ 0 (computed from σ_zz = 0 condition)
- Constitutive: 3×3 reduced stiffness matrix
**Plane Strain (thick section):**
- Strain state: ε_zz = γ_xz = γ_yz = 0
- Stress: σ_zz ≠ 0 (computed from ε_zz = 0 condition)
- Constitutive: 3×3 reduced stiffness matrix (different from plane stress!)
**Axisymmetric:**
- Cylindrical coordinates (r, θ, z)
- No ∂/∂θ terms (axial symmetry)
- 4 stress components: σ_rr, σ_θθ, σ_zz, σ_rz
- Hoop stress σ_θθ from radial displacement
# See Also
- [`ContinuumFormulation`](@ref) - Formulation struct using these theories
"""
abstract type AbstractContinuumTheory end
"""
FullThreeD <: AbstractContinuumTheory
Full 3D analysis with no simplifications.
All six stress components: σ_xx, σ_yy, σ_zz, σ_xy, σ_yz, σ_xz
"""
struct FullThreeD <: AbstractContinuumTheory end
"""
PlaneStress <: AbstractContinuumTheory
2D plane stress assumption (σ_zz = 0).
Applicable to thin plates and membranes where thickness << in-plane dimensions.
"""
struct PlaneStress <: AbstractContinuumTheory end
"""
PlaneStrain <: AbstractContinuumTheory
2D plane strain assumption (ε_zz = 0).
Applicable to thick sections with no variation in z-direction.
"""
struct PlaneStrain <: AbstractContinuumTheory end
"""
Axisymmetric <: AbstractContinuumTheory
Axisymmetric analysis (rotation around z-axis).
Geometry and loading symmetric about z-axis with no circumferential variations.
"""
struct Axisymmetric <: AbstractContinuumTheory end
"""
ContinuumFormulation{Theory} <: AbstractFormulation
Standard continuum mechanics formulation with theory variant.
This is the fundamental FEM formulation for solid mechanics, heat transfer,
and other continuum physics problems.
# Type Parameter
- `Theory <: AbstractContinuumTheory` - Dimensionality/simplification theory
# Examples
```julia
# 3D elasticity
physics = Physics(
formulation = ContinuumFormulation{FullThreeD}(),
field = Displacement{3}(),
mesh = mesh,
material = steel
)
# 2D plane stress (thin plate)
physics_2d = Physics(
formulation = ContinuumFormulation{PlaneStress}(),
field = Displacement{2}(),
mesh = mesh_2d,
material = aluminum
)
# 2D plane strain (thick section)
physics_2d = Physics(
formulation = ContinuumFormulation{PlaneStrain}(),
field = Displacement{2}(),
mesh = mesh_2d,
material = concrete
)
# Axisymmetric (cylinder)
physics_axisym = Physics(
formulation = ContinuumFormulation{Axisymmetric}(),
field = Displacement{2}(), # (r, z) displacements
mesh = mesh_2d,
material = steel
)
```
# Assembly Dispatch
Assembly methods specialize on theory × field combinations:
```julia
# 3D solid mechanics
function assemble!(physics::Physics{ContinuumFormulation{FullThreeD}, Displacement{3}, M, Mat})
# Standard 3D displacement-based assembly
# Full 6×6 strain-displacement matrix (Bε)
# 6×6 constitutive matrix (Dε)
end
# 2D plane stress
function assemble!(physics::Physics{ContinuumFormulation{PlaneStress}, Displacement{2}, M, Mat})
# 2D assembly with plane stress assumptions
# 3×3 reduced strain-displacement matrix
# 3×3 plane stress constitutive matrix
end
# Heat transfer (same formulation, different field!)
function assemble!(physics::Physics{ContinuumFormulation{FullThreeD}, Temperature, M, Mat})
# Thermal assembly (scalar field)
# Thermal conductivity matrix
end
```
# Implementation Location
Concrete assembly implementations are in:
- `src/assembly/continuum_3d.jl` - 3D continuum mechanics
- `src/assembly/continuum_2d.jl` - 2D plane stress/strain
- `src/assembly/axisymmetric.jl` - Axisymmetric problems
# See Also
- [`AbstractContinuumTheory`](@ref) - Theory variants
- Field types: src/fields/api.jl (Displacement, Temperature)
- Physics coupling: src/physics/api.jl (AbstractPhysics)
- Assembly: src/assembly/continuum_*.jl
"""
struct ContinuumFormulation{Theory<:AbstractContinuumTheory} <: AbstractFormulation end
src/gpu_elasticity.jl
-476
View File
@@ -1,476 +0,0 @@
"""
Main solver: Elasticity on GPU
Solves linear elasticity using:
- Two-phase nodal assembly (no atomics)
- Matrix-free conjugate gradient
- GPU-resident throughout
"""
function solve_elasticity_gpu(physics::ElasticityPhysics; tol=1e-6, max_iter=1000)
module GPUElasticity
export solve_elasticity_gpu, ElasticityPhysics, ElasticMaterial
using CUDA
using Tensors
using LinearAlgebra
using Printf
# Re-export mesh reader
include("gmsh_reader.jl")
using .GmshReader
export read_gmsh_mesh, GmshMesh, get_surface_nodes
"""
Elastic material properties
"""
struct ElasticMaterial
E::Float64 # Young's modulus [Pa]
ν::Float64 # Poisson's ratio [-]
end
"""
Elasticity physics definition
"""
struct ElasticityPhysics
mesh::GmshMesh
material::ElasticMaterial
fixed_nodes::Vector{Int} # Dirichlet BC (fixed displacement)
pressure_nodes::Vector{Int} # Neumann BC (pressure load)
pressure_value::Float64 # Pressure magnitude [Pa]
end
"""
Node-to-elements connectivity (CSR format)
"""
struct NodeToElementsMap
ptr::CuArray{Int32,1}
data::CuArray{Int32,1}
end
"""
Build CSR map: which elements touch each node?
"""
function build_node_to_elems_gpu(elements::Matrix{Int}, n_nodes::Int)
# Count connections per node
counts = zeros(Int, n_nodes)
for elem_idx in 1:size(elements, 2)
for i in 1:4
node = elements[i, elem_idx]
counts[node] += 1
end
end
# Build CSR structure
ptr = cumsum([1; counts])
data = Vector{Int32}(undef, sum(counts))
# Fill data array
offset = copy(ptr[1:end-1])
for elem_idx in 1:size(elements, 2)
for i in 1:4
node = elements[i, elem_idx]
data[offset[node]] = elem_idx
offset[node] += 1
end
end
return NodeToElementsMap(CuArray(Int32.(ptr)), CuArray(data))
end
"""
PHASE 1 GPU KERNEL: Compute element stiffness contributions at integration points
For LINEAR ELASTICITY (no plasticity), we don't need state variables.
Just compute stresses from strains using Hooke's law.
"""
function compute_element_stresses_kernel!(
σ_gp::CuDeviceArray{SymmetricTensor{2,3,Float64,6},1},
u::CuDeviceArray{Float64,1},
nodes::CuDeviceArray{Float64,2},
elements::CuDeviceArray{Int32,2},
E, ν
)
gp_idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
if gp_idx <= length(σ_gp)
# Map GP to element
elem_idx = (gp_idx - 1) ÷ 4 + 1 # 4 GPs per Tet4
# Extract element nodes
n1 = elements[1, elem_idx]
n2 = elements[2, elem_idx]
n3 = elements[3, elem_idx]
n4 = elements[4, elem_idx]
# Node coordinates
X1 = Vec{3}((nodes[1, n1], nodes[2, n1], nodes[3, n1]))
X2 = Vec{3}((nodes[1, n2], nodes[2, n2], nodes[3, n2]))
X3 = Vec{3}((nodes[1, n3], nodes[2, n3], nodes[3, n3]))
X4 = Vec{3}((nodes[1, n4], nodes[2, n4], nodes[3, n4]))
# Displacements
u1 = Vec{3}((u[3*n1-2], u[3*n1-1], u[3*n1]))
u2 = Vec{3}((u[3*n2-2], u[3*n2-1], u[3*n2]))
u3 = Vec{3}((u[3*n3-2], u[3*n3-1], u[3*n3]))
u4 = Vec{3}((u[3*n4-2], u[3*n4-1], u[3*n4]))
# Shape derivatives (constant for Tet4)
dN1_dxi = Vec{3}((-1.0, -1.0, -1.0))
dN2_dxi = Vec{3}((1.0, 0.0, 0.0))
dN3_dxi = Vec{3}((0.0, 1.0, 0.0))
dN4_dxi = Vec{3}((0.0, 0.0, 1.0))
# Jacobian
J = dN1_dxi ⊗ X1 + dN2_dxi ⊗ X2 + dN3_dxi ⊗ X3 + dN4_dxi ⊗ X4
invJ = inv(J)
# Physical derivatives
dN1_dx = invJ ⋅ dN1_dxi
dN2_dx = invJ ⋅ dN2_dxi
dN3_dx = invJ ⋅ dN3_dxi
dN4_dx = invJ ⋅ dN4_dxi
# Strain (small strain assumption)
ε = symmetric(dN1_dx ⊗ u1 + dN2_dx ⊗ u2 + dN3_dx ⊗ u3 + dN4_dx ⊗ u4)
# Stress (Hooke's law)
λ = E * ν / ((1 + ν) * (1 - 2ν))
μ = E / (2(1 + ν))
I = one(ε)
σ = λ * tr(ε) * I + 2μ * ε
# Store result
σ_gp[gp_idx] = σ
end
return nothing
end
"""
PHASE 2 GPU KERNEL: Nodal assembly (matrix-free, no atomics!)
"""
function nodal_assembly_kernel!(
r::CuDeviceArray{Float64,1},
σ_gp::CuDeviceArray{SymmetricTensor{2,3,Float64,6},1},
nodes::CuDeviceArray{Float64,2},
elements::CuDeviceArray{Int32,2},
node_to_elems_ptr::CuDeviceArray{Int32,1},
node_to_elems_data::CuDeviceArray{Int32,1}
)
node_idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
if node_idx <= size(nodes, 2)
# Accumulate forces
f_node = zero(Vec{3,Float64})
# Gauss weight for Tet4
gauss_weight = 1.0 / 24.0
# Shape derivatives
dN_dxi = (
Vec{3}((-1.0, -1.0, -1.0)),
Vec{3}((1.0, 0.0, 0.0)),
Vec{3}((0.0, 1.0, 0.0)),
Vec{3}((0.0, 0.0, 1.0))
)
# Get element range for this node
elem_start = node_to_elems_ptr[node_idx]
elem_end = node_to_elems_ptr[node_idx+1] - 1
# Loop over touching elements
for elem_offset in elem_start:elem_end
elem_idx = node_to_elems_data[elem_offset]
# Extract element nodes
n1 = elements[1, elem_idx]
n2 = elements[2, elem_idx]
n3 = elements[3, elem_idx]
n4 = elements[4, elem_idx]
# Find local node index
local_node = 1
if node_idx == n2
local_node = 2
elseif node_idx == n3
local_node = 3
elseif node_idx == n4
local_node = 4
end
# Recompute geometry (matrix-free!)
X1 = Vec{3}((nodes[1, n1], nodes[2, n1], nodes[3, n1]))
X2 = Vec{3}((nodes[1, n2], nodes[2, n2], nodes[3, n2]))
X3 = Vec{3}((nodes[1, n3], nodes[2, n3], nodes[3, n3]))
X4 = Vec{3}((nodes[1, n4], nodes[2, n4], nodes[3, n4]))
J = dN_dxi[1] ⊗ X1 + dN_dxi[2] ⊗ X2 + dN_dxi[3] ⊗ X3 + dN_dxi[4] ⊗ X4
detJ = det(J)
invJ = inv(J)
# Physical derivative for this node
dN_dx = invJ ⋅ dN_dxi[local_node]
# Loop over Gauss points (4 per Tet4)
for local_gp in 1:4
gp_idx = (elem_idx - 1) * 4 + local_gp
σ = σ_gp[gp_idx]
# Accumulate force
f_node += (dN_dx ⋅ σ) * (gauss_weight * detJ)
end
end
# Write result (no atomics!)
r[3*node_idx-2] = f_node[1]
r[3*node_idx-1] = f_node[2]
r[3*node_idx] = f_node[3]
end
return nothing
end
"""
Compute residual on GPU (internal forces)
"""
function compute_residual_gpu!(
r::CuArray{Float64,1},
u::CuArray{Float64,1},
nodes::CuArray{Float64,2},
elements::CuArray{Int32,2},
node_to_elems::NodeToElementsMap,
E, ν
)
n_gp = size(elements, 2) * 4
n_nodes = size(nodes, 2)
# Phase 1: Compute stresses at GPs
σ_gp = CuArray{SymmetricTensor{2,3,Float64,6}}(undef, n_gp)
threads = 256
blocks = cld(n_gp, threads)
@cuda threads = threads blocks = blocks compute_element_stresses_kernel!(
σ_gp, u, nodes, elements, E, ν
)
# Phase 2: Nodal assembly
fill!(r, 0.0)
threads = 256
blocks = cld(n_nodes, threads)
@cuda threads = threads blocks = blocks nodal_assembly_kernel!(
r, σ_gp, nodes, elements,
node_to_elems.ptr, node_to_elems.data
)
return r
end
"""
Apply pressure load to top surface (Neumann BC)
"""
function apply_pressure_load!(
f::CuArray{Float64,1},
pressure_nodes::Vector{Int},
mesh::GmshMesh,
pressure::Float64
)
# Simple uniform distribution (should integrate properly over surface)
# For now, divide pressure equally among nodes
f_cpu = Array(f)
n_pressure_nodes = length(pressure_nodes)
# Estimate surface area (assuming uniform Z = height)
surface_area = (maximum(mesh.nodes[1, :]) - minimum(mesh.nodes[1, :])) *
(maximum(mesh.nodes[2, :]) - minimum(mesh.nodes[2, :]))
# Total force
total_force = pressure * surface_area
force_per_node = total_force / n_pressure_nodes
# Apply in Z direction (negative, pointing down)
for node in pressure_nodes
f_cpu[3*node] += -force_per_node # Z component
end
copyto!(f, f_cpu)
return f
end
"""
Apply Dirichlet boundary conditions (fixed nodes)
"""
function apply_dirichlet_bc!(
K_op::Function,
f::CuArray{Float64,1},
fixed_nodes::Vector{Int}
)
# Zero out DOFs
f_cpu = Array(f)
for node in fixed_nodes
f_cpu[3*node-2] = 0.0 # X
f_cpu[3*node-1] = 0.0 # Y
f_cpu[3*node] = 0.0 # Z
end
copyto!(f, f_cpu)
# Return modified operator that zeros fixed DOFs
function K_bc(u)
r = K_op(u)
r_cpu = Array(r)
for node in fixed_nodes
r_cpu[3*node-2] = 0.0
r_cpu[3*node-1] = 0.0
r_cpu[3*node] = 0.0
end
copyto!(r, r_cpu)
return r
end
return K_bc
end
"""
Conjugate Gradient solver (GPU)
"""
function cg_solve_gpu!(
x::CuArray{Float64,1},
A_op::Function,
b::CuArray{Float64,1};
tol=1e-6,
max_iter=1000
)
n = length(x)
# Initial residual
r = b - A_op(x)
p = copy(r)
rsold = dot(r, r)
println("\nConjugate Gradient solver:")
println(" Initial residual: $(sqrt(rsold))")
for iter in 1:max_iter
Ap = A_op(p)
alpha = rsold / dot(p, Ap)
x .+= alpha .* p
r .-= alpha .* Ap
rsnew = dot(r, r)
if iter % 10 == 0 || iter == 1
@printf(" Iter %4d: ||r|| = %.6e\n", iter, sqrt(rsnew))
end
if sqrt(rsnew) < tol
println(" ✅ Converged in $iter iterations")
return x, iter
end
beta = rsnew / rsold
p .= r .+ beta .* p
rsold = rsnew
end
println(" ❌ Did not converge in $max_iter iterations")
return x, max_iter
end
"""
Solve linear elasticity problem on GPU
"""
function solve_elasticity_gpu(problem::ElasticityProblem; tol=1e-6, max_iter=1000)
println("\n" * "="^70)
println("GPU Linear Elasticity Solver")
println("="^70)
# Check CUDA
if !CUDA.functional()
error("CUDA not available!")
end
println("GPU: ", CUDA.name(CUDA.device()))
# Extract mesh data
mesh = physics.mesh
n_nodes = size(mesh.nodes, 2)
n_elems = size(mesh.elements, 2)
n_dofs = 3 * n_nodes
println("\nMesh:")
println(" Nodes: $n_nodes")
println(" Elements: $n_elems")
println(" DOFs: $n_dofs")
println("\nBoundary conditions:")
println(" Fixed nodes: $(length(physics.fixed_nodes))")
println(" Pressure nodes: $(length(physics.pressure_nodes))")
println(" Pressure value: $(physics.pressure_value) Pa")
println("\nMaterial:")
println(" Young's modulus: $(physics.material.E) Pa")
println(" Poisson's ratio: $(physics.material.ν)")
# Transfer to GPU
println("\nTransferring data to GPU...")
nodes_gpu = CuArray(mesh.nodes)
elements_gpu = CuArray(Int32.(mesh.elements))
# Build CSR map
println("Building node-to-elements map...")
node_to_elems = build_node_to_elems_gpu(mesh.elements, n_nodes)
# Initial guess
u_gpu = CUDA.zeros(Float64, n_dofs)
# External force (pressure load)
f_gpu = CUDA.zeros(Float64, n_dofs)
apply_pressure_load!(f_gpu, physics.pressure_nodes, mesh, physics.pressure_value)
println("External force norm: $(norm(Array(f_gpu)))")
# Define stiffness operator K(u) = internal forces
E = physics.material.E
ν = physics.material.ν
function K_op(u)
r = CUDA.zeros(Float64, n_dofs)
compute_residual_gpu!(r, u, nodes_gpu, elements_gpu, node_to_elems, E, ν)
return r
end
# Apply Dirichlet BC
K_bc = apply_dirichlet_bc!(K_op, f_gpu, physics.fixed_nodes)
# Solve: K * u = f
println("\n" * "-"^70)
println("Solving linear system...")
println("-"^70)
u_gpu, n_iter = cg_solve_gpu!(u_gpu, K_bc, f_gpu, tol=tol, max_iter=max_iter)
# Transfer back to CPU
u_cpu = Array(u_gpu)
println("\n" * "="^70)
println("Solution statistics:")
println("="^70)
println(" Max displacement: $(maximum(abs.(u_cpu))) m")
println(" CG iterations: $n_iter")
# Compute final residual
r_final = K_bc(u_gpu) - f_gpu
println(" Final residual: $(norm(Array(r_final)))")
println("\n" * "="^70)
println("✅ GPU elasticity solver complete!")
println("="^70)
return u_cpu
end
end # module
src/io.jl
-518
View File
@@ -1,518 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using HDF5
using LightXML
mutable struct Xdmf <: AbstractResultsWriter
name::String
xml::XMLElement
hdf::HDF5File
hdf_counter::Int
format::String
end
function Xdmf()
return Xdmf(tempname())
end
function h5file(xdmf::Xdmf)
return xdmf.name * ".h5"
end
function xmffile(xdmf::Xdmf)
return xdmf.name * ".xmf"
end
"""
Xdmf(name, version="3.0", overwrite=false)
Initialize a new Xdmf object.
"""
function Xdmf(name::String; version="3.0", overwrite=false)
xdmf = new_element("Xdmf")
h5file = "$name.h5"
xmlfile = "$name.xmf"
if isfile(h5file)
if overwrite
@debug("Result file $h5file exists, removing old file.")
rm(h5file)
else
error("Result file $h5file exists, use Xdmf($name; overwrite=true) to rewrite results")
end
end
if isfile(xmlfile)
if overwrite
@debug("Result file $xmlfile exists, removing old file.")
rm(xmlfile)
else
error("Result file $xmlfile exists, use Xdmf($name; overwrite=true) to rewrite results")
end
end
set_attribute(xdmf, "xmlns:xi", "http://www.w3.org/2001/XInclude")
set_attribute(xdmf, "Version", version)
flag = isfile(h5file) ? "r+" : "w"
hdf = h5open(h5file, flag)
return Xdmf(name, xdmf, hdf, 1, "HDF")
end
"""
get_temporal_collection(xdmf)
Return the basic structure of Xdmf document.
Creates a new TemporalCollection if not found.
Basic structure for XML part of Xdmf file is
<?xml version="1.0" encoding="utf-8"?>
<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="2.1">
<Domain>
<Grid CollectionType="Temporal" GridType="Collection">
</Grid>
</Domain>
</Xdmf>
"""
function get_temporal_collection(xdmf::Xdmf)
domain = find_element(xdmf.xml, "Domain")
grid = nothing
if domain == nothing
domain = new_child(xdmf.xml, "Domain")
grid = new_child(domain, "Grid")
set_attribute(grid, "CollectionType", "Temporal")
set_attribute(grid, "GridType", "Collection")
end
grid = find_element(domain, "Grid")
return grid
end
"""
xdmf_filter(child_elements, child_name)
Returns some spesific child xml element from an array of XMLElement based on,
"Xdmf extensions" see [1] for details.
Parameters
----------
child_elements :: Vector{XMLElement}
A vector of XMLElements where to perform filtering.
child_name :: String
Child element name, maybe containing Xdmf instructions
Returns
-------
nothing if nothing is found, otherwise XMLElement matching to filtering
#Examples
julia> grid1 = new_element("Grid")
julia> add_text(grid1, "I am first grid")
julia> grid2 = new_element("Grid")
julia> add_text(grid2, "I am second grid")
julia> set_attribute(grid2, "Name", "Frame 2")
julia> grid3 = new_element("Grid")
julia> add_text(grid3, "I am third grid")
julia> grids = [grid1, grid2, grid3]
To return second Grid element, one can use
julia> xdmf_filter(grids, "Grid[2]")
To return Grid which has attribute Name="Frame 2", use
julia> xdmf_filter(grids, "Grid[@name=Frame 2]")
To pick last Grid, use [end], e.g.
julia> xdmf_filter(grids, "Grid[end]").
References
----------
[1] http://www.xdmf.org/index.php/XDMF_Model_and_Format
"""
function xdmf_filter(child_elements, child_name)
if '/' in child_name # needs path traversal
return nothing
end
# filter children elements using syntax child[X] -> rename child_name
m = match(r"(\w+)\[(.+)\]", child_name)
if m != nothing
child_name = m[1]
end
# first find any relevant child elements (has same tag)
childs = []
for child in child_elements
if LightXML.name(child) == child_name
push!(childs, child)
end
end
# childs not found at all
length(childs) == 0 && return nothing
# by default return first
m == nothing && return first(childs)
# if [end] return last
m[2] == "end" && return childs[end]
# otherwise try parse int and return nth children from list
parsed_int = tryparse(Int, m[2])
if !isnull(parsed_int)
idx = get(parsed_int)
if (idx > 0) && (idx <= length(childs))
return childs[idx]
else
# wrong index
return nothing
end
end
# [X] is something else than integer, filter children elements using syntax child[@attr=value]
m2 = match(r"@(.+)=(.+)", m[2])
m2 == nothing && throw("Unable to parse: $(m[2])")
attr_name = convert(String, m2[1])
attr_value = convert(String, m2[2])
for child in childs
has_attribute(child, attr_name) || continue
if attribute(child, attr_name) == attr_value
return child
end
end
# nothing found
return nothing
end
"""
traverse(xdmf, x, attr_name)
Traverse XML path. Xdmf filtering can be used, so it's possible to find
data from xml using syntax e.g.
#Example
julia> traverse(xdmf, x, "/Domain/Grid[2]/Grid[@Name=Frame 1]/DataItem")
"""
function traverse(xdmf::Xdmf, x::XMLElement, attr_name::String)
attr_name = strip(attr_name, '/')
if has_attribute(x, attr_name)
return attribute(x, attr_name)
end
childs = child_elements(x)
if '/' in attr_name
items = split(attr_name, '/')
new_item = xdmf_filter(childs, first(items))
if new_item == nothing
@debug("traverse: childs:")
for child in childs
@debug(LightXML.name(child))
end
error("traverse: failed, items = $items, xdmf_filter not find child")
end
new_path = join(items[2:end], '/')
return traverse(xdmf, new_item, new_path)
end
child = xdmf_filter(childs, attr_name)
return child
end
"""
read(xdmf, path)
Read data from Xdmf file.
#Example
Traversing is supported, so one can easily traverse XML tree e.g.
julia> read(xdmf, "/Domain/Grid/Grid[2]/Geometry")
"""
function read(xdmf::Xdmf, path::String)
result = traverse(xdmf, xdmf.xml, path)
if endswith(path, "DataItem")
format = attribute(result, "Format"; required=true)
if format == "HDF"
h5file, path = map(String, split(content(result), ':'))
h5file = dirname(xdmf.name) * "/" * h5file
isfile(h5file) || throw("Xdmf: h5 file $h5file not found!")
return read(xdmf.hdf, path)
else
error("Read from Xdmf, reading from $format not implemented")
end
else
return result
end
end
"""
save!(xdmf)
Save the xdmf file.
"""
function save!(xdmf::Xdmf)
doc = XMLDocument()
set_root(doc, xdmf.xml)
save_file(doc, xmffile(xdmf))
end
function Base.close(xdmf::Xdmf)
close(xdmf.hdf)
end
function new_dataitem(xdmf::Xdmf, path::String, data::Array{T,N}) where {T,N}
dataitem = new_element("DataItem")
datatype = replace("$T", "64" => "")
dimensions = join(reverse(size(data)), " ")
set_attribute(dataitem, "DataType", datatype)
set_attribute(dataitem, "Dimensions", dimensions)
set_attribute(dataitem, "Format", xdmf.format)
if xdmf.format == "HDF"
hdf = basename(h5file(xdmf))
if exists(xdmf.hdf, path)
@debug("Xdmf: $path already existing in h5 file, not overwriting.")
else
write(xdmf.hdf, path, data)
end
add_text(dataitem, "$hdf:$path")
elseif xdmf.format == "XML"
text_data = string(data')
text_data = strip(text_data, ['[', ']'])
text_data = replace(text_data, ';', '\n')
text_data = "\n" * text_data * "\n"
add_text(dataitem, text_data)
else
error("Unsupported Xdmf big data format $(xdmf.format)")
end
return dataitem
end
"""
new_dataitem(xdmf, data)
Create a new DataItem element, hdf path automatically determined.
"""
function new_dataitem(xdmf::Xdmf, data::Array{T,N}) where {T,N}
if xdmf.format == "XML"
# Path can be whatever as XML format does not store to HDF at all
return new_dataitem(xdmf, "/whatever", data)
else
path = "/DataItem_$(xdmf.hdf_counter)"
while exists(xdmf.hdf, path)
xdmf.hdf_counter += 1
path = "/DataItem_$(xdmf.hdf_counter)"
end
return new_dataitem(xdmf, path, data)
end
end
const global xdmf_element_mapping = Dict(
"Poi1" => "Polyvertex",
"Seg2" => "Polyline",
"Tri3" => "Triangle",
"Quad4" => "Quadrilateral",
"Tet4" => "Tetrahedron",
"Pyramid5" => "Pyramid",
"Wedge6" => "Wedge",
"Hex8" => "Hexahedron",
"Seg3" => "Edge_3",
"Tri6" => "Tri_6",
"Quad8" => "Quad_8",
"Tet10" => "Tet_10",
"Pyramid13" => "Pyramid_13",
"Wedge15" => "Wedge_15",
"Hex20" => "Hex_20")
get_xdmf_element_code(::Element{M,Poi1}) where M = 1
get_xdmf_element_code(::Element{M,Seg2}) where M = 2
# get_xdmf_element_code(::Element{Polygon}) = 3
get_xdmf_element_code(::Element{M,Tri3}) where M = 4
get_xdmf_element_code(::Element{M,Quad4}) where M = 5
get_xdmf_element_code(::Element{M,Tet4}) where M = 6
get_xdmf_element_code(::Element{M,Pyr5}) where M = 7
get_xdmf_element_code(::Element{M,Wedge6}) where M = 8
get_xdmf_element_code(::Element{M,Hex8}) where M = 9
# get_xdmf_element_code(::Element{Polyhedron}) = 16
get_xdmf_element_code(::Element{M,Seg3}) where M = 34
get_xdmf_element_code(::Element{M,Quad9}) where M = 35
get_xdmf_element_code(::Element{M,Tri6}) where M = 36
get_xdmf_element_code(::Element{M,Quad8}) where M = 37
get_xdmf_element_code(::Element{M,Tet10}) where M = 38
# get_xdmf_element_code(::Element{Pyr13}) = 39
get_xdmf_element_code(::Element{M,Wedge15}) where M = 40
# get_xdmf_element_code(::Element{Wedge18}) = 41
get_xdmf_element_code(::Element{M,Hex20}) where M = 48
# get_xdmf_element_code(::Element{Hex24}) = 49
get_xdmf_element_code(::Element{M,Hex27}) where M = 50
"""
get_spatial_collection()
Return a SpatialCollection at given time either by creating new one or returning
existing one.
"""
function get_spatial_collection(temporal_collection, time)
for spatial_collection in get_elements_by_tagname(temporal_collection, "Grid")
time_element = find_element(spatial_collection, "Time")
time_value = Meta.parse(attribute(time_element, "Value"; required=true))
isapprox(time_value, time) && return spatial_collection
end
# did not find, create new one
spatial_collection = new_child(temporal_collection, "Grid")
set_attribute(spatial_collection, "GridType", "Collection")
set_attribute(spatial_collection, "Name", "Problems")
set_attribute(spatial_collection, "CollectionType", "Spatial")
time_element = new_child(spatial_collection, "Time")
set_attribute(time_element, "Value", time)
return spatial_collection
end
"""
update_xdmf!(xdmf, problem, time, fields)
Write new fields to Xdmf file.
#Example
To write displacement and temperature fields from p1 at time t=0.0:
julia> update_xdmf!(p1, 0.0, ["displacement", "temperature"])
"""
function update_xdmf!(xdmf::Xdmf, problem::Problem, time::Float64, fields::Vector)
@debug("Xdmf: storing fields $fields of problem $(problem.name) at time $time")
# 1. find domain
xml = xdmf.xml
domain = find_element(xml, "Domain")
if domain == nothing
@debug("Xdmf: Domain not found, creating.")
domain = new_child(xml, "Domain")
end
# 2. find for TemporalCollection
temporal_collection = find_element(domain, "Grid")
if temporal_collection == nothing
@debug("Xdmf: Temporal collection not found, creating.")
temporal_collection = new_child(domain, "Grid")
set_attribute(temporal_collection, "GridType", "Collection")
set_attribute(temporal_collection, "Name", "Time")
set_attribute(temporal_collection, "CollectionType", "Temporal")
end
# 2.1 make sure that Grid element we found really is TemporalCollection
collection_type = attribute(temporal_collection, "CollectionType"; required=true)
@assert collection_type == "Temporal"
spatial_collection = get_spatial_collection(temporal_collection, time)
for frame in get_elements_by_tagname(spatial_collection, "Grid")
frame_name = attribute(frame, "Name")
if frame_name == problem.name
@warn("Xdmf: Already found Grid with name $frame_name for time $time, skipping.")
return
end
end
frame_name = problem.name
@debug("Xdmf: Creating Grid for problem $frame_name")
frame = new_child(spatial_collection, "Grid")
set_attribute(frame, "Name", frame_name)
# 4. save geometry
X_dict = problem("geometry", time)
node_ids = sort(collect(keys(X_dict)))
node_mapping = Dict(j => i for (i, j) in enumerate(node_ids))
X_array = hcat([X_dict[nid] for nid in node_ids]...)
ndim, nnodes = size(X_array)
geom_type = (ndim == 2 ? "XY" : "XYZ")
@debug("Xdmf: Creating geometry, type = $geom_type, number of nodes = $nnodes")
X_dataitem = new_dataitem(xdmf, X_array)
geometry = new_child(frame, "Geometry")
set_attribute(geometry, "Type", geom_type)
add_child(geometry, X_dataitem)
# 5. save topology
mesh_type = "unstructured"
if mesh_type == "unstructured"
element_conn = Int64[]
for element in get_elements(problem)
xdmf_element_code = get_xdmf_element_code(element)
xdmf_element_code > 0 || continue
push!(element_conn, xdmf_element_code)
if xdmf_element_code == 2
push!(element_conn, length(element))
end
for j in get_connectivity(element)
push!(element_conn, node_mapping[j] - 1)
end
end
topology_dataitem = new_dataitem(xdmf, element_conn)
topology = new_child(frame, "Topology")
set_attribute(topology, "TopologyType", "Mixed")
add_child(topology, topology_dataitem)
else
all_elements = get_elements(problem)
nelements = length(all_elements)
element_types = unique(map(get_element_type, all_elements))
nelement_types = length(element_types)
@debug("Xdmf: Saving topology of $nelements elements total, $nelement_types different element types.")
if nelement_types != 1
error("Xdmf: only single type of element supported by structured grid type!")
end
for element_type in element_types
elements = collect(filter_by_element_type(element_type, all_elements))
nelements = length(elements)
@debug("Xdmf: $nelements elements of type $element_type")
sort!(elements, by=get_element_id)
element_ids = map(get_element_id, elements)
element_conn = map(element -> [node_mapping[j] - 1 for j in get_connectivity(element)], elements)
element_conn = hcat(element_conn...)
element_code = split(string(element_type), ".")[end]
topology_dataitem = new_dataitem(xdmf, element_conn)
topology = new_child(frame, "Topology")
set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code])
set_attribute(topology, "NumberOfElements", length(elements))
add_child(topology, topology_dataitem)
end
end
# 6. save requested fields
for field_name in fields
field_dict = problem(field_name, time)
field_center = "Node"
field_node_ids = sort(collect(keys(field_dict)))
if node_ids != field_node_ids
@error("geom node ids = $node_ids")
@error("field node ids = $field_node_ids")
error("!=, geometry does not match with field.")
end
field_dim = length(field_dict[first(field_node_ids)])
if field_dim == 2
@debug("Xdmf: Field dimension = 2, extending to 3")
for nid in field_node_ids
field_dict[nid] = [field_dict[nid]; 0.0]
end
field_dim = 3
end
field_type = Dict(1 => "Scalar", 3 => "Vector", 6 => "Tensor6")[field_dim]
@debug("Xdmf: Saving field $field_name, type = $field_type, dimension = $field_dim, center = $field_center")
field_array = hcat([field_dict[nid] for nid in field_node_ids]...)
field_dataitem = new_dataitem(xdmf, field_array)
attribute = new_child(frame, "Attribute")
set_attribute(attribute, "Name", uppercasefirst(field_name))
set_attribute(attribute, "Center", field_center)
set_attribute(attribute, "AttributeType", field_type)
add_child(attribute, field_dataitem)
end
save!(xdmf)
@debug("Xdmf: all done.")
end
src/materials_plasticity.jl
-123
View File
@@ -1,123 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using ForwardDiff
"""
Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
"""
function find_root!(f, df, x; max_iter=50, norm_acc=1e-9)
converged = false
for i=1:max_iter
dx = -df(x) \ f(x)
x += dx
norm(dx) < norm_acc && (converged = true; break)
end
converged || error("No convergence in radial return!")
return x
end
"""
Equivalent tensile stress.
More info can be found from: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion
Section: Reduced von Mises equation for different stress conditions
"""
function equivalent_stress(stress, ::Type{Val{:type_3d}})
stress_ten = [stress[1] stress[6] stress[5];
stress[6] stress[2] stress[4];
stress[5] stress[4] stress[3]]
stress_dev = stress_ten - 1/3 * tr(stress_ten) * eye(3)
s = vec(stress_dev)
return sqrt(3/2 * dot(s, s))
end
"""
http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm
von mises: plane stress
https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf
"""
function equivalent_stress(stress, ::Type{Val{:type_2d}})
s1, s2, t12 = stress
# Calculating principal stresses
# http://www.engineersedge.com/material_science/principal_vonmises_stress__13418.htm
se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2)
se2 = (s1 + s2)/2 - sqrt(((s1 - s2)/2)^2 + t12^2)
return sqrt(se1^2 -se1*se2 + se2^2)
end
"""
https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf
"""
function yield_function(stress, stress_y, ::Type{Val{:von_mises}}, type_)
equivalent_stress(stress, type_) - stress_y
end
function radial_return(params, dstrain, D, stress_y, stress_base, yield_surface_, type_)
# Creating wrapper for gradient
vm_wrap(stress_) = yield_function(stress_, stress_y, yield_surface_, type_)
dfds = x -> ForwardDiff.gradient(vm_wrap, x)
# Stress rate and total strain
dstress = params[1:end-1]
stress_tot = stress_base + dstress
# Calculating plastic strain rate
dstrain_p = params[end] * dfds(stress_tot)
# Calculating equations
function_1 = dstress - D * (dstrain - dstrain_p)
function_2 = vm_wrap(stress_tot)
[vec(function_1); function_2]
end
function ideal_plasticity!(stress_new, stress_last, dstrain_vec, pstrain, D, params, Dtan, yield_surface_, time, dt, type_)
# Test stress
dstress = vec(D * dstrain_vec)
stress_trial = stress_last + dstress
stress_y = params["yield_stress"]
yield_curr = x -> yield_function(x, stress_y, yield_surface_, type_)
# Calculating and checking for yield
yield = yield_curr(stress_trial)
if isless(yield, 0.0)
stress_new[:] = stress_trial[:]
Dtan[:,:] = D[:,:]
else
# Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ \ f and initial values
f = stress_ -> radial_return(stress_, dstrain_vec, D, stress_y, stress_last, yield_surface_, type_)
df = x -> ForwardDiff.jacobian(f, x)
# Calculating root (two options)
vals = [vec(stress_trial - stress_last); 0.0]
#results = nlsolve(not_in_place(f), vals).zero
results = find_root!(f, df, vals)
# extracting results
dstress = results[1:end-1]
plastic_multiplier = results[end]
# Updating stress
stress_new[:] = stress_last + dstress
# Calculating plastic strain
dfds_ = x -> ForwardDiff.gradient(yield_curr, x)
dep = plastic_multiplier * dfds_(vec(stress_new))
# Equations for consistent tangent matrix can be found from:
# http://homes.civil.aau.dk/lda/continuum/plast.pdf
# equations: 152 & 153
D2g = x -> ForwardDiff.hessian(yield_curr, x)
Dc = (D^-1 + plastic_multiplier * D2g(stress_new))^-1
dfds = dfds_(stress_new)
Dtan[:,:] = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1]
pstrain[:] = plastic_multiplier * dfds
end
end
src/nodal_assembly_structures.jl
-234
View File
@@ -1,234 +0,0 @@
# Nodal Assembly Data Structures
#
# This module provides the inverse mapping needed for efficient nodal assembly:
# Given a node, find all elements touching it and the local node index within each element.
using Tensors
"""
ElementNodeInfo
Information about how a node appears in an element.
# Fields
- `element_id::Int`: Global element ID
- `local_node_idx::Int`: Local node index within the element (1-based)
"""
struct ElementNodeInfo
element_id::Int
local_node_idx::Int
end
"""
NodeToElementsMap
Inverse connectivity mapping: for each node, lists all elements touching it.
# Fields
- `node_to_elements::Vector{Vector{ElementNodeInfo}}`: For node j, gives all elements touching it
- `nnodes::Int`: Total number of nodes in mesh
- `nelements::Int`: Total number of elements in mesh
# Example
```julia
map = NodeToElementsMap(connectivity)
# Get all elements touching node 5
elements_touching_5 = map.node_to_elements[5]
for info in elements_touching_5
println("Node 5 is local node ", info.local_node_idx, " in element ", info.element_id)
end
```
"""
struct NodeToElementsMap
node_to_elements::Vector{Vector{ElementNodeInfo}}
nnodes::Int
nelements::Int
end
"""
NodeToElementsMap(connectivity::Vector{NTuple{N,Int}}) where N
Build inverse mapping from element connectivity.
# Arguments
- `connectivity`: Vector of element connectivity tuples, e.g., [(1,2,3,4), (2,3,5,6), ...]
# Returns
- `NodeToElementsMap`: Inverse mapping structure
# Example
```julia
# Tet4 mesh with 2 elements
connectivity = [(1,2,3,4), (2,3,4,5)]
map = NodeToElementsMap(connectivity)
# Node 2 appears in both elements
@assert length(map.node_to_elements[2]) == 2
```
"""
function NodeToElementsMap(connectivity::Vector{NTuple{N,Int}}) where N
nelements = length(connectivity)
# Find maximum node ID to determine array size
nnodes = maximum(maximum(conn) for conn in connectivity)
# Pre-allocate vectors for each node
node_to_elements = [Vector{ElementNodeInfo}() for _ in 1:nnodes]
# Build inverse mapping
for (elem_id, conn) in enumerate(connectivity)
for (local_idx, global_node_id) in enumerate(conn)
push!(node_to_elements[global_node_id],
ElementNodeInfo(elem_id, local_idx))
end
end
return NodeToElementsMap(node_to_elements, nnodes, nelements)
end
"""
get_node_spider(map::NodeToElementsMap, node_id::Int) -> Vector{Int}
Get the "spider" of a node - all nodes that couple with it (including itself).
This is the union of all nodes in elements touching `node_id`. These are exactly
the nodes for which we need to compute 3×3 stiffness blocks.
# Arguments
- `map`: Node-to-elements mapping
- `node_id`: Node for which to find the spider
# Returns
- `spider_nodes::Vector{Int}`: Sorted unique list of node IDs in the spider
# Example
```julia
# For node j, find all nodes it couples with
spider = get_node_spider(map, j)
# Now compute K_blocks[k] for each k in spider
```
"""
function get_node_spider(map::NodeToElementsMap, node_id::Int,
connectivity::Vector{NTuple{N,Int}}) where N
spider = Set{Int}()
# For each element touching this node
for elem_info in map.node_to_elements[node_id]
# Add all nodes in that element
for node in connectivity[elem_info.element_id]
push!(spider, node)
end
end
return sort(collect(spider))
end
"""
NodalStiffnessContribution{T}
Storage for nodal assembly contribution at a single node.
# Fields
- `node_id::Int`: Global node ID
- `spider_nodes::Vector{Int}`: Node IDs that couple with this node
- `K_blocks::Vector{Tensor{2,3,T}}`: 3×3 stiffness blocks for each spider node
- `f_int::Vec{3,T}`: Internal force at this node
- `f_ext::Vec{3,T}`: External force at this node
# Notes
- `K_blocks[k]` corresponds to `spider_nodes[k]`
- Diagonal block (self-coupling) is included in spider
- All quantities use Tensors.jl types (zero-allocation)
"""
struct NodalStiffnessContribution{T}
node_id::Int
spider_nodes::Vector{Int}
K_blocks::Vector{Tensor{2,3,T,9}}
f_int::Vec{3,T}
f_ext::Vec{3,T}
end
"""
NodalStiffnessContribution(node_id::Int, spider_nodes::Vector{Int}, ::Type{T}=Float64)
Allocate storage for nodal assembly contribution.
# Example
```julia
spider = get_node_spider(map, 5, connectivity)
contrib = NodalStiffnessContribution(5, spider, Float64)
# Now fill in K_blocks, f_int, f_ext during assembly
```
"""
function NodalStiffnessContribution(node_id::Int, spider_nodes::Vector{Int},
::Type{T}=Float64) where T
nspider = length(spider_nodes)
K_blocks = [zero(Tensor{2,3,T}) for _ in 1:nspider]
f_int = zero(Vec{3,T})
f_ext = zero(Vec{3,T})
return NodalStiffnessContribution{T}(node_id, spider_nodes, K_blocks, f_int, f_ext)
end
"""
matrix_vector_product_nodal(contrib::NodalStiffnessContribution,
u::Vector{Vec{3,T}}) -> Vec{3,T}
Compute the matrix-vector product for one node using nodal assembly.
This computes: w_i = sum_j K_ij * u_j for node i
# Arguments
- `contrib`: Nodal stiffness contribution (contains K_blocks for all j in spider)
- `u`: Displacement field at all nodes (Vec{3} per node)
# Returns
- `w_i::Vec{3}`: Result of K_i * u at this node
# Example
```julia
# Assemble contribution for node i
contrib = assemble_nodal_contribution(element_set, node_i, u, time)
# Matrix-free matvec: w_i = K_i * u
w_i = matrix_vector_product_nodal(contrib, u)
```
"""
function matrix_vector_product_nodal(contrib::NodalStiffnessContribution{T},
u::Vector{Vec{3,T}}) where T
w = zero(Vec{3,T})
# Loop over spider nodes (only non-zero columns)
for (k, node_j) in enumerate(contrib.spider_nodes)
K_ij = contrib.K_blocks[k] # 3×3 block
u_j = u[node_j] # 3×1 displacement
# Block matrix-vector product: K_ij is Tensor{2,3}, u_j is Vec{3}
# Use regular matrix-vector multiplication (single contraction)
w += K_ij ⋅ u_j # Tensor{2,3} ⋅ Vec{3} → Vec{3}
end
return w
end
"""
print_spider_info(map::NodeToElementsMap, node_id::Int,
connectivity::Vector{NTuple{N,Int}}) where N
Print diagnostic information about a node's spider for debugging.
"""
function print_spider_info(map::NodeToElementsMap, node_id::Int,
connectivity::Vector{NTuple{N,Int}}) where N
println("Node $node_id Spider Analysis:")
println(" Touches $(length(map.node_to_elements[node_id])) elements")
for elem_info in map.node_to_elements[node_id]
println(" Element $(elem_info.element_id): local node $(elem_info.local_node_idx)")
println(" Connectivity: $(connectivity[elem_info.element_id])")
end
spider = get_node_spider(map, node_id, connectivity)
println(" Spider has $(length(spider)) nodes: $spider")
println(" → Need to compute $(length(spider)) 3×3 blocks")
println(" → Diagonal block at node $node_id")
end
src/postprocess_utils.jl
-180
View File
@@ -1,180 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
Calculate field values to nodal points from Gauss points using least-squares fitting.
"""
function calc_nodal_values!(elements::Vector, field_name, field_dim, time;
F=nothing, nz=nothing, b=nothing, return_F_and_nz=false)
if F == nothing
A = SparseMatrixCOO()
for element in elements
gdofs = get_connectivity(element)
for ip in get_integration_points(element)
detJ = element(ip, time, Val{:detJ})
w = ip.weight*detJ
N = element(ip, time)
add!(A, gdofs, gdofs, w*kron(N', N))
end
end
A = sparse(A)
nz = get_nonzero_rows(A)
A = 1/2*(A + A')
F = ldlt(A[nz,nz])
end
if b == nothing
b = SparseMatrixCOO()
for element in elements
gdofs = get_connectivity(element)
for ip in get_integration_points(element)
if !haskey(ip, field_name)
@warn("integration point does not have field $field_name")
continue
end
detJ = element(ip, time, Val{:detJ})
w = ip.weight*detJ
f = ip(field_name, time)
N = element(ip, time)
for dim=1:field_dim
add!(b, gdofs, [dim], w*f[dim]*N')
end
end
end
b = sparse(b)
end
x = zeros(size(b)...)
x[nz, :] = F \ b[nz, :]
nodal_values = Dict()
for i=1:size(x,1)
nodal_values[i] = vec(x[i,:])
end
update!(elements, field_name, time => nodal_values)
if return_F_and_nz
return F, nz
end
end
"""
Return node ids + vector of values
"""
function get_nodal_vector(elements::Vector, field_name::AbstractString, time::Float64)
f = Dict()
for element in elements
for (c, v) in zip(get_connectivity(element), element(field_name, time))
if haskey(f, c)
@assert isapprox(f[c], v)
end
f[c] = v
end
end
node_ids = sort(collect(keys(f)))
field = [f[nid] for nid in node_ids]
return node_ids, field
end
"""
problem(field_name, X, time)
Interpolate field from a set of elements defined in problem. Here, `X` is the
location inside domain described by elements.
Internally, function loops through all the elements, finding the one containing
the point `X`. After that, using inverse isoparametric mapping, first find
dimensionless coordinates (ξ,η,ζ) of that element corresponding to the location
of point `X` and after that interpolate the values of field under investigation.
Algorithm can be expected to be somewhat slow for big models, but for tests
models the performance is good.
# Examples
Having a problem called `body`, one can query the field `displacement` at
position `X = (1.0, 2.0, 3.0)` and time `t = 1.0`, with the command
```julia
X = (1.0, 2.0, 3.0)
time = 1.0
u = body("displacement", X, time)
```
"""
function (problem::Problem)(field_name, X, time; fillna=NaN)
for element in get_elements(problem)
if inside(element, X, time)
xi = get_local_coordinates(element, X, time)
return element(field_name, xi, time)
end
end
return fillna
end
function (problem::Problem)(field_name, X, time, ::Type{Val{:Grad}}; fillna=NaN)
for element in get_elements(problem)
if inside(element, X, time)
xi = get_local_coordinates(element, X, time)
return element(field_name, xi, time, Val{:Grad})
end
end
return fillna
end
function (solver::Solver)(field_name::AbstractString, X::Vector, time::Float64; fillna=NaN)
for problem in get_problems(solver)
for element in get_elements(problem)
if inside(element, X, time)
xi = get_local_coordinates(element, X, time)
return element(field_name, xi, time)
end
end
end
return fillna
end
""" Calculate area of cross-section. """
function calculate_area(problem::Problem, X=[0.0, 0.0], time=0.0)
A = 0.0
for element in get_elements(problem)
elsize = size(element)
elsize[1] == 2 || error("wrong dimension of problem for area calculation, element size = $elsize")
for ip in get_integration_points(element)
w = ip.weight*element(ip, time, Val{:detJ})
A += w
end
end
return A
end
""" Calculate center of mass of body with respect to X.
https://en.wikipedia.org/wiki/Center_of_mass
"""
function calculate_center_of_mass(problem::Problem, X=[0.0, 0.0, 0.0], time=0.0)
M = 0.0
Xc = zero(X)
for element in get_elements(problem)
for ip in get_integration_points(element)
w = ip.weight*element(ip, time, Val{:detJ})
M += w
rho = haskey(element, "density") ? element("density", ip, time) : 1.0
Xp = element("geometry", ip, time)
Xc += w*rho*(Xp-X)
end
end
return 1.0/M * Xc
end
""" Calculate second moment of mass with respect to X.
https://en.wikipedia.org/wiki/Second_moment_of_area
"""
function calculate_second_moment_of_mass(problem::Problem, X=[0.0, 0.0, 0.0], time=0.0)
n = length(X)
I = zeros(n, n)
for element in get_elements(problem)
for ip in get_integration_points(element)
w = ip.weight*element(ip, time, Val{:detJ})
rho = haskey(element, "density") ? element("density", ip, time) : 1.0
Xp = element("geometry", ip, time) - X
I += w*rho*Xp*Xp'
end
end
return I
end
src/preprocess.jl
-398
View File
@@ -1,398 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
#=
- read meshes from different formats
- reorder connectivity, create element sets, node sets, ...
- create partitions for parallel runs
- renumber elements / nodes
- maybe precheck for bad elements
- check surface normal direction in boundary elements
- orientation of 2d elements
- etc only topology related stuff
=#
mutable struct Mesh
nodes::Dict{Int,Vector{Float64}}
node_sets::Dict{Symbol,Set{Int}}
elements::Dict{Int,Vector{Int}}
element_types::Dict{Int,Symbol}
element_codes::Dict{Int,Symbol}
element_sets::Dict{Symbol,Set{Int}}
surface_sets::Dict{Symbol,Vector{Tuple{Int,Symbol}}}
surface_types::Dict{Symbol,Symbol}
end
function Mesh()
return Mesh(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict())
end
"""
Mesh(m::Dict)
Create a new `Mesh` using data `m`. It is assumed that `m` is in format what
`abaqus_read_mesh` in `AbaqusReader.jl` is returning.
"""
function Mesh(m::Dict)
mesh = Mesh()
mesh.nodes = m["nodes"]
mesh.elements = m["elements"]
mesh.element_types = m["element_types"]
for (k, v) in m["surface_types"]
mesh.surface_types[Symbol(k)] = v
end
for (nset_name, node_ids) in m["node_sets"]
mesh.node_sets[Symbol(nset_name)] = Set(node_ids)
end
for (elset_name, element_ids) in m["element_sets"]
mesh.element_sets[Symbol(elset_name)] = Set(element_ids)
end
for (surfset_name, surfaces) in m["surface_sets"]
mesh.surface_sets[Symbol(surfset_name)] = surfaces
end
return mesh
end
"""
add_node!(mesh, nid, ncoords)
Add node into the mesh. `nid` is node id and `ncoords` are the node
coordinates.
"""
function add_node!(mesh::Mesh, nid::Int, ncoords::Vector{Float64})
mesh.nodes[nid] = ncoords
end
"""
add_nodes!(mesh, nodes)
Add nodes into the mesh.
"""
function add_nodes!(mesh::Mesh, nodes::Dict{Int,Vector{Float64}})
for (nid, ncoords) in nodes
add_node!(mesh, nid, ncoords)
end
end
"""
add_node_to_node_set!(mesh, nid, ncoords)
Add nodes into a node set. `set_name` is the name of the set and `nids...`
are all the node id:s that wants to be added.
"""
function add_node_to_node_set!(mesh::Mesh, set_name, nids...)
if !haskey(mesh.node_sets, set_name)
mesh.node_sets[set_name] = Set{Int}()
end
push!(mesh.node_sets[set_name], nids...)
return
end
"""
create_node_set_from_element_set!(mesh, set_names...)
Create a new node set from the nodes in an element set. ´set_names...´ are all
the set names to be inserted in the function.
"""
function create_node_set_from_element_set!(mesh::Mesh, set_names::String...)
for set_name in set_names
set_name = Symbol(set_name)
@info("Creating node set $set_name from element set")
node_ids = Set{Int}()
for elid in mesh.element_sets[set_name]
push!(node_ids, mesh.elements[elid]...)
end
mesh.node_sets[set_name] = node_ids
end
return
end
"""
create_node_set_from_element_set!(mesh, set_name)
Create a new node set from an element set.
"""
function create_node_set_from_element_set!(mesh::Mesh, set_name::Symbol)
create_node_set_from_element_set!(mesh, string(set_name))
end
"""
add_element!(mesh, elid, eltype, connectivity)
Add an element into the mesh. ´elid´ is the element id, ´eltype´ is the type of
the element and ´connectivity´ is the connectivity of the element.
"""
function add_element!(mesh::Mesh, elid, eltype, connectivity)
mesh.elements[elid] = connectivity
mesh.element_types[elid] = eltype
return nothing
end
"""
add_elements!(mesh, elements)
Add elements into the mesh.
"""
function add_elements!(mesh::Mesh, elements::Dict{Int,Tuple{Symbol,Vector{Int}}})
for (elid, (eltype, elcon)) in elements
add_element!(mesh, elid, eltype, elcon)
end
return nothing
end
"""
add_element_to_element_set!(mesh, set_name, elids...)
Add elements into the mesh. ´set_name´ is the name of the element set and
´elids..´ are id:s of all the elements that wants to be added.
"""
function add_element_to_element_set!(mesh::Mesh, set_name, elids...)
if !haskey(mesh.element_sets, set_name)
mesh.element_sets[set_name] = Set{Int}()
end
push!(mesh.element_sets[set_name], elids...)
end
"""
copy(mesh)
Return a copy of the mesh.
"""
function Base.copy(mesh::Mesh)
mesh2 = Mesh()
mesh2.nodes = copy(mesh.nodes)
mesh2.node_sets = copy(mesh.node_sets)
mesh2.elements = copy(mesh.elements)
mesh2.element_types = copy(mesh.element_types)
mesh2.element_sets = copy(mesh.element_sets)
return mesh2
end
"""
filter_by_element_id(mesh, element_ids)
Filter elements by their id's.
"""
function filter_by_element_id(mesh::Mesh, element_ids::Vector{Int})
mesh2 = copy(mesh)
mesh2.elements = Dict()
for elid in element_ids
if haskey(mesh.elements, elid)
mesh2.elements[elid] = mesh.elements[elid]
end
end
return mesh2
end
"""
filter_by_element_set(mesh, set_name)
Filter elements by an element set.
"""
function filter_by_element_set(mesh::Mesh, set_name)
filter_by_element_id(mesh::Mesh, collect(mesh.element_sets[set_name]))
end
"""
create_element(mesh, id)
Create an element from the mesh by it's id.
"""
function create_element(mesh::Mesh, id::Int)
connectivity = mesh.elements[id]
element_type = getfield(JuliaFEM, mesh.element_types[id])
element = Element(element_type, connectivity)
element.id = id
update!(element, "geometry", mesh.nodes)
return element
end
function create_elements(mesh::Mesh; element_type=nothing)
element_ids = collect(keys(mesh.elements))
if element_type != nothing
filter!(id -> mesh.element_types[id] == element_type, element_ids)
end
elements = [create_element(mesh, id) for id in element_ids]
return elements
end
function create_elements(mesh::Mesh, element_sets::Symbol...; element_type=nothing)
if isempty(element_sets)
element_ids = collect(keys(mesh.elements))
else
element_ids = Set{Int}()
for set_name in element_sets
element_ids = union(element_ids, mesh.element_sets[set_name])
end
end
if element_type != nothing
filter!(id -> mesh.element_types[id] == element_type, element_ids)
end
elements = [create_element(mesh, id) for id in element_ids]
nelements = length(elements)
content = Dict{Symbol,Int}()
for elid in element_ids
eltype = mesh.element_types[elid]
content[eltype] = get(content, eltype, 0) + 1
end
s = join(("$v x $k" for (k, v) in content), ", ")
v = join(element_sets, ", ")
@info("Created $nelements elements ($s) from element set: $v.")
return elements
end
"""
create_elements(mesh::Mesh, element_set::String)
# Examples
Suppose that there is a `mesh` with element set `Body_1`. Creating elements
based on that element set is done then
```julia
create_elements(mesh, "Body_1")
```
"""
function create_elements(mesh::Mesh, element_sets::String...)
return create_elements(mesh, map(Symbol, element_sets)...)
end
"""
find_nearest_nodes(mesh, coords, npts=1; node_set=nothing)
find npts nearest nodes from the mesh and return their id numbers as a list.
"""
function find_nearest_nodes(mesh::Mesh, coords::Vector{Float64}, npts::Int=1; node_set=nothing)
dist = Dict{Int,Float64}()
for (nid, c) in mesh.nodes
if node_set != nothing && !(nid in mesh.node_sets[Symbol(node_set)])
continue
end
dist[nid] = norm(coords - c)
end
s = sort(collect(dist), by=x -> x[2])
nd = s[1:npts] # [(id1, dist1), (id2, dist2), ..., (id_npts, dist_npts)]
node_ids = [n[1] for n in nd]
return node_ids
end
function find_nearest_node(mesh::Mesh, coords::Vector{Float64}; node_set=nothing)
return first(find_nearest_nodes(mesh, coords, 1; node_set=node_set))
end
"""
reorder_element_connectivity!(mesh, mapping)
Apply a new node ordering to elements. JuliaFEM uses the same node ordering as
ABAQUS. If the mesh is parsed from FEM format with some other node ordering,
this function can be used to reorder the nodes.
Parameters
----------
mapping :: Dict{Symbol, Vector{Int}}
e.g. :Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
"""
function reorder_element_connectivity!(mesh::Mesh, mapping::Dict{Symbol,Vector{Int}})
for (elid, eltype) in mesh.element_types
haskey(mapping, eltype) || continue
new_order = mapping[eltype]
element_connectivity = mesh.elements[elid]
new_element_connectivity = element_connectivity[new_order]
mesh.elements[elid] = new_element_connectivity
end
end
function JuliaFEM.Problem(mesh::Mesh, ::Type{P}, name::AbstractString, dimension::Int) where P<:FieldProblem
problem = Problem(P, name, dimension)
problem.elements = create_elements(mesh, name)
return problem
end
function JuliaFEM.Problem(mesh::Mesh, ::Type{P}, name, dimension, parent_field_name) where P<:BoundaryProblem
problem = Problem(P, name, dimension, parent_field_name)
problem.elements = create_elements(mesh, name)
return problem
end
"""
create_coloring!(mesh::Mesh) -> Dict{Int, Int}
Greedy algorithm for coloring a grid such that no two cells with the same node
have the same color.
The returned value is a mapping between an element id and its color.
It is safe to assemble elements with the same color in parallel
"""
function create_coloring(mesh::Mesh)
# Contains the elements that each node contain
cell_containing_node = Dict{Int,Set{Int}}()
for (cellid, nodes) in mesh.elements
for v in nodes
if !haskey(cell_containing_node, v)
cell_containing_node[v] = Set{Int}()
end
push!(cell_containing_node[v], cellid)
end
end
I, J, V = Int[], Int[], Bool[]
for (node, cells) in cell_containing_node
for cell1 in cells # All these cells have a neighboring node
for cell2 in cells
if cell1 != cell2
push!(I, cell1)
push!(J, cell2)
push!(V, true)
end
end
end
end
incidence_matrix = sparse(I, J, V)
# cell -> color of cell
cell_colors = Dict{Int,Int}()
# color -> list of cells
final_colors = Set{Int}[]
occupied_colors = Set{Int}()
# Zero represents no color set yet
for (cellid, _) in mesh.elements
cell_colors[cellid] = 0
end
total_colors = 0
for (cellid, _) in mesh.elements
empty!(occupied_colors)
# loop over neighbors
for r in nzrange(incidence_matrix, cellid)
cell_neighbour = incidence_matrix.rowval[r]
color = cell_colors[cell_neighbour]
if color != 0
push!(occupied_colors, color)
end
end
# occupied colors now contains all the colors we are not allowed to use
free_color = 0
for attempt_color in 1:total_colors
if attempt_color ∉ occupied_colors
free_color = attempt_color
break
end
end
if free_color == 0 # no free color found, need to bump max colors
total_colors += 1
free_color = total_colors
push!(final_colors, Set{Int}())
end
cell_colors[cellid] = free_color
push!(final_colors[free_color], cellid)
end
return cell_colors
end
src/quadrature.jl
-65
View File
@@ -1,65 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE
#
# Gaussian-Legendre quadrature rules consolidated from FEMQuad.jl
# Core API types and abstract interfaces
include("quadrature/api.jl")
# Quadrature data
include("quadrature/quaddata.jl")
# Gauss-Legendre quadrature rules by element topology
include("quadrature/gl_tensor_product.jl") # Segments, quadrilaterals, hexahedra (tensor products)
include("quadrature/gl_triangles.jl") # 2D triangular elements
include("quadrature/gl_tetrahedra.jl") # 3D tetrahedral elements
include("quadrature/gl_wedges.jl") # 3D wedge/prism elements
include("quadrature/gl_pyramids.jl") # 3D pyramid elements
"""
get_rule(order::Int, rules::Symbol...)
Get the first quadrature rule that meets the required order.
"""
function get_rule(order::Int, rules::Vararg{Symbol})
for rule in rules
if get_order(Val{rule}) >= order
return rule
end
end
@warn("No accurate rule enough found, picking last.", order, rules)
return rules[end]
end
"""
integrate_1d(f::Function, rule::Symbol)
Integrate a 1D function using the specified quadrature rule.
"""
function integrate_1d(f::Function, rule::Symbol)
points = get_quadrature_points(Val{rule})
result = sum(w * f(ip) for (w, ip) in points)
return result
end
"""
integrate_2d(f::Function, rule::Symbol)
Integrate a 2D function using the specified quadrature rule.
"""
function integrate_2d(f::Function, rule::Symbol)
points = get_quadrature_points(Val{rule})
result = sum(w * f(ip) for (w, ip) in points)
return result
end
"""
integrate_3d(f::Function, rule::Symbol)
Integrate a 3D function using the specified quadrature rule.
"""
function integrate_3d(f::Function, rule::Symbol)
points = get_quadrature_points(Val{rule})
result = sum(w * f(ip) for (w, ip) in points)
return result
end
src/quadrature/gauss_points.jl
-300
View File
@@ -1,300 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
get_gauss_points!(::Type{T}, ::Type{S}) where {T<:AbstractTopology, S<:Gauss}
-> NTuple{N, Tuple{Float64, Vec{D}}}
Return Gauss quadrature points for topology T with scheme S.
**Zero allocation:** Returns compile-time tuple of (weight, coordinates) pairs.
Coordinates are `Vec{D}` from Tensors.jl for efficient FEM operations.
# Type Parameters
- `T`: Topology type (Triangle, Tetrahedron, Segment, etc.)
- `S`: Gauss quadrature scheme (Gauss{1}, Gauss{2}, etc.)
# Returns
Tuple of `(weight, Vec{D}(ξ))` pairs where:
- `weight`: Integration weight (Float64)
- `Vec{D}(ξ)`: Parametric coordinates as Tensors.jl Vec
# Examples
```julia
# 1-point Gauss for triangle
ips = get_gauss_points!(Triangle, Gauss{1})
# Returns: ((0.5, Vec{2}((1/3, 1/3))),)
# 4-point Gauss for tetrahedron
ips = get_gauss_points!(Tetrahedron, Gauss{1})
# Returns: ((1/24, Vec{3}((0.25, 0.25, 0.25))),)
# Usage in assembly loop (zero allocation):
for (w, ξ) in get_gauss_points!(Triangle, Gauss{2})
# Get basis functions and derivatives using NEW API
N = get_basis_functions(Triangle(), Lagrange{1}(), ξ)
dN = get_basis_derivatives(Triangle(), Lagrange{1}(), ξ)
# Compute Jacobian
detJ = compute_jacobian(ξ)
# Accumulate element matrix
for i in 1:3, j in 1:3
K[i,j] += w * detJ * dot(dN[i], dN[j])
end
end
```
# Performance
- Zero allocations (fully inlined)
- Type-stable (all types known at compile time)
- ~50× faster than runtime dispatch
- Matches golden standard architecture
See also: [`Gauss`](@ref), [`integration_points`](@ref)
"""
function get_gauss_points! end
# ============================================================================
# 1D: Segment
# ============================================================================
# Gauss{1}: 1-point (exact for linear)
@inline function get_gauss_points!(::Type{Segment}, ::Type{Gauss{1}})
return (
(2.0, Vec{1}((0.0,))),
)
end
# Gauss{2}: 2-point (exact for cubic)
@inline function get_gauss_points!(::Type{Segment}, ::Type{Gauss{2}})
a = 1.0 / sqrt(3.0)
return (
(1.0, Vec{1}((-a,))),
(1.0, Vec{1}((a,))),
)
end
# Gauss{3}: 3-point (exact for quintic)
@inline function get_gauss_points!(::Type{Segment}, ::Type{Gauss{3}})
a = sqrt(3.0 / 5.0)
return (
(5.0 / 9.0, Vec{1}((-a,))),
(8.0 / 9.0, Vec{1}((0.0,))),
(5.0 / 9.0, Vec{1}((a,))),
)
end
# ============================================================================
# 2D: Triangle
# ============================================================================
# Gauss{1}: 1-point (exact for linear)
@inline function get_gauss_points!(::Type{Triangle}, ::Type{Gauss{1}})
return (
(0.5, Vec{2}((1 / 3, 1 / 3))),
)
end
# Gauss{2}: 3-point (exact for quadratic)
@inline function get_gauss_points!(::Type{Triangle}, ::Type{Gauss{2}})
return (
(1 / 6, Vec{2}((1 / 6, 1 / 6))),
(1 / 6, Vec{2}((2 / 3, 1 / 6))),
(1 / 6, Vec{2}((1 / 6, 2 / 3))),
)
end
# Gauss{3}: 4-point (exact for cubic)
@inline function get_gauss_points!(::Type{Triangle}, ::Type{Gauss{3}})
a = 1 / 3
b = 1 / 5
c = 3 / 5
return (
(-27 / 96, Vec{2}((a, a))),
(25 / 96, Vec{2}((b, b))),
(25 / 96, Vec{2}((c, b))),
(25 / 96, Vec{2}((b, c))),
)
end
# ============================================================================
# 2D: Quadrilateral (tensor product)
# ============================================================================
# Gauss{1}: 1×1 = 1-point
@inline function get_gauss_points!(::Type{Quadrilateral}, ::Type{Gauss{1}})
return (
(4.0, Vec{2}((0.0, 0.0))),
)
end
# Gauss{2}: 2×2 = 4-point (standard Q1)
@inline function get_gauss_points!(::Type{Quadrilateral}, ::Type{Gauss{2}})
a = 1.0 / sqrt(3.0)
return (
(1.0, Vec{2}((-a, -a))),
(1.0, Vec{2}((a, -a))),
(1.0, Vec{2}((-a, a))),
(1.0, Vec{2}((a, a))),
)
end
# Gauss{3}: 3×3 = 9-point
@inline function get_gauss_points!(::Type{Quadrilateral}, ::Type{Gauss{3}})
a = sqrt(3.0 / 5.0)
w1 = 5.0 / 9.0
w2 = 8.0 / 9.0
return (
(w1 * w1, Vec{2}((-a, -a))),
(w1 * w2, Vec{2}((0.0, -a))),
(w1 * w1, Vec{2}((a, -a))),
(w2 * w1, Vec{2}((-a, 0.0))),
(w2 * w2, Vec{2}((0.0, 0.0))),
(w2 * w1, Vec{2}((a, 0.0))),
(w1 * w1, Vec{2}((-a, a))),
(w1 * w2, Vec{2}((0.0, a))),
(w1 * w1, Vec{2}((a, a))),
)
end
# ============================================================================
# 3D: Tetrahedron
# ============================================================================
# Gauss{1}: 1-point (exact for linear)
@inline function get_gauss_points!(::Type{Tetrahedron}, ::Type{Gauss{1}})
return (
(1 / 6, Vec{3}((0.25, 0.25, 0.25))),
)
end
# Gauss{2}: 4-point (exact for quadratic)
@inline function get_gauss_points!(::Type{Tetrahedron}, ::Type{Gauss{2}})
a = 0.585410196624968
b = 0.138196601125011
return (
(1 / 24, Vec{3}((a, b, b))),
(1 / 24, Vec{3}((b, a, b))),
(1 / 24, Vec{3}((b, b, a))),
(1 / 24, Vec{3}((b, b, b))),
)
end
# Gauss{3}: 5-point (exact for cubic)
@inline function get_gauss_points!(::Type{Tetrahedron}, ::Type{Gauss{3}})
return (
(-4 / 30, Vec{3}((0.25, 0.25, 0.25))),
(9 / 120, Vec{3}((1 / 6, 1 / 6, 1 / 6))),
(9 / 120, Vec{3}((1 / 2, 1 / 6, 1 / 6))),
(9 / 120, Vec{3}((1 / 6, 1 / 2, 1 / 6))),
(9 / 120, Vec{3}((1 / 6, 1 / 6, 1 / 2))),
)
end
# ============================================================================
# 3D: Hexahedron (tensor product)
# ============================================================================
# Gauss{1}: 1×1×1 = 1-point
@inline function get_gauss_points!(::Type{Hexahedron}, ::Type{Gauss{1}})
return (
(8.0, Vec{3}((0.0, 0.0, 0.0))),
)
end
# Gauss{2}: 2×2×2 = 8-point (standard Hex8)
@inline function get_gauss_points!(::Type{Hexahedron}, ::Type{Gauss{2}})
a = 1.0 / sqrt(3.0)
return (
(1.0, Vec{3}((-a, -a, -a))),
(1.0, Vec{3}((a, -a, -a))),
(1.0, Vec{3}((-a, a, -a))),
(1.0, Vec{3}((a, a, -a))),
(1.0, Vec{3}((-a, -a, a))),
(1.0, Vec{3}((a, -a, a))),
(1.0, Vec{3}((-a, a, a))),
(1.0, Vec{3}((a, a, a))),
)
end
# Gauss{3}: 3×3×3 = 27-point
@inline function get_gauss_points!(::Type{Hexahedron}, ::Type{Gauss{3}})
a = sqrt(3.0 / 5.0)
w1 = 5.0 / 9.0
w2 = 8.0 / 9.0
# Generate all 27 combinations
coords_1d = ((-a, w1), (0.0, w2), (a, w1))
result = ntuple(27) do i
ix = (i - 1) % 3 + 1
iy = div(i - 1, 3) % 3 + 1
iz = div(i - 1, 9) + 1
x, wx = coords_1d[ix]
y, wy = coords_1d[iy]
z, wz = coords_1d[iz]
(wx * wy * wz, Vec{3}((x, y, z)))
end
return result
end
# ============================================================================
# 3D: Wedge (Prism) - tensor product of triangle × segment
# ============================================================================
# Gauss{1}: Triangle(1) × Segment(1) = 1-point
@inline function get_gauss_points!(::Type{Wedge}, ::Type{Gauss{1}})
return (
(1.0, Vec{3}((1 / 3, 1 / 3, 0.0))),
)
end
# Gauss{2}: Triangle(3) × Segment(2) = 6-point
@inline function get_gauss_points!(::Type{Wedge}, ::Type{Gauss{2}})
# Triangle points
tri_pts = ((1 / 6, 1 / 6), (2 / 3, 1 / 6), (1 / 6, 2 / 3))
tri_w = 1 / 6
# Segment points
a = 1.0 / sqrt(3.0)
seg_pts = ((-a,), (a,))
seg_w = 1.0
return (
(tri_w * seg_w, Vec{3}((tri_pts[1]..., seg_pts[1][1]))),
(tri_w * seg_w, Vec{3}((tri_pts[1]..., seg_pts[2][1]))),
(tri_w * seg_w, Vec{3}((tri_pts[2]..., seg_pts[1][1]))),
(tri_w * seg_w, Vec{3}((tri_pts[2]..., seg_pts[2][1]))),
(tri_w * seg_w, Vec{3}((tri_pts[3]..., seg_pts[1][1]))),
(tri_w * seg_w, Vec{3}((tri_pts[3]..., seg_pts[2][1]))),
)
end
# ============================================================================
# 3D: Pyramid - special quadrature (not tensor product)
# ============================================================================
# Gauss{1}: 1-point (centroid)
@inline function get_gauss_points!(::Type{Pyramid}, ::Type{Gauss{1}})
return (
(4 / 3, Vec{3}((0.0, 0.0, 0.25))),
)
end
# Gauss{2}: 5-point
@inline function get_gauss_points!(::Type{Pyramid}, ::Type{Gauss{2}})
# Pyramid quadrature is non-trivial due to singularity at apex
a = 0.584237394672177
b = 0.138196601125011
return (
(0.2378, Vec{3}((0.0, 0.0, 0.5))),
(0.2378, Vec{3}((a, 0.0, b))),
(0.2378, Vec{3}((-a, 0.0, b))),
(0.2378, Vec{3}((0.0, a, b))),
(0.2378, Vec{3}((0.0, -a, b))),
)
end
src/quadrature/integration.jl
-143
View File
@@ -1,143 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
# DEPRECATED: This file contains the old integration API.
# New code should use the types from api.jl:
# - AbstractQuadratureRule (replaces AbstractIntegration)
# - QuadraturePoint (replaces IntegrationPoint)
# - GaussLegendre{N} (replaces Gauss{N})
"""
AbstractIntegration
DEPRECATED: Use `AbstractQuadratureRule` instead.
Abstract base type for all numerical integration (quadrature) schemes.
An integration scheme defines how to numerically integrate over a reference element
by specifying integration point locations and weights. Integration schemes are
independent of element topology and interpolation schemes (though the number of
points needed may depend on polynomial order).
# Key Properties
- Integration points (locations in parametric space)
- Weights
- Accuracy order
# Examples
```julia
Gauss{2}() # 2-point Gauss quadrature
Gauss{3}() # 3-point Gauss quadrature
Lobatto{3}() # 3-point Gauss-Lobatto quadrature
Reduced() # Reduced integration (element-dependent)
```
See also: [`Gauss`](@ref), [`Lobatto`](@ref), [`IntegrationPoint`](@ref)
"""
abstract type AbstractIntegration end
"""
IntegrationPoint{D}
DEPRECATED: Use `QuadraturePoint{D,T}` instead.
Represents a single integration point in D-dimensional parametric space.
# Fields
- `ξ::Vec{D, Float64}`: Location in parametric coordinates
- `weight::Float64`: Integration weight
# Migration
```julia
# Old:
ip = IntegrationPoint(Vec(0.0, 0.0), 1.0)
# New:
qp = QuadraturePoint(SVector(0.0, 0.0), 1.0)
# Access: qp.coords instead of ip.ξ
```
"""
struct IntegrationPoint{D}
ξ::Vec{D,Float64}
weight::Float64
end
"""
integration_points(scheme::AbstractIntegration, topology::AbstractTopology)
-> NTuple{N, IntegrationPoint{D}}
Return the integration points and weights for the given integration scheme
applied to the reference element topology.
**Zero allocation:** Returns compile-time sized tuple of IntegrationPoints for
known quadrature rules. Falls back to Vector for dynamic rules.
# Arguments
- `scheme`: Integration scheme (e.g., `Gauss{3}()`)
- `topology`: Reference element topology (e.g., `Tri3()`)
# Returns
Tuple of `IntegrationPoint` with locations ξ and weights.
# Examples
```julia
julia> ips = integration_points(Gauss{1}(), Tri3())
(IntegrationPoint{2}((0.333..., 0.333...), 0.5),)
julia> typeof(ips)
Tuple{IntegrationPoint{2}}
```
"""
function integration_points end
"""
npoints(scheme::AbstractIntegration, topology::AbstractTopology) -> Int
Return the number of integration points for the given scheme and topology.
# Examples
```julia
julia> npoints(Gauss{2}(), Tri3())
3
julia> npoints(Gauss{2}(), Quad4())
4
```
"""
function npoints end
"""
default_integration(topology::Type{<:AbstractTopology{N}}) where N
-> AbstractIntegration
Return the default (recommended) integration scheme for a given topology type.
# Default Rules
- Linear elements (P1): Use minimal integration that's exact for linear basis
- Quadratic elements (P2): Use integration exact for quadratic basis
# Examples
```julia
julia> default_integration(Hexahedron{8})
Gauss{2}() # 2×2×2 = 8 points (exact for trilinear)
julia> default_integration(Tetrahedron{4})
Gauss{1}() # 1 point (exact for linear)
julia> default_integration(Hexahedron{27})
Gauss{3}() # 3×3×3 = 27 points (exact for triquadratic)
```
"""
function default_integration end
# Default integration rules for common topologies
default_integration(::Type{Tetrahedron{4}}) = Gauss{1}()
default_integration(::Type{Tetrahedron{10}}) = Gauss{2}()
default_integration(::Type{Hexahedron{8}}) = Gauss{2}()
default_integration(::Type{Hexahedron{20}}) = Gauss{3}()
default_integration(::Type{Hexahedron{27}}) = Gauss{3}()
default_integration(::Type{Triangle{3}}) = Gauss{1}()
default_integration(::Type{Triangle{6}}) = Gauss{2}()
default_integration(::Type{Quadrilateral{4}}) = Gauss{2}()
default_integration(::Type{Quadrilateral{8}}) = Gauss{3}()
default_integration(::Type{Quadrilateral{9}}) = Gauss{3}()
src/readers.jl
-15
View File
@@ -1,15 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE
#
# Mesh readers consolidated from AbaqusReader.jl and AsterReader.jl
# AbaqusReader - ABAQUS .inp file format
include("readers/keyword_register.jl")
include("readers/parse_mesh.jl")
include("readers/parse_model.jl")
include("readers/create_surface_elements.jl")
include("readers/abaqus_download.jl")
# AsterReader - Code Aster .med file format (requires HDF5)
# include("readers/read_aster_mesh.jl")
# include("readers/read_aster_results.jl")
src/shells/api.jl
-100
View File
@@ -1,100 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
Shell formulation API definitions.
This file defines shell-specific abstract types and formulation theories.
Must be included after core api.jl.
"""
# ============================================================================
# SHELL FORMULATION THEORIES
# ============================================================================
"""
AbstractShellTheory
Abstract type for shell theory variants.
Shell theories differ in how they model transverse shear deformation and thickness effects.
# Concrete Theories
- `ReissnerMindlin`: Thick shells (includes transverse shear)
- `KirchhoffLove`: Thin shells (no transverse shear)
# See Also
- [`ShellFormulation`](@ref)
"""
abstract type AbstractShellTheory end
"""
ReissnerMindlin <: AbstractShellTheory
Reissner-Mindlin shell theory (thick shells, includes shear).
Assumptions:
- Normals to mid-surface remain straight but NOT perpendicular
- Transverse shear deformation included
- Valid for thick shells (h/L > 1/20)
- 5 DOFs per node: 3 displacements + 2 rotations
# Usage
```julia
formulation = ShellFormulation{ReissnerMindlin}()
physics = Physics(
formulation=formulation,
field=DisplacementRotation{3}(),
mesh=shell_mesh,
material=steel
)
```
"""
struct ReissnerMindlin <: AbstractShellTheory end
"""
KirchhoffLove <: AbstractShellTheory
Kirchhoff-Love shell theory (thin shells, no shear).
Assumptions:
- Normals to mid-surface remain straight and perpendicular
- No transverse shear deformation
- Valid for thin shells (h/L < 1/20)
- 3 DOFs per node: 3 displacements (rotations computed from displacements)
# Usage
```julia
formulation = ShellFormulation{KirchhoffLove}()
physics = Physics(
formulation=formulation,
field=Displacement{3}(), # Only displacements, rotations implicit
mesh=shell_mesh,
material=aluminum
)
```
"""
struct KirchhoffLove <: AbstractShellTheory end
"""
ShellFormulation{Theory<:AbstractShellTheory} <: AbstractFormulation
Shell element formulation with theory variant.
# Type Parameter
- `Theory`: Shell theory type (ReissnerMindlin or KirchhoffLove)
# Examples
```julia
# Thick shell (includes shear)
ShellFormulation{ReissnerMindlin}()
# Thin shell (classical theory)
ShellFormulation{KirchhoffLove}()
```
# Fields per Node
- Reissner-Mindlin: 5 DOFs (ux, uy, uz, θx, θy)
- Kirchhoff-Love: 3 DOFs (ux, uy, uz) - rotations implicit
"""
struct ShellFormulation{Theory<:AbstractShellTheory} <: AbstractFormulation end
src/trusses/api.jl
-79
View File
@@ -1,79 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
Truss formulation API definitions.
This file defines truss-specific abstract types and formulation theories.
Must be included after core api.jl.
"""
# ============================================================================
# TRUSS FORMULATION THEORIES
# ============================================================================
"""
AbstractTrussTheory
Abstract type for truss theory variants.
Truss elements carry only axial forces (tension/compression), no bending.
# Concrete Theories
- `SimpleTruss`: Standard 2-node truss (axial force only)
# Future Extensions
- `CableTruss`: Cable elements (tension only, no compression)
- `PretensionedTruss`: Trusses with initial stress
# See Also
- [`TrussFormulation`](@ref)
"""
abstract type AbstractTrussTheory end
"""
SimpleTruss <: AbstractTrussTheory
Simple truss theory (axial force only).
Assumptions:
- Only axial forces (tension/compression)
- No bending moments
- Pin-jointed connections
- 3 DOFs per node in 3D: (ux, uy, uz)
- 2 DOFs per node in 2D: (ux, uy)
# Usage
```julia
formulation = TrussFormulation{SimpleTruss}()
physics = Physics(
formulation=formulation,
field=Displacement{3}(), # 3D truss
mesh=truss_mesh,
material=steel
)
```
"""
struct SimpleTruss <: AbstractTrussTheory end
"""
TrussFormulation{Theory<:AbstractTrussTheory} <: AbstractFormulation
Truss element formulation with theory variant.
# Type Parameter
- `Theory`: Truss theory type (SimpleTruss, CableTruss, etc.)
# Examples
```julia
# Standard truss
TrussFormulation{SimpleTruss}()
```
# Fields per Node
- 3D: 3 DOFs (ux, uy, uz)
- 2D: 2 DOFs (ux, uy)
Use with `Displacement{Dim}` field type.
"""
struct TrussFormulation{Theory<:AbstractTrussTheory} <: AbstractFormulation end
test/runtests.jl.old
-163
View File
@@ -1,163 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM, Test
@testset "JuliaFEM.jl" begin
@testset "test_dirichlet.jl" begin
include("test_dirichlet.jl")
end
@testset "test_elasticity_1d.jl" begin
include("test_elasticity_1d.jl")
end
@testset "test_elasticity_2d_linear_with_surface_load.jl" begin
include("test_elasticity_2d_linear_with_surface_load.jl")
end
@testset "test_elasticity_2d_nonhomogeneous_boundary_conditions.jl" begin
include("test_elasticity_2d_nonhomogeneous_boundary_conditions.jl")
end
@testset "test_elasticity_2d_nonlinear_with_surface_load.jl" begin
include("test_elasticity_2d_nonlinear_with_surface_load.jl")
end
@testset "test_elasticity_2d_plane_stress_stiffness_matrix.jl" begin
include("test_elasticity_2d_plane_stress_stiffness_matrix.jl")
end
@testset "test_elasticity_2d_residual.jl" begin
include("test_elasticity_2d_residual.jl")
end
@testset "test_elasticity_3d_linear_with_surface_load.jl" begin
include("test_elasticity_3d_linear_with_surface_load.jl")
end
@testset "test_elasticity_3d_nonlinear_with_surface_load.jl" begin
include("test_elasticity_3d_nonlinear_with_surface_load.jl")
end
@testset "test_elasticity_3d_unit_block.jl" begin
include("test_elasticity_3d_unit_block.jl")
end
@testset "test_elasticity_forwarddiff.jl" begin
include("test_elasticity_forwarddiff.jl")
end
@testset "test_elasticity_hollow_sphere_with_surface_pressure.jl" begin
include("test_elasticity_hollow_sphere_with_surface_pressure.jl")
end
@testset "test_elasticity_med_pyr5_point_load.jl" begin
include("test_elasticity_med_pyr5_point_load.jl")
end
@testset "test_elasticity_plane_strain.jl" begin
include("test_elasticity_plane_strain.jl")
end
@testset "test_elasticity_pyr5_point_load.jl" begin
include("test_elasticity_pyr5_point_load.jl")
end
@testset "test_elasticity_tet4_stiffness_matrix.jl" begin
include("test_elasticity_tet4_stiffness_matrix.jl")
end
@testset "test_elasticity_tet10_mass_matrix.jl" begin
include("test_elasticity_tet10_mass_matrix.jl")
end
@testset "test_elasticity_tet10_stiffness_matrix.jl" begin
include("test_elasticity_tet10_stiffness_matrix.jl")
end
@testset "test_elasticity_tetra.jl" begin
include("test_elasticity_tetra.jl")
end
@testset "test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl" begin
include("test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl")
end
@testset "test_elasticplastic_3d_linear_with_surface_load.jl" begin
include("test_elasticplastic_3d_linear_with_surface_load.jl")
end
@testset "test_heat_2d_one_element.jl" begin
include("test_heat_2d_one_element.jl")
end
@testset "test_heat_3d.jl" begin
include("test_heat_3d.jl")
end
@testset "test_heat_tet10_convection.jl" begin
include("test_heat_tet10_convection.jl")
end
@testset "test_heat_3d_2.jl" begin
include("test_heat_3d_2.jl")
end
@testset "test_heat.jl" begin
include("test_heat.jl")
end
@testset "test_heat_2.jl" begin
include("test_heat_2.jl")
end
@testset "test_heat_3.jl" begin
include("test_heat_3.jl")
end
@testset "test_heat_3d_two_rings.jl" begin
include("test_heat_3d_two_rings.jl")
end
@testset "test_heat_4.jl" begin
include("test_heat_4.jl")
end
@testset "test_modal_analysis.jl" begin
include("test_modal_analysis.jl")
end
@testset "test_modal_analysis_elasticity.jl" begin
include("test_modal_analysis_elasticity.jl")
end
@testset "test_modal_analysis_elasticity_2.jl" begin
include("test_modal_analysis_elasticity_2.jl")
end
@testset "test_modal_analysis_zero_eigenmodes.jl" begin
include("test_modal_analysis_zero_eigenmodes.jl")
end
@testset "test_mortar.jl" begin
include("test_mortar.jl")
end
@testset "test_mortar_2d.jl" begin
include("test_mortar_2d.jl")
end
@testset "test_mortar_2d_assembly.jl" begin
include("test_mortar_2d.jl")
end
@testset "test_mortar_2d_contact.jl" begin
include("test_mortar_2d_contact.jl")
end
@testset "test_mortar_2d_mesh_tie.jl" begin
include("test_mortar_2d_mesh_tie.jl")
end
@testset "test_mortar_2d_weighted_gap.jl" begin
include("test_mortar_2d_weighted_gap.jl")
end
@testset "test_mortar_3d_mesh_tie_modal.jl" begin
include("test_mortar_3d_mesh_tie_modal.jl")
end
@testset "test_mortar_3d_mesh_tie_two_rings.jl" begin
include("test_mortar_3d_mesh_tie_two_rings.jl")
end
@testset "test_mortar_3d_polygon_clip.jl" begin
include("test_mortar_3d_polygon_clip.jl")
end
@testset "test_postprocess.jl" begin
include("test_postprocess.jl")
end
@testset "test_potential_energy.jl" begin
include("test_potential_energy.jl")
end
@testset "test_problems_contact_3d.jl" begin
include("test_problems_contact_3d.jl")
end
@testset "test_problems_elasticity.jl" begin
include("test_problems_elasticity.jl")
end
@testset "test_problems_mortar_3d.jl" begin
include("test_problems_mortar_3d.jl")
end
@testset "test_problems_mortar_3d_lowlevel.jl" begin
include("test_problems_mortar_3d_lowlevel.jl")
end
@testset "test_solvers_postprocess.jl" begin
include("test_solvers_postprocess.jl")
end
@testset "test_virtual_work.jl" begin
include("test_virtual_work.jl")
end
@testset "test_von_mises_material.jl" begin
include("test_von_mises_material.jl")
end
end
test/runtests_new.jl
-47
View File
@@ -1,47 +0,0 @@
# JuliaFEM Test Suite - New Structure
# Educational testing with Literate.jl
using Test, JuliaFEM
# Configuration
const RUN_TUTORIALS = get(ENV, "JULIAFEM_TEST_TUTORIALS", "true") == "true"
const RUN_UNIT = get(ENV, "JULIAFEM_TEST_UNIT", "false") == "true"
const RUN_OLD = get(ENV, "JULIAFEM_TEST_OLD", "false") == "true"
println("="^70)
println("JuliaFEM Test Suite (New Structure)")
println("="^70)
println("Tutorials: ", RUN_TUTORIALS ? "✓" : "✗")
println("Unit tests: ", RUN_UNIT ? "✓" : "✗")
println("Old tests: ", RUN_OLD ? "✓" : "✗")
println("="^70)
# Tutorial Tests
if RUN_TUTORIALS
@testset "Tutorials" begin
@testset "01_Fundamentals" begin
include("tutorials/01_fundamentals/creating_elements.jl")
include("tutorials/01_fundamentals/reading_gmsh_meshes.jl")
include("tutorials/01_fundamentals/basis_functions.jl")
include("tutorials/01_fundamentals/validation_1element_quad4.jl")
end
end
end
# Unit Tests
if RUN_UNIT
@testset "Unit Tests" begin
@info "No unit tests yet"
end
end
# Old Tests
if RUN_OLD
@warn "Old tests have 49+ failures - see docs/TEST_FIXES_NEEDED.md"
include("runtests.jl") # Original test suite
end
println()
println("="^70)
println("Complete - See docs/TESTING_PHILOSOPHY.md")
println("="^70)
test/test_problems_elasticity_assemble_3d_seg3.jl
-12
View File
@@ -1,12 +0,0 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM, Test
# Test that if 3d continuum problem have some elements we don't know how to
# deal with, raise error with clear message
elements = [Element(Seg3, (1, 2, 3))]
problem = Problem(Elasticity, "test seg3", 3)
add_elements!(problem, elements)
@test_throws ErrorException assemble!(problem, 0.0)