From 7bd681d7f5f63d489b216a3a663505aee5aeed38 Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Thu, 7 Jul 2016 18:02:47 +0300
Subject: [PATCH] dirichlet boundary condition as nodal collocation by default

---
 src/problems_dirichlet.jl | 81 +++++++++++++++++++++++++++++++++++----
 1 file changed, 74 insertions(+), 7 deletions(-)

diff --git a/src/problems_dirichlet.jl b/src/problems_dirichlet.jl
index 83cafe5..6dda4a4 100644
--- a/src/problems_dirichlet.jl
+++ b/src/problems_dirichlet.jl
@@ -1,18 +1,18 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-""" Here formulation is :total or :incremental meaning that we either give
-constraint for total quantity u or it's increment Δu. For elasticity we are
-using incremental formulation.
+"""
+Problem u(X) = u₀ in Γ(d)
 """
 type Dirichlet <: BoundaryProblem
     formulation :: Symbol
     variational :: Bool
     dual_basis :: Bool
+    order :: Int
 end
 
 function Dirichlet()
-    Dirichlet(:incremental, true, false)
+    Dirichlet(:incremental, false, false, 1)
 end
 
 function get_unknown_field_name(::Type{Dirichlet})
@@ -23,20 +23,87 @@ function get_formulation_type(problem::Problem{Dirichlet})
     return problem.properties.formulation
 end
 
-function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Element, time)
+function assemble!(problem::Problem{Dirichlet}, time::Float64=0.0;
+                   auto_initialize=true)
+    # FIXME: boilerplate
+    if !isempty(problem.assembly)
+        warn("Assemble problem $(problem.name): problem.assembly is not empty and assembling, are you sure you know what are you doing?")
+    end
+    if isempty(problem.elements)
+        warn("Assemble problem $(problem.name): problem.elements is empty, no elements in problem?")
+    else
+        first_element = first(problem.elements)
+        unknown_field_name = get_unknown_field_name(problem)
+        if !haskey(first_element, unknown_field_name)
+            warn("Assemble problem $(problem.name): seems that problem is uninitialized.")
+            if auto_initialize
+                info("Initializing problem $(problem.name) at time $time automatically.")
+                initialize!(problem, time)
+            end
+        end
+    end
+
+    if method_exists(assemble_prehook!, Tuple{typeof(problem), Float64})
+        assemble_prehook!(problem, time)
+    end
+
+    if problem.properties.variational
+        for element in get_elements(problem)
+            assemble!(problem.assembly, problem, element, time)
+        end
+    else # nodal collocation
+        field_vals = Dict{Int64, Float64}()
+        field_name = get_parent_field_name(problem)
+        field_dim = get_unknown_field_dimension(problem)
+        for element in get_elements(problem)
+            gdofs = get_gdofs(problem, element)
+            for i=1:field_dim
+                haskey(element, field_name*" $i") || continue
+                ldofs = gdofs[i:field_dim:end]
+                xis = get_reference_coordinates(typeof(element.properties))
+                vals = Float64[]
+                for xi in xis
+                    g = element(field_name*" $i", xi, time)
+                    # u = u_prev + Δu ⇒ Δu = u - u_prev
+                    if haskey(element, field_name)
+                        g_prev = element(field_name, xi, time)
+                        g -= g_prev[i]
+                    end
+                    push!(vals, g)
+                end
+                for (dof, g) in zip(ldofs, vals)
+                    field_vals[dof] = g
+                end
+            end
+        end
+        for (k, v) in field_vals
+            push!(problem.assembly.C1, k, k, 1.0)
+            push!(problem.assembly.C2, k, k, 1.0)
+            push!(problem.assembly.g, k, 1, v)
+        end
+    end
+
+    if method_exists(assemble_posthook!, Tuple{typeof(problem), Float64})
+        assemble_posthook!(problem, time)
+    end
+end
+
+function assemble!(assembly::Assembly, problem::Problem{Dirichlet},
+                   element::Element, time::Float64)
 
     # get dimension and name of PARENT field
     nnodes = length(element)
     field_dim = get_unknown_field_dimension(problem)
     field_name = get_parent_field_name(problem)
     gdofs = get_gdofs(element, field_dim)
+    props = problem.properties
 
     if problem.properties.dual_basis
         De, Me, Ae = get_dualbasis(element, time)
     else
         Ae = eye(nnodes)
         De = zeros(nnodes, nnodes)
-        for ip in get_integration_points(element, 1)
+        for ip in get_integration_points(element, props.order)
             N = element(ip, time)
             detJ = element(ip, time, Val{:detJ})
             De += ip.weight*N'*N*detJ
@@ -53,7 +120,7 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele
     end
 
     # right hand side
-    for ip in get_integration_points(element, 1)
+    for ip in get_integration_points(element, props.order)
         detJ = element(ip, time, Val{:detJ})
         w = ip.weight*detJ
         N = element(ip, time)