some review fixes

docs/src/literate-tutorials/gpu_qp_heat_equation.jl

@@ -5,31 +5,20 @@ using CUDA
 
 
 left = Tensor{1, 2, Float32}((0, -0)) # define the left bottom corner of the grid.
-
-right = Tensor{1, 2, Float32}((1000.0, 1000.0)) # define the right top corner of the grid.
-
-
+right = Tensor{1, 2, Float32}((1.0, 1.0)) # define the right top corner of the grid.
 grid = generate_grid(Quadrilateral, (1000, 1000), left, right)
 
 
 ip = Lagrange{RefQuadrilateral, 2}() # define the interpolation function (i.e. Bilinear lagrange)
-
-
 qr = QuadratureRule{RefQuadrilateral}(Float32, 3)
-
-
 cellvalues = CellValues(Float32, qr, ip)
 
 
-dh = DofHandler(grid)
-
-add!(dh, :u, ip)
-
-close!(dh);
+dh = DofHandler(grid)  
+add!(dh, :u, ip)  
+close!(dh); 
 
 
-dh |> get_grid
-
 # Standard assembly of the element.
 function assemble_element_std!(Ke::Matrix, fe::Vector, cellvalues::CellValues)
     n_basefuncs = getnbasefunctions(cellvalues)
@@ -55,36 +44,60 @@ function assemble_element_std!(Ke::Matrix, fe::Vector, cellvalues::CellValues)
 end
 
 
-function create_buffers(cellvalues, dh)
-    f = zeros(ndofs(dh))
-    K = allocate_matrix(dh)
-    assembler = start_assemble(K, f)
-    ## Local quantities
+function assemble_element_qpiter!(Ke::Matrix, fe::Vector, cellvalues, cell_coords::AbstractVector)
     n_basefuncs = getnbasefunctions(cellvalues)
-    Ke = zeros(n_basefuncs, n_basefuncs)
-    fe = zeros(n_basefuncs)
-    return (; f, K, assembler, Ke, fe)
+    ## Loop over quadrature points
+    for qv in Ferrite.QuadratureValuesIterator(cellvalues, cell_coords)
+        ## Get the quadrature weight
+        dΩ = getdetJdV(qv)
+        ## Loop over test shape functions
+        for i in 1:n_basefuncs
+            δu  = shape_value(qv, i)
+            ∇δu = shape_gradient(qv, i)
+            ## Add contribution to fe
+            fe[i] += δu * dΩ
+            ## Loop over trial shape functions
+            for j in 1:n_basefuncs
+                ∇u = shape_gradient(qv, j)
+                ## Add contribution to Ke
+                Ke[i, j] += (∇δu ⋅ ∇u) * dΩ
+            end
+        end
+    end
+    return Ke, fe
 end
 
 
 # Standard global assembly
+function assemble_global_qp!(cellvalues, dh::DofHandler, K, f)
+    n_basefuncs = getnbasefunctions(cellvalues)
+    Ke = zeros(eltype(K), n_basefuncs, n_basefuncs)
+    fe = zeros(eltype(f), n_basefuncs)
+    assembler = start_assemble(K, f)
+    ## Loop over all cels
+    for cell in CellIterator(dh)
+        fill!(Ke, 0)
+        fill!(fe, 0)
+        assemble_element_qpiter!(Ke, fe, cellvalues, getcoordinates(cell))
+        ## Assemble Ke and fe into K and f
+        assemble!(assembler, celldofs(cell), Ke, fe)
+    end
+    return K, f
+end
 
-
-function assemble_global!(cellvalues, dh::DofHandler, qp_iter::Val{QPiter}) where {QPiter}
-    (; f, K, assembler, Ke, fe) = create_buffers(cellvalues, dh)
+function assemble_global_std!(cellvalues, dh::DofHandler, K, f)
+    n_basefuncs = getnbasefunctions(cellvalues)
+    Ke = zeros(eltype(K), n_basefuncs, n_basefuncs)
+    fe = zeros(eltype(f), n_basefuncs)
+    assembler = start_assemble(K, f)
     ## Loop over all cels
     for cell in CellIterator(dh)
         fill!(Ke, 0)
         fill!(fe, 0)
-        if QPiter
-            ## reinit!(cellvalues, getcoordinates(cell))
-            assemble_element_qpiter!(Ke, fe, cellvalues, getcoordinates(cell))
-        else
-            ## Reinitialize cellvalues for this cell
-            reinit!(cellvalues, cell)
-            ## Compute element contribution
-            assemble_element_std!(Ke, fe, cellvalues)
-        end
+        ## Reinitialize cellvalues for this cell
+        reinit!(cellvalues, cell)
+        ## Compute element contribution
+        assemble_element_std!(Ke, fe, cellvalues)
         ## Assemble Ke and fe into K and f
         assemble!(assembler, celldofs(cell), Ke, fe)
     end
@@ -93,8 +106,6 @@ end
 
 
 ## gpu version of element assembly
-
-
 function assemble_element!(Ke, fe, cv, cell)
     n_basefuncs = getnbasefunctions(cv)
     for qv in Ferrite.QuadratureValuesIterator(cv, getcoordinates(cell))
@@ -137,14 +148,14 @@ n_basefuncs = getnbasefunctions(cellvalues)
 
 ## Allocate CPU matrix
 K = allocate_matrix(SparseMatrixCSC{Float32, Int32}, dh);
-
 #K = allocate_matrix(SparseMatrixCSC{Float64, Int64}, dh);
-f = zeros(ndofs(dh));
+#f = zeros(eltype(K), ndofs(dh));
+
 
 # Allocate GPU matrix
 ## commented to pass the test
 ## Kgpu = CUSPARSE.CuSparseMatrixCSC(K);
-## fgpu = CUDA.zeros(ndofs(dh));
+## fgpu = CUDA.zeros(Float32, ndofs(dh));
 
 n_cells = dh |> get_grid |> getncells