From 2b326afcc792256ee52c668661ae43fa4a8e875d Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Thu, 31 Dec 2015 07:40:38 +0200
Subject: [PATCH] calculate nodal field in least squares sense

---
 src/problems.jl | 44 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 44 insertions(+)

diff --git a/src/problems.jl b/src/problems.jl
index 4d0e2bd..6990f77 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -39,3 +39,47 @@ function Base.push!(problem::AllProblems, element::Element)
     push!(problem.elements, element)
 end
 
+# TODO: better place for utility functions?
+
+""" Calculate "nodal" vector from set of elements.
+
+For example element 1 with dofs [1, 2, 3, 4] has [1, 1, 1, 1] and
+element 2 with dofs [3, 4, 5, 6] has [2, 2, 2, 2] the result will
+be sparse matrix with values [1, 1, 3, 3, 2, 2].
+
+Parameters
+----------
+field_name
+    name of field, e.g. "geometry"
+field_dim
+    degrees of freedom / node
+elements
+    elements used to calculate vector
+time
+"""
+function calculate_nodal_vector{T}(field_name::ASCIIString, field_dim::Int, elements::Vector{Element{T}}, time::Real)
+    A = SparseMatrixCOO()
+    b = SparseMatrixCOO()
+    for element in elements
+        haskey(element, field_name) || continue
+        gdofs = get_gdofs(element, field_dim)
+#       info("gdofs = $gdofs")
+        for ip in get_integration_points(element, Val{2})
+            J = get_jacobian(element, ip, time)
+            w = ip.weight*norm(J)
+            f = element(field_name, ip, time)
+            N = element(ip, time)
+            for i=1:field_dim
+                ldofs = gdofs[i:field_dim:end]
+                add!(A, ldofs, ldofs, w*kron(N', N))
+            end
+            add!(b, gdofs, w*f*N)
+        end
+    end
+    A = sparse(A)
+    dim = size(A, 1)
+    b = sparse(b, dim, 1)
+    x = A \ full(b)
+    return x
+end
+