From ef667e8f34a903e1c2a708e4d20c9fe791a835d1 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Fri, 27 Nov 2015 10:10:00 +0200 Subject: [PATCH] removed Equation type from code --- src/JuliaFEM.jl | 12 +- src/assembly.jl | 12 +- src/basis.jl | 263 --------------------- src/directsolver.jl | 59 ++--- src/dirichlet.jl | 65 +----- src/elasticity.jl | 145 +++--------- src/elements.jl | 171 ++++---------- src/equations.jl | 95 +++----- src/fields.jl | 429 +++++++++++++++++----------------- src/fields2.jl | 306 ------------------------ src/heat.jl | 88 ++----- src/integrate.jl | 12 +- src/lagrange.jl | 31 ++- src/mortar_elements.jl | 2 +- src/problems.jl | 85 ++----- src/solvers.jl | 142 ++++------- test/test_directsolver.jl | 132 +++++++++++ test/test_dirichlet.jl | 11 +- test/test_elasticity.jl | 15 +- test/test_heat.jl | 22 +- test/test_potential_energy.jl | 100 +++----- test/test_solver.jl | 284 ++-------------------- test/test_virtual_work.jl | 66 +++--- 23 files changed, 730 insertions(+), 1817 deletions(-) create mode 100644 test/test_directsolver.jl diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 78268ba..6ab0373 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -56,7 +56,7 @@ function Base.linspace{T<:Array}(X1::T, X2::T, n) end # fields, see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb -include("fields2.jl") +include("fields.jl") #include("basis.jl") # interpolation of discrete fields include("symbolic.jl") # a thin symbolic layer for fields include("types.jl") # type definitions @@ -65,24 +65,26 @@ include("types.jl") # type definitions include("elements.jl") include("lagrange.jl") # Lagrange elements #include("hierarchical.jl") # P-elements -include("mortar_elements.jl") # Mortar elements +#include("mortar_elements.jl") # Mortar elements ### EQUATIONS ### include("integrate.jl") # default integration points for elements include("sparse.jl") -include("equations.jl") include("problems.jl") +include("equations.jl") ### FORMULATIION ### include("dirichlet.jl") -include("mortar.jl") # mortar projection include("heat.jl") include("elasticity.jl") ### ASSEMBLY + SOLVE ### include("assembly.jl") include("solvers.jl") -include("directsolver.jl") # parallel sparse direct solver for non-lniear problems +include("directsolver.jl") # parallel sparse direct solver for non-linear problems + +### MORTAR STUFF ### +#include("mortar.jl") # mortar projection # PRE AND POSTPROCESS include("xdmf.jl") diff --git a/src/assembly.jl b/src/assembly.jl index 793a082..8222759 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -3,19 +3,19 @@ # Functions to handle global assembly of problem -function assemble!(assembly::Assembly, problem::Problem, time::Number=0.0, empty_assembly::Bool=true) +function assemble!(assembly::Assembly, problem::AllProblems, time::Float64, empty_assembly::Bool=true) if empty_assembly empty!(assembly) end - for equation in get_equations(problem) - assemble!(assembly, equation, time, problem) + for element in get_elements(problem) + assemble!(assembly, problem, element, time) end end -function assemble(problem::Problem, time::Number=0.0) +function assemble(problem::AllProblems, time::Float64) assembly = Assembly() - for equation in get_equations(problem) - assemble!(assembly, equation, time, problem) + for element in get_elements(problem) + assemble!(assembly, problem, element, time) end return assembly end diff --git a/src/basis.jl b/src/basis.jl index 62cacec..e49c6c3 100644 --- a/src/basis.jl +++ b/src/basis.jl @@ -1,266 +1,3 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -type Basis <: ContinuousField - basis :: Function - dbasisdxi :: Function -end - -""" Evaluate basis. """ -function Base.call(basis::Basis, xi::Vector, time::Number=0.0) - basis.basis(xi) # passing time does not make much sense actually for this... -end - -""" Evaluate gradient of basis. This need geometry information to calculate Jacobian. """ -function Base.call(basis::Basis, geometry::Increment, xi::Vector, - ::Type{Val{:grad}}) - dbasis = basis.dbasisdxi(xi) - J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)]) - grad = inv(J)*dbasis - return grad -end - -### INTERPOLATION IN SPATIAL DOMAIN ### - -""" Interpolate increment in spatial domain using Basis. """ -function Base.call(basis::Basis, increment::Increment, xi::Vector) - basis = basis.basis(xi) - sum([basis[i]*increment[i] for i=1:length(increment)]) -end - -""" Return gradient of increment in spatial domain using Basis.. """ -function Base.call(basis::Basis, geometry::Increment, field::Increment, - xi::Vector, ::Type{Val{:grad}}) - grad = basis(geometry, xi, Val{:grad}) - gradf = sum([grad[:,i]*field[i]' for i=1:length(field)])' - return gradf -end - -### INTERPOLATION IN TIME DOMAIN ### - -""" Interpolate discrete field in time domain. Return Increment. """ -function Base.call(field::DiscreteField, time::Number, - time_extrapolation::Symbol=:linear, - time_interpolation::Symbol=:linear) - - # special cases, only 1 timestep defined or time = -Inf -> return first ts - if (length(field) == 1) || (time == -Inf) - #return field[1][end] - return first(field) - end - - # special case, time = +Inf -> return last ts - if time == +Inf - #return field[end][end] - return last(field) - end - - # very likely we are always near some defined timestep, usually field - # defined only on t = 0.0, test neighbourhood for timesteps - for i=1:length(field) - if isapprox(field[i].time, time) - return field[i][end] - end - end - - # special case: out of time domain in positive direction, very likely - # to happen in incremental constitutive models - if time > field[end].time - if time_extrapolation == :constant - # constant time extrapolation, return last field - return field[end][end] - else - # multiple fields, pick last and second last and do linear interpolation - f1 = field[end-1] - f2 = field[end] - dt = abs(f2.time - f1.time) - i1 = f1[end] - i2 = f2[end] - di = i2 - i1 - increment = Increment(i2 + di./dt * (time-f2.time)) - return increment - end - end - - # special case: out of time domain in negative direction - if time < field[1].time - if time_extrapolation == :constant - # constant time extrapolation, return first field - return field[1][end] - else - # multiple fields, pick first and second and do linear interpolation - f1 = field[1] - f2 = field[2] - dt = abs(f2.time - f1.time) - i1 = f1[end] - i2 = f2[end] - di = i2 - i1 - increment = Increment(i1 - di./dt * (f1.time - time)) - return increment - end - end - - # find correct bin and perform interpolation - i = length(field) - while field[i].time >= time - i -= 1 - end - - if time_interpolation == :linear - t1 = field[i].time - t2 = field[i+1].time - inc1 = field[i][end] - inc2 = field[i+1][end] - dt = abs(t2 - t1) - t = (time-t1)/dt - increment = Increment((1-t)*inc1 + t*inc2) - return increment - end - - if time_interpolation == :constant - # nearest neightbour interpolation, i.e. pick nearest defined field - t1 = field[i].time - t2 = field[i+1].time - dt1 = abs(t1-time) - dt2 = abs(t2-time) - if dt1 < dt2 - return field[i][end] - else - return field[i+1][end] - end - end - -end - -""" Interpolate time derivative of field in some time t. This assumes linear -interpolation in time which is then differentiated. - -Parameters ----------- -field - Discrete field to interpolate. Must have timesteps and increments defined -time - Time to interpolate. -derivative - set Val{:diff} to activate this function - -""" -function Base.call(field::DiscreteField, time::Number, ::Type{Val{:diff}}) - - # FieldSet -> Field -> TimeStep -> Increment -> data - - if length(field) == 1 - # just one timestep, time derivative cannot be evaluated. - error("Field length = $(length(field)), cannot evaluate time derivative") - end - - function eval_field(i, j) - t1 = field[i] - t2 = field[j] - J = abs(t2.time - t1.time) - t = (time-t1.time)/J - result = 1/J*(-1*t1[end] + 1*t2[end]) - return Increment(result) - end - - # special cases, +Inf, -Inf, ~0.0 - if (time > field[end].time) || isapprox(time, field[end].time) - return eval_field(endof(field)-1, endof(field)) - end - - if (time < field[1].time) || isapprox(time, field[1].time) - return eval_field(1, 2) - end - - # search for a correct "bin" between time steps - i = length(field) - while (field[i].time > time) && !isapprox(field[i].time, time) - i -= 1 - end - - if isapprox(field[i].time, time) - # This is the hard case, maybe discontinuous time - # derivative if linear approximation. - # we are on the "mid node" in time axis - field1 = eval_field(i-1,i) - field2 = eval_field(i,i+1) - return 1/2*(field1 + field2) - end - - return eval_field(i, i+1) - -end - - - -### ELEMENT FIELD BASIS = ELEMENT BASIS + FIELD -#= -""" Here we add field we are wanting to interpolate with ElementBasis. """ -type ElementFieldBasis <: Basis - element_basis :: ElementBasis - field :: DiscreteField - time_extrapolation :: Symbol - time_interpolation :: Symbol -end - -function Basis(basis::Function, dbasisdxi::Function, field::DiscreteField, - time_extrapolation=:linear, time_interpolation=:linear) - element_basis = ElementBasis(basis, dbasisdxi) - return ElementFieldBasis(element_basis, field, time_extrapolation, - time_interpolation) -end - -function Base.call(basis::ElementFieldBasis, xi::Vector, time::Number) - increment = basis.field(time, basis.time_extrapolation, basis.time_interpolation) - return basis.element_basis(increment, xi) -end -=# - - -### ELEMENT GRADIENT BASIS = ELEMENT BASIS + GEOMETRY -#= -""" Gradient of ElementBasis, needs geometry information. """ -type ElementGradientBasis <: Basis - element_basis :: ElementBasis - geometry :: DiscreteField - time_extrapolation :: Symbol - time_interpolation :: Symbol -end - -function grad(N::ElementBasis, f::ElementFieldBasis, X::ElementFieldBasis) - f.time_extrapolation == X.time_extrapolation || error("interpolation mismatch") - f.time_interpolation == X.time_interpolation || error("interpolation mismatch") - dN = ElementGradientBasis(N, X.field, f.time_extrapolation, f.time_interpolation) - dfdX = ElementFieldGradientBasis(dN, f.field, f.time_extrapolation, f.time_interpolation) - return dfdX -end - -function grad(N::ElementBasis, f::DiscreteField, X::DiscreteField) - dN = ElementGradientBasis(N, X) - dfdX = ElementFieldGradientBasis(dN, f) -end - -function ElementGradientBasis(element_basis::ElementBasis, geometry::DiscreteField) - return ElementGradientBasis(element_basis, geometry, :linear, :linear) -end -=# - -### ELEMENT FIELD GRADIENT BASIS = ELEMENT GRADIENT BASIS + FIELD -#= -""" Gradient of ElementFieldBasis, needs field to interpolate. """ -type ElementFieldGradientBasis <: Basis - element_gradient_basis :: ElementGradientBasis - field :: DiscreteField - time_extrapolation :: Symbol - time_interpolation :: Symbol -end - -function ElementFieldGradientBasis(element_gradient_basis::ElementGradientBasis, - field::DiscreteField) - return ElementFieldGradientBasis(element_gradient_basis, field, :linear, :linear) -end -=# - -### INTERPOLATION IN TIME DOMAIN ### - - diff --git a/src/directsolver.jl b/src/directsolver.jl index 91aa2f9..2815827 100644 --- a/src/directsolver.jl +++ b/src/directsolver.jl @@ -4,7 +4,7 @@ ## Direct solver type DirectSolver <: Solver - field_problems :: Vector{FieldProblem} + field_problems :: Vector{Problem} boundary_problems :: Vector{BoundaryProblem} parallel :: Bool nonlinear_problem :: Bool @@ -12,7 +12,7 @@ type DirectSolver <: Solver tol :: Float64 end -function push!(solver::DirectSolver, problem::FieldProblem) +function push!(solver::DirectSolver, problem::Problem) push!(solver.field_problems, problem) end @@ -45,24 +45,30 @@ function call(solver::DirectSolver, time::Number=0.0) # for this increment for field_problem in solver.field_problems - for equation in get_equations(field_problem) - element = get_element(equation) - gdofs = get_gdofs(field_problem, equation) - if !isapprox(last(element[field_name]).time, time) - last_data = copy(last(element[field_name]).data) - push!(element[field_name], time => last_data) + for element in get_elements(field_problem) + gdofs = get_gdofs(element, field_dim) + if haskey(element, field_name) + if !isapprox(last(element[field_name]).time, time) + last_data = copy(last(element[field_name]).data) + push!(element[field_name], time => last_data) + end + else + data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)] + element[field_name] = (time => data) end end end for boundary_problem in solver.boundary_problems - for equation in get_equations(boundary_problem) - element = get_element(equation) - gdofs = get_gdofs(boundary_problem, equation) - eqdim = size(equation)[2] - data = Vector{Float64}[zeros(field_dim) for i in 1:eqdim] - if !isapprox(last(element["reaction force"]).time, time) - push!(element["reaction force"], time => data) + for element in get_elements(boundary_problem) + gdofs = get_gdofs(element, field_dim) + data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)] + if haskey(element, "reaction force") + if !isapprox(last(element["reaction force"]).time, time) + push!(element["reaction force"], time => data) + end + else + element["reaction force"] = (time => data) end end end @@ -105,32 +111,27 @@ function call(solver::DirectSolver, time::Number=0.0) # update elements in field problems for field_problem in solver.field_problems - for equation in get_equations(field_problem) - element = get_element(equation) - gdofs = get_gdofs(field_problem, equation) - eqsize = size(equation) + for element in get_elements(field_problem) + gdofs = get_gdofs(element, field_dim) local_sol = vec(full(sol[gdofs])) # incremental data for element - local_sol = reshape(local_sol, eqsize) - local_sol = Vector{Float64}[local_sol[:,i] for i=1:size(local_sol,2)] + local_sol = reshape(local_sol, field_dim, length(element)) + local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)] last(element[field_name]).data += local_sol # <-- added end end # update elements in boundary problems for boundary_problem in solver.boundary_problems - for equation in get_equations(boundary_problem) - element = get_element(equation) - gdofs = get_gdofs(boundary_problem, equation) + dim - eqsize = size(equation) + for element in get_elements(boundary_problem) + gdofs = get_gdofs(element, field_dim) + dim local_sol = vec(full(sol[gdofs])) - #info("local sol = $local_sol") - local_sol = reshape(local_sol, field_dim, eqsize[2]) - local_sol = Vector{Float64}[local_sol[:,i] for i=1:size(local_sol,2)] + local_sol = reshape(local_sol, field_dim, length(element)) + local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)] last(element["reaction force"]).data = local_sol # <-- replaced end end - info("Iteration took $(toq()) seconds") + info("Non-linear iteration took $(toq()) seconds") if norm(sol[1:dim]) < solver.tol return (iter, true) diff --git a/src/dirichlet.jl b/src/dirichlet.jl index 2ebada9..1c2312c 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -1,67 +1,26 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -# Dirichlet boundary conditions in weak form +abstract DirichletProblem <: AbstractProblem -abstract DirichletEquation <: BoundaryEquation - -### Dirichlet problem + equations - -type DirichletProblem <: BoundaryProblem - unknown_field_name :: ASCIIString - unknown_field_dimension :: Int - equations :: Vector{DirichletEquation} +function DirichletProblem(parent_field_name, parent_field_dim, dim=1, elements=[]) + return BoundaryProblem{DirichletProblem}(parent_field_name, parent_field_dim, dim, elements) end -""" Initialize new Dirichlet boundary condition. +function assemble!{E}(assembly::Assembly, problem::BoundaryProblem{DirichletProblem}, element::Element{E}, time::Number) -Parameters ----------- -dimension - dimension of unknown field (scalar, vector, ...) + # get dimension and name of PARENT field + field_dim = problem.parent_field_dim + field_name = problem.parent_field_name -Examples --------- - -""" -function DirichletProblem(unknown_field_name::ASCIIString, dimension::Int=1) - DirichletProblem(unknown_field_name, dimension, []) -end - -""" Dirichlet boundary condition element for 2 node line segment """ -type DBC2D2 <: DirichletEquation - element :: Seg2 - integration_points :: Vector{IntegrationPoint} -end - -function Base.size(equation::DBC2D2) - return (1, 2) -end - -function Base.convert(::Type{DirichletEquation}, element::Seg2) - integration_points = get_integration_points(element, Val{3}) - if !haskey(element, "reaction force") - element["reaction force"] = (0.0 => Vector{Float64}[]) - end - DBC2D2(element, integration_points) -end - - -function assemble!(assembly::Assembly, equation::DirichletEquation, time::Number=0.0, problem=nothing) - isa(problem, Void) && error("Dicihlet boundary condition needs problem defined") - field_dim = problem.unknown_field_dimension - field_name = problem.unknown_field_name - element = get_element(equation) gdofs = get_gdofs(element, field_dim) - basis = get_basis(element) - detJ = det(basis) - for ip in get_integration_points(equation) - w = ip.weight * detJ(ip) - N = basis(ip, time) + for ip in get_integration_points(element) + w = ip.weight * det(element, ip, time) + N = element(ip, time) A = w*N'*N if haskey(element, field_name) - # add all dimensions at once + # add all dimensions at once if defined element["blaa"] = 0.0 for i=1:field_dim g = element(field_name, ip, time) ldofs = gdofs[i:field_dim:end] @@ -71,8 +30,8 @@ function assemble!(assembly::Assembly, equation::DirichletEquation, time::Number end for i=1:field_dim + # add per dof if defined element["blaa 1"] = 1.0, element["blaa 2"] = 0.0 etc. if haskey(element, field_name*" $i") - # add single component g = element(field_name*" $i", ip, time) ldofs = gdofs[i:field_dim:end] add!(assembly.stiffness_matrix, ldofs, ldofs, A) diff --git a/src/elasticity.jl b/src/elasticity.jl index ed90847..5ba92a2 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -3,42 +3,20 @@ # Elasticity problems -abstract ElasticityProblem <: FieldProblem -abstract ElasticityEquation <: Equation +abstract ElasticityProblem <: AbstractProblem -function get_unknown_field_name(equation::ElasticityEquation) +abstract PlaneStressElasticityProblem <: ElasticityProblem + +function PlaneStressElasticityProblem(dim::Int=2, elements=[]) + return Problem{PlaneStressElasticityProblem}(dim, elements) +end + +function get_unknown_field_name{P<:ElasticityProblem}(::Type{P}) return "displacement" end -### Formulation ### - -""" Calculate internal energy for elasticity equation. - -Override this to define your own material model. By default we use -Saint Venant-Kirchhoff material model, which is simply - - S(E) = λtr(E) + 2μE -""" -function get_internal_energy(equation::ElasticityEquation, ip::IntegrationPoint, time::Number, F::Matrix) - element = get_element(equation) - basis = get_basis(element) - dbasis = grad(basis) - - # material parameters - young = basis("youngs modulus", ip, time) - poisson = basis("poissons ratio", ip, time) - mu = young/(2*(1+poisson)) - lambda = young*poisson/((1+poisson)*(1-2*poisson)) - if isa(equation, PlaneStressElasticityEquation) - lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d - end - - # material model - E = 1/2*(F'*F - I) # strain - S = lambda*trace(E)*I + 2*mu*E - P = F*S - - return P*dbasis(ip, time) +function get_unknown_field_type{P<:ElasticityProblem}(::Type{P}) + return Vector{Float64} end """ Elasticity equations. @@ -71,97 +49,44 @@ https://en.wikipedia.org/wiki/Plane_stress https://en.wikipedia.org/wiki/Hooke's_law """ -function get_residual_vector(equation::ElasticityEquation, ip::IntegrationPoint, time::Number; variation=nothing) +function get_residual_vector{EL<:CG}(problem::Problem{PlaneStressElasticityProblem}, element::Element{EL}, ip::IntegrationPoint, time::Number; variation=nothing) - element = get_element(equation) - basis = get_basis(element) - dbasis = grad(basis) + basis = element(ip, time) + dbasis = element(ip, time, Val{:grad}) - u = basis("displacement", ip, time, variation) - gradu = dbasis("displacement", ip, time, variation) + u = element("displacement", ip, time, variation) + gradu = element("displacement", ip, time, Val{:grad}, variation) F = I + gradu # deformation gradient - #info("Deformation gradient: $F") - # residual vector - internal energy - r = get_internal_energy(equation, ip, time, F) - #info("boundary element") + + # internal forces + young = element("youngs modulus", ip, time) + poisson = element("poissons ratio", ip, time) + mu = young/(2*(1+poisson)) + lambda = young*poisson/((1+poisson)*(1-2*poisson)) + lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d + E = 1/2*(F'*F - I) # strain + S = lambda*trace(E)*I + 2*mu*E + + r = F*S*dbasis # external forces - volume load if haskey(element, "displacement load") - b = basis("displacement load", ip, time) - r -= b*basis(ip, time) + b = element("displacement load", ip, time) + r -= b*basis end return vec(r) end +""" Surface load for plane stress model. """ +function get_residual_vector(problem::Problem{PlaneStressElasticityProblem}, element::Element{Seg2}, ip::IntegrationPoint, time::Number; variation=nothing) - - - -### Plane stress elasticity ### - -abstract PlaneElasticityProblem <: ElasticityProblem -abstract PlaneStressElasticityEquation <: ElasticityEquation - -type PlaneStressElasticityProblem <: PlaneElasticityProblem - unknown_field_name :: ASCIIString - unknown_field_dimension :: Int - equations :: Vector{PlaneStressElasticityEquation} -end - -function PlaneStressElasticityProblem(equations=[]) - return PlaneStressElasticityProblem("displacement", 2, equations) -end - -### Equations ### - -""" 4-node plane stress element. """ -type CPS4 <: PlaneStressElasticityEquation - element :: Quad4 - integration_points :: Vector{IntegrationPoint} -end - -function Base.size(equation::CPS4) - return (2, 4) -end - -function Base.convert(::Type{PlaneStressElasticityEquation}, element::Quad4) - integration_points = get_integration_points(element) - if !haskey(element, "displacement") - element["displacement"] = 0.0 => [zeros(2) for i=1:4] - end - CPS4(element, integration_points) -end - -""" Boundary element for plane stress problem for surface loads. """ -type CPS2 <: PlaneStressElasticityEquation - element :: Seg2 - integration_points :: Vector{IntegrationPoint} -end - -function Base.size(equation::CPS2) - return (2, 2) -end - -function Base.convert(::Type{PlaneStressElasticityEquation}, element::Seg2) - integration_points = get_integration_points(element) - if !haskey(element, "displacement") - element["displacement"] = 0.0 => [zeros(2) for i=1:2] - end - CPS2(element, integration_points) -end - -function get_residual_vector(equation::CPS2, ip::IntegrationPoint, time::Number; variation=nothing) - - element = get_element(equation) - basis = get_basis(element) - - u = basis("displacement", ip, time, variation) - r = zeros(size(equation)) + u = element("displacement", ip, time, variation) + r = zeros(problem.dim, length(element)) if haskey(element, "displacement traction force") - T = basis("displacement traction force", ip, time) - r -= T*basis(ip, time) + T = element("displacement traction force", ip, time) + r -= T*element(ip, time) end return vec(r) @@ -237,3 +162,5 @@ function get_residual_vector(equation::CPS2, ip::IntegrationPoint, time::Number; end =# + + diff --git a/src/elements.jl b/src/elements.jl index f273812..fdb1028 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -1,10 +1,23 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -using FactCheck -using ForwardDiff +abstract AbstractElement -abstract Element +type Element{E<:AbstractElement} + connectivity :: Vector{Int} +# integration_points :: Vector{IntegrationPoint} + fields :: Dict{ASCIIString, Field} +end + +function convert{E}(::Type{Element{E}}, connectivity::Vector{Int}) +# return Element{E}(connectivity, get_integration_points(E), Dict()) + return Element{E}(connectivity, Dict()) +end + +function get_integration_points{E}(element::Element{E}) +# return element.integration_points + return get_integration_points(E) +end """ Test routine for element. If this passes, element interface is properly @@ -67,9 +80,9 @@ function Base.getindex(element::Element, field_name) return element.fields[field_name] end -#function get_integration_points(element) -# return get_default_integration_points(element) -#end +function Base.length{E}(element::Element{E}) + size(E)[2] +end """Add new Field to element. @@ -89,151 +102,69 @@ function Base.setindex!(element::Element, data::Tuple, name::ASCIIString) element.fields[name] = Field(data...) end -#function Base.setindex!(element::Element, field_data::Tuple, field_name) -# field = Field() -# for (time, data) in field_data -# ts = TimeStep(time, Increment[Increment(data)]) -# push!(field, ts) -# end -# element[field_name] = field -#end - function get_connectivity(el::Element) return el.connectivity end -abstract AbstractFunctionSpace - -type FunctionSpace <: AbstractFunctionSpace - basis :: CVTI - fields :: FieldSet -end - -type GradientFunctionSpace <: AbstractFunctionSpace - basis :: CVTI - fields :: FieldSet -end - -function get_basis(element::Element) - return FunctionSpace(element.basis, element.fields) -end - -function get_dbasis(element::Element) - return GradientFunctionSpace(element.basis, element.fields) -end - -function grad(u::FunctionSpace) - return GradientFunctionSpace(u.basis, u.fields) -end - -""" If basis is called without a field, return basis functions evaluated at that point. """ -function call(u::FunctionSpace, xi::Union{Vector, IntegrationPoint}, t::Number=0.0) - return u.basis(xi) -end - -""" If gradient of basis is called without a field, return "empty" gradient evaluated at that point. """ -function call(gradu::GradientFunctionSpace, xi::Union{Vector, IntegrationPoint}, t::Number=0.0) - geometry = gradu.fields["geometry"](t) - gradu.basis(geometry, xi, Val{:grad}) -end - -""" Evaluate field on element function space. """ -function call(u::FunctionSpace, field_name, xi::Union{Vector, IntegrationPoint}, t::Number=0.0, variation=nothing) - field = !isa(variation, Void) ? variation : u.fields[field_name](t) - u.basis(field, xi) -end - -""" Evaluate gradient of field on element function space. """ -function call(gradu::GradientFunctionSpace, field_name, xi::Union{Vector, IntegrationPoint}, t::Number=0.0, variation=nothing) - field = !isa(variation, Void) ? variation : gradu.fields[field_name](t) - geometry = gradu.fields["geometry"](t) - gradu.basis(geometry, field, xi, Val{:grad}) -end - typealias VecOrIP Union{Vector, IntegrationPoint} -function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number) - return element.basis(element[field_name](time), xi) +function call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number, variation=nothing) + field = isa(variation, Void) ? element[field_name](time) : variation + basis = get_basis(element) + return basis(field, xi) end -function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number, ::Type{Val{:grad}}) - return element.basis(element["geometry"](time), element[field_name](time), xi, Val{:grad}) +function call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number, ::Type{Val{:grad}}, variation=nothing) + field = isa(variation, Void) ? element[field_name](time) : variation + basis = get_basis(element) + geom = element["geometry"](time) + return basis(geom, field, xi, Val{:grad}) end -function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP) +function call(element::Element, field_name::ASCIIString, xi::VecOrIP) return element.basis(element[field_name], xi) end -function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP, ::Type{Val{:grad}}) +function call(element::Element, field_name::ASCIIString, xi::VecOrIP, ::Type{Val{:grad}}) return element.basis(element["geometry"], element[field_name], xi, Val{:grad}) end -function Base.call(element::Element, field_name::ASCIIString, time::Number) +function call(element::Element, field_name::ASCIIString, time::Number) return element[field_name](time) end -function Base.call(element::Element, xi::VecOrIP) - element.basis(xi) +function get_basis{E}(element::Element{E}, ip::IntegrationPoint) + return get_basis(E, ip.xi) end -function Base.call(element::Element, xi::VecOrIP, ::Type{Val{:grad}}) - element.basis(element["geometry"], xi, Val{:grad}) +function call{E}(element::Element{E}, xi::VecOrIP, time::Float64=0) + return get_basis(element, xi) end -function Base.call(element::Element, field_name::ASCIIString) +function get_basis{E}(element::Element{E}) + basis = CVTI( + (xi::Vector) -> get_basis(E, xi), + (xi::Vector) -> get_dbasis(E, xi)) + return basis +end + +function call{E}(element::Element{E}, xi::VecOrIP, time::Float64, ::Type{Val{:grad}}) + basis = get_basis(element) + return basis(element["geometry"], xi, Val{:grad}) +end + +function call(element::Element, field_name::ASCIIString) return element[field_name] end -# on-line functions to get api more easy to use, ip -> xi.ip -#call(u::FunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t) -#call(u::GradientFunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t) -# i think these will be the most called functions. -#call(u::FunctionSpace, field_name, ip::IntegrationPoint, t::Number=0.0, variation=nothing) = call(u, field_name, ip.xi, t, variation) -#call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number=0.0, variation=nothing) = call(u, field_name, ip.xi, t, variation) -#call(u::FunctionSpace, field_name) = (args...) -> call(u, field_name, args...) -#call(u::GradientFunctionSpace, field_name) = (args...) -> call(u, field_name, args...) - -""" Return a field from function space. """ -function get_field(u::FunctionSpace, field_name, time::Number=0.0) - return u.fields[field_name](time) -end - -""" Return a field from function space. """ -function get_field(u::FunctionSpace, field_name, time::Number=0.0, variation=nothing) - return !isa(variation, Void) ? variation : u.fields[field_name](time) -end - -""" Return a field from function space. """ -function get_fieldset(u::FunctionSpace, field_name) - return u.fields[field_name] -end - -""" Get a determinant of element in point ξ. """ -function LinAlg.det(u::FunctionSpace, xi::Vector, time::Number=0.0) - X = u.fields["geometry"](time) - dN = u.basis(xi, Val{:grad}) +function LinAlg.det{E<:AbstractElement}(element::Element{E}, ip::IntegrationPoint, time::Number=0.0) + X = element("geometry", time) + dN = get_dbasis(E, ip.xi) J = sum([dN[:,i]*X[i]' for i=1:length(X)]) m, n = size(J) return m == n ? det(J) : norm(J) end -function LinAlg.det(u::FunctionSpace, ip::IntegrationPoint, time::Number=0.0) - LinAlg.det(u, ip.xi, time) -end - -function LinAlg.det(u::FunctionSpace) - return (args...) -> det(u, args...) -end - -function LinAlg.det(element::Element) - return det(get_basis(element)) -end - -#Base.(:+)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) + v(args...) -#Base.(:-)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) - v(args...) -#Base.(:+)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) + v(args...) -#Base.(:-)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) - v(args...) - """ Check does field exist. """ function Base.haskey(element::Element, what) haskey(element.fields, what) diff --git a/src/equations.jl b/src/equations.jl index 34a4f2e..4672bad 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -3,10 +3,6 @@ # Functions to handle element level things -- integration, assembly, ... -abstract Equation -abstract FieldEquation <: Equation -abstract BoundaryEquation <: Equation - type Assembly mass_matrix :: SparseMatrixIJV stiffness_matrix :: SparseMatrixIJV @@ -41,58 +37,31 @@ end function get_residual_vector end -function has_mass_matrix(equation::Equation) - default_args = Tuple{typeof(equation), IntegrationPoint, Float64} +function has_mass_matrix(problem::Problem, element::Element) + default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64} return method_exists(get_mass_matrix, default_args) end -function has_stiffness_matrix(equation::Equation) - default_args = Tuple{typeof(equation), IntegrationPoint, Float64} +function has_stiffness_matrix(problem::Problem, element::Element) + default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64} return method_exists(get_stiffness_matrix, default_args) end -function has_force_vector(equation::Equation) - default_args = Tuple{typeof(equation), IntegrationPoint, Float64} +function has_force_vector(problem::Problem, element::Element) + default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64} return method_exists(get_force_vector, default_args) end -function has_potential_energy(equation::Equation) - default_args = Tuple{typeof(equation), IntegrationPoint, Float64} +function has_potential_energy(problem::Problem, element::Element) + default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64} return method_exists(get_potential_energy, default_args) end -function has_residual_vector(equation::Equation) - default_args = Tuple{typeof(equation), IntegrationPoint, Float64} +function has_residual_vector(problem::Problem, element::Element) + default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64} return method_exists(get_residual_vector, default_args) end -function get_element(equation::Equation) - return equation.element -end - -function get_integration_points(equation::Equation) - return equation.integration_points -end - -function Base.size(equation::Equation, i::Int) - return size(equation)[i] -end - -""" Return global degrees of freedom of element in matrix level. - -Notes ------ -This is calculated from connectivity and equation dimension. -""" - -function get_gdofs(equation::Equation) - element = get_element(equation) - conn = get_connectivity(element) - dim = size(equation, 1) - gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...)) - return gdofs -end - function get_gdofs(element::Element, dim::Int) conn = get_connectivity(element) gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...)) @@ -100,26 +69,23 @@ function get_gdofs(element::Element, dim::Int) end """ Assemble element. """ -function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, problem=nothing) +function assemble!(assembly::Assembly, problem::Problem, element::Element, time::Number) - element = get_element(equation) - gdofs = get_gdofs(equation) - basis = get_basis(element) - detJ = det(basis) - unknown_field_name = get_unknown_field_name(equation) + gdofs = get_gdofs(element, problem.dim) + unknown_field_name = get_unknown_field_name(problem) # 1. if equations are defined we just integrate them, without caring how they are done - if has_mass_matrix(equation) || has_stiffness_matrix(equation) || has_force_vector(equation) - for ip in get_integration_points(equation) - s = ip.weight*detJ(ip) - if has_mass_matrix(equation) - add!(assembly.mass_matrix, gdofs, gdofs, s*get_mass_matrix(equation, ip, time)) + if has_mass_matrix(problem, element) || has_stiffness_matrix(problem, element) || has_force_vector(problem, element) + for ip in get_integration_points(element) + s = ip.weight*det(element, ip, time) + if has_mass_matrix(element) + add!(assembly.mass_matrix, gdofs, gdofs, s*get_mass_matrix(problem, element, ip, time)) end - if has_stiffness_matrix(equation) - add!(assembly.stiffness_matrix, gdofs, gdofs, s*get_stiffness_matrix(equation, ip, time)) + if has_stiffness_matrix(element) + add!(assembly.stiffness_matrix, gdofs, gdofs, s*get_stiffness_matrix(problem, element, ip, time)) end - if has_force_vector(equation) - add!(assembly.force_vector, gdofs, s*get_force_vector(equation, ip, time)) + if has_force_vector(element) + add!(assembly.force_vector, gdofs, s*get_force_vector(problem, element, ip, time)) end end # external loads -- if any nodal loads is defined add to force vector @@ -129,7 +95,7 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro end # 2. energy form -- user has defined potential energy W -> min! - if has_potential_energy(equation) + if has_potential_energy(problem, element) field = element[unknown_field_name](time) """ Wrapper for potential energy for ForwardDiff. """ @@ -137,9 +103,9 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro W = 0.0 df = similar(field, data) # integrate potential energy - for ip in get_integration_points(equation) - s = ip.weight*detJ(ip) - dw = get_potential_energy(equation, ip, time; variation=df) + for ip in get_integration_points(element) + s = ip.weight*det(element, ip, time) + dw = get_potential_energy(problem, element, ip, time; variation=df) W += s*dw end # external energy -- if any nodal loads is defined, decrease from potential energy @@ -156,7 +122,7 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro end # 3. virtual work -- user has defined some residual r = p - f = 0 - if has_residual_vector(equation) + if has_residual_vector(problem, element) field = DVTI(last(element[unknown_field_name]).data) @@ -165,9 +131,9 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro R = zeros(length(data)) df = similar(field, data) # integrate residual vector - for ip in get_integration_points(equation) - s = ip.weight*detJ(ip) - dr = get_residual_vector(equation, ip, time; variation=df) + for ip in get_integration_points(element) + s = ip.weight*det(element, ip, time) + dr = get_residual_vector(problem, element, ip, time; variation=df) R += s*dr end # external loads -- if any nodal loads is defined, decrease from residual @@ -183,3 +149,4 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro add!(assembly.force_vector, gdofs, -ForwardDiff.value(allresults)) end end + diff --git a/src/fields.jl b/src/fields.jl index e1d26e4..58b53f7 100644 --- a/src/fields.jl +++ b/src/fields.jl @@ -3,301 +3,310 @@ # https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb -abstract Field +abstract AbstractField -abstract DiscreteField <: Field -abstract ContinuousField <: Field +abstract Discrete <: AbstractField +abstract Continuous <: AbstractField +abstract Constant <: AbstractField +abstract Variable <: AbstractField +abstract TimeVariant <: AbstractField +abstract TimeInvariant <: AbstractField -### DEFAULT DISCRETE FIELD ### -# 1. Increment - -type Increment{T} <: AbstractVector{T} - data :: Vector{T} +type Field{A<:Union{Discrete,Continuous}, B<:Union{Constant,Variable}, C<:Union{TimeVariant,TimeInvariant}} + data end -function Base.size(increment::Increment) - return size(increment.data) +### Basic data structure for discrete field + +type Increment{T} + time :: Float64 + data :: T end -function Base.linearindexing(::Type{Increment}) - return LinearFast() +function Base.convert{T}(::Type{Increment{T}}, data::Pair{Float64,T}) + return Increment{T}(data[1], data[2]) end -function Base.getindex(increment::Increment, i::Int) +function Base.convert{T}(::Type{Increment{Vector{Vector{T}}}}, data::Pair{Float64, Matrix{T}}) + time = data[1] + content = data[2] + return Increment(time, Vector{T}[content[:,i] for i=1:size(content,2)]) +end + +function Base.getindex{T}(increment::Increment{Vector{T}}, i::Int64) return increment.data[i] end -function Base.setindex!(increment::Increment, v, i::Int) - increment.data[i] = v +function Base.(:*)(d, increment::Increment) + return d*increment.data end -function Base.dot(k::Number, increment::Increment) - return k*increment +### Basic data structure for continuous field + +type Basis + basis :: Function + dbasis :: Function end -function Base.convert(::Type{Increment}, data::Number) - return Increment([data]) +function Base.call(basis::Basis, xi::Vector) + basis.basis(xi) end -function Base.convert{T}(::Type{Increment}, data::Array{T, 2}) - return Increment([data[:,i] for i=1:size(data, 2)]) +function Base.call(basis::Basis, xi::Vector, ::Type{Val{:grad}}) + basis.dbasis(xi) end -function Base.convert{T}(::Type{Increment}, data::Array{T, 3}) - return Increment([data[:,:,i] for i=1:size(data, 3)]) -end +### Different field combinations and other typealiases -function Base.convert{T}(::Type{Increment}, data::Array{T, 4}) - return Increment([data[:,:,:,i] for i=1:size(data, 4)]) -end +typealias DCTI Field{Discrete, Constant, TimeInvariant} +typealias DVTI Field{Discrete, Variable, TimeInvariant} +typealias DCTV Field{Discrete, Constant, TimeVariant} +typealias DVTV Field{Discrete, Variable, TimeVariant} +typealias CCTI Field{Continuous, Constant, TimeInvariant} +typealias CVTI Field{Continuous, Variable, TimeInvariant} # can be used to interpolate in spatial dimension +typealias CCTV Field{Continuous, Constant, TimeVariant} # can be used to interpolate in time +typealias CVTV Field{Continuous, Variable, TimeVariant} -function Base.convert{T}(::Type{Increment}, data::Array{T, 5}) - return Increment([data[:,:,:,:,i] for i=1:size(data, 5)]) -end +typealias ScalarIncrement{T} Increment{T} +typealias VectorIncrement{T} Increment{Vector{T}} +typealias TensorIncrement{T} Increment{Matrix{T}} -function Base.zeros(::Type{Increment}, T, dims...) - return Increment(zeros(T, dims...)) -end +typealias DiscreteField Union{DCTI, DVTI, DCTV, DVTV} +typealias ContinuousField Union{CCTI, CVTI, CCTV, CVTV} +typealias ConstantField Union{DCTI, DCTV, CCTI, CCTV} +typealias VariableField Union{DVTI, DVTV, CVTI, CVTV} +typealias TimeInvariantField Union{DCTI, DVTI, CCTI, CVTI} +typealias TimeVariantField Union{DCTV, DVTV, CCTV, CVTV} -""" Flatten increment to Vector. -Examples --------- +### Convenient functions to create fields ->>> inc = ones(Increment, 2, 4) ->>> vec(inc) -[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] - -""" -function Base.vec(increment::Increment) - return [increment...;] -end - -function Base.similar{T}(increment::Increment, data::Vector{T}) - return Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment))) -end - -function Base.convert{T}(::Type{Vector{T}}, increment::Increment) - return Increment[increment] -end - -# 2. TimeStep - -type TimeStep - time :: Float64 - increments :: Vector{Increment} -end - -function Base.size(timestep::TimeStep) - return size(timestep.increments) -end - -function Base.endof(timestep::TimeStep) - return endof(timestep.increments) -end - -function Base.length(timestep::TimeStep) - return length(timestep.increments) -end - -function Base.linearindexing(::Type{TimeStep}) - return Base.LinearFast() -end - -function Base.getindex(timestep::TimeStep, i::Int) - return timestep.increments[i] -end - -#function TimeStep(data::Union{Number, Array}...) -# return TimeStep(0.0, Increment[Increment(d) for d in data]) +#function Base.convert(::Type{Field}, data) +# return Field(data) #end -function TimeStep() - return TimeStep(0.0, []) +function Field(data) + return DCTI(data) end -function TimeStep{T}(data::T...) - return TimeStep(0.0, Increment[Increment(d) for d in data]) +function Field(data::Vector) + return DVTI(data) end -function Base.convert(::Type{TimeStep}, value::Number) - return TimeStep(0.0, Increment[Increment(value)]) +function Field{T}(data::Pair{Float64, T}...) + return DCTV([Increment{T}(d[1], d[2]) for d in data]) end -function Base.push!(timestep::TimeStep, increment::Increment) - push!(timestep.increments, increment) +function Field{T}(data::Pair{Float64, Vector{T}}...) + return DVTV([Increment{Vector{T}}(d[1], d[2]) for d in data]) end -# FIXME: having some serious problems here to get tuple form working. - -# 3. DefaultDiscreteField -immutable DefaultDiscreteField <: DiscreteField - timesteps :: Vector{TimeStep} -#= - function DefaultDiscreteField(data::Array) - if (typeof(data) == Vector{Int64}) || (typeof(data) == Vector{Float64}) - new(TimeStep[TimeStep(data)]) - else - new(data) - end - end -=# +function Base.convert{T}(::Type{DCTV}, data::Pair{Float64, Vector{T}}...) + return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data]) end -#= -type DefaultDiscreteField <: DiscreteField - timesteps :: Vector{TimeStep} - function DefaultDiscreteField(data...) - timesteps = TimeStep[] - for (i, d) in enumerate(data) - @debug("i = $i, d = $d") - if isa(d, Tuple) - # contains time vector - increments = Increment[Increment(d[2])] - push!(timesteps, TimeStep(d[1], increments)) - else - increments = Increment[Increment(d)] - push!(timesteps, TimeStep(i-1.0, increments)) - end - end - new(timesteps) +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, Number}) + return CVTV(func) + else + error("no proper definition found for function: check methods.") end end -=# - -function Base.size(field::DefaultDiscreteField) - return size(field.timesteps) +function CVTI(basis::Function, dbasis::Function) + return CVTI(Basis(basis, dbasis)) end -function Base.length(field::DefaultDiscreteField) - return length(field.timesteps) +function Field(basis::Function, dbasis::Function) + return CVTI(basis, dbasis) end -function Base.start(::DefaultDiscreteField) +### Accessing and manipulating discrete fields + +function Base.getindex(field::DVTV, i::Int64) + return field.data[i] +end + +function Base.push!(field::DCTV, data::Pair) + push!(field.data, data) +end + +function Base.push!(field::DVTV, data::Pair) +# info("field.data = \n$(field.data)") +# info("data = \n$data") + push!(field.data, data) +end + +function Base.getindex(field::DVTV, i::Int64) + return field.data[i] +end + +function Base.getindex(field::DVTI, i::Int64) + return field.data[i] +end + +function Base.getindex(field::DCTV, i::Int64) + return field.data[i] +end + +function Base.getindex(field::Field, i::Int64) + return field.data[i] +end + +function Base.length(field::DVTI) + return length(field.data) +end + +function Base.length(field::DCTI) return 1 end -function Base.next(field::DefaultDiscreteField, state) - return (field[state+1], state+1) +function Base.length(field::DVTV) + return length(field.data) end -function Base.done(field::DefaultDiscreteField, state) - return state > length(field) +function Base.length(field::DCTV) + return length(field.data) end -function eltype(::Type{DefaultDiscreteField}) - return TimeStep +for op = (:+, :*, :/, :-) + @eval ($op)(increment::Increment, field::DCTI) = ($op)(increment.data, field.data) + @eval ($op)(field::DCTI, increment::Increment) = ($op)(increment.data, field.data) + @eval ($op)(field1::DCTI, field2::DCTI) = ($op)(field1.data, field2.data) + @eval ($op)(field::DCTI, k) = ($op)(field.data, k) + @eval ($op)(k, field::DCTI) = ($op)(field.data, k) end -function Base.linearindexing(::Type{DefaultDiscreteField}) - return LinearFast() +function Base.vec(field::DVTI) + return [field.data...;] end -function Base.getindex(field::DefaultDiscreteField, i::Int) - return field.timesteps[i] +function Base.vec(field::DCTV) + info("trying to vectorize $field") + error("does not make sense") end -function Base.endof(field::DefaultDiscreteField) - return endof(field.timesteps) +function Base.endof(field::Field) + return endof(field.data) end -function Base.first(field::DefaultDiscreteField) - return field[1][end] +#function Base.similar{T}(field::DVTI, data::Vector{T}) +# return Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment))) +#end + +function Base.similar{T}(field::DVTI, data::Vector{T}) + n = length(field.data) + data = reshape(data, round(Int, length(data)/n), n) + newdata = Vector[data[:,i] for i=1:n] + return typeof(field)(newdata) end -function Base.last(field::DefaultDiscreteField) - return field[end][end] +function Base.start(::DVTI) + return 1 end -function Base.push!(field::DefaultDiscreteField, timestep::TimeStep) - push!(field.timesteps, timestep) +function Base.next(f::DVTI, state) + return f.data[state], state+1 end -function Base.push!(field::DefaultDiscreteField, data::Union{Vector, Matrix}) - push!(field[end], Increment(data)) +function Base.done(f::DVTI, s) + return s > length(f.data) end -function Base.push!(field::DefaultDiscreteField, data::Pair) - ts = TimeStep(data[1], Increment(data[2])) - push!(field, ts) +### Accessing continuous fields + +function Base.call(field::CVTI, xi::Vector) + field.data(xi) end -"""Quickly create fields. +function Base.call(field::CVTI, xi::Vector, ::Type{Val{:grad}}) + field.data(xi, Val{:grad}) +end -Examples --------- ->>> Field([1, 2]) # creates field with one timestep and vector value [1, 2] ->>> Field(1, 2) # creates field with two timesteps, each having scalar value ->>> Field([1, 2], [3, 4]) # creates field with two timesteps, each having vector value ->>> Field( (0.0, [1, 2]), (0.5, [3, 4]) ) # like above, but give time also -""" -function Base.convert(::Type{DefaultDiscreteField}, data...) - timesteps = TimeStep[] - for (i, d) in enumerate(data) - if isa(d, Tuple) - @debug("is tuple, has time, d = $d") - # contains time vector - increments = Increment[Increment(d[2])] - push!(timesteps, TimeStep(d[1], increments)) - else -# @debug("array without time, d = $d") -# @debug(typeof(d)) - increments = Increment[Increment(d)] - push!(timesteps, TimeStep(i-1.0, increments)) +function Base.convert(::Type{Basis}, field::CVTI) + return field.data +end + +function Base.call(field::CCTV, time::Number) + return field.data(time) +end + +### Interpolation + +""" Interpolate time-invariant field in time direction. """ +function Base.call(field::DVTI, time::Float64) + return field +end +function Base.call(field::DCTI, time::Float64) + return field +end +function Base.call(field::CVTI, time::Float64) + return field.data() +end +function Base.call(field::CCTI, time::Float64) + return field.data() +end + +""" Interpolate time-variant field in time direction. """ +function Base.call(field::DCTV, time::Float64) + for i=reverse(1:length(field)) + if isapprox(field[i].time, time) + return DCTI(field[i].data) end end - field = DefaultDiscreteField(timesteps) - return field + info(field.data) + info(time) + error("interpolate DCTV: not implemented yet") end -function Base.convert(::Type{DefaultDiscreteField}, data::Vector{TimeStep}) - field = DefaultDiscreteField(data) -# @debug(field) - return field +function Base.call(field::DVTV, time::Float64, time_extrapolation::Symbol=:linear) + for i=reverse(1:length(field)) + if isapprox(field[i].time, time) + return DVTI(field[i].data) + end + end + info(field.data) + info(time) + error("interpolate DVTV: not implemented yet") end -### CONTINUOUS FIELDS ### - -type DefaultContinuousField <: ContinuousField - field :: Function +""" Interpolate constant field in spatial dimension. """ +function Base.call(basis::CVTI, field::DCTI, xi::Vector) + return field.data end -function Base.call(field::DefaultContinuousField, xi::Vector, time::Number) - return field.field(xi, time) +""" Interpolate variable field in spatial dimension. """ +function Base.call(basis::CVTI, values::DVTI, xi::Vector) + N = basis(xi) + return sum([N[i]*values[i] for i=1:length(N)]) end -function Base.convert(::Type{DefaultContinuousField}, f::Function) - return DefaultContinuousField(f) +function Base.call(basis::CVTI, geometry::DVTI, xi::Vector, ::Type{Val{:grad}}) + dbasis = basis(xi, Val{:grad}) + J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)]) + invJ = isa(J, Vector) ? inv(J[1]) : inv(J) + grad = invJ * dbasis + return grad end +function Base.call(basis::CVTI, geometry::DVTI, values::DVTI, xi::Vector, ::Type{Val{:grad}}) + grad = call(basis, geometry, xi, Val{:grad}) + gradf = sum([grad[:,i]*values[i]' for i=1:length(geometry)])' + return length(gradf) == 1 ? gradf[1] : gradf +end + +function Base.call(basis::CVTI, xi::Vector, time::Number) + call(basis, xi) +end ### FIELDSET ### typealias FieldSet Dict{ASCIIString, Field} -# 1. given numbers, arrays or tuples -> discrete field - -function Base.convert(::Type{Field}, data::Union{Number, Array, Tuple}...) - return DiscreteField(data...) -end - -function Base.convert(::Type{DiscreteField}, data::Union{Number, Array, Tuple}...) - return convert(DefaultDiscreteField, data...) -end - -# 2. given function -> continuous field - -function Base.convert(::Type{Field}, data::Function) - return ContinuousField(data) -end - -function Base.convert(::Type{ContinuousField}, data::Function) - return convert(DefaultContinuousField, data) -end - -function Base.length(::Field) - return 1 -end - diff --git a/src/fields2.jl b/src/fields2.jl index 58b53f7..ff01a75 100644 --- a/src/fields2.jl +++ b/src/fields2.jl @@ -3,310 +3,4 @@ # https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb -abstract AbstractField - -abstract Discrete <: AbstractField -abstract Continuous <: AbstractField -abstract Constant <: AbstractField -abstract Variable <: AbstractField -abstract TimeVariant <: AbstractField -abstract TimeInvariant <: AbstractField - - -type Field{A<:Union{Discrete,Continuous}, B<:Union{Constant,Variable}, C<:Union{TimeVariant,TimeInvariant}} - data -end - -### Basic data structure for discrete field - -type Increment{T} - time :: Float64 - data :: T -end - -function Base.convert{T}(::Type{Increment{T}}, data::Pair{Float64,T}) - return Increment{T}(data[1], data[2]) -end - -function Base.convert{T}(::Type{Increment{Vector{Vector{T}}}}, data::Pair{Float64, Matrix{T}}) - time = data[1] - content = data[2] - return Increment(time, Vector{T}[content[:,i] for i=1:size(content,2)]) -end - -function Base.getindex{T}(increment::Increment{Vector{T}}, i::Int64) - return increment.data[i] -end - -function Base.(:*)(d, increment::Increment) - return d*increment.data -end - -### Basic data structure for continuous field - -type Basis - basis :: Function - dbasis :: Function -end - -function Base.call(basis::Basis, xi::Vector) - basis.basis(xi) -end - -function Base.call(basis::Basis, xi::Vector, ::Type{Val{:grad}}) - basis.dbasis(xi) -end - -### Different field combinations and other typealiases - -typealias DCTI Field{Discrete, Constant, TimeInvariant} -typealias DVTI Field{Discrete, Variable, TimeInvariant} -typealias DCTV Field{Discrete, Constant, TimeVariant} -typealias DVTV Field{Discrete, Variable, TimeVariant} -typealias CCTI Field{Continuous, Constant, TimeInvariant} -typealias CVTI Field{Continuous, Variable, TimeInvariant} # can be used to interpolate in spatial dimension -typealias CCTV Field{Continuous, Constant, TimeVariant} # can be used to interpolate in time -typealias CVTV Field{Continuous, Variable, TimeVariant} - -typealias ScalarIncrement{T} Increment{T} -typealias VectorIncrement{T} Increment{Vector{T}} -typealias TensorIncrement{T} Increment{Matrix{T}} - -typealias DiscreteField Union{DCTI, DVTI, DCTV, DVTV} -typealias ContinuousField Union{CCTI, CVTI, CCTV, CVTV} -typealias ConstantField Union{DCTI, DCTV, CCTI, CCTV} -typealias VariableField Union{DVTI, DVTV, CVTI, CVTV} -typealias TimeInvariantField Union{DCTI, DVTI, CCTI, CVTI} -typealias TimeVariantField Union{DCTV, DVTV, CCTV, CVTV} - - -### Convenient functions to create fields - -#function Base.convert(::Type{Field}, data) -# return Field(data) -#end - -function Field(data) - return DCTI(data) -end - -function Field(data::Vector) - return DVTI(data) -end - -function Field{T}(data::Pair{Float64, T}...) - return DCTV([Increment{T}(d[1], d[2]) for d in data]) -end - -function Field{T}(data::Pair{Float64, Vector{T}}...) - return DVTV([Increment{Vector{T}}(d[1], d[2]) for d in data]) -end - -function Base.convert{T}(::Type{DCTV}, data::Pair{Float64, Vector{T}}...) - return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data]) -end - -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, Number}) - return CVTV(func) - else - error("no proper definition found for function: check methods.") - end -end - -function CVTI(basis::Function, dbasis::Function) - return CVTI(Basis(basis, dbasis)) -end - -function Field(basis::Function, dbasis::Function) - return CVTI(basis, dbasis) -end - -### Accessing and manipulating discrete fields - -function Base.getindex(field::DVTV, i::Int64) - return field.data[i] -end - -function Base.push!(field::DCTV, data::Pair) - push!(field.data, data) -end - -function Base.push!(field::DVTV, data::Pair) -# info("field.data = \n$(field.data)") -# info("data = \n$data") - push!(field.data, data) -end - -function Base.getindex(field::DVTV, i::Int64) - return field.data[i] -end - -function Base.getindex(field::DVTI, i::Int64) - return field.data[i] -end - -function Base.getindex(field::DCTV, i::Int64) - return field.data[i] -end - -function Base.getindex(field::Field, i::Int64) - return field.data[i] -end - -function Base.length(field::DVTI) - return length(field.data) -end - -function Base.length(field::DCTI) - return 1 -end - -function Base.length(field::DVTV) - return length(field.data) -end - -function Base.length(field::DCTV) - return length(field.data) -end - -for op = (:+, :*, :/, :-) - @eval ($op)(increment::Increment, field::DCTI) = ($op)(increment.data, field.data) - @eval ($op)(field::DCTI, increment::Increment) = ($op)(increment.data, field.data) - @eval ($op)(field1::DCTI, field2::DCTI) = ($op)(field1.data, field2.data) - @eval ($op)(field::DCTI, k) = ($op)(field.data, k) - @eval ($op)(k, field::DCTI) = ($op)(field.data, k) -end - -function Base.vec(field::DVTI) - return [field.data...;] -end - -function Base.vec(field::DCTV) - info("trying to vectorize $field") - error("does not make sense") -end - -function Base.endof(field::Field) - return endof(field.data) -end - -#function Base.similar{T}(field::DVTI, data::Vector{T}) -# return Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment))) -#end - -function Base.similar{T}(field::DVTI, data::Vector{T}) - n = length(field.data) - data = reshape(data, round(Int, length(data)/n), n) - newdata = Vector[data[:,i] for i=1:n] - return typeof(field)(newdata) -end - -function Base.start(::DVTI) - return 1 -end - -function Base.next(f::DVTI, state) - return f.data[state], state+1 -end - -function Base.done(f::DVTI, s) - return s > length(f.data) -end - -### Accessing continuous fields - -function Base.call(field::CVTI, xi::Vector) - field.data(xi) -end - -function Base.call(field::CVTI, xi::Vector, ::Type{Val{:grad}}) - field.data(xi, Val{:grad}) -end - -function Base.convert(::Type{Basis}, field::CVTI) - return field.data -end - -function Base.call(field::CCTV, time::Number) - return field.data(time) -end - -### Interpolation - -""" Interpolate time-invariant field in time direction. """ -function Base.call(field::DVTI, time::Float64) - return field -end -function Base.call(field::DCTI, time::Float64) - return field -end -function Base.call(field::CVTI, time::Float64) - return field.data() -end -function Base.call(field::CCTI, time::Float64) - return field.data() -end - -""" Interpolate time-variant field in time direction. """ -function Base.call(field::DCTV, time::Float64) - for i=reverse(1:length(field)) - if isapprox(field[i].time, time) - return DCTI(field[i].data) - end - end - info(field.data) - info(time) - error("interpolate DCTV: not implemented yet") -end - -function Base.call(field::DVTV, time::Float64, time_extrapolation::Symbol=:linear) - for i=reverse(1:length(field)) - if isapprox(field[i].time, time) - return DVTI(field[i].data) - end - end - info(field.data) - info(time) - error("interpolate DVTV: not implemented yet") -end - -""" Interpolate constant field in spatial dimension. """ -function Base.call(basis::CVTI, field::DCTI, xi::Vector) - return field.data -end - -""" Interpolate variable field in spatial dimension. """ -function Base.call(basis::CVTI, values::DVTI, xi::Vector) - N = basis(xi) - return sum([N[i]*values[i] for i=1:length(N)]) -end - -function Base.call(basis::CVTI, geometry::DVTI, xi::Vector, ::Type{Val{:grad}}) - dbasis = basis(xi, Val{:grad}) - J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)]) - invJ = isa(J, Vector) ? inv(J[1]) : inv(J) - grad = invJ * dbasis - return grad -end - -function Base.call(basis::CVTI, geometry::DVTI, values::DVTI, xi::Vector, ::Type{Val{:grad}}) - grad = call(basis, geometry, xi, Val{:grad}) - gradf = sum([grad[:,i]*values[i]' for i=1:length(geometry)])' - return length(gradf) == 1 ? gradf[1] : gradf -end - -function Base.call(basis::CVTI, xi::Vector, time::Number) - call(basis, xi) -end - -### FIELDSET ### - -typealias FieldSet Dict{ASCIIString, Field} diff --git a/src/heat.jl b/src/heat.jl index 70caa93..b72cc73 100644 --- a/src/heat.jl +++ b/src/heat.jl @@ -3,14 +3,20 @@ # Heat problems -abstract HeatProblem <: Problem -abstract HeatEquation <: Equation +abstract HeatProblem <: AbstractProblem -function get_unknown_field_name(equation::HeatEquation) +function HeatProblem(dim::Int=1, elements=[]) + return Problem{HeatProblem}(dim, elements) +end + +function get_unknown_field_name{P<:HeatProblem}(::Type{P}) return "temperature" end -### Formulation ### +function get_unknown_field_type{P<:HeatProblem}(::Type{P}) + # scalar field + return Float64 +end """ Heat equations. @@ -36,83 +42,29 @@ References https://en.wikipedia.org/wiki/Heat_equation """ -function assemble!(assembly::Assembly, equation::HeatEquation, time::Number=0.0, problem=nothing) +function assemble!{E<:CG}(assembly::Assembly, problem::Problem{HeatProblem}, element::Element{E}, time::Number) - element = get_element(equation) - gdofs = get_gdofs(equation) - basis = get_basis(element) - dbasis = grad(basis) - detJ = det(basis) - for ip in get_integration_points(equation) - w = ip.weight*detJ(ip) - N = basis(ip, time) + gdofs = get_gdofs(element, problem.dim) + for ip in get_integration_points(element) + w = ip.weight*det(element, ip, time) + N = element(ip, time) if haskey(element, "density") - rho = basis("density", ip, time) + rho = element("density", ip, time) add!(assembly.mass_matrix, gdofs, gdofs, w*rho*N'*N) end if haskey(element, "temperature thermal conductivity") - dN = dbasis(ip, time) - k = basis("temperature thermal conductivity", ip, time) + dN = element(ip, time, Val{:grad}) + k = element("temperature thermal conductivity", ip, time) add!(assembly.stiffness_matrix, gdofs, gdofs, w*k*dN'*dN) end if haskey(element, "temperature load") - f = basis("temperature load", ip, time) + f = element("temperature load", ip, time) add!(assembly.force_vector, gdofs, w*N'*f) end if haskey(element, "temperature flux") - info("assemble boundary flux") - g = basis("temperature flux", ip, time) + g = element("temperature flux", ip, time) add!(assembly.force_vector, gdofs, w*N'*g) end end end - -### Equations ### - -""" Diffusive heat transfer for 4-node bilinear element. """ -type DC2D4 <: HeatEquation - element :: Quad4 - integration_points :: Vector{IntegrationPoint} -end - -function Base.size(equation::DC2D4) - return (1, 4) -end - -""" Diffusive heat transfer for 2-node linear segment. """ -type DC2D2 <: HeatEquation - element :: Seg2 - integration_points :: Vector{IntegrationPoint} -end - -function Base.size(equation::DC2D2) - return (1, 2) -end - -# Conversions element -> equation - -function Base.convert(::Type{HeatEquation}, element::Quad4) - integration_points = get_integration_points(element) - haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(4)) - DC2D4(element, integration_points) -end - -function Base.convert(::Type{HeatEquation}, element::Seg2) - integration_points = get_integration_points(element) - haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(2)) - DC2D2(element, integration_points) -end - -### Problems ### - -type PlaneHeatProblem <: HeatProblem - unknown_field_name :: ASCIIString - unknown_field_dimension :: Int - equations :: Vector{HeatEquation} -end - -""" Default constructor for problem takes no arguments. """ -function PlaneHeatProblem(equations=[]) - return PlaneHeatProblem("temperature", 1, equations) -end diff --git a/src/integrate.jl b/src/integrate.jl index 97fed45..4dfcd0c 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -3,7 +3,7 @@ # Let's drop here all integration schemes and some defaults for different element types -function get_integration_points(Quad4::Element) +function get_integration_points(::Type{Quad4}) [ IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0), IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0), @@ -14,13 +14,13 @@ end typealias LineElement Union{Seg2, Seg3} -function get_integration_points(element::LineElement, ::Type{Val{1}}) +function get_integration_points(::Type{Seg2}, ::Type{Val{1}}) [ IntegrationPoint([0.0], 2.0) ] end -function get_integration_points(element::LineElement, ::Type{Val{2}}) +function get_integration_points(::Type{Seg2}, ::Type{Val{2}}) [ IntegrationPoint([-sqrt(1/3)], 1) IntegrationPoint([+sqrt(1/3)], 1) @@ -54,12 +54,12 @@ function get_integration_points(element::LineElement, ::Type{Val{5}}) ] end -function get_integration_points(element::Seg2) - return get_integration_points(element, Val{1}) +function get_integration_points(::Type{Seg2}) + return get_integration_points(Seg2, Val{2}) end function get_integration_points(element::Seg3) - return get_integration_points(element, Val{2}) + return get_integration_points(element, Val{3}) end ### 3D elements diff --git a/src/lagrange.jl b/src/lagrange.jl index c97d630..3dd2161 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -3,7 +3,7 @@ # Lagrange (Continous Galerkin) finite elements -abstract CG <: Element +abstract CG <: AbstractElement """ Given polynomial P and coordinates of reference element, calculate @@ -29,29 +29,26 @@ Examples >>> @create_lagrange_element(Seg2, "2 node linear segment", X, P) """ macro create_lagrange_element(element_name, element_description, X, P) -# Logging.debug("Creating element ", element_name, ": ", element_description, "\n") eltype = esc(element_name) quote - global get_element_description - basis, dbasisdxi = calculate_lagrange_basis($P, $X) - type $eltype <: CG - connectivity :: Array{Int, 1} - basis :: CVTI - fields :: FieldSet + global get_basis, get_dbasis + basis, dbasis = calculate_lagrange_basis($P, $X) + abstract $eltype <: CG + function get_basis(::Type{$eltype}, xi::Vector{Float64}) + return basis(xi) end - function $eltype(connectivity, args...) - $eltype(connectivity, CVTI(basis, dbasisdxi), FieldSet()) + function get_dbasis(::Type{$eltype}, xi::Vector{Float64}) + return dbasis(xi) + end + function $eltype(args...) + return Element{$eltype}(args...) + end + function Base.size(::Type{$eltype}) + return Base.size($X) end - get_element_description(el::Type{$eltype}) = $element_description - Base.size(element::Type{$eltype}) = Base.size($X) - Base.size(element::$eltype) = Base.size($X) end end -# 0d Lagrange element - -#@create_element(Point1, CG, "1 node point element") - # 1d Lagrange elements @create_lagrange_element(Seg2, "2 node linear line element", diff --git a/src/mortar_elements.jl b/src/mortar_elements.jl index d10643c..4de6950 100644 --- a/src/mortar_elements.jl +++ b/src/mortar_elements.jl @@ -6,7 +6,7 @@ # slave element = non-mortar element where integration happens # master element = mortar element projected to non-mortar side -abstract MortarElement <: Element +abstract MortarElement <: AbstractElement type MSeg2 <: MortarElement connectivity :: Vector{Int} diff --git a/src/problems.jl b/src/problems.jl index 6d2f4a7..033e7bc 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -1,77 +1,36 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -abstract Problem -abstract BoundaryProblem <: Problem -abstract FieldProblem <: Problem +abstract AbstractProblem -""" Return all equations beloging to this problem. """ -function get_equations(problem::Problem) - problem.equations +type Problem{T<:AbstractProblem} + dim :: Int + elements :: Vector{Element} +end + +type BoundaryProblem{T<:AbstractProblem} + parent_field_name :: ASCIIString + parent_field_dim :: Int + dim :: Int + elements :: Vector{Element} +end + +typealias AllProblems Union{Problem, BoundaryProblem} + +function get_elements(problem::AllProblems) + return problem.elements end """ Return the dimension of the unknown field of this problem. """ function get_unknown_field_dimension(problem::Problem) - problem.unknown_field_dimension + return problem.dim end """ Return the name of the unknown field of this problem. """ -function get_unknown_field_name(problem::Problem) - problem.unknown_field_name +function get_unknown_field_name{P<:AbstractProblem}(problem::Problem{P}) + return get_unknown_field_name(P) end -""" -Add new equation to problem. - -Parameters ----------- -problem -element - -Notes ------ -Equation is automatically created during process based on problem -element -> equation mapping and element type. -""" -function Base.push!(problem::Problem, element::Element, args...) -# element_type = typeof(element) -# equation_type = problem.element_mapping[element_type] -# push!(problem.equations, equation_type(element, args...)) - push!(problem.equations, element) +function Base.push!(problem::AllProblems, element::Element) + push!(problem.elements, element) end - -""" -Return the size of the problem, i.e. -maximum number of connectivity × unknown field dimension. -""" -function Base.size(problem::Problem) - mc = 0 - for equation in get_equations(problem) - element = get_element(equation) - mc = max(mc, get_connectivity(element)...) - end - dim = get_unknown_field_dimension(problem) - return (dim, dim*mc) -end - -""" Assign new equation mapping to problem for some element. - -Examples --------- ->>> p = PlaneHeatProblem() ->>> p[Seg2] = DC2D2 - -""" -function Base.setindex!(problem::Problem, equation, element) - problem.element_mapping[element] = equation -end - -""" Return global degrees of freedom of element in matrix level. """ -function get_gdofs(problem::Problem, equation::Equation) - dim = get_unknown_field_dimension(problem) - element = get_element(equation) - conn = get_connectivity(element) - gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...)) - return gdofs -end - diff --git a/src/solvers.jl b/src/solvers.jl index d1992ef..0b06e05 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -5,39 +5,6 @@ abstract Solver -""" -Solve field equations for single element with some dofs fixed. This can be used -to test nonlinear element formulations. -""" -function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false, callback=nothing) - unknown_field_name = get_unknown_field_name(equation) - element = get_element(equation) - x0 = element[unknown_field_name](0.0) - x = zeros(prod(size(equation))) - dx = fill!(similar(x), 0.0) - ass = Assembly() - for i=1:max_iterations - empty!(ass) - assemble!(ass, equation) - A = full(ass.stiffness_matrix)[free_dofs, free_dofs] - b = full(ass.force_vector)[free_dofs] - if dump_matrices - dump(full(A)) - dump(full(b)') - end - dx[free_dofs] = A \ b - x += dx - eqsize = size(equation) - data = eqsize[1] != 1 ? reshape(x, eqsize) : x - push!(element[unknown_field_name], time => data) - norm(dx) < tolerance && return - if !isa(callback, Void) - callback(x) - end - end - error("Did not converge in $max_iterations iterations") -end - """ Solve field equations for a single problem with some dofs fixed. This can be used to test nonlinear element formulations. Dirichlet boundary is assumed to be homogeneous @@ -72,13 +39,12 @@ function solve!(problem::Problem, free_dofs::Vector{Int}, time::Float64; max_ite if !(isa(callback, Void)) callback(x) end - for equation in get_equations(problem) - element = get_element(equation) - gdofs = get_gdofs(equation) + for element in get_elements(problem) + gdofs = get_gdofs(element, problem.dim) data = full(x[gdofs]) - eqsize = size(equation) - if eqsize[1] != 1 - data = reshape(data, eqsize) + if length(data) != length(element) + data = reshape(data, problem.dim, length(element)) + data = [data[:,i] for i=1:size(data,2)] end push!(element[field_name], time => data) end @@ -87,24 +53,11 @@ function solve!(problem::Problem, free_dofs::Vector{Int}, time::Float64; max_ite error("Did not converge in $max_iterations iterations") end -function Base.push!(solver::Solver, problem::Problem) - push!(solver.problems, problem) -end -""" Get all problems assigned to solver. """ -function get_problems(solver::Solver) - return solver.problems -end - -## SimpleSolver -- tiny direct demo solver -""" Simple solver for educational purposes. """ -type SimpleSolver <: Solver - problems :: Vector{Problem} -end - -""" Default initializer. """ -function SimpleSolver() - SimpleSolver([]) +""" Simple linear solver for educational purposes. """ +type LinearSolver <: Solver + field_problem :: Problem + boundary_problem :: BoundaryProblem end """ @@ -113,28 +66,30 @@ Call solver to solve a set of problems. This is a simple direct solver for demonstration purposes. It handles the common situation, i.e., some main field problem and it's Dirichlet boundary. - Au + C'λ = f + Ku + C'λ = f Cu = g """ -function call(solver::SimpleSolver, time::Number=0.0) - problem1, problem2 = get_problems(solver) +function call(solver::LinearSolver, time::Float64) - assembly1 = Assembly() - assemble!(assembly1, problem1, time) - assembly2 = Assembly() - assemble!(assembly2, problem2, time) + field_name = get_unknown_field_name(solver.field_problem) + field_dim = get_unknown_field_dimension(solver.field_problem) + info("solving $field_name problem, $field_dim dofs / nodes") + + field_assembly = assemble(solver.field_problem, time) + boundary_assembly = assemble(solver.boundary_problem, time) + + info("Creating sparse matrices") + K = sparse(field_assembly.stiffness_matrix) + dim = size(K, 1) + f = sparse(field_assembly.force_vector, dim, 1) + + C = sparse(boundary_assembly.stiffness_matrix, dim, dim) + g = sparse(boundary_assembly.force_vector, dim, 1) -# info("Creating sparse matrices") - A1 = sparse(assembly1.stiffness_matrix) - dims = size(A1) - b1 = sparse(assembly1.force_vector, dims[1], 1) - A2 = sparse(assembly2.stiffness_matrix, dims[1], dims[2]) - b2 = sparse(assembly2.force_vector, dims[1], 1) - # create a saddle point problem - A = [A1 A2; A2' zeros(A2)] - b = [b1; b2] + A = [K C'; C' zeros(C)] + b = [f; g] # solve problem nz = unique(rowvals(A)) # take only non-zero rows @@ -142,37 +97,26 @@ function call(solver::SimpleSolver, time::Number=0.0) x[nz] = lufact(A[nz,nz]) \ full(b[nz]) # get "problem-wise" solution vectors - x1 = x[1:length(b1)] - x2 = x[length(b1)+1:end] + u = x[1:dim] + la = x[dim+1:end] # update field for elements in problem 1 - for equation in get_equations(problem1) - element = get_element(equation) - field_name = get_unknown_field_name(problem1) - gdofs = get_gdofs(problem1, equation) - local_sol = vec(full(x1[gdofs])) - eqsize = size(equation) - if eqsize[1] != 1 - local_sol = reshape(local_sol, eqsize) + for element in get_elements(solver.field_problem) + gdofs = get_gdofs(element, field_dim) + local_sol = vec(full(u[gdofs])) + # if solving vector field, modify local solution vector + # to array of vectors + if field_dim != 1 + local_sol = reshape(local_sol, field_dim, length(element)) + local_sol = [local_sol[:,i] for i=1:size(local_sol,2)] + end + if haskey(element, field_name) + push!(element[field_name], time => local_sol) + else + element[field_name] = (time => local_sol) end - #info("problem1: pushing to $field_name") - push!(element[field_name], time => local_sol) end - # update field for elements in problem 2 (Dirichlet boundary) - for equation in get_equations(problem2) - element = get_element(equation) - field_name = "reaction force" #get_unknown_field_name(problem2) - gdofs = get_gdofs(problem2, equation) - local_sol = vec(full(x1[gdofs])) - eqsize = size(equation) - if eqsize[1] != 1 - local_sol = reshape(local_sol, eqsize) - end - #info("problem2: pushing to $field_name") - #push!(element[field_name], time => local_sol) - end - - return norm(x1) + return norm(u) end diff --git a/test/test_directsolver.jl b/test/test_directsolver.jl new file mode 100644 index 0000000..5c6b160 --- /dev/null +++ b/test/test_directsolver.jl @@ -0,0 +1,132 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +module DirectSolverTests + +using JuliaFEM.Test +using JuliaFEM + +using JuliaFEM: Seg2, Quad4 +using JuliaFEM: PlaneStressElasticityProblem, DirichletProblem +using JuliaFEM: DirectSolver + +function test_solver_multiple_dirichlet_bc() + N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]] + + e1 = Quad4([1, 2, 4, 3]) + e1["geometry"] = Vector[N[1], N[2], N[4], N[3]] + e1["youngs modulus"] = 900.0 + e1["poissons ratio"] = 0.25 + b1 = Seg2([3, 4]) + b1["geometry"] = Vector[N[3], N[4]] + b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]] + + problem = PlaneStressElasticityProblem() + push!(problem, e1) + push!(problem, b1) + + # manually solve problem 1 + # free_dofs = [3, 5, 6, 8] + # free_dofs = [3, 6, 7, 8] + #solve!(problem, free_dofs, 0.0; max_iterations=10) + #disp = e1("displacement", [1.0, 1.0], 0.0) + #info("displacement at tip: $disp") + #@test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01]) + + # boundary elements for dirichlet dx=0 + dx = Seg2([1, 3]) + dx["geometry"] = Vector[N[1], N[3]] + dx["displacement 1"] = 0.0 + + # boundary elements for dirichlet dy=0 + dy = Seg2([1, 2]) + dy["geometry"] = Vector[N[1], N[2]] + dy["displacement 2"] = 0.0 + + problem2 = DirichletProblem("displacement", 2) + push!(problem2, dx) + + problem3 = DirichletProblem("displacement", 2) + push!(problem3, dy) + + solver = DirectSolver() + push!(solver, problem) + push!(solver, problem2) + push!(solver, problem3) + + # launch solver + norm = solver(0.0) + + disp = e1("displacement", [1.0, 1.0], 0.0) + info("displacement at tip: $disp") + @test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01]) + +end + +function test_solver_multiple_bodies_multiple_dirichlet_bc() + N = Vector[ + [0.0, 0.0], [1.0, 0.0], + [0.0, 1.0], [1.0, 1.0], + [0.0, 2.0], [1.0, 2.0]] + + e1 = Quad4([1, 2, 4, 3]) + e1["geometry"] = Vector[N[1], N[2], N[4], N[3]] + e2 = Quad4([3, 4, 6, 5]) + e2["geometry"] = Vector[N[3], N[4], N[6], N[5]] + for el in [e1, e2] + el["youngs modulus"] = 900.0 + el["poissons ratio"] = 0.25 + end + b1 = Seg2([5, 6]) + b1["geometry"] = Vector[N[5], N[6]] + b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]] + + body1 = PlaneStressElasticityProblem() + push!(body1, e1) + + body2 = PlaneStressElasticityProblem() + push!(body2, e2) + push!(body2, b1) + + # boundary elements for dirichlet dx=0 + dx1 = Seg2([1, 3]) + dx1["geometry"] = Vector[N[1], N[3]] + dx2 = Seg2([3, 5]) + dx2["geometry"] = Vector[N[3], N[5]] + for dx in [dx1, dx2] + dx["displacement 1"] = 0.0 + end + + boundary1 = DirichletProblem("displacement", 2) + push!(boundary1, dx1) + push!(boundary1, dx2) + + # boundary elements for dirichlet dy=0 + dy1 = Seg2([1, 2]) + dy1["geometry"] = Vector[N[1], N[2]] + dy1["displacement 2"] = 0.0 + + boundary2 = DirichletProblem("displacement", 2) + push!(boundary2, dy1) + + + solver = DirectSolver() + push!(solver, body1) + push!(solver, body2) + push!(solver, boundary1) + push!(solver, boundary2) + + # launch solver + norm = solver(0.0) + + disp = e2("displacement", [1.0, 1.0], 0.0) + info("displacement at tip: $disp") + # code aster verification, two_elements.comm + @test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01]) + +end + +# test_solver_multiple_bodies_multiple_dirichlet_bc() + + +end diff --git a/test/test_dirichlet.jl b/test/test_dirichlet.jl index ab448ea..faba641 100644 --- a/test/test_dirichlet.jl +++ b/test/test_dirichlet.jl @@ -5,7 +5,7 @@ module TestDirichletBoundaryCondition using JuliaFEM.Test using JuliaFEM -using JuliaFEM: Seg2, DirichletProblem, Assembly, assemble! +using JuliaFEM: Seg2, DirichletProblem, Assembly, assemble function test_dirichlet_problem_1_dim() element = Seg2([1, 2]) @@ -13,8 +13,7 @@ function test_dirichlet_problem_1_dim() element["temperature"] = 0.0 problem = DirichletProblem("temperature", 1) push!(problem, element) - assembly = Assembly() - assemble!(assembly, problem) + assembly = assemble(problem, 0.0) A = full(assembly.stiffness_matrix) b = full(assembly.force_vector) @test isapprox(A, 1/6*[2 1; 1 2]) @@ -27,8 +26,7 @@ function test_dirichlet_problem_2_dim() element["displacement"] = 0.0 problem = DirichletProblem("displacement", 2) push!(problem, element) - assembly = Assembly() - assemble!(assembly, problem) + assembly = assemble(problem, 0.0) A = full(assembly.stiffness_matrix) b = full(assembly.force_vector) A_expected = 1/6*[2 0 1 0; 0 2 0 1; 1 0 2 0; 0 1 0 2] @@ -42,8 +40,7 @@ function test_dirichlet_problem_2_dim_single_dof_fixed() element["displacement 2"] = 0.0 problem = DirichletProblem("displacement", 2) push!(problem, element) - assembly = Assembly() - assemble!(assembly, problem) + assembly = assemble(problem, 0.0) A = full(assembly.stiffness_matrix) b = full(assembly.force_vector) info(b) diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl index 3d6a837..cad3d69 100644 --- a/test/test_elasticity.jl +++ b/test/test_elasticity.jl @@ -4,9 +4,7 @@ module ElasticityTests using JuliaFEM.Test -using JuliaFEM: Seg2, Quad4, Field, FieldSet, CPS4, - get_basis, solve!, - PlaneStressElasticityProblem +using JuliaFEM: Seg2, Quad4, PlaneStressElasticityProblem, solve! function test_elasticity_volume_load() element = Quad4([1, 2, 3, 4]) @@ -14,11 +12,13 @@ function test_elasticity_volume_load() element["youngs modulus"] = 500.0 element["poissons ratio"] = 0.3 element["displacement load"] = Vector[[0.0, -10.0], [0.0, -10.0], [0.0, -10.0], [0.0, -10.0]] - free_dofs = [3, 4, 5, 6] + element["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]) problem = PlaneStressElasticityProblem() push!(problem, element) + + free_dofs = [3, 4, 5, 6] solve!(problem, free_dofs, 0.0; max_iterations=10) - disp = get_basis(element)("displacement", [1.0, 1.0], 0.0) + disp = element("displacement", [1.0, 1.0], 0.0) info("displacement at tip: $disp") # verified using Code Aster. @test isapprox(disp[2], -8.77303119819776) @@ -31,9 +31,12 @@ function test_elasticity_surface_load() element1["geometry"] = Vector[N[1], N[2], N[4], N[3]] element1["youngs modulus"] = 900.0 element1["poissons ratio"] = 0.25 + element1["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]) + element2 = Seg2([3, 4]) element2["geometry"] = Vector[N[3], N[4]] element2["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]] + element2["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]]) #free_dofs = [3, 5, 6, 8] free_dofs = [3, 6, 7, 8] @@ -41,8 +44,6 @@ function test_elasticity_surface_load() push!(problem, element1) push!(problem, element2) solve!(problem, free_dofs, 0.0; max_iterations=10) - #disp = get_basis(element1)("displacement", [1.0, 1.0], 1.0)[2] - info(last(element1["displacement"])) disp = element1("displacement", [1.0, 1.0], 0.0) info("displacement at tip: $disp") # verified using Code Aster. diff --git a/test/test_heat.jl b/test/test_heat.jl index 8dd4408..61a60c7 100644 --- a/test/test_heat.jl +++ b/test/test_heat.jl @@ -7,8 +7,7 @@ module HeatTests # always wrap tests to module ending with "Tests" using JuliaFEM.Test # always use JuliaFEM.Test, not Base.Test -using JuliaFEM: HeatEquation -using JuliaFEM: Seg2, Quad4, DC2D4, DC2D2, Assembly, assemble! +using JuliaFEM: Seg2, Quad4, HeatProblem, assemble function test_one_element() # always start test function with name test_ @@ -26,11 +25,13 @@ function test_one_element() # always start test function with name test_ # linear ramp from 0 to 6 in time 0 to 1 boundary_element["temperature flux"] = (0.0 => 0.0, 1.0 => 6.0) + problem = HeatProblem() + push!(problem, element) + push!(problem, boundary_element) + # Set constant source f=12 with k=6. Accurate solution is # T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x) - equation = convert(HeatEquation, element) - assembly = Assembly() - assemble!(assembly, equation) + assembly = assemble(problem, 0.0) fdofs = [1, 2] A = full(assembly.stiffness_matrix) b = full(assembly.force_vector) @@ -38,18 +39,15 @@ function test_one_element() # always start test function with name test_ # Set constant flux g=6 on boundary. Accurate solution is # u(x,y) = x which equals T=1 on boundary. - boundary_equation = convert(HeatEquation, boundary_element) - empty!(assembly) - - time = 1.0 - assemble!(assembly, equation, time) - assemble!(assembly, boundary_equation, time) + # at time t=1.0 all loads should be on. + assembly = assemble(problem, 1.0) A = full(assembly.stiffness_matrix) b = full(assembly.force_vector) T = A[fdofs, fdofs] \ b[fdofs] info("T = $T") - @test isapprox(T, [2.0, 2.0]) # always use @test to test things. + @test isapprox(T, [2.0, 2.0]) end end + diff --git a/test/test_potential_energy.jl b/test/test_potential_energy.jl index 3890767..8144ff0 100644 --- a/test/test_potential_energy.jl +++ b/test/test_potential_energy.jl @@ -3,73 +3,43 @@ module ElementTests +using JuliaFEM using JuliaFEM.Test -using JuliaFEM -using JuliaFEM: Equation, Quad4, IntegrationPoint, Assembly, assemble!, - get_element, get_basis, grad, get_unknown_field_name, - PlaneHeatProblem, Seg2, Problem, solve!, - get_default_integration_points, Equation +using JuliaFEM: AbstractProblem, Problem +using JuliaFEM: Element, Seg2, Quad4 +using JuliaFEM: IntegrationPoint, solve! -abstract MyEquation <: Equation +abstract HeatProblem <: AbstractProblem -function JuliaFEM.get_unknown_field_name(equation::MyEquation) +function HeatProblem(dim::Int=1, elements=[]) + return Problem{HeatProblem}(dim, elements) +end + +function JuliaFEM.get_unknown_field_name{P<:HeatProblem}(::Type{P}) return "temperature" end -""" Diffusive heat transfer for 4-node bilinear element, with a nonlinear source term. """ -type DC2D4NL <: MyEquation - element :: Quad4 - integration_points :: Vector{IntegrationPoint} +function JuliaFEM.get_unknown_field_type{P<:HeatProblem}(::Type{P}) + return Float64 end -function Base.size(equation::DC2D4NL) - return (1, 4) -end - -""" Nonlinear flux term. """ -type DC2D2NL <: MyEquation - element :: Seg2 - integration_points :: Vector{IntegrationPoint} -end - -function Base.size(equation::DC2D2NL) - return (1, 2) -end - -function Base.convert(::Type{MyEquation}, element::Quad4) - integration_points = get_default_integration_points(element) - haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(4)) - DC2D4NL(element, integration_points) -end - -function Base.convert(::Type{MyEquation}, element::Seg2) - integration_points = JuliaFEM.line5() - haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(2)) - DC2D2NL(element, integration_points) -end - - """ Calculate a potential Π = Wint - Wext of system. """ -function JuliaFEM.get_potential_energy(equation::DC2D4NL, ip, time; variation=nothing) - element = get_element(equation) - basis = get_basis(element) - k = basis("temperature thermal conductivity", ip, time) - f = basis("temperature load", ip, time) - T = basis("temperature", ip, time, variation) - c = basis("temperature nonlinearity coefficient", ip, time) - gradT = grad(basis)("temperature", ip, time, variation) +function JuliaFEM.get_potential_energy(problem::Problem{HeatProblem}, element::Element{Quad4}, ip::IntegrationPoint, time::Number; variation=nothing) + k = element("temperature thermal conductivity", ip, time) + f = element("temperature load", ip, time) + T = element("temperature", ip, time, variation) + c = element("temperature nonlinearity coefficient", ip, time) + gradT = element("temperature", ip, time, Val{:grad}, variation) Wint = (k + c*T) * 1/2*vecdot(gradT, gradT) Wext = f*T return Wint - Wext end -function JuliaFEM.get_potential_energy(equation::DC2D2NL, ip, time; variation=nothing) - element = get_element(equation) - basis = get_basis(element) - T = basis("temperature", ip, time, variation)[1] - T_ext = basis("temperature external", ip, time)[1] - coeff = basis("temperature coefficient", ip, time)[1] +function JuliaFEM.get_potential_energy(problem::Problem{HeatProblem}, element::Element{Seg2}, ip::IntegrationPoint, time::Number; variation=nothing) + T = element("temperature", ip, time, variation)[1] + T_ext = element("temperature external", ip, time)[1] + coeff = element("temperature coefficient", ip, time)[1] q0 = coeff*(T_ext^4 - T^4) Wint = 0.0 Wext = q0*T @@ -86,28 +56,19 @@ function test_potential_energy_method() element["temperature load"] = [0.0, 0.0, 0.0, 0.0] element["temperature nodal load"] = [3.0, 3.0, 0.0, 0.0] element["temperature nonlinearity coefficient"] = 6.0 - equation = convert(MyEquation, element) + element["temperature"] = (0.0 => zeros(Float64, 4)) + problem = HeatProblem() + push!(problem, element) # create model -- end - solve!(equation, [1, 2], 0.0) - basis = get_basis(element) - temp = basis("temperature", [0.0, -1.0], 0.0) + solve!(problem, [1, 2], 0.0) + temp = element("temperature", [0.0, -1.0], 0.0) err = temp - 2/3 info("error: $err") @test isapprox(err, 0.0) end -type TestProblem <: Problem - unknown_field_name :: ASCIIString - unknown_field_dimension :: Int - equations :: Vector{MyEquation} -end - -function TestProblem(equations=[]) - TestProblem("temperature", 1, equations) -end - function test_potential_energy_method_2() # create model -- start @@ -117,20 +78,21 @@ function test_potential_energy_method_2() element1["temperature thermal conductivity"] = 6.0 element1["temperature load"] = [0.0, 0.0, 0.0, 0.0] element1["temperature nonlinearity coefficient"] = [0.0, 0.0, 0.0, 0.0] + element1["temperature"] = (0.0 => zeros(Float64, 4)) element2 = Seg2([1, 2]) element2["geometry"] = Vector[N[1], N[2]] element2["temperature coefficient"] = 3.0e-8 # ~ 5.7e-8 * 0.5 element2["temperature external"] = 100.0 + element2["temperature"] = (0.0 => zeros(Float64, 2)) # create model -- end - problem = TestProblem() + problem = HeatProblem() push!(problem, element1) push!(problem, element2) solve!(problem, [1, 2], 0.0) - basis = get_basis(element1) - temp = basis("temperature", [0.0, -1.0], 0.0) + temp = element1("temperature", [0.0, -1.0], 0.0) err = temp - 0.5 info("error: $err") @test isapprox(err, 0.0, atol=1.0e-6) diff --git a/test/test_solver.jl b/test/test_solver.jl index beeba9d..f19a68e 100644 --- a/test/test_solver.jl +++ b/test/test_solver.jl @@ -6,14 +6,11 @@ module SolverTests using JuliaFEM.Test using JuliaFEM -using JuliaFEM: DirichletProblem, Seg2, PlaneHeatProblem, Quad4, SimpleSolver, get_element, get_basis, MortarElement, MortarProblem, PlaneStressElasticityProblem, solve!, DirectSolver +using JuliaFEM: Seg2, Quad4 +using JuliaFEM: DirichletProblem, HeatProblem +using JuliaFEM: LinearSolver -""" Define Problem 1: - -- Field function: Laplace equation Δu=0 in Ω={u∈R²|(x,y)∈[0,1]×[0,1]} -- Neumann boundary on Γ₁={0<=x<=1, y=0}, ∂u/∂n=600 on Γ₁ -""" -function get_heatproblem() +function test_linearsolver() el1 = Quad4([1, 2, 3, 4]) el1["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]] el1["temperature thermal conductivity"] = 6.0 @@ -26,278 +23,31 @@ function get_heatproblem() (1.0 => 600.0) ) - problem1 = PlaneHeatProblem() - push!(problem1, el1) - push!(problem1, el2) - return problem1 -end + field_problem = HeatProblem() + push!(field_problem, el1) + push!(field_problem, el2) -""" Define Problem 2: - - Dirichlet boundary Γ₂={0<=x<=1, y=1}, u=0 on Γ₂ -""" -function get_boundaryproblem() el3 = Seg2([3, 4]) el3["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]] el3["temperature"] = 0.0 - problem2 = DirichletProblem("temperature", 1) - push!(problem2, el3) - return problem2 -end -function test_simplesolver() - info("construct heat problem") - problem1 = get_heatproblem() - info("construct boundary problem") - problem2 = get_boundaryproblem() + boundary_problem = DirichletProblem("temperature", 1) + + push!(boundary_problem, el3) + # Create a solver for a set of problems - info("create SimpleSolver with problems.") - solver = SimpleSolver() - push!(solver, problem1) - push!(solver, problem2) - info("solve!") + solver = LinearSolver(field_problem, boundary_problem) + # Solve problem at time t=1.0 and update fields - call(solver, 1.0) + solver(1.0) + # Postprocess. # Interpolate temperature field along boundary of Γ₁ at time t=1.0 xi = [0.0, -1.0] - el2 = get_element(problem1.equations[2]) - basis = get_basis(el2) - X = basis("geometry", xi, 1.0) - T = basis("temperature", xi, 1.0) + X = el2("geometry", xi, 1.0) + T = el2("temperature", xi, 1.0) info("Temperature at point X = $X is T = $T") @test isapprox(T, 100.0) end -#test_simplesolver() - -function atest_direct_solver() - - N = Dict{Int, Vector}( - 1 => [0.0, 0.0], - 2 => [2.0, 0.0], - 3 => [4.0, 0.0], - 4 => [0.0, 1.0], - 5 => [2.0, 1.0], - 6 => [4.0, 1.0], - 7 => [0.0, 1.0], - 8 => [1.0, 1.0], - 9 => [3.0, 1.0], - 10 => [4.0, 1.0], - 11 => [0.0, 2.0], - 12 => [1.0, 2.0], - 13 => [3.0, 2.0], - 13 => [4.0, 1.0]) - - # volume elements - e1 = Quad4([1, 2, 5, 4]) - e1["geometry"] = Vector[N[1], N[2], N[5], N[4]] - e2 = Quad4([2, 3, 6, 5]) - e2["geometry"] = Vector[N[2], N[3], N[6], N[5]] - e3 = Quad4([7, 8, 12, 11]) - e3["geometry"] = Vector[N[7], N[8], N[12], N[11]] - e4 = Quad4([8, 9, 13, 12]) - e4["geometry"] = Vector[N[8], N[9], N[13], N[12]] - e5 = Quad4([9, 10, 14, 13]) - e5["geometry"] = Vector[N[9], N[10], N[14], N[13]] - - # boundary elements for boundary load - b1 = Seg2([11, 12]) - b1["geometry"] = Vector[N[11], N[12]] - b1["displacement traction force"] = Vector[[0.0, -10.0], [0.0, -10.0]] - b2 = Seg2([12, 13]) - b2["geometry"] = Vector[N[12], N[13]] - b2["displacement traction force"] = Vector[[0.0, -10.0], [0.0, -10.0]] - b3 = Seg3([13, 14]) - b3["geometry"] = Vector[N[13], N[14]] - b3["displacement traction force"] = Vector[[0.0, -10.0], [0.0, -10.0]] - - # boundary elements for dirichlet dy=0 - d1 = Seg2([1, 2]) - d1["geometry"] = Vector[N[1], N[2]] - d1["displacement 2"] = 0.0 - d2 = Seg2([2, 3]) - d2["geometry"] = Vector[N[2], N[3]] - d2["displacement 2"] = 0.0 - - # boundary elements for dirichlet dx=0 - d3 = Seg2([1, 4]) - d3["geometry"] = Vector[N[1], N[4]] - d3["displacement 1"] = 0.0 - d4 = Seg2([4, 11]) - d4["geometry"] = Vector[N[4], N[11]] - d4["displacmeent 1"] = 0.0 - - # mortar elements to tie meshes -- masters - m1 = MSeg2([4, 5]) - m1["geometry"] = Vector[N[4], N[5]] - m2 = MSeg2([5, 6]) - m2["geometry"] = Vector[N[5], N[6]] - - # mortar elements to tie meshes -- slaves - rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)] - phi = rotation_matrix(-pi/2) - m3 = MSeg2([7, 8]) - m3["geometry"] = Vector[N[7], N[8]] - m3["nodal ntsys"] = Matrix[phi, phi] - m3["master elements"] = MortarElement[m1, m2] - m4 = MSeg2([8, 9]) - m4["geometry"] = Vector[N[8], N[9]] - m4["nodal ntsys"] = Matrix[phi, phi] - m4["master elements"] = MortarElement[m1, m2] - m5 = MSeg2([9, 10]) - m5["geometry"] = Vector[N[9], N[10]] - m5["nodal ntsys"] = Matrix[phi, phi] - m5["master elements"] = MortarElement[m1, m2] - - problem1 = PlaneStressElasticityProblem() - push!(problem1, e1) - push!(problem1, e2) - push!(problem1, e3) - push!(problem1, e4) - push!(problem1, e5) - push!(problem1, b1) - push!(problem1, b2) - push!(problem1, b3) - - problem2 = DirichletProblem() - push!(problem2, d1) - push!(problem2, d2) - push!(problem2, d3) - push!(problem2, d4) - - problem3 = MortarProblem() - push!(problem3, m1) - push!(problem3, m2) - push!(problem3, m3) - push!(problem3, m4) - push!(problem3, m5) - - solver = DirectSolver() - push!(solver, problem1) - push!(solver, problem2) - push!(solver, problem3) - - call(solver) - -end - -function test_solver_multiple_dirichlet_bc() - N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]] - - e1 = Quad4([1, 2, 4, 3]) - e1["geometry"] = Vector[N[1], N[2], N[4], N[3]] - e1["youngs modulus"] = 900.0 - e1["poissons ratio"] = 0.25 - b1 = Seg2([3, 4]) - b1["geometry"] = Vector[N[3], N[4]] - b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]] - - problem = PlaneStressElasticityProblem() - push!(problem, e1) - push!(problem, b1) - - # manually solve problem 1 - # free_dofs = [3, 5, 6, 8] - # free_dofs = [3, 6, 7, 8] - #solve!(problem, free_dofs, 0.0; max_iterations=10) - #disp = e1("displacement", [1.0, 1.0], 0.0) - #info("displacement at tip: $disp") - #@test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01]) - - # boundary elements for dirichlet dx=0 - dx = Seg2([1, 3]) - dx["geometry"] = Vector[N[1], N[3]] - dx["displacement 1"] = 0.0 - - # boundary elements for dirichlet dy=0 - dy = Seg2([1, 2]) - dy["geometry"] = Vector[N[1], N[2]] - dy["displacement 2"] = 0.0 - - problem2 = DirichletProblem("displacement", 2) - push!(problem2, dx) - - problem3 = DirichletProblem("displacement", 2) - push!(problem3, dy) - - solver = DirectSolver() - push!(solver, problem) - push!(solver, problem2) - push!(solver, problem3) - - # launch solver - norm = solver(0.0) - -# info(e1("displacement")) -# info(last(e1["displacement"])) - disp = e1("displacement", [1.0, 1.0], 0.0) - info("displacement at tip: $disp") - @test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01]) - -end - -function test_solver_multiple_bodies_multiple_dirichlet_bc() - N = Vector[ - [0.0, 0.0], [1.0, 0.0], - [0.0, 1.0], [1.0, 1.0], - [0.0, 2.0], [1.0, 2.0]] - - e1 = Quad4([1, 2, 4, 3]) - e1["geometry"] = Vector[N[1], N[2], N[4], N[3]] - e2 = Quad4([3, 4, 6, 5]) - e2["geometry"] = Vector[N[3], N[4], N[6], N[5]] - for el in [e1, e2] - el["youngs modulus"] = 900.0 - el["poissons ratio"] = 0.25 - end - b1 = Seg2([5, 6]) - b1["geometry"] = Vector[N[5], N[6]] - b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]] - - body1 = PlaneStressElasticityProblem() - push!(body1, e1) - - body2 = PlaneStressElasticityProblem() - push!(body2, e2) - push!(body2, b1) - - # boundary elements for dirichlet dx=0 - dx1 = Seg2([1, 3]) - dx1["geometry"] = Vector[N[1], N[3]] - dx2 = Seg2([3, 5]) - dx2["geometry"] = Vector[N[3], N[5]] - for dx in [dx1, dx2] - dx["displacement 1"] = 0.0 - end - - boundary1 = DirichletProblem("displacement", 2) - push!(boundary1, dx1) - push!(boundary1, dx2) - - # boundary elements for dirichlet dy=0 - dy1 = Seg2([1, 2]) - dy1["geometry"] = Vector[N[1], N[2]] - dy1["displacement 2"] = 0.0 - - boundary2 = DirichletProblem("displacement", 2) - push!(boundary2, dy1) - - - solver = DirectSolver() - push!(solver, body1) - push!(solver, body2) - push!(solver, boundary1) - push!(solver, boundary2) - - # launch solver - norm = solver(0.0) - - disp = e2("displacement", [1.0, 1.0], 0.0) - info("displacement at tip: $disp") - # code aster verification, two_elements.comm - @test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01]) - -end - -# test_solver_multiple_bodies_multiple_dirichlet_bc() end diff --git a/test/test_virtual_work.jl b/test/test_virtual_work.jl index 0713004..e6dd270 100644 --- a/test/test_virtual_work.jl +++ b/test/test_virtual_work.jl @@ -5,57 +5,48 @@ module TestAutoDiffWeakForm using JuliaFEM.Test using JuliaFEM -using JuliaFEM: Quad4, Equation, IntegrationPoint, assemble!, Assembly, - solve!, get_field, get_element, get_basis, - grad, get_default_integration_points -""" Plane stress formulation for 4-node bilinear element. """ -type CPS4 <: Equation - element :: Quad4 - integration_points :: Vector{IntegrationPoint} +using JuliaFEM: Problem, AbstractProblem, CG, Element, IntegrationPoint, Quad4, solve! + +abstract PlaneStressElasticityProblem <: AbstractProblem + +function PlaneStressElasticityProblem(dim::Int=2, elements=[]) + return Problem{PlaneStressElasticityProblem}(dim, elements) end -function JuliaFEM.get_unknown_field_name(equation::CPS4) +function JuliaFEM.get_unknown_field_name{P<:PlaneStressElasticityProblem}(::Type{P}) return "displacement" end -function CPS4(element::Quad4) - integration_points = get_default_integration_points(element) - if !haskey(element, "displacement") - element["displacement"] = 0.0 => Vector{Float64}[[0.0,0.0], [0.0,0.0], [0.0,0.0], [0.0,0.0]] - end - CPS4(element, integration_points) +function JuliaFEM.get_unknown_field_type{P<:PlaneStressElasticityProblem}(::Type{P}) + return Vector{Float64} end -function Base.size(eq::CPS4) - return (2, 4) -end +function JuliaFEM.get_residual_vector{EL<:CG}(problem::Problem{PlaneStressElasticityProblem}, element::Element{EL}, ip::IntegrationPoint, time::Number; variation=nothing) -function JuliaFEM.get_residual_vector(equation::CPS4, ip, time; variation=nothing) - element = get_element(equation) - basis = get_basis(element) - dbasis = grad(basis) + + basis = element(ip, time) + dbasis = element(ip, time, Val{:grad}) # material parameters - E = basis("youngs modulus", ip, time) - nu = basis("poissons ratio", ip, time) + E = element("youngs modulus", ip, time) + nu = element("poissons ratio", ip, time) mu = E/(2*(1+nu)) la = E*nu/((1+nu)*(1-2*nu)) la = 2*la*mu/(la + 2*mu) # <- correction for 2d # elasticity formulation - u = basis("displacement", ip, time, variation) - gradu = dbasis("displacement", ip, time, variation) + u = element("displacement", ip, time, variation) + gradu = element("displacement", ip, time, Val{:grad}, variation) F = I + gradu - b = basis("displacement volume load", ip, time) + E = 1/2*(F'*F - I) S = la*trace(E)*I + 2*mu*E - P = F*S + r = F*S*dbasis + + b = element("displacement volume load", ip, time) + r -= b*basis - # residual vector - r_int = P*dbasis(ip,time) - r_ext = b*basis(ip,time) - r = r_int - r_ext return vec(r) end @@ -66,15 +57,18 @@ function test_residual_form() element["youngs modulus"] = 500.0 element["poissons ratio"] = 0.3 element["displacement volume load"] = Vector[[0.0,-10.0], [0.0,-10.0], [0.0,-10.0], [0.0,-10.0]] - equation = CPS4(element) + element["displacement"] = (0.0 => Vector{Float64}[zeros(2) for i=1:length(element)]) + problem = PlaneStressElasticityProblem() + push!(problem, element) # create model -- end free_dofs = [3, 4, 5, 6] - solve!(equation, free_dofs, 0.0) # launch a newton solver for single element - disp = get_basis(element)("displacement", [1.0, 1.0], 0.0)[2] - println("displacement at tip: $disp") + solve!(problem, free_dofs, 0.0) # launch a newton solver for single element + disp = element("displacement", [1.0, 1.0], 0.0) + info("displacement at tip: $disp") + # verified using Code Aster. - @test isapprox(disp, -8.77303119819776E+00) + @test isapprox(disp[2], -8.77303119819776E+00) end end