diff --git a/examples/2d_hertz_contact.jl b/examples/2d_hertz_contact.jl deleted file mode 100644 index 62daf37..0000000 --- a/examples/2d_hertz_contact.jl +++ /dev/null @@ -1,214 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -# # 2D Hertz contact problem - -# ![](2d_hertz_contact/model.png) - -# In the example, a cylinder is pressed agains block with a force of 35 kN. -# A similar example can be found from NAFEMS report FENET D3613 (advanced -# finite element contact benchmarks). -# -# Solution for maximum pressure ``p_0`` and contact radius ``a`` is -# ```math -# p_{0} = \sqrt{\frac{FE}{2\pi R}}, \\ -# a = \sqrt{\frac{8FR}{\pi E}}, -# ``` -# where -# -# ```math -# E = \frac{2E_{1}E_{2}}{E_{2}\left(1-\nu_{1}^{2}\right)+E_{1}\left(1-\nu_{2}^{2}\right)}. -# ``` -# -# Substituting values, one gets accurate solution to be ``p_0 = 3585 \;\mathrm{MPa}`` and -# ``a = 6.21 \;\mathrm{mm}``. - -using JuliaFEM, LinearAlgebra - -# Simulation starts by reading the mesh. Model is constructed and meshed using -# SALOME, thus mesh format is .med. Mesh type is quite simple structure, -# containing things like `mesh.nodes`, `mesh.elements` and so on. Keep on mind, -# that Mesh contains only standard Julia types and we think it as a structure -# helping us to construct elements needed in simulation. In principle, we don't -# need to use `Mesh` in simulation anyway if we figure some other way to define -# the geometry for elements. - -datadir = abspath(joinpath(pathof(JuliaFEM), "..", "..", "examples", "2d_hertz_contact")) -meshfile = joinpath(datadir, "hertz_2d_full.med") -mesh = aster_read_mesh(meshfile) -for (elset_name, element_ids) in mesh.element_sets - nel = length(element_ids) - println("Element set $elset_name contains $nel elements.") -end -for (nset_name, node_ids) in mesh.node_sets - nno = length(node_ids) - println("Node set $nset_name contains $nno nodes.") -end -nnodes = length(mesh.nodes) -println("Total number of nodes in mesh: $nnodes") -nelements = length(mesh.elements) -println("Total number of elements in mesh: $nelements") - -# Next, define two bodies. Technically, we could have only one problem and add -# elements from both bodies to the same problem, but defining two different -# problems is recommended for clarity. Plain strain assumption is used. -# To make clear what is happening here: we first create a set of elements -# (elements are in vector called `upper_elements`), then we define new -# problem which type is `Elasticity`, give it some meaningful name (this time -# `cylinder`), and last value 2 means that problems does have two degrees of -# freedom per node. - -upper_elements = create_elements(mesh, "CYLINDER") -update!(upper_elements, "youngs modulus", 70.0e3) -update!(upper_elements, "poissons ratio", 0.3) -upper = Problem(Elasticity, "cylinder", 2) -upper.properties.formulation = :plane_strain -add_elements!(upper, upper_elements) - -lower_elements = create_elements(mesh, "BLOCK") -update!(lower_elements, "youngs modulus", 210.0e3) -update!(lower_elements, "poissons ratio", 0.3) -lower = Problem(Elasticity, "block", 2) -lower.properties.formulation = :plane_strain -add_elements!(lower, lower_elements) - -# Next we define some boundary conditions: creating "boundary" problems goes -# in the same way than defining "field" problems, the only difference is that -# we add extra argument giving what field are we tring to fix. This time, -# we have 2 dofs / node and we fix displacement in direction 2. - -bc_fixed_elements = create_elements(mesh, "FIXED") -update!(bc_fixed_elements, "displacement 2", 0.0) -bc_fixed = Problem(Dirichlet, "fixed", 2, "displacement") -add_elements!(bc_fixed, bc_fixed_elements) - -# Defining symmetry boundary condition goes with the same idea - -bc_sym_23_elements = create_elements(mesh, "SYM23") -update!(bc_sym_23_elements, "displacement 1", 0.0) -bc_sym_23 = Problem(Dirichlet, "symmetry line 23", 2, "displacement") -add_elements!(bc_sym_23, bc_sym_23_elements) - -# Next we define point load. To define that, we first need to find some node -# near the top of cylinder, using function `find_nearest_node`. Then we create -# a new problem, again of type Elasticity. Like told already, we don't need to -# use `Mesh` if we have some other procedure to define the geometry of the -# element (and it's connectivity, of course). So we can directly create an -# element of type `Poi1`, meaning 1-node point element, update it's geometry -# and apply 35.0e3 kN load in negative y-direction: - -nid = find_nearest_node(mesh, [0.0, 100.0]) -load = Problem(Elasticity, "point load", 2) -load.properties.formulation = :plane_strain -load.elements = [Element(Poi1, [nid])] -update!(load.elements, "geometry", mesh.nodes) -update!(load.elements, "displacement traction force 2", -35.0e3) - -# Next, we define another boudary problem, this time the type of problem is -# Contact2D, which is a mortar contact formulation for two dimensions. -# Elements are added using `add_slave_elements!` and `add_master_elements!`. -# Problems, in general, can have some properties defined, like the formulation -# in `Elasticity` (we also have `:plane_stress`). For contact, we need to swap -# normal direction for meshes created by SALOME because in Code Aster, element -# orientation is defined opposite to what is used in ABAQUS, and in JuliaFEM in -# general we follow the same conventions what are used in ABAQUS. - -contact = Problem(Contact2D, "contact", 2, "displacement") -contact.properties.rotate_normals = true -contact_slave_elements = create_elements(mesh, "BLOCK_TO_CYLINDER") -contact_master_elements = create_elements(mesh, "CYLINDER_TO_BLOCK") -add_master_elements!(contact, contact_master_elements) -add_slave_elements!(contact, contact_slave_elements) - -# After all problems are defined, we define some `Analysis`, which can be e.g. -# static analysis, dynamic analysis, modal analysis, linear perturbation -# analysis and so on. Here, the analysis type is `Nonlinear`, which is nonlinear -# quasistatic analysis. In the same manner as we do `add_elements!` to add -# elements to `Problem`, we use `add_problems!` to add problems to analysis. -# Because we are not restricted to some particular input and output formats, -# we "connect" a `ResultsWriter` to our analysis, this time we want to visualize -# results using ParaView, thus we write our results to Xdmf format, which uses -# well defined standards XML and HDF to store model data. - -analysis = Analysis(Nonlinear) -add_problems!(analysis, upper, lower, bc_fixed, bc_sym_23, load, contact) -xdmf = Xdmf("2d_hertz_results"; overwrite=true) -add_results_writer!(analysis, xdmf) - -# In last part, we run the analysis. - -run!(analysis) -close(xdmf) - -# # Results - -# Results are stored in `2d_hertz_results.xmf` and `2d_hertz_results.h5` for -# visual inspection. We can also postprocess results programmatically because -# we are inside a real scripting / programming environment all the time. For -# example, we can integrate the resultant force in normal and tangential direction -# in contact surface to validate our result. - -Rn = 0.0 -Rt = 0.0 -time = 0.0 -for sel in contact_slave_elements - for ip in get_integration_points(sel) - global Rn, Rt - w = ip.weight*sel(ip, time, Val{:detJ}) - n = sel("normal", ip, time) - t = sel("tangent", ip, time) - la = sel("lambda", ip, time) - Rn += w*dot(n, la) - Rt += w*dot(t, la) - end -end - -println("2d hertz contact resultant forces: Rn = $Rn, Rt = $Rt") - -using Test -@test isapprox(Rn, 35.0e3) -@test isapprox(Rt, 0.0) - -# Visualization of the results can be done using ParaView: -# ![](2d_hertz_contact/results_displacement.png) - -# For optimization loops, we want to programmatically find, for example, maximum -# contact pressure. We can, for example, get all the values in nodes: - -lambda = contact("lambda", time) -normal = contact("normal", time) -p0 = 0.0 -p0_acc = 3585.0 -for (nid, n) in normal - lan = dot(n, lambda[nid]) - println("$nid => $lan") - global p0 - p0 = max(p0, lan) -end -p0 = round(p0, digits=2) -rtol = round(norm(p0-p0_acc)/max(p0,p0_acc)*100, digits=2) -println("Maximum contact pressure p0 = $p0, p0_acc = $p0_acc, rtol = $rtol %") - -# To get rough approximation where does the contact open, we can find the element -# from slave contact surface, where contact pressure is zero in the other node -# and something nonzero in the other node. - -a_rad = 0.0 -for element in contact_slave_elements - la1, la2 = element("lambda", time) - p1, p2 = norm(la1), norm(la2) - a, b = isapprox(p1, 0.0), isapprox(p2, 0.0) - if (a && !b) || (b && !a) - X1, X2 = element("geometry", time) - println("Contact opening element geometry: X1 = $X1, X2 = $X2") - println("Contact opening element lambda: la1 = $la1, la2 = $la2") - x11, y11 = X1 - x12, y12 = X2 - global a_rad - a_rad = 1/2*abs(x11+x12) - break - end -end -println("Contact radius: $a_rad") - -# This example briefly described some of the core features of JuliaFEM. diff --git a/examples/2d_hertz_contact/hertz_2d_full.med b/examples/2d_hertz_contact/hertz_2d_full.med deleted file mode 100644 index 2733f0b..0000000 Binary files a/examples/2d_hertz_contact/hertz_2d_full.med and /dev/null differ diff --git a/examples/2d_hertz_contact/model.png b/examples/2d_hertz_contact/model.png deleted file mode 100644 index e7ecd68..0000000 Binary files a/examples/2d_hertz_contact/model.png and /dev/null differ diff --git a/examples/2d_hertz_contact/results_displacement.png b/examples/2d_hertz_contact/results_displacement.png deleted file mode 100644 index 6c6f6f7..0000000 Binary files a/examples/2d_hertz_contact/results_displacement.png and /dev/null differ diff --git a/examples/3d_frame.jl b/examples/3d_frame.jl deleted file mode 100644 index 8bcfb38..0000000 --- a/examples/3d_frame.jl +++ /dev/null @@ -1,90 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -# # Natural frequency analysis of 3d frame structure - -# For general information about Euler-Bernoulli beam theory, see -# [this](https://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory) -# wikipedia page. - -# The model is a 3d frame, shown in picture. - -# ![](3d_frame/model.png) - -using JuliaFEM, LinearAlgebra - -# Reading mesh - -datadir = Pkg.dir("JuliaFEM", "examples", "3d_frame") -mesh = aster_read_mesh(joinpath(datadir, "model.med")) -println("Number of nodes in a model: ", length(mesh.nodes)) - -# Create beam elements. For 3d model, we need to define at least -# [Young's modulus](https://en.wikipedia.org/wiki/Young%27s_modulus), -# [shear modulus](https://en.wikipedia.org/wiki/Shear_modulus), -# [density](https://en.wikipedia.org/wiki/Density) -# cross-section area, moment of inertia in local coordinate -# system and polar moment of inertia. - -beam_elements = create_elements(mesh, "FRAME") -@info("Number of elements: ", length(beam_elements)) -update!(beam_elements, "youngs modulus", 210.0e6) -update!(beam_elements, "shear modulus", 84.0e6) -update!(beam_elements, "density", 7850.0e-3) -update!(beam_elements, "cross-section area", 20.0e-2) -update!(beam_elements, "torsional moment of inertia 1", 10.0e-5) -update!(beam_elements, "torsional moment of inertia 2", 10.0e-5) -update!(beam_elements, "polar moment of inertia", 30.0e-5) - -# The direction of beam is defined in same way than in ABAQUS. -# That is, we have a tangent direction and one normal direction. -# The third direction is then cross product of tangent and normal. -# Because the second area moment is same in both directions, we can -# choose normal direction freely. - -for element in beam_elements - X1, X2 = element("geometry", 0.0) - t = (X2-X1)/norm(X2-X1) - I = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] - k = indmax([norm(cross(t, I[:,k])) for k in 1:3]) - n = cross(t, I[:,k])/norm(cross(t, I[:,k])) - update!(element, "normal", n) -end - -# Create boundary conditions: fix all degrees of freedom for nodes in -# a set FIXED. Here we first create elements of type `Poi1` for each -# node j in set FIXED, update geometry field and then create new fields -# `fixed displacmeent 1`, `fixed displacement 2`, and so on, where the -# displacement / rotation is prescribed. - -bc_elements = [Element(Poi1, [j]) for j in mesh.node_sets[:FIXED]] -update!(bc_elements, "geometry", mesh.nodes) -for i=1:3 - update!(bc_elements, "fixed displacement $i", 0.0) - update!(bc_elements, "fixed rotation $i", 0.0) -end - -# Create a problem, containing beam elements and boundary conditions: - -frame = Problem(Beam, "3d frame", 6) -add_elements!(frame, beam_elements) -add_elements!(frame, bc_elements) - -# Perform modal analysis - -analysis = Analysis(Modal) -xdmf = Xdmf(joinpath(datadir, "3d_frame_results"); overwrite=true) -add_results_writer!(analysis, xdmf) -add_problems!(analysis, frame) -run!(analysis) -close(xdmf) - -# Each `Analysis` can have properties, e.g. time, maximum number of iterations, -# convergence tolerance and so on. Eigenvalues of calculation are stored as a -# properties of analysis: - -freqs = sqrt.(step.properties.eigvals) / (2*pi) -println("Natural frequencies [Hz]: $(round.(freqs, 2))") - -# [![mode5](3d_frame/natfreq.png)](https://www.youtube.com/watch?v=GzktCqeASmo) - diff --git a/examples/3d_frame/model.med b/examples/3d_frame/model.med deleted file mode 100644 index a3f9c50..0000000 Binary files a/examples/3d_frame/model.med and /dev/null differ diff --git a/examples/3d_frame/model.png b/examples/3d_frame/model.png deleted file mode 100644 index 616fb37..0000000 Binary files a/examples/3d_frame/model.png and /dev/null differ diff --git a/examples/3d_frame/natfreq.png b/examples/3d_frame/natfreq.png deleted file mode 100644 index e91024b..0000000 Binary files a/examples/3d_frame/natfreq.png and /dev/null differ diff --git a/examples/academic_matrix_extraction/README.md b/examples/academic_matrix_extraction/README.md deleted file mode 100644 index 4c78e89..0000000 --- a/examples/academic_matrix_extraction/README.md +++ /dev/null @@ -1,123 +0,0 @@ -# Academic Example: Matrix Extraction for External Solvers - -**Addresses Issue #183**: Demonstrates using JuliaFEM for spatial discretization only, extracting matrices for external solvers. - -## The Three Requirements - -This example demonstrates exactly what was requested in Issue #183: - -### a) Discretize space (into a mesh) -- Shows programmatic mesh generation -- Element connectivity accessible -- Compatible with Gmsh mesh files - -### b) Assemble the stiffness matrix -- Assembly framework demonstrated -- Currently works with Dirichlet BC -- Heat/Elasticity coming in Phase 2 (2-4 months) - -### c) Get back vectors and matrices -- Extract `K` (stiffness), `M` (mass), `f` (force) as standard Julia types -- `SparseMatrixCSC{Float64,Int64}` and `Vector{Float64}` -- Direct compatibility with entire Julia ecosystem - -## Quick Start - -```bash -cd examples/academic_matrix_extraction -julia --project=../.. academic_example.jl -``` - -## What You Get - -After assembly, matrices are extracted as: - -```julia -K = problem.assembly.K # Stiffness matrix (sparse) -M = problem.assembly.M # Mass matrix (sparse) -f = problem.assembly.f # Force vector -``` - -These are standard Julia types that work with: - -### DifferentialEquations.jl (Transient Problems) - -```julia -using DifferentialEquations - -function fem_ode!(du, u, p, t) - K, M, f = p - du .= M \ (-K * u .+ f) -end - -u0 = zeros(N) -prob = ODEProblem(fem_ode!, u0, (0.0, 1.0), (K, M, f)) -sol = solve(prob, Tsit5()) -``` - -### LinearSolve.jl (Steady-State) - -```julia -using LinearSolve - -prob = LinearProblem(K, f) -sol = solve(prob, KrylovJL_GMRES()) -``` - -### Krylov.jl (Iterative Methods) - -```julia -using Krylov - -u, stats = gmres(K, f; atol=1e-10, rtol=1e-8) -``` - -### Custom Research Solvers - -```julia -using SparseArrays, LinearAlgebra - -u = K \ f # Direct solve -L = cholesky(K) # Factorization -λ, v = eigs(K, M) # Eigenvalue analysis -``` - -## Current Status - -**What Works NOW:** -- ✅ Mesh generation and element connectivity -- ✅ Matrix extraction API (`problem.assembly.K`, `.M`, `.f`) -- ✅ Dirichlet boundary conditions -- ✅ Integration with Julia solver ecosystem - -**Coming in Phase 2 (2-4 months):** -- ⏳ Heat equation problem type -- ⏳ Elasticity problem type -- ⏳ Full assembly for physics problems -- ⏳ 40-130x performance improvement - -## Why Phase 2? - -JuliaFEM is undergoing architecture refactoring (Nov 2025): -- Replacing Dict-based fields (100x performance penalty) with type-stable system -- New immutable element architecture (already 40-130x faster) -- Heat/Elasticity problem types depend on old system -- Being restored with new architecture - -## See Also - -- **Issue #183**: Original request from Chris Rackauckas (2017) -- **examples/gmsh_heat_equation/**: Full workflow with Gmsh mesh files -- **docs/book/gmsh_tutorial.md**: Comprehensive step-by-step tutorial -- **llm/ARCHITECTURE.md**: Architecture design and roadmap -- **docs/blog/immutability_performance.md**: Performance analysis - -## Academic Use Case - -Perfect for research where you need: -1. Spatial discretization (FEM assembly) -2. Custom time integration schemes -3. Novel solver algorithms -4. Integration with other Julia packages - -JuliaFEM handles the messy FEM assembly; you control the solving. diff --git a/examples/academic_matrix_extraction/academic_example.jl b/examples/academic_matrix_extraction/academic_example.jl deleted file mode 100644 index 998f52b..0000000 --- a/examples/academic_matrix_extraction/academic_example.jl +++ /dev/null @@ -1,255 +0,0 @@ -#!/usr/bin/env julia -# Academic Example: Matrix Extraction for External Solvers -# Addresses Issue #183 - Demonstrates a), b), and c) -# -# Shows how to: -# a) Discretize space (tetrahedral/triangular mesh) -# b) Assemble stiffness matrix -# c) Get back vectors and matrices for external solvers -# -# This is a WORKING example using Dirichlet BC (which is currently available) - -using JuliaFEM -using LinearAlgebra -using SparseArrays - -println("="^80) -println("Academic Example: FEM Matrix Extraction (Issue #183)") -println("="^80) -println() -println("This demonstrates the three requirements:") -println(" a) Discretize space into mesh") -println(" b) Assemble stiffness matrix") -println(" c) Extract vectors/matrices for external solvers") -println() -println("-"^80) -println() - -# ============================================================================= -# Step (a): Discretize Space - Create Mesh -# ============================================================================= - -println("Step (a): Spatial Discretization") -println("-"^80) - -# Create a simple 2D triangular mesh programmatically -# Unit square divided into triangles -# -# 4 ------- 3 -# | \ / | -# | \ / | -# | / \ | -# | / \ | -# 1 ------- 2 - -nodes = Dict{Int64,Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [1.0, 1.0], - 4 => [0.0, 1.0], - 5 => [0.5, 0.5] # Center node -) - -# Element connectivity (node IDs for each triangle) -elements = [ - ("Tri3", [1, 2, 5]), - ("Tri3", [2, 3, 5]), - ("Tri3", [3, 4, 5]), - ("Tri3", [4, 1, 5]) -] - -# Boundary nodes (for BC application) -left_boundary_nodes = [1, 4] - -println("✓ Mesh created:") -println(" Nodes: $(length(nodes))") -println(" Elements: $(length(elements)) triangles") -println(" Boundary nodes: $(length(left_boundary_nodes)) (left edge)") -println() -println(" Mesh topology:") -println(" Element 1: nodes $(elements[1][2])") -println(" Element 2: nodes $(elements[2][2])") -println(" Element 3: nodes $(elements[3][2])") -println(" Element 4: nodes $(elements[4][2])") -println() - -# ============================================================================= -# Step (b): Assemble Stiffness Matrix - Create Problem -# ============================================================================= - -println("Step (b): Stiffness Matrix Assembly") -println("-"^80) - -# Create a simple Laplacian problem: -∇²u = f -# We'll construct the stiffness matrix K and force vector f directly -# to demonstrate matrix extraction without needing Heat problem type - -N = length(nodes) # 5 nodes -println("Assembling $(N)×$(N) Laplacian system...") - -# For this simple example, construct a basic 1D Laplacian-like system -# This represents a discretized -∇²u = f problem - -# Simple tridiagonal stiffness matrix (like 1D Laplacian) -# K = [-2 1 0 0 0] -# [ 1 -2 1 0 0] -# [ 0 1 -2 1 0] -# [ 0 0 1 -2 1] -# [ 0 0 0 1 -2] -K = spdiagm(0 => -2.0 * ones(N), - 1 => ones(N - 1), - -1 => ones(N - 1)) - -# Force vector (right-hand side) -f = ones(N) # Uniform source term - -println("✓ System assembled:") -println(" K: $(N)×$(N) sparse matrix ($(nnz(K)) non-zeros)") -println(" f: $(N)-element force vector") -println() -println(" Matrix K (Laplacian-like stiffness):") -println(" $(Matrix(K))") -println() -println(" Force vector f:") -println(" $f") -println() - -# ============================================================================= -# Step (c): Extract Matrices and Solve -# ============================================================================= - -println("Step (c): Matrix Extraction and Solution") -println("-"^80) -println() -println("The assembled system is K * u = f") -println() -println("Solving using direct method: u = K \\ f") -println() - -# Solve the system -u = K \ f - -println("✓ Solution computed!") -println() -println("Solution vector u:") -for i in 1:N - println(" u[$i] = $(u[i])") -end -println() - -# Verify solution -residual = K * u - f -residual_norm = norm(residual) -println("Verification:") -println(" Residual ||K*u - f|| = $residual_norm") -println(" $(residual_norm < 1e-10 ? "✓" : "✗") Solution is $(residual_norm < 1e-10 ? "correct" : "incorrect")") -println() - -println("This demonstrates Issue #183 requirement (c):") -println(" ✓ Extracted K (stiffness matrix) as SparseMatrixCSC{Float64,Int64}") -println(" ✓ Extracted f (force vector) as Vector{Float64}") -println(" ✓ Solved K * u = f to get solution vector u") -println(" ✓ Solution available for further analysis or time integration") -println() - -# ============================================================================= -# Step (d): Integration with External Solvers -# ============================================================================= - -println("Step (d): Using Matrices with External Solvers") -println("-"^80) -println() -println("The matrices K and f are standard Julia types compatible with:") -println() -println("1. DifferentialEquations.jl (for transient problems):") -println(" ------------------------------------------------------") -println(" using DifferentialEquations") -println(" ") -println(" # Define ODE system: M * du/dt = -K * u + f") -println(" function fem_ode!(du, u, p, t)") -println(" K, M, f = p") -println(" du .= M \\ (-K * u .+ f)") -println(" end") -println(" ") -println(" u0 = zeros(N) # Initial condition") -println(" tspan = (0.0, 1.0)") -println(" prob = ODEProblem(fem_ode!, u0, tspan, (K, M, f))") -println(" sol = solve(prob, Tsit5())") -println() -println("2. LinearSolve.jl (for steady-state problems):") -println(" ---------------------------------------------") -println(" using LinearSolve") -println(" ") -println(" # Solve K * u = f") -println(" prob = LinearProblem(K, f)") -println(" sol = solve(prob, KrylovJL_GMRES())") -println(" u_solution = sol.u") -println() -println("3. Krylov.jl (for iterative methods):") -println(" ------------------------------------") -println(" using Krylov") -println(" ") -println(" # Direct iterative solve") -println(" u, stats = gmres(K, f; atol=1e-10, rtol=1e-8)") -println(" ") -println(" # With preconditioner") -println(" using IncompleteLU") -println(" P = ilu(K, τ=0.01)") -println(" u, stats = gmres(K, f; M=P, atol=1e-10)") -println() -println("4. Custom research solvers:") -println(" -------------------------") -println(" # Matrices are standard SparseArrays, so any Julia") -println(" # linear algebra works:") -println(" ") -println(" using SparseArrays, LinearAlgebra") -println(" u = K \\ f # Direct solve (for small systems)") -println(" L = cholesky(K) # Factorization (if K is SPD)") -println(" λ, v = eigs(K, M) # Eigenvalue analysis") -println() - -# ============================================================================= -# Summary -# ============================================================================= - -println("="^80) -println("Summary: Issue #183 Requirements - ALL DEMONSTRATED") -println("="^80) -println() -println("✓ (a) Discretize space:") -println(" • Mesh created: 5 nodes, 4 triangular elements") -println(" • Element connectivity accessible") -println(" • Gmsh .msh file import available (see examples/gmsh_heat_equation/)") -println() -println("✓ (b) Assemble stiffness matrix:") -println(" • System assembled: K (5×5 sparse), f (5 elements)") -println(" • Matrix structure: Laplacian-like (tridiagonal)") -println(" • 9 non-zero entries in K") -println() -println("✓ (c) Extract vectors/matrices:") -println(" • K extracted as SparseMatrixCSC{Float64,Int64}") -println(" • f extracted as Vector{Float64}") -println(" • Solution computed: u = K \\ f") -println(" • Residual verified: ||K*u - f|| = $residual_norm") -println() -println("Solution:") -println(" u = $u") -println() -println("Current Status:") -println(" [WORKING] Matrix extraction API and data structures") -println(" [WORKING] Mesh generation and element creation") -println(" [WORKING] Dirichlet boundary conditions") -println(" [PENDING] Heat/Elasticity problem types (Phase 2)") -println() -println("Next Steps:") -println(" 1. See examples/gmsh_heat_equation/ for workflow with Gmsh") -println(" 2. See docs/book/gmsh_tutorial.md for comprehensive tutorial") -println(" 3. Architecture refactoring underway (40-130x performance improvement)") -println(" 4. Heat equation example will be fully functional in Phase 2") -println() -println("Reference:") -println(" • Issue: https://github.com/JuliaFEM/JuliaFEM.jl/issues/183") -println(" • Architecture: llm/ARCHITECTURE.md") -println(" • Performance: docs/blog/immutability_performance.md") -println() -println("="^80) diff --git a/examples/generate_stiffness_matrices.jl b/examples/generate_stiffness_matrices.jl deleted file mode 100644 index dc9e9c4..0000000 --- a/examples/generate_stiffness_matrices.jl +++ /dev/null @@ -1,19 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -# # Generating local matrices for problems - -using JuliaFEM - -# Plane stress Quad4 element with linear material model: - -X = Dict(1 => [0.0, 0.0], 2 => [2.0, 0.0], 3 => [2.0, 2.0], 4 => [0.0, 2.0]) -element = Element(Quad4, [1, 2, 3, 4]) -update!(element, "geometry", X) -update!(element, "youngs modulus", 288.0) -update!(element, "poissons ratio", 1/3) -problem = Problem(Elasticity, "test problem", 2) -problem.properties.formulation = :plane_stress -add_elements!(problem, [element]) -assemble!(problem, 0.0) -K = round.(Matrix(problem.assembly.K), 5) diff --git a/examples/gmsh_heat_equation/QUICK_START.md b/examples/gmsh_heat_equation/QUICK_START.md deleted file mode 100644 index bcb7fe6..0000000 --- a/examples/gmsh_heat_equation/QUICK_START.md +++ /dev/null @@ -1,35 +0,0 @@ -# Gmsh Heat Equation Example - -Complete workflow demonstration addressing [Issue #183](https://github.com/JuliaFEM/JuliaFEM.jl/issues/183). - -## Quick Links - -- **Working Example:** `gmsh_heat_equation.jl` -- **Gmsh Geometry:** `unit_square.geo` -- **Comprehensive Tutorial:** `../../docs/book/gmsh_tutorial.md` - -## What's Here - -This example shows: - -1. Mesh generation with Gmsh -2. Loading mesh into JuliaFEM -3. FEM assembly (K and M matrices) -4. Extracting matrices for external solvers -5. Integration with DifferentialEquations.jl - -## Academic Usage (Issue #183) - -If you want JuliaFEM for **discretization only** (not built-in physics): - -```julia -# After assembly: -K = problem.assembly.K # Stiffness matrix -M = problem.assembly.M # Mass matrix -f = problem.assembly.f # Force vector - -# ODE system: M * du/dt = -K * u + f -# Now use with your own solver! -``` - -See tutorial for complete explanation. diff --git a/examples/gmsh_heat_equation/README.md b/examples/gmsh_heat_equation/README.md deleted file mode 100644 index e4a18e3..0000000 --- a/examples/gmsh_heat_equation/README.md +++ /dev/null @@ -1,74 +0,0 @@ -# Heat Equation Example: From Gmsh to Physics -# Addresses Issue #183: Academic usage without built-in physics - -This example demonstrates the complete workflow: -1. Generate mesh using Gmsh -2. Load mesh into JuliaFEM -3. Assemble stiffness matrix and mass matrix -4. Extract matrices for external solvers (e.g., DifferentialEquations.jl) -5. Solve the heat equation - -## Problem Statement - -Solve the transient heat equation on a unit square: - -``` -∂u/∂t = α∇²u + f(x,y,t) -``` - -with boundary conditions: -- u = 0 on left edge (Dirichlet) -- ∂u/∂n = 0 on other edges (Neumann, natural BC) - -Initial condition: u(x,y,0) = sin(πx)sin(πy) - -## Quick Start - -### 1. Generate mesh - -```bash -gmsh -2 unit_square.geo -o unit_square.msh -``` - -This creates a triangular mesh of the unit square. - -### 2. Run the example - -```bash -julia --project gmsh_heat_equation.jl -``` - -## What You Get - -The example shows how to: -- Load Gmsh mesh files -- Create FEM elements with material properties -- Assemble global stiffness matrix K and mass matrix M -- Apply Dirichlet boundary conditions -- Extract the resulting ODE system: M du/dt = -K u + f -- Solve using your own time integrator - -## For Academic Users (Issue #183) - -If you want to use JuliaFEM just for discretization (not the built-in physics): - -```julia -# After assembly, extract the matrices: -K = problem.assembly.K # Stiffness matrix (SparseMatrixCSC) -M = problem.assembly.M # Mass matrix (SparseMatrixCSC) -f = problem.assembly.f # Force vector - -# Now use these with DifferentialEquations.jl, Krylov.jl, etc. -# The ODE system is: M * du/dt = -K * u + f -``` - -## Files - -- `unit_square.geo` - Gmsh geometry definition -- `gmsh_heat_equation.jl` - Complete working example -- `README.md` - This file - -## See Also - -- Tutorial: `docs/book/gmsh_tutorial.md` (comprehensive step-by-step) -- Issue #183: https://github.com/JuliaFEM/JuliaFEM.jl/issues/183 diff --git a/examples/gmsh_heat_equation/gmsh_heat_equation.jl b/examples/gmsh_heat_equation/gmsh_heat_equation.jl deleted file mode 100644 index 9720d34..0000000 --- a/examples/gmsh_heat_equation/gmsh_heat_equation.jl +++ /dev/null @@ -1,280 +0,0 @@ -#!/usr/bin/env julia -# Heat Equation Example: From Gmsh Mesh to FEM Assembly -# Addresses Issue #183 - Academic usage for spatial discretization -# -# NOTE: This is a DEMONSTRATION of the workflow. The Heat problem type -# is currently disabled in JuliaFEM pending architecture refactoring. -# This shows the STRUCTURE of how to go from mesh → assembly → matrices. -# -# Problem: ∂u/∂t = α∇²u + f(x,y,t) on unit square -# Boundary: u = 0 on left edge, ∂u/∂n = 0 elsewhere -# Initial: u(x,y,0) = sin(πx)sin(πy) - -# using JuliaFEM # Commented out - Heat problem not yet available -using LinearAlgebra -using SparseArrays - -println("="^80) -println("Heat Equation: Mesh Generation Demonstration") -println("Issue #183: Workflow for Academic Usage") -println("="^80) -println() -println("NOTE: This demonstrates the WORKFLOW structure.") -println(" Full Heat problem assembly coming in Phase 2 refactoring.") -println() -println("What this shows:") -println(" 1. Mesh generation (programmatic or from Gmsh)") -println(" 2. Element connectivity structure") -println(" 3. Boundary identification") -println(" 4. How matrices WOULD be assembled (K, M, f)") -println(" 5. Integration with external solvers (DifferentialEquations.jl)") -println() -println("-"^80) -println() - -# ============================================================================= -# Step 1: Generate and Load Mesh -# ============================================================================= - -println("Step 1: Mesh Generation") -println("-"^80) - -# Check if mesh exists, otherwise generate it -mesh_file = "unit_square.msh" -geo_file = "unit_square.geo" - -if !isfile(mesh_file) - if !isfile(geo_file) - error("Geometry file $geo_file not found. Please create it first.") - end - - println("Generating mesh with Gmsh...") - run(`gmsh -2 $geo_file -o $mesh_file`) - println("✓ Mesh generated: $mesh_file") -else - println("✓ Using existing mesh: $mesh_file") -end - -# Load mesh (Note: Gmsh reader needs to be implemented or use existing) -# For now, we'll create a simple unit square mesh programmatically -println("\nCreating unit square mesh...") - -# Simple structured mesh: 10×10 grid -n = 10 # divisions per side -nodes = Dict{Int,Vector{Float64}}() -node_id = 1 -for j in 0:n - for i in 0:n - x = i / n - y = j / n - nodes[node_id] = [x, y, 0.0] - node_id += 1 - end -end - -# Create triangular elements (two triangles per square) -elements = Vector{Tuple{Symbol,Vector{Int}}}() -element_sets = Dict{String,Vector{Int}}() -body_elements = Int[] -left_elements = Int[] -right_elements = Int[] -bottom_elements = Int[] -top_elements = Int[] - -elem_id = 1 -for j in 1:n - for i in 1:n - # Node indices for square [i,j] - n1 = (j - 1) * (n + 1) + i # bottom-left - n2 = (j - 1) * (n + 1) + i + 1 # bottom-right - n3 = j * (n + 1) + i + 1 # top-right - n4 = j * (n + 1) + i # top-left - - # Triangle 1: [n1, n2, n3] - push!(elements, (:Tri3, [n1, n2, n3])) - push!(body_elements, elem_id) - elem_id += 1 - - # Triangle 2: [n1, n3, n4] - push!(elements, (:Tri3, [n1, n3, n4])) - push!(body_elements, elem_id) - elem_id += 1 - end -end - -# Boundary edges (1D line elements for visualization/BC) -# Left edge: x = 0 -for j in 1:n - n1 = (j - 1) * (n + 1) + 1 - n2 = j * (n + 1) + 1 - push!(elements, (:Seg2, [n1, n2])) - push!(left_elements, elem_id) - elem_id += 1 -end - -element_sets["body"] = body_elements -element_sets["left"] = left_elements - -println("✓ Mesh created: $(length(nodes)) nodes, $(length(body_elements)) triangles") -println() - -# ============================================================================= -# Step 2: Element Connectivity and Boundary Identification -# ============================================================================= - -println("Step 2: Element Connectivity") -println("-"^80) - -println("✓ Body elements: $(length(body_elements)) triangles") -println("✓ Boundary elements: $(length(left_elements)) edges on left") -println() -println("Element connectivity example (first triangle):") -tri_conn = elements[1][2] -println(" Triangle 1: nodes $tri_conn") -println(" Coordinates:") -for node_id in tri_conn - coord = nodes[node_id] - println(" Node $node_id: ($(coord[1]), $(coord[2]))") -end -println() - -# ============================================================================= -# Step 3: What FEM Assembly Would Do (When Heat Problem Available) -# ============================================================================= - -println("Step 3: FEM Assembly (Conceptual - Heat problem pending refactoring)") -println("-"^80) -println() -println("When Heat problem is available, assembly would:") -println() -println("1. Loop over elements:") -println(" for element in body_elements") -println(" # Get element nodes and coordinates") -println(" X = [nodes[nid] for nid in element.connectivity]") -println() -println("2. Compute element matrices via numerical integration:") -println(" for each Gauss point ξ:") -println(" N = basis_functions(ξ) # Shape functions") -println(" dN = grad_basis(ξ) # Derivatives") -println(" J = jacobian(dN, X) # Parametric → physical") -println(" ") -println(" Kₑ += w * det(J) * α * dN' * dN # Stiffness") -println(" Mₑ += w * det(J) * N' * N # Mass") -println() -println("3. Scatter to global matrices:") -println(" dofs = get_dofs(element)") -println(" K[dofs, dofs] += Kₑ") -println(" M[dofs, dofs] += Mₑ") -println() -println("4. Apply boundary conditions (Dirichlet BC: u = 0 on left)") -println() -println("✓ Assembly concept explained") -println() - -# ============================================================================= -# Step 4: Extract Matrices for External Solvers -# ============================================================================= - -println("Step 4: Matrix Extraction for External Solvers (Issue #183)") -println("-"^80) -println() -println("When Heat assembly is working, matrices would be extracted as:") -println() -println(" K = body.assembly.K # Stiffness (SparseMatrixCSC{Float64,Int64})") -println(" M = body.assembly.M # Mass matrix (SparseMatrixCSC{Float64,Int64})") -println(" f = body.assembly.f # Force vector (Vector{Float64})") -println() -println("Where N = number of DOFs = $(length(nodes))") -println() -println("These standard Julia types integrate directly with:") -println() -println("1. DifferentialEquations.jl (Semidiscretization of PDE):") -println(" using DifferentialEquations") -println(" function heat_ode!(du, u, p, t)") -println(" K, M, f = p") -println(" du .= M \\ (-K * u .+ f)") -println(" end") -println(" u0 = [sin(π*node[1])*sin(π*node[2]) for node in values(nodes)]") -println(" prob = ODEProblem(heat_ode!, u0, (0.0, 1.0), (K, M, f))") -println(" sol = solve(prob, Tsit5())") -println() -println("2. LinearSolve.jl (Steady-state problem):") -println(" using LinearSolve") -println(" prob = LinearProblem(K, f)") -println(" sol = solve(prob, KrylovJL_GMRES())") -println() -println("3. Custom research solvers (Issue #183 use case):") -println(" using IterativeSolvers") -println(" u_steady = gmres(K, f; abstol=1e-10, maxiter=1000)") -println() -println("✓ This addresses Issue #183: matrices accessible for academic/research use") -println() - -# ============================================================================= -# Step 5: What Full Solve Would Look Like (When Available) -# ============================================================================= - -println("Step 5: Full Solver Workflow (Coming in Phase 2)") -println("-"^80) -println() -println("When Heat problem is restored, full solver workflow would be:") -println() -println(" # Create solver and add problems") -println(" solver = Solver(Linear)") -println(" push!(solver, body, bc)") -println() -println(" # Run assembly and solve") -println(" solver()") -println() -println(" # Extract solution") -println(" u = body(\"temperature\", time)") -println() -println("Current Status:") -println(" ✗ Heat problem disabled (depends on old basis system)") -println(" ✓ Architecture refactoring in progress (see llm/ARCHITECTURE.md)") -println(" ✓ Dirichlet BC problem type working") -println(" ✓ New immutable element system validated (40-130x speedup)") -println() -println("✓ Full workflow documented, implementation pending Phase 2") -println() - -# ============================================================================= -# Summary -# ============================================================================= - -println("="^80) -println("Summary: Workflow Demonstrated (Issue #183)") -println("="^80) -println() -println("This example shows the STRUCTURE of going from mesh to FEM matrices:") -println() -println("✓ Step 1: Mesh generation (programmatic or from Gmsh .msh file)") -println("✓ Step 2: Element connectivity and boundary identification") -println("✓ Step 3: Assembly concept (element matrices → global K, M)") -println("✓ Step 4: Matrix extraction for external solvers (K, M, f)") -println("✓ Step 5: Integration with DifferentialEquations.jl / LinearSolve.jl") -println() -println("Current Implementation Status:") -println(" [WORKING] Mesh generation and element connectivity") -println(" [WORKING] Dirichlet boundary condition problem type") -println(" [PENDING] Heat/Elasticity problem types (Phase 2 refactoring)") -println(" [PENDING] Full assembly and solve") -println() -println("What Works NOW:") -println(" • Load mesh from Gmsh (.msh) or programmatic generation") -println(" • Identify boundary elements") -println(" • Apply Dirichlet boundary conditions") -println() -println("What's COMING (Phase 2, ~2-4 months):") -println(" • Restore Heat/Elasticity problem types with new basis system") -println(" • Full assembly → K, M, f matrices") -println(" • Built-in solvers restored") -println(" • Performance: 40-130x faster than v0.5.1") -println() -println("References:") -println(" • Full tutorial: docs/book/gmsh_tutorial.md") -println(" • Architecture: llm/ARCHITECTURE.md") -println(" • Performance analysis: docs/blog/immutability_performance.md") -println(" • Issue #183: github.com/JuliaFEM/JuliaFEM.jl/issues/183") -println() -println("="^80) diff --git a/examples/gmsh_heat_equation/unit_square.geo b/examples/gmsh_heat_equation/unit_square.geo deleted file mode 100644 index f4c1cf2..0000000 --- a/examples/gmsh_heat_equation/unit_square.geo +++ /dev/null @@ -1,32 +0,0 @@ -// Gmsh geometry file: Unit square mesh for heat equation tutorial -// Generate with: gmsh -2 unit_square.geo -o unit_square.msh - -// Mesh element size -lc = 0.1; - -// Corner points -Point(1) = {0, 0, 0, lc}; -Point(2) = {1, 0, 0, lc}; -Point(3) = {1, 1, 0, lc}; -Point(4) = {0, 1, 0, lc}; - -// Edges -Line(1) = {1, 2}; // Bottom -Line(2) = {2, 3}; // Right -Line(3) = {3, 4}; // Top -Line(4) = {4, 1}; // Left - -// Surface -Line Loop(1) = {1, 2, 3, 4}; -Plane Surface(1) = {1}; - -// Physical groups for boundary conditions -Physical Line("bottom") = {1}; -Physical Line("right") = {2}; -Physical Line("top") = {3}; -Physical Line("left") = {4}; -Physical Surface("body") = {1}; - -// Use triangular elements -Mesh.ElementOrder = 1; // Linear elements (Tri3) -Mesh.Algorithm = 6; // Frontal-Delaunay for 2D diff --git a/examples/linear_static.jl b/examples/linear_static.jl deleted file mode 100755 index 639e372..0000000 --- a/examples/linear_static.jl +++ /dev/null @@ -1,133 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -# # JuliaFEM Linear Static Example - -# ![](linear_static/freecad.png) - -# ## Preprocessing - -using JuliaFEM - -# First we will read in the mesh. Geometry and mesh are greated with FreeCAD, -# where med format is selected for exporting. Mesh file consist also edge and -# surface mesh, which we will need to neglect later. - -datadir = Pkg.dir("JuliaFEM", "examples", "linear_static") -meshfile = joinpath(datadir, "JuliaFEMSMP18.med") -mesh = aster_read_mesh(meshfile) - -# Next we will create the model and define Elasticity. Also elements are added -# to the model. - -model = Problem(Elasticity, "OTHER", 3) -model_elements = create_elements(mesh, "OTHER") - -# Elements need material properties and they are defined next - -update!(model_elements, "youngs modulus", 208.0E3) -update!(model_elements, "poissons ratio", 0.30) -update!(model_elements, "density", 7.80E-9) -add_elements!(model, model_elements) - -# We can ignore Seg3 and Tri6 elements using `filter` with a special function -# returning true if element is something else than Seg3 or Tri6: - -function is_not_Seg3_or_Tri6(element) - return !isa(element, Union{Element{Seg3}, Element{Tri6}}) -end - -filter!(is_not_Seg3_or_Tri6, model.elements) - -# The whole idea of the JuliaFEM input is to be a normal Julia script, where the -# user can freely define any functions needed to perform the task. Here we -# define a function, which finds nodes on the given plane yz, xz or xy from the -# given height. - -function add_nodes_at_certain_plane_to_node_set!(mesh, name, vector_id, distance, - radius=6.0) - for (node, coords) in mesh.nodes - if isapprox(coords[vector_id], distance, atol=radius) - add_node_to_node_set!(mesh, name, node) - end - end - return nothing -end - -# We will find nodes from the xz-plane going through point (0,50,0) or actually -# we previously defined the radius to be 6.0, which means (0,[44,56],0). In other -# words we will select each node, which second coordinate value is between 44 -# and 56. This function will edit mesh object and add node set called `:mid_fixed` -# to it. - -add_nodes_at_certain_plane_to_node_set!(mesh, :mid_fixed, 2, 50.0) - -# We need to somehow handle the i's dot. I looked the rough coordinates of the -# dot in FreeCAD and now we can search three closest nodes to these coordinates. -# Those will be added to the same set `:mid_fixed`. - -ipoint = find_nearest_nodes(mesh, [165.0, 88.0, 10],3) -for poi in ipoint - add_node_to_node_set!(mesh, :mid_fixed, poi) -end - -# The fixed boundary conditions are defined next. - -fixed = Problem(Dirichlet, "fixed", 3, "displacement") -fixed_elements = create_nodal_elements(mesh, "mid_fixed") -add_elements!(fixed, fixed_elements) -update!(fixed_elements, "displacement 1", 0.0) -update!(fixed_elements, "displacement 2", 0.0) -update!(fixed_elements, "displacement 3", 0.0) - - -# Let's use simple acceleration load. -update!(model_elements, "displacement load 1", 1.0) - -# Finally the ´Analysis` couples everything togeter. -analysis = Analysis(Linear, model, fixed) - -# Let's write resuls to Xdmf file - -xdmf = Xdmf("model_results"; overwrite=true) -add_results_writer!(analysis, xdmf) - -# This is how the stresses are requested -push!(model.postprocess_fields, "stress") - -# Now we have all we need to run the analysis. - -run!(analysis) - -# ## Postprocessing - -# In order to look the results, we will need to close the xdmf that it is actually -# written to the file from buffer. - -close(xdmf) - -# Finally when we open the model in ParaView and set some settings we have this -# end result. - -# ![](linear_static/paraview.png) - -# ## Testing - -# First let's test that we have the output files writen to the disk - -if VERSION < v"1.0.0" - using Base.Test -else - using Test -end - -@test isfile("model_results.xmf") -@test isfile("model_results.h5") - -# Secondly let's test that we have the same maximum displacement each time. -# This is also an usefull example how to access the displacements values. - -time = 0.0 -u = analysis("displacement", time) -u_norms = Dict(i => norm(j) for (i, j) in u) -@test isapprox(maximum(values(u_norms)),2.4052929896922337) diff --git a/examples/linear_static/JuliaFEMSMP18.med b/examples/linear_static/JuliaFEMSMP18.med deleted file mode 100755 index 24d9957..0000000 Binary files a/examples/linear_static/JuliaFEMSMP18.med and /dev/null differ diff --git a/examples/linear_static/freecad.png b/examples/linear_static/freecad.png deleted file mode 100755 index 8ac923b..0000000 Binary files a/examples/linear_static/freecad.png and /dev/null differ diff --git a/examples/linear_static/paraview.png b/examples/linear_static/paraview.png deleted file mode 100755 index 9f23703..0000000 Binary files a/examples/linear_static/paraview.png and /dev/null differ