mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
Moved core functionality to FEMBase
Core functionality is moved to base package called FEMBase.jl. The aim is that when developing new elements, solvers, materials and so on, user only imports FEMBase.jl and uses the functionality there. JuliaFEM.jl is a sort of "metapackage" collecting together all the packages and features can be programmed in smaller packages focusing only on one thing. This structure makes it attractive to contribute smaller amount of code e.g. in the form of thesis. Moreover, FEMBase.jl is under 2000 lines of code, which will be very clearly documented thus everyone can understand the basic concepts behing JuliaFEM easily.
This commit is contained in:
Jukka Aho
+48
-77
@@ -8,91 +8,50 @@ This is JuliaFEM -- Finite Element Package
|
||||
"""
|
||||
module JuliaFEM
|
||||
|
||||
using FEMBase
|
||||
using FEMBase: SparseMatrixCOO, SparseVectorCOO, Node, BasisInfo,
|
||||
Discrete, Variable, TimeVariant, TimeInvariant, Field,
|
||||
DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment,
|
||||
IP, AbstractProblem, IntegrationPoint
|
||||
using FEMBase: is_field_problem, is_boundary_problem, get_elements,
|
||||
get_connectivity, assemble_prehook!, assemble_posthook!,
|
||||
get_parent_field_name, get_reference_coordinates,
|
||||
get_assembly, get_nonzero_rows, get_nonzero_columns,
|
||||
eval_basis!, get_basis, get_dbasis, grad!, get_dualbasis,
|
||||
assemble_mass_matrix!, get_local_coordinates, inside,
|
||||
get_element_type, filter_by_element_type, get_element_id,
|
||||
optimize!, resize_sparse, resize_sparsevec
|
||||
import FEMBase: get_unknown_field_name, get_unknown_field_dimension,
|
||||
assemble!, update!, initialize!
|
||||
|
||||
|
||||
# from other packages TimerOutputs.jl and Logging.jl
|
||||
using TimerOutputs
|
||||
export @timeit, print_timer
|
||||
|
||||
import Base: getindex, setindex!, convert, length, size, isapprox, similar,
|
||||
start, first, next, done, last, endof, vec, ==, +, -, *, /, haskey, copy,
|
||||
push!, isempty, empty!, append!, sparse, full, read
|
||||
|
||||
using FEMBasis
|
||||
using FEMBasis: AbstractBasis
|
||||
using FEMQuad
|
||||
using AbaqusReader
|
||||
using AsterReader
|
||||
|
||||
using Logging
|
||||
|
||||
Logging.configure(level=INFO)
|
||||
|
||||
if haskey(ENV, "JULIAFEM_LOGLEVEL")
|
||||
Logging.configure(level=LogLevel(ENV["JULIAFEM_LOGLEVEL"]))
|
||||
end
|
||||
|
||||
export info, debug
|
||||
|
||||
import Base: getindex, setindex!, convert, length, size, isapprox,
|
||||
similar, start, first, next, done, last, endof, vec,
|
||||
==, +, -, *, /, haskey, copy, push!, isempty, empty!,
|
||||
append!, sparse, full, read
|
||||
|
||||
module Testing
|
||||
|
||||
using Base.Test
|
||||
|
||||
export @test, @testset, @test_throws
|
||||
|
||||
using Base.Test
|
||||
export @test, @testset, @test_throws
|
||||
end
|
||||
|
||||
include("fields.jl")
|
||||
export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment
|
||||
|
||||
include("types.jl") # data types: Point, IntegrationPoint, ...
|
||||
export AbstractPoint, Point, IntegrationPoint, IP, Node
|
||||
|
||||
### ELEMENTS ###
|
||||
include("elements.jl") # common element routines
|
||||
export Node, Element, update!, get_connectivity, get_basis,
|
||||
get_dbasis, inside, get_local_coordinates, get_element_type,
|
||||
filter_by_element_type, get_element_id
|
||||
|
||||
include("elements_lagrange.jl") # Continuous Galerkin (Lagrange) elements
|
||||
export get_reference_coordinates,
|
||||
get_interpolation_polynomial,
|
||||
description
|
||||
export Poi1,
|
||||
Seg2, Seg3,
|
||||
Tri3, Tri6, Tri7,
|
||||
Quad4, Quad8, Quad9,
|
||||
Tet4, Tet10,
|
||||
Pyr5,
|
||||
Wedge6,
|
||||
Hex8, Hex20, Hex27
|
||||
|
||||
include("integrate.jl") # default integration points for elements
|
||||
export get_integration_points
|
||||
|
||||
include("sparse.jl")
|
||||
export add!, SparseMatrixCOO, SparseVectorCOO, get_nonzero_rows, get_nonzero_columns, optimize!, resize_sparse, resize_sparsevec
|
||||
|
||||
include("problems.jl") # common problem routines
|
||||
export Problem, AbstractProblem, FieldProblem, BoundaryProblem,
|
||||
get_unknown_field_dimension, get_gdofs, Assembly,
|
||||
get_parent_field_name, get_elements, add_elements!
|
||||
|
||||
using AbaqusReader
|
||||
using AsterReader
|
||||
include("problems_elasticity.jl")
|
||||
export Elasticity
|
||||
|
||||
include("materials_plasticity.jl")
|
||||
export plastic_von_mises
|
||||
|
||||
include("problems_dirichlet.jl")
|
||||
export Dirichlet
|
||||
|
||||
include("problems_heat.jl")
|
||||
export Heat
|
||||
|
||||
export assemble!, postprocess!
|
||||
|
||||
function assemble!(problem::Problem, element::Element, time=0.0)
|
||||
assemble!(problem.assembly, problem, element, time)
|
||||
end
|
||||
|
||||
### Mortar methods ###
|
||||
include("problems_mortar.jl")
|
||||
include("problems_mortar_2d.jl")
|
||||
@@ -101,12 +60,8 @@ include("problems_mortar_2d_autodiff.jl")
|
||||
export calculate_normals, calculate_normals!, project_from_slave_to_master,
|
||||
project_from_master_to_slave, Mortar, get_slave_elements,
|
||||
get_polygon_clip
|
||||
|
||||
include("io.jl")
|
||||
export Xdmf, h5file, xmffile, xdmf_filter, new_dataitem, update_xdmf!, save!
|
||||
|
||||
### ASSEMBLY + SOLVE ###
|
||||
include("assembly.jl")
|
||||
include("solvers.jl")
|
||||
export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver,
|
||||
get_unknown_field_name, get_formulation_type, get_problems,
|
||||
@@ -116,8 +71,6 @@ export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver,
|
||||
is_field_problem, is_boundary_problem
|
||||
include("solvers_modal.jl")
|
||||
export Modal
|
||||
|
||||
### Mortar methods, contact mechanics extension ###
|
||||
include("problems_contact.jl")
|
||||
include("problems_contact_2d.jl")
|
||||
include("problems_contact_3d.jl")
|
||||
@@ -125,22 +78,26 @@ include("problems_contact_2d_autodiff.jl")
|
||||
#include("problems_contact_3d_autodiff.jl")
|
||||
export Contact
|
||||
|
||||
# Preprocess module
|
||||
|
||||
module Preprocess
|
||||
using FEMBase
|
||||
include("preprocess.jl")
|
||||
export create_elements, Mesh, add_node!, add_nodes!, add_element!,
|
||||
add_elements!, add_element_to_element_set!, add_node_to_node_set!,
|
||||
find_nearest_nodes, find_nearest_node, reorder_element_connectivity!,
|
||||
create_node_set_from_element_set!, filter_by_element_set
|
||||
|
||||
include("preprocess_abaqus_reader.jl")
|
||||
export abaqus_read_mesh, create_surface_elements, create_nodal_elements
|
||||
|
||||
include("preprocess_aster_reader.jl")
|
||||
export aster_read_mesh
|
||||
|
||||
end
|
||||
|
||||
# Postprocess module
|
||||
|
||||
module Postprocess
|
||||
using FEMBase
|
||||
using FEMBase: get_elements
|
||||
include("postprocess_utils.jl")
|
||||
export calc_nodal_values!, get_nodal_vector, get_nodal_dict, copy_field!,
|
||||
calculate_area, calculate_center_of_mass,
|
||||
@@ -149,4 +106,18 @@ end
|
||||
|
||||
include("deprecations.jl")
|
||||
|
||||
export SparseMatrixCOO, SparseVectorCOO, optimize!, resize_sparse
|
||||
export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment
|
||||
export FieldProblem, BoundaryProblem, Problem, Node, Element, Assembly
|
||||
export Poi1, Seg2, Seg3, Tri3, Tri6, Tri7, Quad4, Quad8, Quad9,
|
||||
Tet4, Tet10, Pyr5, Wedge6, Wedge15, Hex8, Hex20, Hex27
|
||||
export update!, add_elements!, get_unknown_field_name, add!,
|
||||
is_field_problem, is_boundary_problem, get_gdofs,
|
||||
initialize!, get_integration_points, group_by_element_type,
|
||||
get_unknown_field_dimension, get_connectivity
|
||||
|
||||
export get_nonzero_rows, get_local_coordinates, inside, IP, get_element_type,
|
||||
get_elements, AbstractProblem, IntegrationPoint, filter_by_element_type,
|
||||
get_element_id, get_nonzero_columns, resize_sparse, resize_sparsevec
|
||||
|
||||
end
|
||||
|
||||
-173
@@ -1,173 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
function isapprox(a1::Assembly, a2::Assembly)
|
||||
T = isapprox(a1.K, a2.K)
|
||||
T &= isapprox(a1.C1, a2.C1)
|
||||
T &= isapprox(a1.C2, a2.C2)
|
||||
T &= isapprox(a1.D, a2.D)
|
||||
T &= isapprox(a1.f, a2.f)
|
||||
T &= isapprox(a1.g, a2.g)
|
||||
return T
|
||||
end
|
||||
|
||||
function assemble_prehook!
|
||||
end
|
||||
|
||||
function assemble_posthook!
|
||||
end
|
||||
|
||||
function assemble!(problem::Problem, time=0.0; auto_initialize=true)
|
||||
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
|
||||
assemble!(get_assembly(problem), problem, get_elements(problem), time)
|
||||
if method_exists(assemble_posthook!, Tuple{typeof(problem), Float64})
|
||||
assemble_posthook!(problem, time)
|
||||
end
|
||||
return true
|
||||
end
|
||||
|
||||
function assemble!(assembly::Assembly, problem::Problem, elements::Vector{Element}, time)
|
||||
warn("assemble!() this is default assemble operation, decreased performance can be expected without preallocation of memory!")
|
||||
for element in elements
|
||||
assemble!(assembly, problem, element, time)
|
||||
end
|
||||
return nothing
|
||||
end
|
||||
|
||||
function assemble_mass_matrix!(problem::Problem, time)
|
||||
if !isempty(problem.assembly.M)
|
||||
info("Mass matrix for $(problem.name) is already assembled, skipping assemble routine")
|
||||
return
|
||||
end
|
||||
elements = get_elements(problem)
|
||||
for (element_type, elements) in group_by_element_type(get_elements(problem))
|
||||
assemble_mass_matrix!(problem::Problem, elements, time)
|
||||
end
|
||||
return
|
||||
end
|
||||
|
||||
function assemble_mass_matrix!{Basis}(problem::Problem, elements::Vector{Element{Basis}}, time)
|
||||
nnodes = length(Basis)
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
M = zeros(nnodes, nnodes)
|
||||
N = zeros(1, nnodes)
|
||||
NtN = zeros(nnodes, nnodes)
|
||||
ldofs = zeros(Int, nnodes)
|
||||
for element in elements
|
||||
fill!(M, 0.0)
|
||||
for ip in get_integration_points(element, 2)
|
||||
detJ = element(ip, time, Val{:detJ})
|
||||
rho = element("density", time)
|
||||
w = ip.weight*rho*detJ
|
||||
eval_basis!(Basis, N, ip)
|
||||
N = element(ip, time)
|
||||
At_mul_B!(NtN, N, N)
|
||||
scale!(NtN, w)
|
||||
for i=1:nnodes^2
|
||||
M[i] += NtN[i]
|
||||
end
|
||||
end
|
||||
for (i, j) in enumerate(get_connectivity(element))
|
||||
@inbounds ldofs[i] = (j-1)*dim
|
||||
end
|
||||
for i=1:dim
|
||||
add!(problem.assembly.M, ldofs+i, ldofs+i, M)
|
||||
end
|
||||
end
|
||||
return
|
||||
end
|
||||
|
||||
"""
|
||||
assemble_mass_matrix!(problem, elements::Vector{Element{Tet10}}, time)
|
||||
|
||||
Assemble Tet10 mass matrices using special method. If Tet10 has constant metric
|
||||
if can be integrated analytically to gain performance.
|
||||
"""
|
||||
function assemble_mass_matrix!(problem::Problem, elements::Vector{Element{Tet10}}, time)
|
||||
nnodes = length(Tet10)
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
M = zeros(nnodes, nnodes)
|
||||
N = zeros(1, nnodes)
|
||||
NtN = zeros(nnodes, nnodes)
|
||||
ldofs = zeros(Int, nnodes)
|
||||
|
||||
M_CM = 1.0/2520.0 * [
|
||||
6 1 1 1 -4 -6 -4 -4 -6 -6
|
||||
1 6 1 1 -4 -4 -6 -6 -4 -6
|
||||
1 1 6 1 -6 -4 -4 -6 -6 -4
|
||||
1 1 1 6 -6 -6 -6 -4 -4 -4
|
||||
-4 -4 -6 -6 32 16 16 16 16 8
|
||||
-6 -4 -4 -6 16 32 16 8 16 16
|
||||
-4 -6 -4 -6 16 16 32 16 8 16
|
||||
-4 -6 -6 -4 16 8 16 32 16 16
|
||||
-6 -4 -6 -4 16 16 8 16 32 16
|
||||
-6 -6 -4 -4 8 16 16 16 16 32]
|
||||
|
||||
function is_CM(element::Element{Tet10}, X; rtol=1.0e-6)
|
||||
isapprox(X[5], 1/2*(X[1]+X[2]); rtol=rtol) || return false
|
||||
isapprox(X[6], 1/2*(X[2]+X[3]); rtol=rtol) || return false
|
||||
isapprox(X[7], 1/2*(X[3]+X[1]); rtol=rtol) || return false
|
||||
isapprox(X[8], 1/2*(X[1]+X[4]); rtol=rtol) || return false
|
||||
isapprox(X[9], 1/2*(X[2]+X[4]); rtol=rtol) || return false
|
||||
isapprox(X[10], 1/2*(X[3]+X[4]); rtol=rtol) || return false
|
||||
return true
|
||||
end
|
||||
|
||||
|
||||
n_CM = 0
|
||||
for element in elements
|
||||
for (i, j) in enumerate(get_connectivity(element))
|
||||
@inbounds ldofs[i] = (j-1)*dim
|
||||
end
|
||||
|
||||
X = element("geometry", time)
|
||||
rho = element("density", time)
|
||||
if is_CM(element, X) && length(rho) == 1
|
||||
ip = (1.0/3.0, 1.0/3.0, 1.0/3.0)
|
||||
detJ = element(ip, time, Val{:detJ})
|
||||
rho = element("density", ip, time)
|
||||
CM_s = detJ*rho
|
||||
n_CM += 1
|
||||
for i=1:dim
|
||||
add!(problem.assembly.M, ldofs+i, ldofs+i, CM_s * M_CM)
|
||||
end
|
||||
else
|
||||
fill!(M, 0.0)
|
||||
for ip in get_integration_points(element, 2)
|
||||
detJ = element(ip, time, Val{:detJ})
|
||||
rho = element("density", ip, time)
|
||||
w = ip.weight*rho*detJ
|
||||
eval_basis!(Tet10, N, ip)
|
||||
N = element(ip, time)
|
||||
At_mul_B!(NtN, N, N)
|
||||
scale!(NtN, w)
|
||||
for i=1:nnodes^2
|
||||
M[i] += NtN[i]
|
||||
end
|
||||
end
|
||||
for i=1:dim
|
||||
add!(problem.assembly.M, ldofs+i, ldofs+i, M)
|
||||
end
|
||||
end
|
||||
end
|
||||
info("$n_CM of $(length(elements)) was constant metric.")
|
||||
return
|
||||
end
|
||||
@@ -5,6 +5,10 @@ function assemble!(problem::Problem, time, ::Type{Val{:mass_matrix}})
|
||||
assemble_mass_matrix!(problem, time)
|
||||
end
|
||||
|
||||
function assemble!(problem::Problem, element::Element, time=0.0)
|
||||
assemble!(problem.assembly, problem, element, time)
|
||||
end
|
||||
|
||||
module Abaqus
|
||||
using JuliaFEM.Preprocess: create_surface_elements
|
||||
end
|
||||
|
||||
-418
@@ -1,418 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
type Element{E<:AbstractBasis}
|
||||
id :: Int
|
||||
connectivity :: Vector{Int}
|
||||
integration_points :: Vector{IP}
|
||||
fields :: Dict{String, Field}
|
||||
properties :: E
|
||||
end
|
||||
|
||||
"""
|
||||
Element(element_type, connectivity_vector)
|
||||
|
||||
Construct a new element where element_type is the type of the element
|
||||
and connectivity_vector is the vector of nodes that the element is connected to.
|
||||
|
||||
Examples
|
||||
--------
|
||||
In the example a new element (E in the figure below) of type Tri3 is created.
|
||||
This spesific element connects to nodes 89, 43, 12 in the finite element mesh.
|
||||
|
||||
```@example
|
||||
element = Element(Tri3, [89, 43, 12])
|
||||
```
|
||||
|
||||
"""
|
||||
function Element{E<:AbstractBasis}(::Type{E}, connectivity::Vector{Int})
|
||||
return Element{E}(-1, connectivity, [], Dict(), E())
|
||||
end
|
||||
|
||||
"""
|
||||
length(element::Element)
|
||||
|
||||
Return the number of nodes in element.
|
||||
"""
|
||||
function length{B}(element::Element{B})
|
||||
return length(B)
|
||||
end
|
||||
|
||||
function size{B}(element::Element{B})
|
||||
return size(B)
|
||||
end
|
||||
|
||||
function getindex(element::Element, field_name::AbstractString)
|
||||
return element.fields[field_name]
|
||||
end
|
||||
|
||||
function setindex!(element::Element, data::Field, field_name)
|
||||
element.fields[field_name] = data
|
||||
end
|
||||
|
||||
function get_element_type{E}(element::Element{E})
|
||||
return E
|
||||
end
|
||||
|
||||
function get_element_id{E}(element::Element{E})
|
||||
return element.id
|
||||
end
|
||||
|
||||
function is_element_type{E}(element::Element{E}, element_type)
|
||||
return E === element_type
|
||||
end
|
||||
|
||||
function filter_by_element_type(element_type, elements)
|
||||
return filter(element -> is_element_type(element, element_type), elements)
|
||||
end
|
||||
|
||||
"""
|
||||
group_by_element_type(elements::Vector{Element})
|
||||
|
||||
Given a vector of elements, group elements by element type to several vectors.
|
||||
Returns a dictionary, where key is the element type and value is a vector
|
||||
containing all elements of type `element_type`.
|
||||
"""
|
||||
function group_by_element_type(elements::Vector{Element})
|
||||
results = Dict{DataType, Any}()
|
||||
basis_types = map(element -> typeof(element.properties), elements)
|
||||
for basis in unique(basis_types)
|
||||
element_type = Element{basis}
|
||||
subset = filter(element -> isa(element, element_type), elements)
|
||||
results[element_type] = convert(Vector{element_type}, subset)
|
||||
end
|
||||
return results
|
||||
end
|
||||
|
||||
function setindex!(element::Element, data::Function, field_name)
|
||||
if method_exists(data, Tuple{Element, Vector, Float64})
|
||||
# create enclosure to pass element as argument
|
||||
function wrapper_(ip, time)
|
||||
return data(element, ip, time)
|
||||
end
|
||||
field = Field(wrapper_)
|
||||
else
|
||||
field = Field(data)
|
||||
end
|
||||
element.fields[field_name] = field
|
||||
end
|
||||
|
||||
function setindex!(element::Element, data, field_name)
|
||||
element.fields[field_name] = Field(data)
|
||||
end
|
||||
|
||||
""" Return a Field object from element.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> element = Element(Seg2, [1, 2])
|
||||
>>> data = Dict(1 => 1.0, 2 => 2.0)
|
||||
>>> update!(element, "my field", data)
|
||||
>>> element("my field")
|
||||
|
||||
"""
|
||||
function (element::Element)(field_name::String)
|
||||
return element[field_name]
|
||||
end
|
||||
|
||||
""" Return a Field object from element and interpolate in time direction.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> element = Element(Seg2, [1, 2])
|
||||
>>> data1 = Dict(1 => 1.0, 2 => 2.0)
|
||||
>>> data2 = Dict(1 => 2.0, 2 => 3.0)
|
||||
>>> update!(element, "my field", 0.0 => data1, 1.0 => data2)
|
||||
>>> element("my field", 0.5)
|
||||
|
||||
"""
|
||||
function (element::Element)(field_name::String, time::Float64)
|
||||
return element[field_name](time)
|
||||
end
|
||||
|
||||
function last(element::Element, field_name::String)
|
||||
return last(element[field_name])
|
||||
end
|
||||
|
||||
function (element::Element)(ip, time::Float64=0.0)
|
||||
return get_basis(element, ip, time)
|
||||
end
|
||||
|
||||
"""
|
||||
Examples
|
||||
|
||||
julia> el = Element(Quad4, [1, 2, 3, 4]);
|
||||
|
||||
julia> el([0.0, 0.0], 0.0, 1)
|
||||
1x4 Array{Float64,2}:
|
||||
0.25 0.25 0.25 0.25
|
||||
|
||||
julia> el([0.0, 0.0], 0.0, 2)
|
||||
2x8 Array{Float64,2}:
|
||||
0.25 0.0 0.25 0.0 0.25 0.0 0.25 0.0
|
||||
0.0 0.25 0.0 0.25 0.0 0.25 0.0 0.25
|
||||
|
||||
"""
|
||||
function (element::Element)(ip, time::Float64, dim::Int)
|
||||
dim == 1 && return get_basis(element, ip, time)
|
||||
Ni = vec(get_basis(element, ip, time))
|
||||
N = zeros(dim, length(element)*dim)
|
||||
for i=1:dim
|
||||
N[i,i:dim:end] += Ni
|
||||
end
|
||||
return N
|
||||
end
|
||||
|
||||
function (element::Element)(ip, time::Float64, ::Type{Val{:Jacobian}})
|
||||
X = element("geometry", time)
|
||||
dN = get_dbasis(element, ip, time)
|
||||
nbasis = length(element)
|
||||
if isa(X.data, Vector)
|
||||
J = sum([kron(dN[:,i], X[i]') for i=1:nbasis])
|
||||
else
|
||||
c = get_connectivity(element)
|
||||
J = sum([kron(dN[:,i], X[c[i]]') for i=1:nbasis])
|
||||
end
|
||||
return J
|
||||
end
|
||||
|
||||
function (element::Element)(ip, time::Float64, ::Type{Val{:detJ}})
|
||||
J = element(ip, time, Val{:Jacobian})
|
||||
n, m = size(J)
|
||||
if n == m # volume element
|
||||
return det(J)
|
||||
end
|
||||
JT = transpose(J)
|
||||
if size(JT, 2) == 1 # boundary of 2d problem, || ∂X/∂ξ ||
|
||||
return norm(JT)
|
||||
else # manifold on 3d problem, || ∂X/∂ξ₁ × ∂X/∂ξ₂ ||
|
||||
return norm(cross(JT[:,1], JT[:,2]))
|
||||
end
|
||||
end
|
||||
|
||||
function (element::Element)(ip, time::Float64, ::Type{Val{:Grad}})
|
||||
J = element(ip, time, Val{:Jacobian})
|
||||
return inv(J)*get_dbasis(element, ip, time)
|
||||
end
|
||||
|
||||
function (element::Element)(field_name::String, ip, time::Float64, ::Type{Val{:Grad}})
|
||||
return element(ip, time, Val{:Grad})*element[field_name](time)
|
||||
end
|
||||
|
||||
function (element::Element)(field_name::String, ip, time::Float64)
|
||||
field = element[field_name]
|
||||
return element(field, ip, time)
|
||||
end
|
||||
|
||||
function (element::Element)(field::DCTI, ip, time::Float64)
|
||||
return field.data
|
||||
end
|
||||
|
||||
function (element::Element)(field::DCTV, ip, time::Float64)
|
||||
return field(time).data
|
||||
end
|
||||
|
||||
function (element::Element)(field::CVTV, ip, time::Float64)
|
||||
return field(ip, time)
|
||||
end
|
||||
|
||||
function (element::Element)(field::Field, ip, time::Float64)
|
||||
field_ = field(time)
|
||||
basis = element(ip, time)
|
||||
n = length(element)
|
||||
if isa(field_.data, Vector)
|
||||
m = length(field_)
|
||||
if n != m
|
||||
error("Error when trying to interpolate field $field at coords $ip and time $time: element length is $n and field length is $m, f = Nᵢfᵢ makes no sense!")
|
||||
end
|
||||
return sum([field_[i]*basis[i] for i=1:n])
|
||||
else
|
||||
c = get_connectivity(element)
|
||||
return sum([field_[c[i]]*basis[i] for i=1:n])
|
||||
end
|
||||
end
|
||||
|
||||
function size(element::Element, dim)
|
||||
return size(element)[dim]
|
||||
end
|
||||
|
||||
""" Update element field based on a dictionary of nodal data and connectivity information.
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> data = Dict(1 => [0.0, 0.0], 2 => [1.0, 2.0])
|
||||
julia> element = Seg2([1, 2])
|
||||
julia> update!(element, "geometry", data)
|
||||
|
||||
As a result element now have time invariant (variable) vector field "geometry" with data ([0.0, 0.0], [1.0, 2.0]).
|
||||
|
||||
"""
|
||||
function update!{E}(element::Element{E}, field_name::AbstractString, data::Dict)
|
||||
#element[field_name] = Field(data)
|
||||
element_id = element.id
|
||||
local_connectivity = get_connectivity(element)
|
||||
for i in local_connectivity
|
||||
if !haskey(data, i)
|
||||
ndata = length(data)
|
||||
critical("Unable to set field data $field_name for element $E with
|
||||
id $element_id and connectivity $local_connectivity: no data for
|
||||
node id $i found. Length of data dict = $ndata")
|
||||
end
|
||||
end
|
||||
local_data = [data[i] for i in local_connectivity]
|
||||
element[field_name] = local_data
|
||||
end
|
||||
|
||||
function update!{K,V}(element::Element, field_name, data::Pair{Float64, Dict{K, V}})
|
||||
time, field_data = data
|
||||
element_data = V[field_data[i] for i in get_connectivity(element)]
|
||||
update!(element, field_name, time => element_data)
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Any}})
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
element[field_name] = data
|
||||
end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Int64}})
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
element[field_name] = data
|
||||
end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Float64}})
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
element[field_name] = data
|
||||
end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Vector{Float64}}})
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
element[field_name] = data
|
||||
end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Float64})
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
element[field_name] = data
|
||||
end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name::AbstractString, data::Union{Float64, Vector})
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
if length(data) != length(element)
|
||||
update!(element, field_name, DCTI(data))
|
||||
else
|
||||
element[field_name] = data
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
function update!(element::Element, datas::Pair...)
|
||||
for (field_name, data) in datas
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
element[field_name] = data
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name::String, data::Function)
|
||||
element[field_name] = data
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name::String, field::Field)
|
||||
element[field_name] = field
|
||||
end
|
||||
|
||||
function update!(elements::Vector, field_name::String, data)
|
||||
for element in elements
|
||||
update!(element, field_name, data)
|
||||
end
|
||||
end
|
||||
|
||||
""" Check existence of field. """
|
||||
function haskey(element::Element, field_name::String)
|
||||
haskey(element.fields, field_name)
|
||||
end
|
||||
|
||||
function get_connectivity(element::Element)
|
||||
return element.connectivity
|
||||
end
|
||||
|
||||
function get_integration_points{E}(element::Element{E})
|
||||
# first time initialize default integration points
|
||||
if length(element.integration_points) == 0
|
||||
ips = get_integration_points(element.properties)
|
||||
if E in (Poi1, Seg2, Seg3, NSeg)
|
||||
element.integration_points = [IP(i, w, (xi,)) for (i, (w, xi)) in enumerate(ips)]
|
||||
else
|
||||
element.integration_points = [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
|
||||
end
|
||||
end
|
||||
return element.integration_points
|
||||
end
|
||||
|
||||
""" This is a special case, temporarily change order
|
||||
of integration scheme mainly for mass matrix.
|
||||
"""
|
||||
function get_integration_points{E}(element::Element{E}, change_order::Int)
|
||||
ips = get_integration_points(element.properties, Val{change_order})
|
||||
if E in (Poi1, Seg2, Seg3, NSeg)
|
||||
return [IP(i, w, (xi,)) for (i, (w, xi)) in enumerate(ips)]
|
||||
else
|
||||
return [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
|
||||
end
|
||||
end
|
||||
|
||||
""" Return dual basis transformation matrix Ae. """
|
||||
function get_dualbasis(element::Element, time::Float64, order=1)
|
||||
nnodes = length(element)
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(element, order)
|
||||
detJ = element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
N = element(ip, time)
|
||||
De += w*diagm(vec(N))
|
||||
Me += w*N'*N
|
||||
end
|
||||
return De, Me, De*inv(Me)
|
||||
end
|
||||
|
||||
""" Find inverse isoparametric mapping of element. """
|
||||
function get_local_coordinates(element::Element, X::Vector, time::Float64; max_iterations=10, tolerance=1.0e-6)
|
||||
haskey(element, "geometry") || error("element geometry not defined, cannot calculate inverse isoparametric mapping")
|
||||
dim = size(element, 1)
|
||||
dim == length(X) || error("manifolds not supported.")
|
||||
xi = zeros(dim)
|
||||
dX = element("geometry", xi, time) - X
|
||||
for i=1:max_iterations
|
||||
J = element(xi, time, Val{:Jacobian})'
|
||||
xi -= J \ dX
|
||||
dX = element("geometry", xi, time) - X
|
||||
norm(dX) < tolerance && return xi
|
||||
end
|
||||
info("X = $X, dX = $dX, xi = $xi")
|
||||
error("Unable to find inverse isoparametric mapping for element $element for X = $X")
|
||||
end
|
||||
|
||||
""" Test is X inside element. """
|
||||
function inside{E}(element::Element{E}, X, time)
|
||||
xi = get_local_coordinates(element, X, time)
|
||||
return inside(E, xi)
|
||||
end
|
||||
@@ -1,70 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using FEMBasis
|
||||
|
||||
# "Poi1" => (0, 1),
|
||||
"1 node discrete point element",
|
||||
type Poi1 <: AbstractBasis
|
||||
end
|
||||
|
||||
function get_basis(element::Element{Poi1}, ip, time)
|
||||
return [1]
|
||||
end
|
||||
|
||||
function get_dbasis(element::Element{Poi1}, ip, time)
|
||||
return [0]
|
||||
end
|
||||
|
||||
function (element::Element{Poi1})(ip, time::Float64, ::Type{Val{:detJ}})
|
||||
return 1.0
|
||||
end
|
||||
|
||||
function get_integration_order(element::Poi1)
|
||||
return 1
|
||||
end
|
||||
|
||||
function get_integration_points(element::Poi1, order::Int64)
|
||||
return [ (1.0, 0.0) ]
|
||||
end
|
||||
|
||||
function size(::Type{Poi1})
|
||||
return (0, 1)
|
||||
end
|
||||
|
||||
function length(::Type{Poi1})
|
||||
return 1
|
||||
end
|
||||
|
||||
function FEMBasis.get_reference_element_coordinates(::Type{Poi1})
|
||||
Vector{Float64}[[0.0]]
|
||||
end
|
||||
|
||||
function get_basis{B}(element::Element{B}, ip, time)
|
||||
T = typeof(first(ip))
|
||||
N = zeros(T, 1, length(B))
|
||||
eval_basis!(B, N, tuple(ip...))
|
||||
return N
|
||||
end
|
||||
|
||||
function get_dbasis{B}(element::Element{B}, ip, time)
|
||||
T = typeof(first(ip))
|
||||
dN = zeros(T, size(B)...)
|
||||
eval_dbasis!(B, dN, tuple(ip...))
|
||||
return dN
|
||||
end
|
||||
|
||||
function inside(::Union{Type{Seg2}, Type{Seg3}, Type{Quad4}, Type{Quad8},
|
||||
Type{Quad9}, Type{Pyr5}, Type{Hex8}, Type{Hex20},
|
||||
Type{Hex27}}, xi)
|
||||
return all(-1.0 .<= xi .<= 1.0)
|
||||
end
|
||||
|
||||
function inside(::Union{Type{Tri3}, Type{Tri6}, Type{Tri7}, Type{Tet4}, Type{Tet10}}, xi)
|
||||
return all(xi .>= 0.0) && (sum(xi) <= 1.0)
|
||||
end
|
||||
|
||||
function get_reference_coordinates{B}(element::Element{B})
|
||||
return get_reference_element_coordinates(B)
|
||||
end
|
||||
|
||||
-457
@@ -1,457 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
abstract type AbstractField end
|
||||
|
||||
abstract type Discrete<:AbstractField end
|
||||
abstract type Continuous<:AbstractField end
|
||||
abstract type Constant<:AbstractField end
|
||||
abstract type Variable<:AbstractField end
|
||||
abstract type TimeVariant<:AbstractField end
|
||||
abstract type TimeInvariant<:AbstractField end
|
||||
|
||||
type Field{A<:Union{Discrete, Continuous},
|
||||
B<:Union{Constant, Variable},
|
||||
C<:Union{TimeVariant, TimeInvariant},
|
||||
T}
|
||||
data :: T
|
||||
end
|
||||
|
||||
const FieldSet = Dict{String, Field}
|
||||
|
||||
### Different field combinations and other typealiases
|
||||
|
||||
const DCTI{T} = Field{Discrete, Constant, TimeInvariant, T}
|
||||
const DVTI{T} = Field{Discrete, Variable, TimeInvariant, T}
|
||||
const DCTV{T} = Field{Discrete, Constant, TimeVariant, T}
|
||||
const DVTV{T} = Field{Discrete, Variable, TimeVariant, T}
|
||||
const CCTI{T} = Field{Continuous, Constant, TimeInvariant, T}
|
||||
const CVTI{T} = Field{Continuous, Variable, TimeInvariant, T}
|
||||
const CCTV{T} = Field{Continuous, Constant, TimeVariant, T}
|
||||
const CVTV{T} = Field{Continuous, Variable, TimeVariant, T}
|
||||
|
||||
# Discrete fields
|
||||
|
||||
|
||||
""" Discrete, constant, time-invariant field. This is constant in both spatial
|
||||
direction and time direction, i.e. df/dX = 0 and df/dt = 0.
|
||||
|
||||
This is the most basic type of field having no anything special functionality.
|
||||
|
||||
Examples
|
||||
--------
|
||||
|
||||
julia> f = DCTI()
|
||||
julia> update!(f, 1.0)
|
||||
|
||||
Multiplying by constant works:
|
||||
|
||||
julia> 2*f
|
||||
2.0
|
||||
|
||||
Interpolation in time direction gives the same constant:
|
||||
|
||||
julia> f(1.0)
|
||||
1.0
|
||||
|
||||
By default, when calling Field with scalar, DCTI is assumed, i.e.
|
||||
|
||||
julia> Field(0.0) == DCTI(0.0)
|
||||
true
|
||||
|
||||
"""
|
||||
function DCTI()
|
||||
return DCTI(nothing)
|
||||
end
|
||||
|
||||
function DCTI{T}(a::T)
|
||||
return DCTI{T}(a)
|
||||
end
|
||||
|
||||
function DVTI{T}(a::T)
|
||||
return DVTI{T}(a)
|
||||
end
|
||||
|
||||
function DCTV{T}(a::T)
|
||||
return DCTV{T}(a)
|
||||
end
|
||||
|
||||
function DVTV{T}(a::T)
|
||||
return DVTV{T}(a)
|
||||
end
|
||||
|
||||
function CCTI{T}(a::T)
|
||||
return CCTI{T}(a)
|
||||
end
|
||||
|
||||
function CVTI{T}(a::T)
|
||||
return CVTI{T}(a)
|
||||
end
|
||||
|
||||
function CCTV{T}(a::T)
|
||||
return CCTV{T}(a)
|
||||
end
|
||||
|
||||
function CVTV{T}(a::T)
|
||||
return CVTV{T}(a)
|
||||
end
|
||||
|
||||
function Field()
|
||||
return DCTI()
|
||||
end
|
||||
|
||||
function Field{T}(data::T)
|
||||
return DCTI{T}(data)
|
||||
end
|
||||
|
||||
function ==(x::DCTI, y::DCTI)
|
||||
return ==(x.data, y.data)
|
||||
end
|
||||
|
||||
function ==(x::DCTI, y)
|
||||
return ==(x.data, y)
|
||||
end
|
||||
|
||||
function isapprox(x::DCTI, y::DCTI)
|
||||
isapprox(x.data, y.data)
|
||||
end
|
||||
|
||||
function isapprox(x::DCTI, y)
|
||||
isapprox(x.data, y)
|
||||
end
|
||||
|
||||
function length(f::DCTI)
|
||||
return 1
|
||||
end
|
||||
|
||||
function *(c::Number, f::DCTI)
|
||||
return c*f.data
|
||||
end
|
||||
|
||||
""" Kind of spatial interpolation of DCTI. """
|
||||
function *(N::Matrix, f::DCTI)
|
||||
@assert length(N) == 1
|
||||
return N[1]*f.data
|
||||
end
|
||||
|
||||
function update!(field::DCTI, data)
|
||||
field.data = data
|
||||
end
|
||||
|
||||
""" Interpolate time-invariant field in time direction. """
|
||||
function (field::DCTI)(time::Float64)
|
||||
return field.data
|
||||
end
|
||||
|
||||
""" Discrete, variable, time-invariant field. This is constant in time direction,
|
||||
but not in spatial direction, i.e. df/dt = 0 but df/dX != 0. The basic structure
|
||||
of data is Vector, and it is implicitly assumed that length of field matches to
|
||||
the number of shape functions, so that interpolation in spatial direction works.
|
||||
|
||||
Examples
|
||||
--------
|
||||
"""
|
||||
function DVTI()
|
||||
return DVTI([])
|
||||
end
|
||||
|
||||
""" For vector data, DVTI is automatically created.
|
||||
|
||||
julia> DVTI([1.0, 2.0]) == Field([1.0, 2.0])
|
||||
true
|
||||
|
||||
"""
|
||||
function Field(data::Vector)
|
||||
return DVTI(data)
|
||||
end
|
||||
|
||||
""" For dictionary data, DVTI is automatically created.
|
||||
|
||||
Define e.g. nodal coordinates in dictionary
|
||||
julia> X = Dict(1 => [1.0, 2.0], 2 => [3.0, 4.0])
|
||||
julia> Field(X) == DVTI(X)
|
||||
|
||||
"""
|
||||
function Field(data::Dict)
|
||||
return DVTI(data)
|
||||
end
|
||||
|
||||
function ==(x::DVTI, y::DVTI)
|
||||
return ==(x.data, y.data)
|
||||
end
|
||||
|
||||
function isapprox(x::DVTI, y)
|
||||
return isapprox(x.data, y)
|
||||
end
|
||||
|
||||
""" Default slicing of field.
|
||||
|
||||
julia> f = DVTI([1.0, 2.0])
|
||||
julia> f[1]
|
||||
1.0
|
||||
|
||||
"""
|
||||
function getindex(field::DVTI, i::Int64)
|
||||
return field.data[i]
|
||||
end
|
||||
|
||||
""" Multi-slicing of field.
|
||||
|
||||
julia> f = DVTI([1.0, 2.0, 3.0])
|
||||
julia> f[[1, 3]]
|
||||
[1.0, 3.0]
|
||||
|
||||
"""
|
||||
function getindex(field::DVTI, I::Array{Int64, 1})
|
||||
return [field.data[i] for i in I]
|
||||
end
|
||||
|
||||
function length(field::DVTI)
|
||||
return length(field.data)
|
||||
end
|
||||
|
||||
function start(field::DVTI)
|
||||
return 1
|
||||
end
|
||||
|
||||
function +(f1::DVTI, f2::DVTI)
|
||||
return DVTI(f1.data + f2.data)
|
||||
end
|
||||
|
||||
function -(f1::DVTI, f2::DVTI)
|
||||
return DVTI(f1.data - f2.data)
|
||||
end
|
||||
|
||||
function update!(field::DVTI, data::Union{Vector, Dict})
|
||||
field.data = data
|
||||
end
|
||||
|
||||
""" Take scalar product of DVTI and constant T. """
|
||||
function *(T::Number, field::DVTI)
|
||||
return DVTI(T*field.data)
|
||||
end
|
||||
|
||||
""" Take dot product of DVTI field and vector T. Vector length must match to the
|
||||
field length and this can be used mainly for interpolation purposes, i.e., u = ∑ Nᵢuᵢ.
|
||||
"""
|
||||
function *(T::Union{Vector, RowVector}, f::DVTI)
|
||||
@assert length(T) <= length(f)
|
||||
return sum([T[i]*f[i] for i=1:length(T)])
|
||||
end
|
||||
|
||||
""" Take outer product of DVTI field and matrix T. """
|
||||
function *(T::Matrix, f::DVTI)
|
||||
n, m = size(T)
|
||||
return sum([kron(T[:,i], f[i]') for i=1:m])'
|
||||
end
|
||||
|
||||
function vec(field::DVTI)
|
||||
return [field.data...;]
|
||||
end
|
||||
|
||||
""" Interpolate time-invariant field in time direction. """
|
||||
function (field::DVTI)(time::Float64)
|
||||
return field
|
||||
end
|
||||
|
||||
""" Create a similar DVTI field from vector data.
|
||||
|
||||
julia> f1 = DVTI(Vector[[1.0, 2.0], [3.0, 4.0]])
|
||||
julia> f2 = similar(f1, [2.0, 3.0, 4.0, 5.0])
|
||||
julia> f2 == DVTI(Vector[[2.0, 3.0], [4.0, 5.0]])
|
||||
true
|
||||
|
||||
"""
|
||||
function similar(field::DVTI, data::Vector)
|
||||
n = length(field)
|
||||
m = length(data)
|
||||
dim = round(Int, m/n)
|
||||
@assert dim*n == m
|
||||
new_data = reshape(data, dim, n)
|
||||
new_field = DVTI()
|
||||
new_field.data = [new_data[:,i] for i=1:n]
|
||||
return new_field
|
||||
end
|
||||
|
||||
function next(f::DVTI, state)
|
||||
return f.data[state], state+1
|
||||
end
|
||||
|
||||
function done(f::DVTI, s)
|
||||
return s > length(f.data)
|
||||
end
|
||||
|
||||
""" Simple time frame / increment to contain both time and data. """
|
||||
type Increment{T}
|
||||
time :: Float64
|
||||
data :: T
|
||||
end
|
||||
|
||||
""" Discrete, constant, time variant field. This is constant in spatial
|
||||
direction but non-constant in time direction, i.e. df/dX = 0 but df/dt != 0.
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> t0 = 0.0; t1=1.0; y0 = 0.0; y1 = 1.0
|
||||
julia> f = DCTV(t0 => y0, t1 => y1)
|
||||
|
||||
"""
|
||||
function DCTV(data::Pair...)
|
||||
return DCTV([Increment(d[1],d[2]) for d in data])
|
||||
end
|
||||
|
||||
function Field{T}(data::Pair{Float64, T}...)
|
||||
return DCTV([Increment{T}(d[1], d[2]) for d in data])
|
||||
end
|
||||
|
||||
function getindex(field::DCTV, i::Int64)
|
||||
return field.data[i]
|
||||
end
|
||||
|
||||
function length(field::DCTV)
|
||||
return length(field.data)
|
||||
end
|
||||
|
||||
function first(field::DCTV)
|
||||
return field[1]
|
||||
end
|
||||
|
||||
""" Interpolate constant time-variant field in time direction. """
|
||||
function (field::DCTV)(time::Number)
|
||||
time < first(field).time && return DCTI(first(field).data)
|
||||
time > last(field).time && return DCTI(last(field).data)
|
||||
for i=reverse(1:length(field))
|
||||
isapprox(field[i].time, time) && return DCTI(field[i].data)
|
||||
end
|
||||
for i=reverse(2:length(field))
|
||||
t0 = field[i-1].time
|
||||
t1 = field[i].time
|
||||
if t0 < time < t1
|
||||
y0 = field[i-1].data
|
||||
y1 = field[i].data
|
||||
dt = t1-t0
|
||||
new_data = y0*(1-(time-t0)/dt) + y1*(1-(t1-time)/dt)
|
||||
return DCTI(new_data)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
function endof(field::DCTV)
|
||||
return endof(field.data)
|
||||
end
|
||||
|
||||
""" Discrete, variable, time variant fields. """
|
||||
function DVTV()
|
||||
return DVTV(Increment[])
|
||||
end
|
||||
|
||||
function DVTV{T<:Union{Vector, Dict}}(data::Pair{Float64, T}...)
|
||||
return DVTV([Increment{T}(d[1], d[2]) for d in data])
|
||||
end
|
||||
|
||||
function Field{T<:Union{Vector, Dict}}(data::Pair{Float64, T}...)
|
||||
return DVTV([Increment{T}(d[1], d[2]) for d in data])
|
||||
end
|
||||
|
||||
function length(field::DVTV)
|
||||
return length(field.data)
|
||||
end
|
||||
|
||||
function getindex(field::DVTV, i::Int64)
|
||||
return field.data[i]
|
||||
end
|
||||
|
||||
function first(field::DVTV)
|
||||
return field[1]
|
||||
end
|
||||
|
||||
function endof(field::DVTV)
|
||||
return endof(field.data)
|
||||
end
|
||||
|
||||
""" Interpolate discrete, variable, time-variant field in time direction. """
|
||||
function (field::DVTV)(time::Float64)
|
||||
time < first(field).time && return DVTI(first(field).data)
|
||||
time > last(field).time && return DVTI(last(field).data)
|
||||
for i=reverse(1:length(field))
|
||||
isapprox(field[i].time, time) && return DVTI(field[i].data)
|
||||
end
|
||||
for i=reverse(2:length(field))
|
||||
t0 = field[i-1].time
|
||||
t1 = field[i].time
|
||||
if t0 < time < t1
|
||||
y0 = field[i-1].data
|
||||
y1 = field[i].data
|
||||
dt = t1-t0
|
||||
new_data = y0*(1-(time-t0)/dt) + y1*(1-(t1-time)/dt)
|
||||
return DVTI(new_data)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
""" Update time-dependent fields with new values.
|
||||
|
||||
Examples
|
||||
--------
|
||||
|
||||
julia> f = Field(0.0 => 1.0)
|
||||
julia> update!(f, 1.0 => 2.0)
|
||||
|
||||
Now field has two (time, value) pairs: (0.0, 1.0) and (1.0, 2.0)
|
||||
|
||||
Notes
|
||||
-----
|
||||
Time vector is assumed to be ordered t_i-1 < t_i < t_i+1. If updating
|
||||
field with already existing time the old value is replaced with new one.
|
||||
|
||||
"""
|
||||
function update!{T}(field::Union{DCTV, DVTV}, val::Pair{Float64, T})
|
||||
time, data = val
|
||||
if isapprox(last(field).time, time)
|
||||
last(field).data = data
|
||||
else
|
||||
push!(field.data, Increment(val...))
|
||||
end
|
||||
end
|
||||
|
||||
### Basic data structure for continuous field
|
||||
|
||||
type Basis
|
||||
basis :: Function
|
||||
dbasis :: Function
|
||||
end
|
||||
|
||||
### Convenient functions to create fields
|
||||
|
||||
function Field(func::Function)
|
||||
if method_exists(func, Tuple{})
|
||||
return CCTI(func)
|
||||
elseif method_exists(func, Tuple{Float64})
|
||||
return CCTV(func)
|
||||
elseif method_exists(func, Tuple{Vector})
|
||||
return CVTI(func)
|
||||
elseif method_exists(func, Tuple{Vector, Float64})
|
||||
return CVTV(func)
|
||||
else
|
||||
error("no proper definition found for function: check methods.")
|
||||
end
|
||||
end
|
||||
|
||||
### Accessing continuous fields
|
||||
|
||||
function (field::CCTI)(xi::Vector, time::Number)
|
||||
return field.data()
|
||||
end
|
||||
|
||||
function (field::CVTI)(xi::Vector, time::Number)
|
||||
return field.data(xi)
|
||||
end
|
||||
|
||||
function (field::CCTV)(xi::Vector, time::Number)
|
||||
return field.data(time)
|
||||
end
|
||||
|
||||
function (field::CVTV)(xi, time)
|
||||
return field.data(xi, time)
|
||||
end
|
||||
|
||||
@@ -1,58 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using FEMQuad: get_quadrature_points
|
||||
|
||||
# Default number of integration points for each element. First rule is the
|
||||
# default integration rule returned by `get_integration_points(element)`.
|
||||
# Sometimes we want to increase integration order, e.g. when integrating mass
|
||||
# matrix or boundary conditions. For that reason, additional rules are provied
|
||||
# in list, so e.g. `get_integration_points(element, 1)` returns the second rule,
|
||||
# `get_integration_points(element, 2)` third rule and so on. Rules should be
|
||||
# ordered so that picking next one integrates more accurately.
|
||||
integration_rule_mapping = (
|
||||
:Seg2 => (:GLSEG1, :GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5),
|
||||
:Seg3 => (:GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5),
|
||||
:NSeg => (:GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5),
|
||||
:Quad4 => (:GLQUAD4, :GLQUAD9, :GLQUAD16, :GLQUAD25),
|
||||
:Quad8 => (:GLQUAD9, :GLQUAD16, :GLQUAD25),
|
||||
:Quad9 => (:GLQUAD9, :GLQUAD16, :GLQUAD25),
|
||||
:NSurf => (:GLQUAD9, :GLQUAD16, :GLQUAD25),
|
||||
:Hex8 => (:GLHEX8, :GLHEX27, :GLHEX81, :GLHEX243),
|
||||
:Hex20 => (:GLHEX27, :GLHEX81, :GLHEX243),
|
||||
:Hex27 => (:GLHEX27, :GLHEX81, :GLHEX243),
|
||||
:NSolid => (:GLHEX27, :GLHEX81, :GLHEX243),
|
||||
:Tri3 => (:GLTRI1, :GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12),
|
||||
:Tri6 => (:GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12),
|
||||
:Tri7 => (:GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12),
|
||||
:Tet4 => (:GLTET1, :GLTET4, :GLTET5, :GLTET15),
|
||||
:Tet10 => (:GLTET4, :GLTET5, :GLTET15),
|
||||
:Pyr5 => (:GLPYR5, ),
|
||||
:Wedge6 => (:GLWED6, :GLWED21),
|
||||
:Wedge15 => (:GLWED21, ))
|
||||
|
||||
for (E, R) in integration_rule_mapping
|
||||
for i in 1:length(R)
|
||||
P = Val{R[i]}
|
||||
order = Val{i-1}
|
||||
if i == 1
|
||||
code = quote
|
||||
function get_integration_points(element::$E)
|
||||
return get_quadrature_points($P)
|
||||
end
|
||||
end
|
||||
else
|
||||
code = quote
|
||||
function get_integration_points(element::$E, ::Type{$order})
|
||||
return get_quadrature_points($P)
|
||||
end
|
||||
end
|
||||
end
|
||||
eval(code)
|
||||
end
|
||||
end
|
||||
|
||||
# All good codes needs a special case. Here we have it: Poi1
|
||||
function get_integration_points(element::Poi1)
|
||||
[ (1.0, 0.0) ]
|
||||
end
|
||||
-442
@@ -1,442 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
abstract type AbstractProblem end
|
||||
abstract type FieldProblem<:AbstractProblem end
|
||||
abstract type BoundaryProblem<:AbstractProblem end
|
||||
abstract type MixedProblem<:AbstractProblem end
|
||||
|
||||
"""
|
||||
General linearized problem to solve
|
||||
(K₁+K₂)Δu + C1'*Δλ = f₁+f₂
|
||||
C2Δu + D*Δλ = g
|
||||
"""
|
||||
type Assembly
|
||||
|
||||
M :: SparseMatrixCOO # mass matrix
|
||||
|
||||
# for field assembly
|
||||
K :: SparseMatrixCOO # stiffness matrix
|
||||
Kg :: SparseMatrixCOO # geometric stiffness matrix
|
||||
f :: SparseMatrixCOO # force vector
|
||||
fg :: SparseMatrixCOO #
|
||||
|
||||
# for boundary assembly
|
||||
C1 :: SparseMatrixCOO
|
||||
C2 :: SparseMatrixCOO
|
||||
D :: SparseMatrixCOO
|
||||
g :: SparseMatrixCOO
|
||||
c :: SparseMatrixCOO
|
||||
|
||||
u :: Vector{Float64} # solution vector u
|
||||
u_prev :: Vector{Float64} # previous solution vector u
|
||||
u_norm_change :: Real # change of norm in u
|
||||
|
||||
la :: Vector{Float64} # solution vector la
|
||||
la_prev :: Vector{Float64} # previous solution vector u
|
||||
la_norm_change :: Real # change of norm in la
|
||||
|
||||
removed_dofs :: Vector{Int64} # manually remove dofs from assembly
|
||||
end
|
||||
|
||||
function Assembly()
|
||||
return Assembly(
|
||||
SparseMatrixCOO(),
|
||||
SparseMatrixCOO(),
|
||||
SparseMatrixCOO(),
|
||||
SparseMatrixCOO(),
|
||||
SparseMatrixCOO(),
|
||||
SparseMatrixCOO(),
|
||||
SparseMatrixCOO(),
|
||||
SparseMatrixCOO(),
|
||||
SparseMatrixCOO(),
|
||||
SparseMatrixCOO(),
|
||||
[], [], Inf,
|
||||
[], [], Inf,
|
||||
[])
|
||||
end
|
||||
|
||||
function empty!(assembly::Assembly)
|
||||
empty!(assembly.K)
|
||||
empty!(assembly.Kg)
|
||||
empty!(assembly.f)
|
||||
empty!(assembly.fg)
|
||||
empty!(assembly.C1)
|
||||
empty!(assembly.C2)
|
||||
empty!(assembly.D)
|
||||
empty!(assembly.g)
|
||||
empty!(assembly.c)
|
||||
end
|
||||
|
||||
function isempty(assembly::Assembly)
|
||||
T = isempty(assembly.K)
|
||||
T &= isempty(assembly.Kg)
|
||||
T &= isempty(assembly.f)
|
||||
T &= isempty(assembly.fg)
|
||||
T &= isempty(assembly.C1)
|
||||
T &= isempty(assembly.C2)
|
||||
T &= isempty(assembly.D)
|
||||
T &= isempty(assembly.g)
|
||||
T &= isempty(assembly.c)
|
||||
return T
|
||||
end
|
||||
|
||||
"""
|
||||
Defines types for Problem variables.
|
||||
|
||||
# Examples
|
||||
|
||||
The type of 'elements' is Vector{Element}
|
||||
|
||||
Add elements into the Problem element list.
|
||||
```@example
|
||||
a = [1, 2, 3]
|
||||
Problem.elements = a
|
||||
```
|
||||
|
||||
"""
|
||||
type Problem{P<:AbstractProblem}
|
||||
name :: AbstractString # descriptive name for the problem
|
||||
dimension :: Int # degrees of freedom per node
|
||||
parent_field_name :: AbstractString # (optional) name of the parent field e.g. "displacement"
|
||||
elements :: Vector{Element}
|
||||
dofmap :: Dict{Element, Vector{Int64}} # connects the element local dofs to the global dofs
|
||||
assembly :: Assembly
|
||||
fields :: Dict{AbstractString, Field}
|
||||
postprocess_fields :: Vector{String}
|
||||
properties :: P
|
||||
end
|
||||
|
||||
"""
|
||||
Problem(problem_type, problem_name::String, problem_dimension)
|
||||
|
||||
Construct a new field problem where `problem_type` is the type of the problem
|
||||
(Elasticity, Dirichlet, etc.), `problem_name` is the name of the problem and
|
||||
`problem_dimension` is the number of DOF:s in one node (2 in a 2D problem, 3
|
||||
in an elastic 3D problem, 6 in a 3D beam problem, etc.).
|
||||
|
||||
# Examples
|
||||
|
||||
Create a vector-valued (dim=3) elasticity problem:
|
||||
|
||||
```@example
|
||||
prob1 = Problem(Elasticity, "this is my problem", 3)
|
||||
```
|
||||
|
||||
"""
|
||||
function Problem{P<:FieldProblem}(::Type{P}, name::AbstractString, dimension::Int64)
|
||||
return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), Dict(), Vector(), P())
|
||||
end
|
||||
|
||||
"""
|
||||
Construct a new boundary problem.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Create a Dirichlet boundary problem for a vector-valued (dim=3) elasticity problem.
|
||||
|
||||
julia> bc1 = Problem(Dirichlet, "support", 3, "displacement")
|
||||
solver.
|
||||
"""
|
||||
function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name)
|
||||
return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), Vector(), P())
|
||||
end
|
||||
|
||||
function get_formulation_type(problem::Problem)
|
||||
return :incremental
|
||||
end
|
||||
|
||||
function get_unknown_field_name{P<:BoundaryProblem}(::Type{P})
|
||||
return "lambda"
|
||||
end
|
||||
|
||||
function get_assembly(problem)
|
||||
return problem.assembly
|
||||
end
|
||||
|
||||
"""
|
||||
update!(problem.properties, attr...)
|
||||
|
||||
Update properties for a problem.
|
||||
|
||||
# Example
|
||||
|
||||
```julia
|
||||
update!(body.properties, "finite_strain" => "false")
|
||||
```
|
||||
"""
|
||||
function update!{P<:AbstractProblem}(problem::P, attr::Pair{String, String}...)
|
||||
for (name, value) in attr
|
||||
debug("$P: set $name to $value")
|
||||
setfield!(problem, parse(name), parse(value))
|
||||
end
|
||||
end
|
||||
|
||||
"""
|
||||
function initialize!(problem_type, element_name, time)
|
||||
|
||||
Initialize the element ready for calculation, where `problem_type` is the type
|
||||
of the problem (Elasticity, Dirichlet, etc.), `element_name` is the name of a
|
||||
constructed element (see Element(element_type, connectivity_vector)) and `time`
|
||||
is the starting time of the initializing process.
|
||||
"""
|
||||
function initialize!(problem::Problem, element::Element, time::Float64)
|
||||
field_name = get_unknown_field_name(problem)
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
nnodes = length(element)
|
||||
|
||||
# initialize primary field
|
||||
if !haskey(element, field_name)
|
||||
if field_dim == 1
|
||||
update!(element, field_name, time => zeros(nnodes))
|
||||
else
|
||||
update!(element, field_name, time => [zeros(field_dim) for i=1:nnodes])
|
||||
end
|
||||
end
|
||||
|
||||
# if a boundary problem, initialize also a field for the main problem
|
||||
is_boundary_problem(problem) || return
|
||||
field_name = get_parent_field_name(problem)
|
||||
if !haskey(element, field_name)
|
||||
if field_dim == 1
|
||||
update!(element, field_name, time => zeros(nnodes))
|
||||
else
|
||||
update!(element, field_name, time => [zeros(field_dim) for i=1:nnodes])
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
function initialize!(problem::Problem, time::Float64=0.0)
|
||||
for element in get_elements(problem)
|
||||
initialize!(problem, element, time)
|
||||
end
|
||||
end
|
||||
|
||||
"""
|
||||
update!(problem, assembly, u, la)
|
||||
|
||||
Update the problem solution vector for assembly.
|
||||
"""
|
||||
function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector)
|
||||
|
||||
# resize & fill with zeros vectors if length mismatch with current solution
|
||||
|
||||
if length(u) != length(assembly.u)
|
||||
info("resizing solution vector u")
|
||||
resize!(assembly.u, length(u))
|
||||
fill!(assembly.u, 0.0)
|
||||
end
|
||||
|
||||
if length(la) != length(assembly.la)
|
||||
info("resizing lagrange multiplier vector la")
|
||||
resize!(assembly.la, length(la))
|
||||
fill!(assembly.la, 0.0)
|
||||
end
|
||||
|
||||
# copy current solutions to previous ones and add/replace new solution
|
||||
# TODO: here we have couple of options and they need to be clarified
|
||||
# for total formulation we are solving total quantity Ku = f while in
|
||||
# incremental formulation we solve KΔu = f and u = u + Δu
|
||||
assembly.u_prev = copy(assembly.u)
|
||||
assembly.la_prev = copy(assembly.la)
|
||||
|
||||
if get_formulation_type(problem) == :total
|
||||
assembly.u = u
|
||||
assembly.la = la
|
||||
elseif get_formulation_type(problem) == :incremental
|
||||
assembly.u += u
|
||||
assembly.la = la
|
||||
elseif get_formulation_type(problem) == :forwarddiff
|
||||
assembly.u += u
|
||||
assembly.la += la
|
||||
else
|
||||
info("$(problem.name): unknown formulation type, don't know what to do with results")
|
||||
error("serious failure with problem formulation: $(get_formulation_type(problem))")
|
||||
end
|
||||
|
||||
# calculate change of norm
|
||||
assembly.u_norm_change = norm(assembly.u - assembly.u_prev)
|
||||
assembly.la_norm_change = norm(assembly.la - assembly.la_prev)
|
||||
return assembly.u, assembly.la
|
||||
end
|
||||
|
||||
"""
|
||||
get_global_solution(problem, assembly)
|
||||
|
||||
Return a global solution (u, la) for a problem.
|
||||
|
||||
Notes
|
||||
-----
|
||||
If the length of solution vector != number of nodes, i.e. the field dimension is
|
||||
something else than 1, reshape vectors so that their length matches to the
|
||||
number of nodes. This helps to get nodal results easily.
|
||||
"""
|
||||
function get_global_solution(problem::Problem, assembly::Assembly)
|
||||
u = assembly.u
|
||||
la = assembly.la
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
if field_dim == 1
|
||||
return u, la
|
||||
else
|
||||
nnodes = round(Int, length(u)/field_dim)
|
||||
u = reshape(u, field_dim, nnodes)
|
||||
u = Vector{Float64}[u[:,i] for i in 1:nnodes]
|
||||
la = reshape(la, field_dim, nnodes)
|
||||
la = Vector{Float64}[la[:,i] for i in 1:nnodes]
|
||||
return u, la
|
||||
end
|
||||
end
|
||||
|
||||
"""
|
||||
update!(problem, assembly, elements, time)
|
||||
|
||||
Update a solution from the assebly to elements.
|
||||
"""
|
||||
function update!{P<:FieldProblem}(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64)
|
||||
u, la = get_global_solution(problem, assembly)
|
||||
field_name = get_unknown_field_name(problem)
|
||||
# update solution u for elements
|
||||
for element in elements
|
||||
connectivity = get_connectivity(element)
|
||||
update!(element, field_name, time => u[connectivity])
|
||||
end
|
||||
end
|
||||
|
||||
function update!{P<:BoundaryProblem}(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64)
|
||||
u, la = get_global_solution(problem, assembly)
|
||||
parent_field_name = get_parent_field_name(problem) # displacement
|
||||
field_name = get_unknown_field_name(problem) # lambda
|
||||
# update solution and lagrange multipliers for boundary elements
|
||||
for element in elements
|
||||
connectivity = get_connectivity(element)
|
||||
update!(element, parent_field_name, time => u[connectivity])
|
||||
update!(element, field_name, time => la[connectivity])
|
||||
end
|
||||
end
|
||||
|
||||
"""
|
||||
add_elements!(problem::Problem, elements)
|
||||
|
||||
Add new elements into the problem.
|
||||
"""
|
||||
function add_elements!(problem::Problem, elements)
|
||||
for element in elements
|
||||
push!(problem.elements, element)
|
||||
end
|
||||
end
|
||||
|
||||
function get_elements(problem::Problem)
|
||||
return problem.elements
|
||||
end
|
||||
|
||||
function get_assembly(problem::Problem)
|
||||
return problem.assembly
|
||||
end
|
||||
|
||||
function length(problem::Problem)
|
||||
return length(problem.elements)
|
||||
end
|
||||
|
||||
function update!(problem::Problem, field_name::AbstractString, data)
|
||||
#if haskey(problem.fields, field_name)
|
||||
# update!(problem.fields[field_name], field_name::AbstractString, data)
|
||||
#else
|
||||
# problem.fields[field_name] = Field(data)
|
||||
#end
|
||||
update!(problem.elements, field_name::AbstractString, data)
|
||||
end
|
||||
|
||||
function haskey(problem::Problem, field_name::AbstractString)
|
||||
return haskey(problem.fields, field_name)
|
||||
end
|
||||
|
||||
function getindex(problem::Problem, field_name::AbstractString)
|
||||
return problem.fields[field_name]
|
||||
end
|
||||
|
||||
""" Return field calculated to nodal points for elements in problem p. """
|
||||
function (problem::Problem)(field_name::AbstractString, time::AbstractFloat=0.0)
|
||||
#if haskey(problem, field_name)
|
||||
# return problem[field_name](time)
|
||||
#end
|
||||
f = nothing
|
||||
for element in get_elements(problem)
|
||||
haskey(element, field_name) || continue
|
||||
for (c, v) in zip(get_connectivity(element), element(field_name, time))
|
||||
if f == nothing
|
||||
f = Dict(c => v)
|
||||
continue
|
||||
end
|
||||
if haskey(f, c)
|
||||
if !isapprox(f[c], v)
|
||||
info("several values for single node when returning field $field_name")
|
||||
info("already have: $(f[c]), and trying to set $v")
|
||||
end
|
||||
else
|
||||
f[c] = v
|
||||
end
|
||||
end
|
||||
end
|
||||
#f == nothing && return f
|
||||
#update!(problem, field_name, time => f)
|
||||
return f
|
||||
end
|
||||
|
||||
""" Return the dimension of the unknown field of this problem. """
|
||||
function get_unknown_field_dimension(problem::Problem)
|
||||
return problem.dimension
|
||||
end
|
||||
|
||||
""" Return the name of the unknown field of this problem. """
|
||||
function get_unknown_field_name{P}(problem::Problem{P})
|
||||
return get_unknown_field_name(P)
|
||||
end
|
||||
|
||||
""" Return the name of the parent field of this (boundary) problem. """
|
||||
function get_parent_field_name{P<:BoundaryProblem}(problem::Problem{P})
|
||||
return problem.parent_field_name
|
||||
end
|
||||
|
||||
function push!(problem::Problem, elements...)
|
||||
push!(problem.elements, elements...)
|
||||
end
|
||||
|
||||
function push!(problem::Problem, elements::Vector)
|
||||
push!(problem.elements, elements...)
|
||||
end
|
||||
|
||||
function push!(problem::Problem, elements_::Vector...)
|
||||
for elements in elements_
|
||||
push!(problem.elements, elements...)
|
||||
end
|
||||
end
|
||||
|
||||
function get_gdofs(element::Element, dim::Int)
|
||||
conn = get_connectivity(element)
|
||||
if length(conn) == 0
|
||||
error("element connectivity not defined, cannot determine global dofs for element: $element")
|
||||
end
|
||||
gdofs = vec([dim*(i-1)+j for j=1:dim, i in conn])
|
||||
return gdofs
|
||||
end
|
||||
|
||||
function empty!(problem::Problem)
|
||||
empty!(problem.assembly)
|
||||
end
|
||||
|
||||
""" Return global degrees of freedom for element.
|
||||
|
||||
Notes
|
||||
-----
|
||||
First look dofs from problem.dofmap, it not found, update dofmap from
|
||||
element.element connectivity using formula gdofs = [dim*(nid-1)+j for j=1:dim]
|
||||
1. look element dofs from problem.dofmap
|
||||
2. if not found, use element.connectivity to update dofmap and 1.
|
||||
"""
|
||||
function get_gdofs(problem::Problem, element::Element)
|
||||
if !haskey(problem.dofmap, element)
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
problem.dofmap[element] = get_gdofs(element, dim)
|
||||
end
|
||||
return problem.dofmap[element]
|
||||
end
|
||||
@@ -52,12 +52,6 @@ function haskey(solver::Solver, field_name::String)
|
||||
return haskey(solver.fields, field_name)
|
||||
end
|
||||
|
||||
# one-liner helpers to identify problem types
|
||||
|
||||
is_field_problem(problem) = false
|
||||
is_field_problem{P<:FieldProblem}(problem::Problem{P}) = true
|
||||
is_boundary_problem(problem) = false
|
||||
is_boundary_problem{P<:BoundaryProblem}(problem::Problem{P}) = true
|
||||
get_field_problems(solver::Solver) = filter(is_field_problem, get_problems(solver))
|
||||
get_boundary_problems(solver::Solver) = filter(is_boundary_problem, get_problems(solver))
|
||||
|
||||
|
||||
-192
@@ -1,192 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# Sparse utils to make assembly of local and global matrices easier.
|
||||
# Unoptimized but should do all necessary stuff for at start.
|
||||
|
||||
type SparseMatrixCOO{T<:Real}
|
||||
I :: Vector{Int}
|
||||
J :: Vector{Int}
|
||||
V :: Vector{T}
|
||||
end
|
||||
|
||||
const SparseVectorCOO = SparseMatrixCOO
|
||||
|
||||
function SparseMatrixCOO()
|
||||
return SparseMatrixCOO{Float64}([], [], [])
|
||||
end
|
||||
|
||||
function SparseVectorCOO(I::Vector, V::Vector)
|
||||
return SparseVectorCOO(I, ones(I), V)
|
||||
end
|
||||
|
||||
function convert(::Type{SparseMatrixCOO}, A::SparseMatrixCSC)
|
||||
return SparseMatrixCOO(findnz(A)...)
|
||||
end
|
||||
|
||||
function convert(::Type{SparseVectorCOO}, A::SparseVector)
|
||||
return SparseVectorCOO(findnz(A)...)
|
||||
end
|
||||
|
||||
function convert(::Type{SparseMatrixCOO}, A::Matrix)
|
||||
return SparseMatrixCOO(findnz(A)...)
|
||||
end
|
||||
|
||||
""" Convert from COO format to CSC.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
tol
|
||||
used to drop near zero values less than tol.
|
||||
"""
|
||||
function sparse(A::SparseMatrixCOO; tol=1.0e-12)
|
||||
B = sparse(A.I, A.J, A.V)
|
||||
SparseArrays.droptol!(B, tol)
|
||||
return B
|
||||
end
|
||||
|
||||
function sparse(A::SparseMatrixCOO, n::Int, m::Int; tol=1.0e-12)
|
||||
B = sparse(A.I, A.J, A.V, n, m)
|
||||
SparseArrays.droptol!(B, tol)
|
||||
return B
|
||||
end
|
||||
|
||||
function sparse(A::SparseMatrixCOO, n::Int, m::Int, f::Function; tol=1.0e-12)
|
||||
B = sparse(A.I, A.J, A.V, n, m, f)
|
||||
SparseArrays.droptol!(B, tol)
|
||||
return B
|
||||
end
|
||||
|
||||
function push!(A::SparseMatrixCOO, I::Int, J::Int, V::Float64)
|
||||
push!(A.I, I)
|
||||
push!(A.J, J)
|
||||
push!(A.V, V)
|
||||
end
|
||||
|
||||
function empty!(A::SparseMatrixCOO)
|
||||
empty!(A.I)
|
||||
empty!(A.J)
|
||||
empty!(A.V)
|
||||
end
|
||||
|
||||
function append!(A::SparseMatrixCOO, B::SparseMatrixCOO)
|
||||
append!(A.I, B.I)
|
||||
append!(A.J, B.J)
|
||||
append!(A.V, B.V)
|
||||
end
|
||||
|
||||
function isempty(A::SparseMatrixCOO)
|
||||
return isempty(A.I) && isempty(A.J) && isempty(A.V)
|
||||
end
|
||||
|
||||
function full(A::SparseMatrixCOO, args...)
|
||||
return full(sparse(A.I, A.J, A.V, args...))
|
||||
end
|
||||
|
||||
""" Add local element matrix to sparse matrix. This basically does:
|
||||
|
||||
>>> A[dofs1, dofs2] = A[dofs1, dofs2] + data
|
||||
|
||||
Example
|
||||
-------
|
||||
|
||||
>>> S = [3, 4]
|
||||
>>> M = [6, 7, 8]
|
||||
>>> data = Float64[5 6 7; 8 9 10]
|
||||
>>> A = SparseMatrixCOO()
|
||||
>>> add!(A, S, M, data)
|
||||
>>> full(A)
|
||||
4x8 Array{Float64,2}:
|
||||
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 5.0 6.0 7.0
|
||||
0.0 0.0 0.0 0.0 0.0 8.0 9.0 10.0
|
||||
|
||||
"""
|
||||
function add!(A::SparseMatrixCOO, dofs1::Vector{Int}, dofs2::Vector{Int}, data::Matrix)
|
||||
n, m = size(data)
|
||||
for j=1:m
|
||||
for i=1:n
|
||||
push!(A.I, dofs1[i])
|
||||
push!(A.J, dofs2[j])
|
||||
end
|
||||
end
|
||||
append!(A.V, vec(data))
|
||||
end
|
||||
|
||||
""" Add sparse matrix of CSC to COO. """
|
||||
function add!(A::SparseMatrixCOO, B::SparseMatrixCSC)
|
||||
I, J, V = findnz(B)
|
||||
C = SparseMatrixCOO(I, J, V)
|
||||
append!(A, C)
|
||||
end
|
||||
|
||||
""" Add new data to COO Sparse vector. """
|
||||
function add!(A::SparseMatrixCOO, dofs::Vector{Int}, data::Array{Float64}, dim::Int=1)
|
||||
if length(dofs) != length(data)
|
||||
info("dofs = $dofs")
|
||||
info("data = $(vec(data))")
|
||||
error("when adding to sparse vector dimension mismatch!")
|
||||
end
|
||||
append!(A.I, dofs)
|
||||
append!(A.J, dim*ones(Int, length(dofs)))
|
||||
append!(A.V, vec(data))
|
||||
end
|
||||
|
||||
""" Add SparseVector to SparseVectorCOO. """
|
||||
function add!(a::SparseVectorCOO, b::SparseVector)
|
||||
I, V = findnz(b)
|
||||
c = SparseVectorCOO(I, V)
|
||||
append!(a, c)
|
||||
return
|
||||
end
|
||||
|
||||
""" Combine (I,J,V) values if possible to reduce memory usage. """
|
||||
function optimize!(A::SparseMatrixCOO)
|
||||
I, J, V = findnz(sparse(A))
|
||||
A.I = I
|
||||
A.J = J
|
||||
A.V = V
|
||||
return
|
||||
end
|
||||
|
||||
""" Find all nonzero rows from sparse matrix.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
||||
Ordered list of row indices.
|
||||
"""
|
||||
function get_nonzero_rows(A::SparseMatrixCSC)
|
||||
return sort(unique(rowvals(A)))
|
||||
end
|
||||
|
||||
function get_nonzero_columns(A::SparseMatrixCSC)
|
||||
return get_nonzero_rows(transpose(A))
|
||||
end
|
||||
|
||||
function size(A::SparseMatrixCOO)
|
||||
isempty(A) && return (0, 0)
|
||||
return maximum(A.I), maximum(A.J)
|
||||
end
|
||||
|
||||
function size(A::SparseMatrixCOO, idx::Int)
|
||||
return size(A)[idx]
|
||||
end
|
||||
|
||||
""" Resize sparse matrix A to (higher) dimension n x m. """
|
||||
function resize_sparse(A, n, m)
|
||||
return sparse(findnz(A)..., n, m)
|
||||
end
|
||||
|
||||
""" Resize sparse vector b to (higher) dimension n. """
|
||||
function resize_sparsevec(b, n)
|
||||
return sparsevec(findnz(b)..., n)
|
||||
end
|
||||
|
||||
""" Approximative comparison of two matricse A and B. """
|
||||
function isapprox(A::SparseMatrixCOO, B::SparseMatrixCOO)
|
||||
A2 = sparse(A)
|
||||
B2 = sparse(B, size(A2)...)
|
||||
return isapprox(A2, B2)
|
||||
end
|
||||
@@ -1,76 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
const Node = Vector{Float64}
|
||||
|
||||
abstract type AbstractPoint end
|
||||
|
||||
type Point{P<:AbstractPoint}
|
||||
id :: Int
|
||||
weight :: Float64
|
||||
coords :: Tuple{Vararg{Float64}}
|
||||
fields :: Dict{AbstractString, Field}
|
||||
properties :: P
|
||||
end
|
||||
|
||||
function setindex!{T}(point::Point, val::Pair{Float64, T}, field_name)
|
||||
point.fields[field_name] = Field(val)
|
||||
end
|
||||
|
||||
function getindex(point::Point, field_name)
|
||||
return point.fields[field_name]
|
||||
end
|
||||
|
||||
function getindex(point::Point, idx::Int)
|
||||
return point.coords[idx]
|
||||
end
|
||||
|
||||
function haskey(point::Point, field_name)
|
||||
return haskey(point.fields, field_name)
|
||||
end
|
||||
|
||||
function (point::Point)(field_name, time=0.0)
|
||||
point.fields[field_name](time).data
|
||||
end
|
||||
|
||||
function start(point::Point)
|
||||
return start(point.coords)
|
||||
end
|
||||
|
||||
function done(point::Point, i)
|
||||
return done(point.coords, i)
|
||||
end
|
||||
|
||||
function next(point::Point, i)
|
||||
return next(point.coords, i)
|
||||
end
|
||||
|
||||
function update!{T}(point::Point, field_name, val::Pair{Float64, T})
|
||||
if haskey(point, field_name)
|
||||
update!(point[field_name], val)
|
||||
else
|
||||
point[field_name] = val
|
||||
end
|
||||
end
|
||||
|
||||
#= TODO: in future
|
||||
type Node <: AbstractPoint
|
||||
end
|
||||
|
||||
type MaterialPoint <: AbstractPoint
|
||||
end
|
||||
=#
|
||||
|
||||
type IntegrationPoint <: AbstractPoint
|
||||
end
|
||||
|
||||
const IP = Point{IntegrationPoint}
|
||||
|
||||
function IP(id, weight, coords::Tuple)
|
||||
return IP(id, weight, coords, Dict(), IntegrationPoint())
|
||||
end
|
||||
|
||||
function IP(id, weight, coords::Vector)
|
||||
warn("Consider giving coords as Tuple.")
|
||||
return IP(id, weight, [c for c in coords], Dict(), IntegrationPoint())
|
||||
end
|
||||
Reference in New Issue
Block a user