From e1540b609bd3676876c500c485ea215c4038a615 Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Fri, 5 Feb 2016 11:32:09 +0200
Subject: [PATCH] matrices for dual formulation for elements

---
 src/dirichlet.jl | 39 +++++++++++++--------------------------
 src/elements.jl  | 36 +++++++++++++++++++++++++++---------
 2 files changed, 40 insertions(+), 35 deletions(-)

diff --git a/src/dirichlet.jl b/src/dirichlet.jl
index e916bed..06c84b1 100644
--- a/src/dirichlet.jl
+++ b/src/dirichlet.jl
@@ -23,27 +23,19 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele
     field_name = get_parent_field_name(problem)
     gdofs = get_gdofs(element, field_dim)
 
-    # calculate bi-orthogonal basis transformation matrix Ae
-    nnodes = size(element, 2)
-    De = zeros(nnodes, nnodes)
-    Me = zeros(nnodes, nnodes)
-    for ip in get_integration_points(element, Val{2})
-        w = ip.weight
-        J = get_jacobian(element, ip, time)
-        JT = transpose(J)
-        if size(JT, 2) == 1  # plane problem
-            w *= norm(JT)
-        else
-            w *= norm(cross(JT[:,1], JT[:,2]))
-        end
-        N = element(ip, time)
-        De += w*diagm(vec(N))
-        Me += w*N'*N
-    end
-    Ae = De*inv(Me)
+    De, Me, Ae = get_dualbasis(element, time)
 
-    # do the actual integration
-    for ip in get_integration_points(element, Val{2})
+    # left hand side
+    for i=1:field_dim
+        ldofs = gdofs[i:field_dim:end]
+        if haskey(element, field_name*" $i") || haskey(element, field_name)
+            add!(assembly.C1, ldofs, ldofs, De)
+            add!(assembly.C2, ldofs, ldofs, De)
+        end
+    end
+
+    # right hand side
+    for ip in get_integration_points(element, Val{3})
         w = ip.weight
         J = get_jacobian(element, ip, time)
         JT = transpose(J)
@@ -54,15 +46,11 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele
         end
         N = element(ip, time)
         Phi = (Ae*N')'
-        A = w*Phi'*N
-        A[abs(A) .< 1.0e-12] = 0
 
         if haskey(element, field_name)
             for i=1:field_dim
                 g = element(field_name, ip, time)
                 ldofs = gdofs[i:field_dim:end]
-                add!(assembly.C1, ldofs, ldofs, A)
-                add!(assembly.C2, ldofs, ldofs, A)
                 add!(assembly.g, ldofs, w*g*Phi')
             end
         else
@@ -70,13 +58,12 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele
                 ldofs = gdofs[i:field_dim:end]
                 if haskey(element, field_name*" $i")
                     g = element(field_name*" $i", ip, time)
-                    add!(assembly.C1, ldofs, ldofs, A)
-                    add!(assembly.C2, ldofs, ldofs, A)
                     add!(assembly.g, ldofs, w*g*Phi')
                 end
             end
         end
     end
+
 end
 
 #=
diff --git a/src/elements.jl b/src/elements.jl
index c69e79d..c7e2ec7 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -6,7 +6,10 @@ abstract AbstractElement
 type Element{E}
     connectivity :: Vector{Int}
     fields :: Dict{ASCIIString, Field}
-    dualbasis :: Matrix{Float64} # coefficients to construct dual basis
+    # matrices to construct dual basis
+    D :: Matrix{Float64}
+    M :: Matrix{Float64}
+    A :: Matrix{Float64}
 end
 
 function Base.size{E}(::Element{E})
@@ -18,8 +21,7 @@ function Base.size{E}(::Element{E}, i::Int64)
 end
 
 function convert{E}(::Type{Element{E}}, connectivity::Vector{Int})
-#   return Element{E}(connectivity, get_integration_points(E), Dict())
-    return Element{E}(connectivity, Dict(), Matrix())
+    return Element{E}(connectivity, Dict(), Matrix(), Matrix(), Matrix())
 end
 
 function get_integration_points{E}(element::Element{E}, args...)
@@ -184,22 +186,38 @@ function call{E}(element::Element{E}, xi::VecOrIP, time::Float64=0.0)
     return get_basis(element, xi)
 end
 
-function call(element::Element, xi::VecOrIP, time::Real, ::Type{Val{:dualbasis}})
+""" Return dual basis transformation matrix Ae. """
+function get_dualbasis(element::Element, time::Real)
     if length(element.dualbasis) == 0
         nnodes = size(element, 2)
-        De = zeros(nnodes, nnodes)
-        Me = zeros(nnodes, nnodes)
+        D = zeros(nnodes, nnodes)
+        M = zeros(nnodes, nnodes)
         for ip in get_integration_points(element, Val{3})
+            w = ip.weight
             J = get_jacobian(element, ip, time)
-            w = ip.weight*norm(J)
+            JT = transpose(J)
+            if size(JT, 2) == 1  # plane problem
+                # || ∂X/∂ξ ||
+                w *= norm(JT)
+            else
+                # || ∂X/∂ξ₁ × ∂X/∂ξ₂ ||
+                w *= norm(cross(JT[:,1], JT[:,2]))
+            end
             N = element(ip, time)
             De += w*diagm(vec(N))
             Me += w*N'*N
         end
-        element.dualbasis = De*inv(Me)
+        element.D = D
+        element.M = M
+        element.A = De*inv(Me)
     end
+    return element.D, element.M, element.A
+end
+
+function call(element::Element, xi::VecOrIP, time::Real, ::Type{Val{:dualbasis}})
+    De, Me, Ae = get_dualbasis(element, time)
     N = get_basis(element, xi)
-    Phi = element.dualbasis*N'
+    Phi = Ae*N'
     return Phi'
 end