stuff

src/Dofs/DirichletBoundaryConditions.jl

@@ -32,7 +32,7 @@ julia> t = 1.0;
 julia> update!(dbc, t)
 ```
 
-The boundary conditions can be applied to a vector with
+The boundary conditions can be applied to a vector:
 
 ```jldoctest
 julia> u = zeros(ndofs(dh))
@@ -89,6 +89,7 @@ function add!(dbcs::DirichletBoundaryConditions, field::Symbol,
     add!(dbcs, field, nodes, f, [component])
 end
 
+# Adds a boundary condition
 function add!(dbcs::DirichletBoundaryConditions, field::Symbol,
                           nodes::Union{Set{Int}, Vector{Int}}, f::Function, components::Vector{Int})
     @assert field in dbcs.dh.field_names || error("Missing: $field")
@@ -118,6 +119,7 @@ function add!(dbcs::DirichletBoundaryConditions, field::Symbol,
 
 end
 
+# Updates the DBC's to the current time `time`
 function update!(dbcs::DirichletBoundaryConditions, time::Float64 = 0.0)
     @assert dbcs.closed[]
     bc_offset = 0
@@ -144,6 +146,8 @@ function _update!(values::Vector{Float64}, f::Function, nodes::Set{Int}, field::
     end
 end
 
+# Saves the dirichlet boundary conditions to a vtkfile.
+# Values will have a 1 where bcs are active and 0 otherwise
 function vtk_point_data(vtkfile, dbcs::DirichletBoundaryConditions)
     unique_fields = []
     for dbc in dbcs.bcs
@@ -175,3 +179,45 @@ function apply!(v::Vector, bc::DirichletBoundaryConditions)
     v[bc.dofs] = bc.values
     return v
 end
+
+function apply!(K::SparseMatrixCSC, bc::DirichletBoundaryConditions)
+    @assert all(iszero, bc.values)
+    _apply(K, eltype(K)[], bc)
+end
+
+function apply!(K::SparseMatrixCSC, f::AbstractVector, bc::DirichletBoundaryConditions)
+    @assert length(f) == size(K, 1)
+    @boundscheck checkbounds(K, bc.dofs, bc.dofs)
+    @boundscheck checkbounds(f, bc.dofs)
+
+    m = meandiag(K) # Use the mean of the diagonal here to not ruin things for iterative solver
+    @inbounds for i in 1:length(bc.values)
+        d = bc.dofs[i]
+        v = bc.values[i]
+
+        if v != 0
+            for j in nzrange(K, d)
+                f[K.rowval[j]] -= v * K.nzval[j]
+            end
+        end
+    end
+    K[:, bc.dofs] = 0
+    K[bc.dofs, :] = 0
+    @inbounds for i in 1:length(bc.values)
+        d = bc.dofs[i]
+        v = bc.values[i]
+        K[d, d] = m
+        # We will only enter here with an empty f vector if we have assured that v == 0 for all dofs
+        if v != 0
+            f[d] = v * m
+        end
+    end
+end
+
+function meandiag(K::AbstractMatrix)
+    z = zero(eltype(K))
+    for i in 1:size(K,1)
+        z += K[i,i]
+    end
+    return z / size(K,1)
+end

src/Dofs/DofHandler.jl

@@ -75,7 +75,7 @@ julia> vtk_point_data(vtkfile
 
 """
 
-# TODO: Make this immutable
+# TODO: Make this immutable?
 type DofHandler{dim, N, T, M}
     dofs_nodes::Matrix{Int}
     dofs_cells::Matrix{Int} # TODO <- Is this needed or just extract from dofs_nodes?
@@ -87,7 +87,7 @@ type DofHandler{dim, N, T, M}
 end
 
 function DofHandler(m::Grid)
-    DofHandler(Matrix{Int}(), Matrix{Int}(), Symbol[], Int[], Ref(false), Int[], m)
+    DofHandler(Matrix{Int}(0, 0), Matrix{Int}(0, 0), Symbol[], Int[], Ref(false), Int[], m)
 end
 
 function show(io::IO, dh::DofHandler)
@@ -109,6 +109,7 @@ ndofs_per_cell(dh::DofHandler) = size(dh.dofs_cells, 1)
 isclosed(dh::DofHandler) = dh.closed[]
 dofs_node(dh::DofHandler, i::Int) = dh.dof_nodes[:, i]
 
+# Stores the dofs for the cell with number `i` into the vector `global_dofs`
 function celldofs!(global_dofs::Vector{Int}, dh::DofHandler, i::Int)
     @assert isclosed(dh)
     @assert length(global_dofs) == ndofs_per_cell(dh)
@@ -118,13 +119,15 @@ function celldofs!(global_dofs::Vector{Int}, dh::DofHandler, i::Int)
     return global_dofs
 end
 
+# Add a collection of fields
 function push!(dh::DofHandler, names::Vector{Symbol}, dims)
     @assert length(names) == length(dims)
     for i in 1:length(names)
         push!(dh, names[i], dims[i])
     end
 end
 
+# Add a field to the dofhandler ex `push!(dh, :u, 3)`
 function push!(dh::DofHandler, name::Symbol, dim::Int)
     @assert !isclosed(dh)
     if name in dh.field_names
@@ -173,6 +176,7 @@ end
 
 getnvertices{dim, N, M}(::Type{JuAFEM.Cell{dim, N, M}}) = N
 
+# Computes the "edof"-matrix
 function add_element_dofs!(dh::DofHandler)
     n_elements = getncells(dh.grid)
     n_vertices = getnvertices(getcelltype(dh.grid))
@@ -192,7 +196,8 @@ function add_element_dofs!(dh::DofHandler)
     dh.dofs_cells = reshape(element_dofs, (ndofs * n_vertices, n_elements))
 end
 
-
+# Creates a sparsity pattern from the dofs in a dofhandler.
+# Returns a sparse matrix with the correct pattern.
 function create_sparsity_pattern(dh::DofHandler)
     grid = dh.grid
     n = ndofs_per_cell(dh)

src/Export/VTK.jl

@@ -94,12 +94,12 @@ getVTKtype(::Lagrange{3, RefTetrahedron, 1}) = VTKCellTypes.VTK_TETRA
 getVTKtype(::Type{Line}) = VTKCellTypes.VTK_LINE
 getVTKtype(::Type{QuadraticLine}) = VTKCellTypes.VTK_QUADRATIC_EDGE
 getVTKtype(::Type{Quadrilateral}) = VTKCellTypes.VTK_QUAD
-# getVTKtype(::Type{QuadraticQuadrilateral}) = VTKCellTypes.VTK_BIQUADRATIC_QUAD
 getVTKtype(::Type{Triangle}) = VTKCellTypes.VTK_TRIANGLE
 getVTKtype(::Type{QuadraticTriangle}) = VTKCellTypes.VTK_QUADRATIC_TRIANGLE
-getVTKtype(::Type{QuadraticQuadrilateral}) = VTKCellTypes.VTK_QUADRATIC_QUAD
+getVTKtype(::Type{QuadraticQuadrilateral}) = VTKCellTypes.VTK_BIQUADRATIC_QUAD
 getVTKtype(::Type{Hexahedron}) = VTKCellTypes.VTK_HEXAHEDRON
 getVTKtype(::Type{Tetrahedron}) = VTKCellTypes.VTK_TETRA
+getVTKtype(::Type{Cell{2,8,4}}) = VTKCellTypes.VTK_QUADRATIC_QUAD
 
 
 function vtk_grid{dim, N, T}(filename::AbstractString, grid::Grid{dim, N, T})

src/Grid/grid.jl

@@ -56,11 +56,9 @@ nfaces{dim, N, M}(::Type{Cell{dim, N, M}}) = M
 @compat const Line = Cell{1, 2, 2}
 @compat const QuadraticLine = Cell{1, 3, 2}
 
-
 @compat const Triangle = Cell{2, 3, 3}
 @compat const QuadraticTriangle = Cell{2, 6, 3}
 
-
 @compat const Quadrilateral = Cell{2, 4, 4}
 @compat const QuadraticQuadrilateral = Cell{2, 9, 4}
 

src/Grid/grid_generators.jl

@@ -210,7 +210,7 @@ function generate_grid{T}(::Type{Hexahedron}, nel::NTuple{3, Int}, left::Vec{3,
                               [(cl, 5) for cl in cell_array[1,:,:][:]];
                               [(cl, 6) for cl in cell_array[:,:,end][:]]]
 
-    _apply_onboundary!(QuadraticQuadrilateral, cells, boundary)
+    _apply_onboundary!(Hexahedron, cells, boundary)
 
     # Cell face sets
     offset = 0

test/test_grid.jl

@@ -17,8 +17,8 @@
     end
 
     # VTK types of different Cells.
-    for (dim, n) in ((1,2), (1,3), (2,4), (2,9), (3,8), (2,3), (2,6), (3,4), (2,8))
-        @test getVTKtype(Cell{dim, n}).nodes == n
+    for (dim, n, m) in ((1,2,2), (1,3,2), (2,4,4), (2,9,4), (3,8,6), (2,3,3), (2,6,3), (3,4,4), (2,8,4))
+        @test getVTKtype(Cell{dim, n, m}).nodes == n
     end
 
 end # of testset