mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-16 16:53:19 +00:00
chore: removed files from examples
Removed files from examples directory: - examples/2d_hertz_contact.jl - examples/2d_hertz_contact/hertz_2d_full.med - examples/2d_hertz_contact/model.png - examples/2d_hertz_contact/results_displacement.png - examples/3d_frame.jl - examples/3d_frame/model.med - examples/3d_frame/model.png - examples/3d_frame/natfreq.png - examples/academic_matrix_extraction/academic_example.jl - examples/academic_matrix_extraction/README.md - examples/generate_stiffness_matrices.jl - examples/gmsh_heat_equation/gmsh_heat_equation.jl - examples/gmsh_heat_equation/QUICK_START.md - examples/gmsh_heat_equation/README.md - examples/gmsh_heat_equation/unit_square.geo - examples/linear_static.jl - examples/linear_static/freecad.png - examples/linear_static/JuliaFEMSMP18.med - examples/linear_static/paraview.png
This commit is contained in:
Jukka Aho
@@ -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
|
||||
|
||||
# 
|
||||
|
||||
# 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:
|
||||
# 
|
||||
|
||||
# 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.
|
||||
Binary file not shown.
Binary file not shown.
|
Before Width: | Height: | Size: 25 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 138 KiB |
@@ -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.
|
||||
|
||||
# 
|
||||
|
||||
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))")
|
||||
|
||||
# [](https://www.youtube.com/watch?v=GzktCqeASmo)
|
||||
|
||||
Binary file not shown.
Binary file not shown.
|
Before Width: | Height: | Size: 14 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 84 KiB |
@@ -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.
|
||||
@@ -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)
|
||||
@@ -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)
|
||||
@@ -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.
|
||||
@@ -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
|
||||
@@ -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)
|
||||
@@ -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
|
||||
@@ -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
|
||||
|
||||
# 
|
||||
|
||||
# ## 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.
|
||||
|
||||
# 
|
||||
|
||||
# ## 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)
|
||||
Binary file not shown.
Binary file not shown.
|
Before Width: | Height: | Size: 278 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 272 KiB |
Reference in New Issue
Block a user