improved eigenvalue solver code

src/solvers_modal.jl

@@ -23,24 +23,15 @@ function Modal(nev=10, which=:SM)
     solver = Modal(false, Vector(), Matrix(), nev, which)
 end
 
-function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible=false)
+function call(solver::Solver{Modal}; show_info=true, debug=false,
+              bc_invertible=false, P=nothing,
+              empty_assemblies_before_solution=true)
     show_info && info(repeat("-", 80))
     show_info && info("Starting natural frequency solver")
     show_info && info("Increment time t=$(round(solver.time, 3))")
     show_info && info(repeat("-", 80))
     initialize!(solver)
-    # assemble all field problems
-    info("Assembling problems ...")
-    tic()
-    for problem in get_field_problems(solver)
-        assemble!(problem, solver.time)
-        assemble!(problem, solver.time, Val{:mass_matrix})
-    end
-    for problem in get_boundary_problems(solver)
-        assemble!(problem, solver.time)
-    end
-    t1 = round(toq(), 2)
-    info("Assembled in $t1 seconds.")
+    assemble!(solver; with_mass_matrix=true)
     M, K, Kg, f = get_field_assembly(solver)
     Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
     K = K + Kb
@@ -49,6 +40,14 @@ function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible=
         K += Kg
     end
 
+    # free up some memory before solution
+    if empty_assemblies_before_solution
+        for problem in get_problems(solver)
+            empty!(problem.assembly)
+        end
+        gc()
+    end
+
     @assert nnz(D) == 0
     @assert C1 == C2
 
@@ -56,26 +55,38 @@ function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible=
 
     nboundary_problems = length(get_boundary_problems(solver))
 
-    if nboundary_problems != 0
+    if !(P == nothing)
+        info("using custom P")
+        K_red = P'*K*P
+        M_red = P'*M*P
+    elseif nboundary_problems != 0
         if bc_invertible
+            info("Invertible C, calculating P")
             P, h = create_projection(C1, g, Val{:invertible})
         else
+            info("Contacts, calculate P")
             P, h = create_projection(C1, g)
         end
+        K_red = P'*K*P
+        M_red = P'*M*P
     else
+        info("No dirichlet boundaryes, P = I")
         P = speye(size(K, 1))
+        K_red = K
+        M_red = M
     end
 
-    K_red = P'*K*P
-    M_red = P'*M*P
-
-    # make sure matrices are symmetric
-    K_red = 1/2*(K_red + K_red')
-    M_red = 1/2*(M_red + M_red')
-
     t1 = round(toq(), 2)
     info("Eliminated dirichlet boundaries in $t1 seconds.")
 
+    # make sure matrices are symmetric
+    info("Making matrices symmetric")
+    tic()
+    K_red = 1/2*(K_red + K_red')
+    M_red = 1/2*(M_red + M_red')
+    t1 = round(toq(), 2)
+    info("Finished in $t1 seconds.")
+
     nz = get_nonzero_rows(K_red)
     ndofs = solver.ndofs
     props = solver.properties
@@ -116,8 +127,10 @@ function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible=
 
         rethrow()
     end
+    t1 = round(toq(), 2)
     info("Eigenvalues computed in $t1 seconds. Eigenvalues: $om2")
 
+    tic()
     props.eigvals = om2
     props.eigvecs = zeros(ndofs, length(om2))
     v = zeros(ndofs)
@@ -128,6 +141,7 @@ function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible=
     end
     t1 = round(toq(), 2)
 
+#=    
     for i=1:length(om2)
         freq = real(sqrt(om2[i])/(2.0*pi))
         u = props.eigvecs[:,i]
@@ -145,7 +159,99 @@ function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible=
             end
         end
     end
+=#
 
     return true
 end
 
+function update_xdmf!(solver::Solver{Modal}; show_info=true)
+    xdmf = get(solver.xdmf)
+    temporal_collection = get_temporal_collection(xdmf)
+
+    frame = new_element("Grid")
+    new_child(frame, "Time", Dict("Value" => solver.time))
+
+    # save geometry
+    X = solver("geometry", solver.time)
+    node_ids = sort(collect(keys(X)))
+    geometry = hcat([X[nid] for nid in node_ids]...)
+    ndim, nnodes = size(geometry)
+    geom_type = ndim == 2 ? "XY" : "XYZ"
+    dataitem = new_dataitem(xdmf, "/Node IDs", node_ids)
+    geom = new_child(frame, "Geometry", Dict("Type" => geom_type))
+    dataitem = new_dataitem(xdmf, "/Geometry", geometry)
+    add_child(geom, dataitem)
+
+    # save topology
+    all_elements = get_all_elements(solver)
+    nelements = length(all_elements)
+    element_types = unique(map(get_element_type, all_elements))
+
+    xdmf_element_mapping = Dict(
+        "Seg2" => "Polyline",
+        "Tri3" => "Triangle",
+        "Quad4" => "Quadrilateral",
+        "Tet4" => "Tetrahedron",
+        "Pyramid5" => "Pyramid",
+        "Wedge6" => "Wedge",
+        "Hex8" => "Hexahedron",
+        "Seg3" => "Edge_3",
+        "Tri6" => "Tri_6",
+        "Quad8" => "Quad_8",
+        "Tet10" => "Tet_10",
+        "Pyramid13" => "Pyramid_13",
+        "Wedge15" => "Wedge_15",
+        "Hex20" => "Hex_20")
+
+    for element_type in element_types
+        elements = filter_by_element_type(element_type, all_elements)
+        sort!(elements, by=get_element_id)
+        element_ids = map(get_element_id, elements)
+        element_conn = map(get_connectivity, elements)
+        element_conn = transpose(hcat(element_conn...)) - 1
+        element_code = split(string(element_type), ".")[end]
+        dataitem = new_dataitem(xdmf, "/Topology/$element_code/Element IDs", element_ids)
+        dataitem = new_dataitem(xdmf, "/Topology/$element_code/Connectivity", element_conn)
+        topology = new_child(frame, "Topology")
+        set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code])
+        set_attribute(topology, "NumberOfElements", length(elements))
+        add_child(topology, dataitem)
+    end
+
+    # save solved fields
+    unknown_field_name = get_unknown_field_name(solver)
+    U = solver(unknown_field_name, solver.time)
+    node_ids2 = sort(collect(keys(U)))
+
+    @assert node_ids == node_ids2
+    ndim = length(U[first(node_ids)])
+    field_type = ndim == 1 ? "Scalar" : "Vector"
+    field_center = "Node"
+    if ndim == 2
+        for nid in node_ids
+            U[nid] = [U[nid]; 0.0]
+        end
+        ndim = 3
+    end
+    U = hcat([U[nid] for nid in node_ids]...)
+    unknown_field_name = ucfirst(unknown_field_name)
+    time = solver.time
+    path = ""
+    if S == Nonlinear
+        iteration = solver.properties.iteration
+        path = "/Results/Time $time/Iteration $iteration/Nodal Fields/$unknown_field_name"
+    elseif S == Linear
+        path = "/Results/Time $time/Nodal Fields/$unknown_field_name"
+    end
+    dataitem = new_dataitem(xdmf, path, U)
+    attribute = new_child(frame, "Attribute")
+    set_attribute(attribute, "Name", unknown_field_name)
+    set_attribute(attribute, "Center", field_center)
+    set_attribute(attribute, "AttributeType", field_type)
+    add_child(attribute, dataitem)
+    if (S == Linear) || ((S == Nonlinear) && has_converged(solver))
+        add_child(temporal_collection, frame)
+    end
+    save!(xdmf)
+end
+