Re-enable iga tutorial (#31)

* Generalize to use Abstract[Cell/Facet]Values in codebase
* Apply renaming following Ferrite-FEM/Ferrite.jl#916
* Change default states to Vector{Nothing} instead of Nothing

.github/workflows/CI.yml

@@ -79,7 +79,7 @@ jobs:
                 PackageSpec(;url="https://github.com/KnutAM/MaterialModelsBase.jl.git"),
                 PackageSpec(;url="https://github.com/KnutAM/Newton.jl.git"),
                 PackageSpec(;url="https://github.com/KnutAM/MechanicalMaterialModels.jl.git"),
-                #PackageSpec(;url="https://github.com/lijas/IGA.jl.git"),
+                PackageSpec(;url="https://github.com/lijas/IGA.jl.git"),
                 PackageSpec(path=pwd())
                 ]
             Pkg.develop(packages)

docs/Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.10.5"
 manifest_format = "2.0"
-project_hash = "1cc08f4c6506ca353b7656297763418ecbb18cd4"
+project_hash = "3cccd2fa35dbbf97d5cd332812d31f9109f498c6"
 
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
@@ -411,11 +411,11 @@ version = "3.3.10+0"
 
 [[deps.Ferrite]]
 deps = ["EnumX", "ForwardDiff", "LinearAlgebra", "NearestNeighbors", "OrderedCollections", "Preferences", "Reexport", "SparseArrays", "StaticArrays", "Tensors", "WriteVTK"]
-git-tree-sha1 = "886da974305982fd13637abc00c1e85c2bcc09cd"
+git-tree-sha1 = "a0442408f40f96d6477d1de0035fbb92507bcf63"
 repo-rev = "master"
 repo-url = "https://github.com/Ferrite-FEM/Ferrite.jl.git"
 uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
-version = "0.3.14"
+version = "1.0.0"
 
     [deps.Ferrite.extensions]
     FerriteBlockArrays = "BlockArrays"
@@ -598,6 +598,14 @@ git-tree-sha1 = "f218fe3736ddf977e0e772bc9a586b2383da2685"
 uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 version = "0.3.23"
 
+[[deps.IGA]]
+deps = ["Ferrite", "InteractiveUtils", "LinearAlgebra", "OrderedCollections", "Reexport", "SparseArrays", "StaticArrays", "Tensors", "WriteVTK"]
+git-tree-sha1 = "8b625d3a5be1e2156997604b978c5c151a48a2fa"
+repo-rev = "master"
+repo-url = "https://github.com/lijas/IGA.jl.git"
+uuid = "e7b8d123-e02a-40ba-ad18-d3943ed54f1c"
+version = "0.2.6"
+
 [[deps.IOCapture]]
 deps = ["Logging", "Random"]
 git-tree-sha1 = "b6d6bfdd7ce25b0f9b2f6b3dd56b2673a66c8770"

docs/Project.toml

@@ -6,6 +6,7 @@ Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
 FerriteAssembly = "fd21fc07-c509-4fe1-9468-19963fd5935d"
 FerriteMeshParser = "0f8c756f-80dd-4a75-85c6-b0a5ab9d4620"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+IGA = "e7b8d123-e02a-40ba-ad18-d3943ed54f1c"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 MaterialModelsBase = "af893363-701d-44dc-8b1e-d9a2c129bfc9"
 MechanicalMaterialModels = "b3282f9b-607f-4337-ab95-e5488ab5652c"

docs/make.jl

@@ -13,7 +13,7 @@ tutorials = [
     "viscoelasticity.jl", 
     "incompressible_elasticity.jl", 
     "mixed_materials.jl", 
-    #"iga.jl", # Doesn't currently support Ferrite master/v1.0
+    "iga.jl",
     ]
 generated_tutorials = build_examples(tutorials; type="tutorials")
 howto = [

docs/src/Convenience/MaterialModelsBase.md

@@ -8,8 +8,8 @@ and [`allocate_cell_cache`](@ref FerriteAssembly.allocate_cell_cache)
 are implemented in `FerriteAssembly.jl`.
 
 ```@docs
-FerriteAssembly.element_routine!(Ke, re, state::Vector{<:MMB.AbstractMaterialState}, ae, material::MMB.AbstractMaterial, cellvalues::CellValues, buffer)
-FerriteAssembly.element_residual!(re, state::Vector{<:MMB.AbstractMaterialState}, ae, material::MMB.AbstractMaterial, cellvalues::CellValues, buffer)
-FerriteAssembly.create_cell_state(m::MMB.AbstractMaterial, cv::CellValues, args...)
+FerriteAssembly.element_routine!(Ke, re, state::Vector{<:MMB.AbstractMaterialState}, ae, material::MMB.AbstractMaterial, cellvalues::AbstractCellValues, buffer)
+FerriteAssembly.element_residual!(re, state::Vector{<:MMB.AbstractMaterialState}, ae, material::MMB.AbstractMaterial, cellvalues::AbstractCellValues, buffer)
+FerriteAssembly.create_cell_state(m::MMB.AbstractMaterial, cv::AbstractCellValues, args...)
 FerriteAssembly.allocate_cell_cache(m::MMB.AbstractMaterial, ::Any)
 ```

docs/src/literate_howto/robin_bc.jl

@@ -24,7 +24,7 @@ struct RobinBC{T}
 end
 
 function FerriteAssembly.facet_routine!(
-        Ke, re, ae, rbc::RobinBC, fv::FacetValues, facetbuffer
+        Ke, re, ae, rbc::RobinBC, fv::AbstractFacetValues, facetbuffer
         )
     for q_point in 1:getnquadpoints(fv)
         dΓ = getdetJdV(fv, q_point)

docs/src/literate_tutorials/iga.jl

@@ -1,4 +1,4 @@
-# # Isogeometric analysis with `IGA.jl`
+# # Using `IGA.jl`
 # This tutorial shows how to integrate FerriteAssembly with the 
 # isogeometric analysis toolbox IGA.jl, directly based on `IGA.jl`'s 
 # [Infinite plate with hole](https://lijas.github.io/IGA.jl/dev/examples/plate_with_hole/)
@@ -15,11 +15,13 @@ using Ferrite, IGA, LinearAlgebra, FerriteAssembly
 import FerriteAssembly.ExampleElements: ElasticPlaneStrain
 
 # ## Setup
-# We begin by generating the mesh by using `IGA.jl`'s built-in mesh generators 
-# In this example, we will generate the patch called "plate with hole". 
-# Note, currently this function can only generate the patch with second order basefunctions. 
-function create_mesh(;nels=(20,10))
-    nurbsmesh = generate_nurbs_patch(:plate_with_hole, nels) 
+# To clarify the differences, we split the setup into `IGA.jl`, `Ferrite.jl`,
+# and `FerriteAssembly.jl` specific parts.
+# ### `IGA.jl` setup 
+# We begin by generating the mesh by using `IGA.jl`'s built-in mesh generators, 
+# specifically a "plate with hole". 
+function create_mesh(; nels = (20,10), order)
+    nurbsmesh = generate_nurbs_patch(:plate_with_hole, nels, order)
     grid = BezierGrid(nurbsmesh)
 
     addfacetset!(grid, "left", (x) -> x[1] ≈ -4.0)
@@ -28,17 +30,19 @@ function create_mesh(;nels=(20,10))
 
     return grid 
 end
-grid = create_mesh();
+order = 2
+grid = create_mesh(; order);
 
-# Create the cellvalues storing the shape function values. 
-orders=(2,2)
-ip = BernsteinBasis{2,orders}()
-qr_cell = QuadratureRule{2,RefCube}(4)
-qr_facet = QuadratureRule{1,RefCube}(3)
+# Create the `IGA.jl` cell and facet values for storing the 
+# the `IGA.jl` shape function values and gradients. 
+ip = IGAInterpolation{RefQuadrilateral, order}()
+qr_cell = QuadratureRule{RefQuadrilateral}(4)
+qr_face = FacetQuadratureRule{RefQuadrilateral}(3)
 
-cv = BezierCellValues( CellValues(qr_cell, ip^2) );
-fv = BezierFacetValues( FacetValues(qr_facet, ip^2) );
+cv = BezierCellValues(qr_cell, ip^2);
+fv = BezierFacetValues(qr_face, ip^2);
 
+# ### `Ferrite.jl` setup
 # Distribute dofs as normal
 dh = DofHandler(grid)
 add!(dh, :u, ip^2)
@@ -49,73 +53,67 @@ a = zeros(ndofs(dh))
 r = zeros(ndofs(dh))
 K = allocate_matrix(dh);
 
-# Starting with Dirichlet conditions: 
-# 1) Bottom facet should only be able to move in x-direction
+# Before adding boundary conditions, starting with Dirichlet: 
+# 1) Bottom face should only be able to move in x-direction
 # 2) Right boundary should only be able to move in y-direction
 ch = ConstraintHandler(dh);
 add!(ch, Dirichlet(:u, getfacetset(grid, "bot"), Returns(0.0), 2))
 add!(ch, Dirichlet(:u, getfacetset(grid, "right"), Returns(0.0), 1))
 close!(ch)
 update!(ch, 0.0);
 
-# Then Neumann conditions:
-# Apply outwards traction on the left surface,
-# taking the negative value since r = fint - fext.
+# ### `FerriteAssembly.jl` setup 
+# Neumann boundary conditions are added using `FerriteAssembly`'s 
+# `LoadHandler`. We apply outwards traction on the left surface,
+# and take the negative value since r = fint - fext.
 traction = Vec((-10.0, 0.0))
 lh = LoadHandler(dh)
 add!(lh, Neumann(:u, fv, getfacetset(grid, "left"), Returns(-traction)));
 
-# FerriteAssembly setup 
 material = ElasticPlaneStrain(;E=100.0, ν=0.3)
-domain = DomainSpec(FerriteAssembly.SubDofHandler(dh, dh.fieldhandlers[1]), material, cv)
+domain = DomainSpec(dh, material, cv)
 buffer = setup_domainbuffer(domain);
 
-# ## Assemble 
+# ## Assembly and solution
+# We first assemble the equation system,
 assembler = start_assemble(K, r)
 apply!(a, ch)
 work!(assembler, buffer; a=a)
 apply!(r, lh, 0.0);
 
-# ## Solve
+# before solving it,
 apply_zero!(K, r, ch)
 a .-= K\r
 apply!(a, ch);
 
 # ## Postprocessing
-# First, the stresses in each integration point are calculated by using the Integrator.  
-struct CellStress{TT}
-    s::Vector{Vector{TT}}
-end
-function CellStress(grid::Ferrite.AbstractGrid)
-    Tb = SymmetricTensor{2,Ferrite.getspatialdim(grid)}
-    TT = Tb{Float64,Tensors.n_components(Tb)}
-    return CellStress([TT[] for _ in 1:getncells(grid)])
-end
+# We use the `QuadPointEvaluator` to calculate stresses in each integration point,
+# The `QuadPointEvaluator` requires a function with the following input arguments
+function calculate_stress(m, u, ∇u, qp_state)
+    ϵ = symmetric(∇u)
+    return m.C ⊡ ϵ
+end;
 
-function FerriteAssembly.integrate_cell!(stress::CellStress, state, ae, material, cv, cb)
-    σ = stress.s[cellid(cb)]
-    for q_point in 1:getnquadpoints(cv)
-        ϵ = function_symmetric_gradient(cv, q_point, ae)
-        push!(σ, material.C⊡ϵ)
-    end
-end
-cellstresses = CellStress(grid)
-integrator = Integrator(cellstresses)
-work!(integrator, buffer; a=a);
+# We can then use this to calculate the stress tensor in each integration point
+qe = QuadPointEvaluator{SymmetricTensor{2,2,Float64,3}}(buffer, calculate_stress);
+work!(qe, buffer; a=a);
 
 # Now we want to export the results to VTK. So we project the stresses at 
 # the quadrature points to the nodes using the L2Projector from Ferrite. 
-projector = L2Projector(ip, grid)
-σ_nodes = IGA.igaproject(projector, cellstresses.s, qr_cell; project_to_nodes=true);
+# Currently, however, IGA doesn't support L2 projection. 
+# ```julia
+# # projector = L2Projector(ip, grid)
+# # σ_nodes = IGA.igaproject(projector, qe.data, qr_cell; project_to_nodes=true);
+# ```
 
 # Output results to VTK
-VTKGridFile("plate_with_hole.vtu", grid) do vtk
+IGA.VTKIGAFile("plate_with_hole.vtu", grid) do vtk
     write_solution(vtk, dh, a)
-    write_node_data(vtk, σ_nodes, "sigma", grid)
 end
 
-using Test                                    #src
-@test sum(norm, σ_nodes) ≈ 3087.2447327126742 #src
+using Test                                      #src
+# @test sum(norm, σ_nodes) ≈ 3087.2447327126742 #src
+@test norm(norm.(qe.data)) ≈ 679.3207411544098  #src
 
 #md # ## [Plain program](@id iga_plain_program)
 #md #

docs/src/literate_tutorials/viscoelasticity.jl

@@ -53,14 +53,14 @@ end;
 # We then define how to the initial state variables should look like, which also defines the structure 
 # of the state variables. In this case, we will just have states being a single tensor (viscous strain)
 # for each integration point
-function FerriteAssembly.create_cell_state(::ZenerMaterial, cv::CellValues, args...)
+function FerriteAssembly.create_cell_state(::ZenerMaterial, cv::AbstractCellValues, args...)
     ϵ_template = shape_symmetric_gradient(cv, 1, 1) # ::SymmetricTensor
     return [zero(ϵ_template) for _ in 1:getnquadpoints(cv)]
 end;
 
 # Following this, we define the `element_residual!` function (we will use automatic differentiation 
 # to calculate the element stiffness).
-function FerriteAssembly.element_residual!(re, state, ae, m::ZenerMaterial, cv::CellValues, buffer)
+function FerriteAssembly.element_residual!(re, state, ae, m::ZenerMaterial, cv::AbstractCellValues, buffer)
     Δt = FerriteAssembly.get_time_increment(buffer)
     old_ϵvs = FerriteAssembly.get_old_state(buffer)
     for q_point in 1:getnquadpoints(cv)

src/ExampleElements/LinearElasticity.jl

@@ -43,32 +43,32 @@ LinearElastic(::Val{:planestrain}; kwargs...) = LinearElastic(Val(2); kwargs...)
 # Type-unstable switch could be made for convenience
 # LinearElastic(type::Symbol; kwargs...) = LinearElastic(Val(type); kwargs...)
 
-function FerriteAssembly.element_routine!(Ke, re, state, ae, material::LinearElastic, cv::CellValues, buffer)
+function FerriteAssembly.element_routine!(Ke, re, state, ae, material::LinearElastic, cv::AbstractCellValues, buffer)
     for q_point in 1:getnquadpoints(cv)
         dΩ = getdetJdV(cv, q_point)
         ϵ = function_symmetric_gradient(cv, q_point, ae)
         σ = material.C ⊡ ϵ
         ## Assemble residual contributions
         for i in 1:getnbasefunctions(cv)
-            ∇δNu = shape_symmetric_gradient(cv, q_point, i)
+            ∇δNu = symmetric(shape_gradient(cv, q_point, i))
             re[i] += (∇δNu ⊡ σ )*dΩ
             ∇δNu_C = ∇δNu ⊡  material.C
             for j in 1:getnbasefunctions(cv)
-                ∇Nu = shape_symmetric_gradient(cv, q_point, j)
+                ∇Nu = symmetric(shape_gradient(cv, q_point, j))
                 Ke[j,i] += (∇δNu_C ⊡ ∇Nu)*dΩ # Since Ke is symmetric, we calculate Ke' to index faster
             end
         end
     end
 end
 
-function FerriteAssembly.element_residual!(re, state, ae, material::LinearElastic, cv::CellValues, buffer)
+function FerriteAssembly.element_residual!(re, state, ae, material::LinearElastic, cv::AbstractCellValues, buffer)
     for q_point in 1:getnquadpoints(cv)
         dΩ = getdetJdV(cv, q_point)
         ϵ = function_symmetric_gradient(cv, q_point, ae)
         σ = material.C ⊡ ϵ
         ## Assemble residual contributions
         for i in 1:getnbasefunctions(cv)
-            ∇δNu = shape_symmetric_gradient(cv, q_point, i)
+            ∇δNu = symmetric(shape_gradient(cv, q_point, i))
             re[i] += (∇δNu ⊡ σ )*dΩ
         end
     end

src/LoadHandler/BodyLoad.jl

@@ -43,7 +43,7 @@ struct BodyLoadMaterial{FUN}
     dr::UnitRange{Int}
 end
 
-function element_residual!(fe::Vector, ::Any, ::Vector, m::BodyLoadMaterial, cv::CellValues, cellbuffer)
+function element_residual!(fe::Vector, ::Any, ::Vector, m::BodyLoadMaterial, cv::AbstractCellValues, cellbuffer)
     checkbounds(fe, m.dr)
     t = get_time_increment(cellbuffer) # Abuse...
     for q_point in 1:getnquadpoints(cv)

src/LoadHandler/Neumann.jl

@@ -40,7 +40,7 @@ struct NeumannMaterial{FUN}
     f::FUN
     dr::UnitRange{Int}
 end
-function facet_residual!(fe::Vector, ::Vector, m::NeumannMaterial, fv::FacetValues, facetbuffer)
+function facet_residual!(fe::Vector, ::Vector, m::NeumannMaterial, fv::AbstractFacetValues, facetbuffer)
     checkbounds(fe, m.dr)
     t = get_time_increment(facetbuffer) # Abuse...
     for q_point in 1:getnquadpoints(fv)

src/LoadHandler/defaultvalues.jl

@@ -1,7 +1,7 @@
 # Create FacetValues or CellValues automatically with minimal input 
 
 """
-    autogenerate_facetvalues(fv::FacetValues, args...)
+    autogenerate_facetvalues(fv::AbstractFacetValues, args...)
 
 Just return the provided FacetValues 
 
@@ -10,13 +10,13 @@ Just return the provided FacetValues
 Using quadrature rule, `fqr = FacetQuadratureRule{RefShape}(order)`, 
 create `FacetValues(fqr, ip_fun, ip_geo)`
 """
-autogenerate_facetvalues(fv::FacetValues, args...) = fv
+autogenerate_facetvalues(fv::AbstractFacetValues, args...) = fv
 function autogenerate_facetvalues(order::Int, ip::Interpolation{RefShape}, ip_geo::Interpolation{RefShape}) where RefShape
     return FacetValues(FacetQuadratureRule{RefShape}(order), ip, ip_geo)
 end
 
 """
-    autogenerate_cellvalues(cv::CellValues, args...)
+    autogenerate_cellvalues(cv::AbstractCellValues, args...)
 
 Just return the provided CellValues 
 
@@ -25,7 +25,7 @@ Just return the provided CellValues
 Using quadrature rule, `qr = QuadratureRule{RefShape}(order)`, 
 return `CellValues(qr, ip_fun, ip_geo)`
 """
-autogenerate_cellvalues(cv::CellValues, args...) = cv
+autogenerate_cellvalues(cv::AbstractCellValues, args...) = cv
 function autogenerate_cellvalues(order::Int, ip::Interpolation{RefShape}, ip_geo::Interpolation{RefShape}) where RefShape
     return CellValues(QuadratureRule{RefShape}(order), ip, ip_geo)
 end

src/Utils/MaterialModelsBase.jl

@@ -3,7 +3,7 @@ import MaterialModelsBase as MMB
 """
     FerriteAssembly.element_routine!(
         Ke, re, state::Vector{<:MMB.AbstractMaterialState}, ae, 
-        m::MMB.AbstractMaterial, cv::CellValues, buffer)
+        m::MMB.AbstractMaterial, cv::AbstractCellValues, buffer)
 
 Solve the weak form 
 ```math
@@ -17,7 +17,7 @@ Note that `create_cell_state` is already implemented for `<:AbstractMaterial`.
 """
 function FerriteAssembly.element_routine!(
     Ke, re, state::Vector{<:MMB.AbstractMaterialState},
-    ae, material::MMB.AbstractMaterial, cellvalues::CellValues, buffer)
+    ae, material::MMB.AbstractMaterial, cellvalues::AbstractCellValues, buffer)
     cache = FerriteAssembly.get_user_cache(buffer)
     Δt = FerriteAssembly.get_time_increment(buffer)
     state_old = FerriteAssembly.get_old_state(buffer)
@@ -43,14 +43,14 @@ end
 """
     FerriteAssembly.element_residual!(
         re, state::Vector{<:MMB.AbstractMaterialState}, ae, 
-        m::MMB.AbstractMaterial, cv::CellValues, buffer)
+        m::MMB.AbstractMaterial, cv::AbstractCellValues, buffer)
 
 The `element_residual!` implementation corresponding to the `element_routine!` implementation
 for a `MaterialModelsBase.AbstractMaterial`
 """
 function FerriteAssembly.element_residual!(
     re, state::Vector{<:MMB.AbstractMaterialState},
-    ae, material::MMB.AbstractMaterial, cellvalues::CellValues, buffer)
+    ae, material::MMB.AbstractMaterial, cellvalues::AbstractCellValues, buffer)
     cache = FerriteAssembly.get_user_cache(buffer)
     Δt = FerriteAssembly.get_time_increment(buffer)
     state_old = FerriteAssembly.get_old_state(buffer)
@@ -68,12 +68,12 @@ function FerriteAssembly.element_residual!(
 end
 
 """
-    FerriteAssembly.create_cell_state(m::MMB.AbstractMaterial, cv::CellValues, args...)
+    FerriteAssembly.create_cell_state(m::MMB.AbstractMaterial, cv::AbstractCellValues, args...)
 
 Create a `Vector{<:MMM.AbstractMaterialState}` where each element is the output from 
 `MMB.initial_material_state(m)` and the length is the number of quadrature points in `cv`.
 """
-function create_cell_state(m::MMB.AbstractMaterial, cv::CellValues, args...)
+function create_cell_state(m::MMB.AbstractMaterial, cv::AbstractCellValues, args...)
     return [MMB.initial_material_state(m) for _ in 1:getnquadpoints(cv)]
 end