From d7bef419ed4946671f53fab17a4e532ce021cf93 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Fri, 19 Jan 2018 21:01:33 +0700 Subject: [PATCH] use FEMBase v0.1.x (#185) Lots of stuff moved from JuliaFEM.jl to FEMBase.jl. --- REQUIRE | 2 +- src/JuliaFEM.jl | 15 +- src/io.jl | 2 +- src/postprocess_utils.jl | 2 +- src/problems_contact_2d.jl | 22 +- src/problems_contact_2d_autodiff.jl | 61 +++--- src/problems_contact_3d.jl | 36 ++-- src/problems_dirichlet.jl | 15 ++ src/problems_heat.jl | 8 +- src/problems_mortar.jl | 8 +- src/problems_mortar_2d.jl | 28 +-- src/problems_mortar_2d_autodiff.jl | 60 +++--- src/problems_mortar_3d.jl | 49 +++-- src/solvers.jl | 8 +- src/solvers_modal.jl | 4 +- test/test_add_elements.jl | 13 -- test/test_common_failures.jl | 15 -- test/test_create_nodal_elements.jl | 16 -- ...ticity_2d_plane_stress_stiffness_matrix.jl | 8 +- ...asticity_3d_nonlinear_with_surface_load.jl | 24 +-- ...ity_hollow_sphere_with_surface_pressure.jl | 4 +- test/test_elasticity_tet10_mass_matrix.jl | 4 +- .../test_elasticity_tet10_stiffness_matrix.jl | 3 +- test/test_elasticity_tet4_stiffness_matrix.jl | 3 +- ...c_2d_nonhomogenious_boundary_conditions.jl | 9 +- ...sticplastic_3d_linear_with_surface_load.jl | 10 +- test/test_elements.jl | 75 ------- test/test_elements_2.jl | 53 ----- test/test_elements_add_fields.jl | 18 -- test/test_extrapolate_to_nodes.jl | 34 ---- test/test_fields.jl | 175 ---------------- test/test_fields_time_interpolation.jl | 23 --- test/test_heat.jl | 2 +- test/test_heat_3.jl | 2 +- test/test_integration_points.jl | 15 -- test/test_io.jl | 188 ------------------ test/test_modal_analysis.jl | 12 +- test/test_mortar_2d_calculate_projection.jl | 27 ++- test/test_nodal_constraints.jl | 89 --------- test/test_node_dof_mapping.jl | 26 --- test/test_preprocess.jl | 89 --------- test/test_problem.jl | 82 -------- test/test_solvers.jl | 78 -------- test/test_solvers_geometry_not_defined.jl | 14 -- test/test_sparse.jl | 47 ----- 45 files changed, 232 insertions(+), 1246 deletions(-) delete mode 100644 test/test_add_elements.jl delete mode 100644 test/test_common_failures.jl delete mode 100644 test/test_create_nodal_elements.jl delete mode 100644 test/test_elements.jl delete mode 100644 test/test_elements_2.jl delete mode 100644 test/test_elements_add_fields.jl delete mode 100644 test/test_extrapolate_to_nodes.jl delete mode 100644 test/test_fields.jl delete mode 100644 test/test_fields_time_interpolation.jl delete mode 100644 test/test_integration_points.jl delete mode 100644 test/test_io.jl delete mode 100644 test/test_nodal_constraints.jl delete mode 100644 test/test_node_dof_mapping.jl delete mode 100644 test/test_preprocess.jl delete mode 100644 test/test_problem.jl delete mode 100644 test/test_solvers.jl delete mode 100644 test/test_solvers_geometry_not_defined.jl delete mode 100644 test/test_sparse.jl diff --git a/REQUIRE b/REQUIRE index 4bae784..d2a2cab 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,5 +1,5 @@ julia 0.6 -FEMBase 0.0 0.1- +FEMBase 0.1 0.2- ForwardDiff LightXML 0.4 HDF5 0.7 diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 3fc2b6e..6686244 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -9,18 +9,6 @@ This is JuliaFEM -- Finite Element Package module JuliaFEM using FEMBase -using FEMBase: SparseMatrixCOO, SparseVectorCOO, Node, BasisInfo, - Discrete, Variable, TimeVariant, TimeInvariant, Field, - DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment, - IP, AbstractProblem, IntegrationPoint -using FEMBase: is_field_problem, is_boundary_problem, get_elements, - get_connectivity, assemble_prehook!, assemble_posthook!, - get_parent_field_name, get_reference_coordinates, - get_assembly, get_nonzero_rows, get_nonzero_columns, - eval_basis!, get_basis, get_dbasis, grad!, get_dualbasis, - assemble_mass_matrix!, get_local_coordinates, inside, - get_element_type, filter_by_element_type, get_element_id, - optimize!, resize_sparse, resize_sparsevec import FEMBase: get_unknown_field_name, get_unknown_field_dimension, assemble!, update!, initialize! @@ -107,7 +95,7 @@ end include("deprecations.jl") export SparseMatrixCOO, SparseVectorCOO, optimize!, resize_sparse -export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment +export DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment export FieldProblem, BoundaryProblem, Problem, Node, Element, Assembly export Poi1, Seg2, Seg3, Tri3, Tri6, Tri7, Quad4, Quad8, Quad9, Tet4, Tet10, Pyr5, Wedge6, Wedge15, Hex8, Hex20, Hex27 @@ -115,7 +103,6 @@ export update!, add_elements!, get_unknown_field_name, add!, is_field_problem, is_boundary_problem, get_gdofs, initialize!, get_integration_points, group_by_element_type, get_unknown_field_dimension, get_connectivity - export get_nonzero_rows, get_local_coordinates, inside, IP, get_element_type, get_elements, AbstractProblem, IntegrationPoint, filter_by_element_type, get_element_id, get_nonzero_columns, resize_sparse, resize_sparsevec diff --git a/src/io.jl b/src/io.jl index 9d16f92..c6ebc92 100644 --- a/src/io.jl +++ b/src/io.jl @@ -429,7 +429,7 @@ function update_xdmf!(xdmf::Xdmf, problem::Problem, time::Float64, fields::Vecto info("Xdmf: Saving topology of $nelements elements total, $nelement_types different element types.") for element_type in element_types - elements = filter_by_element_type(element_type, all_elements) + elements = collect(filter_by_element_type(element_type, all_elements)) nelements = length(elements) info("Xdmf: $nelements elements of type $element_type") sort!(elements, by=get_element_id) diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl index 140b4ed..ac909f8 100644 --- a/src/postprocess_utils.jl +++ b/src/postprocess_utils.jl @@ -69,7 +69,7 @@ Return node ids + vector of values function get_nodal_vector(elements::Vector, field_name::AbstractString, time::Float64) f = Dict() for element in elements - for (c, v) in zip(get_connectivity(element), element[field_name](time)) + for (c, v) in zip(get_connectivity(element), element(field_name, time)) if haskey(f, c) @assert isapprox(f[c], v) end diff --git a/src/problems_contact_2d.jl b/src/problems_contact_2d.jl index 722414b..c7b7d1c 100644 --- a/src/problems_contact_2d.jl +++ b/src/problems_contact_2d.jl @@ -80,7 +80,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::T la1 = slave_element("lambda", time) n1 = slave_element("normal", time) t1 = slave_element("tangent", time) - x1 = X1 + u1 + x1 = map(+, X1, u1) contact_area = 0.0 contact_error = 0.0 @@ -116,7 +116,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::T nm = length(master_element) X2 = master_element("geometry", time) u2 = master_element("displacement", time) - x2 = X2 + u2 + x2 = map(+, X2, u2) # 3.3. loop integration points of one integration segment and calculate # local mortar matrices @@ -136,20 +136,20 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::T Phi = Ae*N1 # project gauss point from slave element to master element in direction n_s - X_s = N1*X1 # coordinate in gauss point - n_s = N1*n1 # normal direction in gauss point - t_s = N1*t1 # tangent condition in gauss point + X_s = interpolate(N1, X1) # coordinate in gauss point + n_s = interpolate(N1, n1) # normal direction in gauss point + t_s = interpolate(N1, t1) # tangent condition in gauss point n_s /= norm(n_s) t_s /= norm(t_s) xi_m = project_from_slave_to_master(master_element, X_s, n_s, time) N2 = vec(get_basis(master_element, xi_m, time)) - X_m = N2*X2 + X_m = interpolate(N2, X2) - u_s = N1*u1 - u_m = N2*u2 - x_s = X_s + u_s - x_m = X_m + u_m - la_s = Phi*la1 + u_s = interpolate(N1, u1) + u_m = interpolate(N2, u2) + x_s = map(+, X_s, u_s) + x_m = map(+, X_m, u_m) + la_s = interpolate(Phi, la1) # virtual work De += w*Phi*N1' diff --git a/src/problems_contact_2d_autodiff.jl b/src/problems_contact_2d_autodiff.jl index 3033197..e4077ea 100644 --- a/src/problems_contact_2d_autodiff.jl +++ b/src/problems_contact_2d_autodiff.jl @@ -19,14 +19,20 @@ xi projected master """ -function project_from_master_to_slave{E<:MortarElements2D}( +function project_from_master_to_slave_ad{E<:MortarElements2D}( slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector; tol=1.0e-10, max_iterations=20, debug=false) - x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_ - dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_ - n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_ - dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_ + """ Multiply basis / dbasis at `xi` with field. """ + function mul(func, xi, field) + B = func(slave_element, [xi], time) + return sum(B[i]*field[i] for i=1:length(B)) + end + + x1(xi1) = mul(get_basis, xi1, x1_) + dx1(xi1) = mul(get_dbasis, xi1, x1_) + n1(xi1) = mul(get_basis, xi1, n1_) + dn1(xi1) = mul(get_dbasis, xi1, n1_) cross2(a, b) = cross([a; 0], [b; 0])[3] R(xi1) = cross2(x1(xi1)-x2, n1(xi1)) dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1)) @@ -62,12 +68,12 @@ function project_from_master_to_slave{E<:MortarElements2D}( end -function project_from_slave_to_master{E<:MortarElements2D}( - master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI; +function project_from_slave_to_master_ad{E<:MortarElements2D}( + master_element::Element{E}, x1, n1, x2_; tol=1.0e-10, max_iterations=20) - x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_ - dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_ + x2(xi2) = interpolate(vec(get_basis(master_element, [xi2], time)), x2_) + dx2(xi2) = interpolate(vec(get_dbasis(master_element, [xi2], time)), x2_) cross2(a, b) = cross([a; 0], [b; 0])[3] R(xi2) = cross2(x2(xi2)-x1, n1) dR(xi2) = cross2(dx2(xi2), n1) @@ -119,8 +125,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, push!(S, conn...) gdofs = get_gdofs(element, field_dim) X_el = element("geometry", time) - u_el = Field(Vector[u[:,i] for i in conn]) - x_el = X_el + u_el + x_el = tuple( (X_el[i] + u[:,j] for (i,j) in enumerate(conn))... ) #= for ip in get_integration_points(element, 3) dN = get_dbasis(element, ip, time) @@ -157,10 +162,10 @@ function assemble!(problem::Problem{Contact}, time::Float64, slave_element_nodes = get_connectivity(slave_element) X1 = slave_element("geometry", time) - u1 = Field(Vector[u[:,i] for i in slave_element_nodes]) - x1 = X1 + u1 - la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) - n1 = Field(Vector[normals[:,i] for i in slave_element_nodes]) + u1 = ((u[:,i] for i in slave_element_nodes)...) + x1 = ((Xi+ui for (Xi,ui) in zip(X1,u1))...) + la1 = ((la[:,i] for i in slave_element_nodes)...) + n1 = ((normals[:,i] for i in slave_element_nodes)...) nnodes = size(slave_element, 2) # construct dual basis @@ -170,12 +175,12 @@ function assemble!(problem::Problem{Contact}, time::Float64, master_element_nodes = get_connectivity(master_element) X2 = master_element("geometry", time) - u2 = Field(Vector[u[:,i] for i in master_element_nodes]) - x2 = X2 + u2 + u2 = ((u[:,i] for i in master_element_nodes)...) + x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...) # calculate segmentation: we care only about endpoints - xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1]) - xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2]) + xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1]) + xi1b = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[2]) xi1 = clamp.([xi1a; xi1b], -1.0, 1.0) l = 1/2*abs(xi1[2]-xi1[1]) isapprox(l, 0.0) && continue # no contribution in this master element @@ -200,8 +205,8 @@ function assemble!(problem::Problem{Contact}, time::Float64, master_element_nodes = get_connectivity(master_element) X2 = master_element("geometry", time) - u2 = Field(Vector[u[:,i] for i in master_element_nodes]) - x2 = X2 + u2 + u2 = ((u[:,i] for i in master_element_nodes)...) + x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...) #x1_midpoint = 1/2*(x1[1]+x1[2]) #x2_midpoint = 1/2*(x2[1]+x2[2]) @@ -209,8 +214,8 @@ function assemble!(problem::Problem{Contact}, time::Float64, #distance > props.maximum_distance && continue # calculate segmentation: we care only about endpoints - xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1]) - xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2]) + xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1]) + xi1b = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[2]) xi1 = clamp.([xi1a; xi1b], -1.0, 1.0) l = 1/2*abs(xi1[2]-xi1[1]) isapprox(l, 0.0) && continue # no contribution in this master element @@ -229,15 +234,15 @@ function assemble!(problem::Problem{Contact}, time::Float64, xi = ip.coords[1] xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1) N1 = vec(get_basis(slave_element, xi_s, time)) - x_s = N1*x1 # coordinate in gauss point - n_s = N1*n1 # normal direction in gauss point + x_s = interpolate(N1, x1) # coordinate in gauss point + n_s = interpolate(N1, n1) # normal direction in gauss point t_s = Q'*n_s # tangent direction in gauss point - xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2) + xi_m = project_from_slave_to_master_ad(master_element, x_s, n_s, x2) N2 = vec(get_basis(master_element, xi_m, time)) - x_m = N2*x2 + x_m = interpolate(N2, x2) Phi = Ae*N1 - la_s = Phi*la1 # traction force in gauss point + la_s = interpolate(Phi, la1) # traction force in gauss point gn = -dot(n_s, x_s - x_m) # normal gap fc[:,slave_element_nodes] += w*la_s*N1' diff --git a/src/problems_contact_3d.jl b/src/problems_contact_3d.jl index a23838c..6ce6028 100644 --- a/src/problems_contact_3d.jl +++ b/src/problems_contact_3d.jl @@ -83,13 +83,13 @@ function create_contact_segmentation(slave_element, master_elements, x0, n0, tim result = [] x1 = slave_element("geometry", time) if deformed - x1 += slave_element("displacement", time) + x1 = map(+, x1, slave_element("displacement", time)) end S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in x1] for master_element in master_elements x2 = master_element("geometry", time) if deformed - x2 += master_element("displacement", time) + x2 = map(+, x2, master_element("displacement", time)) end M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in x2] P = get_polygon_clip(S, M, n0) @@ -115,7 +115,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time nsl = length(slave_element) X1 = slave_element("geometry", time) u1 = slave_element("displacement", time) - x1 = X1 + u1 + x1 = map(+, X1, u1) n1 = slave_element("normal", time) la = slave_element("lambda", time) @@ -124,8 +124,8 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time # project slave nodes to auxiliary plane (x0, Q) xi = get_mean_xi(slave_element) N = vec(get_basis(slave_element, xi, time)) - x0 = N*X1 - n0 = N*n1 + x0 = interpolate(N, X1) + n0 = interpolate(N, n1) # create contact segmentation segmentation = create_contact_segmentation(slave_element, slave_element("master elements", time), x0, n0, time) @@ -147,7 +147,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time # loop integration cells for cell in get_cells(P, C0) virtual_element = Element(Tri3, Int[]) - update!(virtual_element, "geometry", cell) + update!(virtual_element, "geometry", tuple(cell...)) for ip in get_integration_points(virtual_element, 3) detJ = virtual_element(ip, time, Val{:detJ}) w = ip.weight*detJ @@ -173,7 +173,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time nm = length(master_element) X2 = master_element("geometry", time) u2 = master_element("displacement", time) - x2 = X2 + u2 + x2 = map(+, X2, u2) De = zeros(nsl, nsl) Me = zeros(nsl, nm) @@ -183,7 +183,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time # loop integration cells for cell in get_cells(P, C0) virtual_element = Element(Tri3, Int[]) - update!(virtual_element, "geometry", cell) + update!(virtual_element, "geometry", tuple(cell...)) # loop integration point of integration cell for ip in get_integration_points(virtual_element, 3) @@ -202,8 +202,8 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time De += w*Phi*N1' Me += w*Phi*N2' - x_s = N1*(X1+u1) - x_m = N2*(X2+u2) + x_s = interpolate(N1, map(+,X1,u1)) + x_m = interpolate(N2, map(+,X2,u2)) ge += w*vec((x_m-x_s)*Phi') end # integration points done @@ -282,8 +282,8 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time # create auxiliary plane xi = get_mean_xi(sub_slave_element) N = vec(get_basis(sub_slave_element, xi, time)) - x0 = N*X1 - n0 = N*n1 + x0 = interpolate(N, X1) + n0 = interpolate(N, n1) # project slave nodes to auxiliary plane S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1] @@ -323,7 +323,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time # 4. loop integration cells for cell in get_cells(P, C0) virtual_element = Element(Tri3, Int[]) - update!(virtual_element, "geometry", cell) + update!(virtual_element, "geometry", tuple(cell...)) for ip in get_integration_points(virtual_element, 3) detJ = virtual_element(ip, time, Val{:detJ}) w = ip.weight*detJ @@ -358,8 +358,8 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time # create auxiliary plane xi = get_mean_xi(sub_slave_element) N = vec(get_basis(sub_slave_element, xi, time)) - x0 = N*X1 - n0 = N*n1 + x0 = interpolate(N, X1) + n0 = interpolate(N, n1) # project slave nodes to auxiliary plane S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1] @@ -406,7 +406,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time # 4. loop integration cells for cell in get_cells(P, C0) virtual_element = Element(Tri3, Int[]) - update!(virtual_element, "geometry", cell) + update!(virtual_element, "geometry", tuple(cell...)) # 5. loop integration point of integration cell for ip in get_integration_points(virtual_element, 3) @@ -429,8 +429,8 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time us = slave_element("displacement", time) um = master_element("displacement", time) - xs = N1*(Xs+us) - xm = N2*(Xs+um) + xs = interpolate(N1, map(+,Xs,us)) + xm = interpolate(N2, map(+,Xs,um)) ge += w*vec((xm-xs)*Phi') end # integration points done diff --git a/src/problems_dirichlet.jl b/src/problems_dirichlet.jl index 1519556..2d8a2fb 100644 --- a/src/problems_dirichlet.jl +++ b/src/problems_dirichlet.jl @@ -15,6 +15,21 @@ function Dirichlet() Dirichlet(:incremental, false, false, 1) end +""" Return dual basis transformation matrix Ae. """ +function get_dualbasis(element::Element, time::Float64, order=1) + nnodes = length(element) + De = zeros(nnodes, nnodes) + Me = zeros(nnodes, nnodes) + for ip in get_integration_points(element, order) + detJ = element(ip, time, Val{:detJ}) + w = ip.weight*detJ + N = element(ip, time) + De += w*diagm(vec(N)) + Me += w*N'*N + end + return De, Me, De*inv(Me) +end + function get_formulation_type(problem::Problem{Dirichlet}) return problem.properties.formulation end diff --git a/src/problems_heat.jl b/src/problems_heat.jl index 7159d82..5198e4e 100644 --- a/src/problems_heat.jl +++ b/src/problems_heat.jl @@ -99,7 +99,7 @@ function assemble!{E<:Heat3DVolumeElements}(assembly::Assembly, problem::Problem fq += w*N'*f end end - T = vec(element[field_name](time)) + T = [interpolate(element[field_name], time)...] fq -= K*T add!(assembly.K, gdofs, gdofs, K) add!(assembly.f, gdofs, fq) @@ -136,7 +136,7 @@ function assemble!{E<:Heat3DSurfaceElements}(assembly::Assembly, problem::Proble fq += w*N'*h*Tu end end - T = vec(element[field_name](time)) + T = collect(element(field_name, time)) fq -= K*T add!(assembly.K, gdofs, gdofs, K) add!(assembly.f, gdofs, fq) @@ -180,7 +180,7 @@ function assemble!{E<:Heat2DVolumeElements}(assembly::Assembly, problem::Problem fq += w*N'*f end end - T = vec(element[field_name](time)) + T = collect(element(field_name, time)) fq -= K*T add!(assembly.K, gdofs, gdofs, K) add!(assembly.f, gdofs, fq) @@ -217,7 +217,7 @@ function assemble!{E<:Heat2DSurfaceElements}(assembly::Assembly, problem::Proble fq += w*N'*h*Tu end end - T = vec(element[field_name](time)) + T = collect(element(field_name, time)) fq -= K*T add!(assembly.K, gdofs, gdofs, K) add!(assembly.f, gdofs, fq) diff --git a/src/problems_mortar.jl b/src/problems_mortar.jl index fc89010..58a9a5d 100644 --- a/src/problems_mortar.jl +++ b/src/problems_mortar.jl @@ -143,13 +143,13 @@ function diagnose_interface(problem::Problem{Mortar}, time::Float64) info("Slave element connectivity = $slave_element_nodes") nsl = length(slave_element) X1 = slave_element("geometry", time) - n1 = Field([normals[j] for j in slave_element_nodes]) + n1 = tuple(collect(normals[j] for j in slave_element_nodes)...) # project slave nodes to auxiliary plane (x0, Q) xi = get_mean_xi(slave_element) N = vec(get_basis(slave_element, xi, time)) - x0 = N*X1 - n0 = N*n1 + x0 = interpolate(N,X1) + n0 = interpolate(N,n1) info("Auxiliary plane x0 = $x0, n0 = $n0") S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl] check_orientation!(S, n0) @@ -211,7 +211,7 @@ function diagnose_interface(problem::Problem{Mortar}, time::Float64) for (cell_id, cell) in enumerate(all_cells) C_area = 0.0 virtual_element = Element(Tri3, Int[]) - update!(virtual_element, "geometry", cell) + update!(virtual_element, "geometry", tuple(cell...)) # 5. loop integration point of integration cell for ip in get_integration_points(virtual_element, 3) diff --git a/src/problems_mortar_2d.jl b/src/problems_mortar_2d.jl index d0d5851..927503d 100644 --- a/src/problems_mortar_2d.jl +++ b/src/problems_mortar_2d.jl @@ -24,12 +24,12 @@ function get_slave_elements(problem::Problem) end function project_from_master_to_slave{E<:MortarElements2D}(slave_element::Element{E}, x2, time) - x1_ = slave_element["geometry"](time) - n1_ = slave_element["normal"](time) - x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_ - dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_ - n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_ - dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_ + x1_ = slave_element("geometry", time) + n1_ = slave_element("normal", time) + x1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), x1_) + dx1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), x1_) + n1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), n1_) + dn1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), n1_) R(xi1) = cross2(x1(xi1)-x2, n1(xi1)) dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1)) xi1 = nothing @@ -55,9 +55,9 @@ function project_from_master_to_slave{E<:MortarElements2D}(slave_element::Elemen end function project_from_slave_to_master{E<:MortarElements2D}(master_element::Element{E}, x1, n1, time) - x2_ = master_element["geometry"](time) - x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_ - dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_ + x2_ = master_element("geometry", time) + x2(xi2) = interpolate(vec(get_basis(master_element, [xi2], time)), x2_) + dx2(xi2) = interpolate(vec(get_dbasis(master_element, [xi2], time)), x2_) cross2(a, b) = cross([a; 0], [b; 0])[3] R(xi2) = cross2(x2(xi2)-x1, n1) dR(xi2) = cross2(dx2(xi2), n1) @@ -173,11 +173,11 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty N1 = vec(get_basis(slave_element, xi_s, time)) Phi = Ae*N1 # project gauss point from slave element to master element in direction n_s - X_s = N1*X1 # coordinate in gauss point - n_s = N1*n1 # normal direction in gauss point + X_s = interpolate(N1, X1) # coordinate in gauss point + n_s = interpolate(N1, n1) # normal direction in gauss point xi_m = project_from_slave_to_master(master_element, X_s, n_s, time) N2 = vec(get_basis(master_element, xi_m, time)) - X_m = N2*X2 + X_m = interpolate(N2, X2) De += w*Phi*N1' Me += w*Phi*N2' if props.adjust @@ -187,8 +187,8 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) < props.distval || continue u1 = slave_element("displacement", time) u2 = master_element("displacement", time) - x_s = X_s + N1*u1 - x_m = X_m + N2*u2 + x_s = X_s + interpolate(N1, u1) + x_m = X_m + interpolate(N2, u2) ge += w*vec((x_m-x_s)*Phi') end end diff --git a/src/problems_mortar_2d_autodiff.jl b/src/problems_mortar_2d_autodiff.jl index b85c762..b753676 100644 --- a/src/problems_mortar_2d_autodiff.jl +++ b/src/problems_mortar_2d_autodiff.jl @@ -5,14 +5,14 @@ using ForwardDiff # forwarddiff version of mesh tying in 2d -function project_from_master_to_slave{E<:MortarElements2D}( - slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector, time::Float64; +function project_from_master_to_slave_ad{E<:MortarElements2D}( + slave_element::Element{E}, x1_, n1_, x2, time; tol=1.0e-10, max_iterations=20) - x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_ - dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_ - n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_ - dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_ + x1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), x1_) + dx1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), x1_) + n1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), n1_) + dn1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), n1_) cross2(a, b) = cross([a; 0], [b; 0])[3] R(xi1) = cross2(x1(xi1)-x2, n1(xi1)) dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1)) @@ -37,12 +37,12 @@ function project_from_master_to_slave{E<:MortarElements2D}( end -function project_from_slave_to_master{E<:MortarElements2D}( - master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI, time::Float64; +function project_from_slave_to_master_ad{E<:MortarElements2D}( + master_element::Element{E}, x1, n1, x2_, time; tol=1.0e-10, max_iterations=20) - x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_ - dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_ + x2(xi2) = interpolate(vec(get_basis(master_element, [xi2], time)), x2_) + dx2(xi2) = interpolate(vec(get_dbasis(master_element, [xi2], time)), x2_) cross2(a, b) = cross([a; 0], [b; 0])[3] R(xi2) = cross2(x2(xi2)-x1, n1) dR(xi2) = cross2(dx2(xi2), n1) @@ -93,8 +93,8 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty conn = get_connectivity(element) push!(S, conn...) X1 = element("geometry", time) - u1 = Field([u[:,i] for i in conn]) - x1 = X1 + u1 + u1 = ((u[:,i] for i in conn)...) + x1 = map(+, X1, u1) dN = get_dbasis(element, [0.0], time) tangent = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) for nid in conn @@ -123,11 +123,11 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty nsl = length(slave_element) slave_element_nodes = get_connectivity(slave_element) - X1 = slave_element["geometry"](time) - u1 = Field(Vector[u[:,i] for i in slave_element_nodes]) - x1 = X1 + u1 - la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) - n1 = Field(Vector[normals[:,i] for i in slave_element_nodes]) + X1 = slave_element("geometry", time) + u1 = ((u[:,i] for i in slave_element_nodes)...) + x1 = map(+, X1, u1) + la1 = ((la[:,i] for i in slave_element_nodes)...) + n1 = ((normals[:,i] for i in slave_element_nodes)...) # 3. loop all master elements @@ -136,12 +136,12 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty nm = length(master_element) master_element_nodes = get_connectivity(master_element) X2 = master_element("geometry", time) - u2 = Field(Vector[u[:,i] for i in master_element_nodes]) - x2 = X2 + u2 + u2 = ((u[:,i] for i in master_element_nodes)...) + x2 = map(+, X2, u2) # 3.1 calculate segmentation - xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1], time) - xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2], time) + xi1a = project_from_master_to_slave_ad(slave_element, x1, n1, x2[1], time) + xi1b = project_from_master_to_slave_ad(slave_element, x1, n1, x2[2], time) # xi1a = project_from_master_to_slave(slave_element, X2[1], time) # xi1b = project_from_master_to_slave(slave_element, X2[2], time) xi1 = clamp.([xi1a; xi1b], -1.0, 1.0) @@ -181,20 +181,20 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty N1 = vec(get_basis(slave_element, xi_s, time)) Phi = Ae*N1 # project gauss point from slave element to master element in direction n_s - x_s = N1*x1 # coordinate in gauss point - n_s = N1*n1 # normal direction in gauss point + x_s = interpolate(N1, x1) # coordinate in gauss point + n_s = interpolate(N1, n1) # normal direction in gauss point #xi_m = project_from_slave_to_master(master_element, X_s, n_s, time) - xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2, time) + xi_m = project_from_slave_to_master_ad(master_element, x_s, n_s, x2, time) N2 = vec(get_basis(master_element, xi_m, time)) - x_m = N2*x2 + x_m = interpolate(N2, x2) - la_s = Phi*la1 + la_s = interpolate(Phi, la1) gn = dot(n_s, x_s-x_m) - u_s = N1*u1 - u_m = N2*u2 - X_s = N1*X1 - X_m = N2*X2 + u_s = interpolate(N1, u1) + u_m = interpolate(N2, u2) + X_s = interpolate(N1, X1) + X_m = interpolate(N2, X2) fc[:,slave_element_nodes] += w*la_s*N1' fc[:,master_element_nodes] -= w*la_s*N2' diff --git a/src/problems_mortar_3d.jl b/src/problems_mortar_3d.jl index 88dde0d..071a6e2 100644 --- a/src/problems_mortar_3d.jl +++ b/src/problems_mortar_3d.jl @@ -154,14 +154,23 @@ function get_polygon_clip{T}(xs::Vector{T}, xm::Vector{T}, n::T) end """ Project some vertex p to surface of element E using Newton's iterations. """ -function project_vertex_to_surface{E}(p::Vector, x0::Vector, n0::Vector, - element::Element{E}, x::DVTI, time::Real; - max_iterations::Int=10, iter_tol::Float64=1.0e-6) +function project_vertex_to_surface{E}(p, x0, n0, + element::Element{E}, x, time; + max_iterations=10, iter_tol=1.0e-6) basis(xi) = get_basis(element, xi, time) - dbasis(xi) = get_dbasis(element, xi, time) + function dbasis(xi) + return get_dbasis(element, xi, time) + end nnodes = length(element) - f(theta) = basis(theta[1:2])*x - theta[3]*n0 - p - L(theta) = inv3([dbasis(theta[1:2])*x -n0]) + mul(a,b) = sum((a[:,i]*b[i]')' for i=1:length(b)) + + function f(theta) + b = [basis(theta[1:2])*collect(x)...;] + b = b - theta[3]*n0 - p + return b + end + + L(theta) = inv3([mul(dbasis(theta[1:2]), x) -n0]) theta = zeros(3) dtheta = zeros(3) for i=1:max_iterations @@ -371,8 +380,8 @@ function assemble!{E<:Union{Tri3, Quad4}}(problem::Problem{Mortar}, slave_elemen xi = get_mean_xi(slave_element) first_slave_element && debug("midpoint xi = $xi") N = vec(get_basis(slave_element, xi, time)) - x0 = N*X1 - n0 = N*n1 + x0 = interpolate(N, X1) + n0 = interpolate(N, n1) S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl] master_elements = slave_element("master elements", time) @@ -412,7 +421,7 @@ function assemble!{E<:Union{Tri3, Quad4}}(problem::Problem{Mortar}, slave_elemen all_cells = get_cells(P, C0) for cell in all_cells virtual_element = Element(Tri3, Int[]) - update!(virtual_element, "geometry", cell) + update!(virtual_element, "geometry", tuple(cell...)) for ip in get_integration_points(virtual_element, 3) detJ = virtual_element(ip, time, Val{:detJ}) w = ip.weight*detJ @@ -476,7 +485,7 @@ function assemble!{E<:Union{Tri3, Quad4}}(problem::Problem{Mortar}, slave_elemen all_cells = get_cells(P, C0) for cell in all_cells virtual_element = Element(Tri3, Int[]) - virtual_element.fields["geometry"] = DVTI(cell) + update!(virtual_element, "geometry", tuple(cell...)) # 5. loop integration point of integration cell for ip in get_integration_points(virtual_element, 3) @@ -499,8 +508,8 @@ function assemble!{E<:Union{Tri3, Quad4}}(problem::Problem{Mortar}, slave_elemen if props.adjust && haskey(slave_element, "displacement") && haskey(master_element, "displacement") u1 = slave_element("displacement", time) u2 = master_element("displacement", time) - x_s = N1*(X1+u1) - x_m = N2*(X2+u2) + x_s = interpolate(N1, map(+,X1,u1)) + x_m = interpolate(N2, map(+,X2,u2)) ge += w*vec((x_m-x_s)*Phi') end area += w @@ -593,8 +602,8 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem xi = get_mean_xi(sub_slave_element) first_slave_element && debug("midpoint xi = $xi") N = vec(get_basis(sub_slave_element, xi, time)) - x0 = N*X1 - n0 = N*n1 + x0 = interpolate(N, X1) + n0 = interpolate(N, n1) # project slave nodes to auxiliary plane S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl] @@ -633,7 +642,7 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem all_cells = get_cells(P, C0) for cell in all_cells virtual_element = Element(Tri3, Int[]) - update!(virtual_element, "geometry", cell) + update!(virtual_element, "geometry", tuple(cell...)) for ip in get_integration_points(virtual_element, 3) x_gauss = virtual_element("geometry", ip, time) xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time) @@ -676,8 +685,8 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem xi = get_mean_xi(sub_slave_element) first_slave_element && debug("midpoint xi = $xi") N = vec(get_basis(sub_slave_element, xi, time)) - x0 = N*X1 - n0 = N*n1 + x0 = interpolate(N, X1) + n0 = interpolate(N, n1) # project slave nodes to auxiliary plane S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl] @@ -736,7 +745,7 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem all_cells = get_cells(P, C0) for cell in all_cells virtual_element = Element(Tri3, Int[]) - update!(virtual_element, "geometry", cell) + update!(virtual_element, "geometry", tuple(cell...)) # 5. loop integration point of integration cell for ip in get_integration_points(virtual_element, 3) @@ -758,8 +767,8 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem if props.adjust && haskey(slave_element, "displacement") && haskey(master_element, "displacement") u1 = slave_element("displacement", time) u2 = master_element("displacement", time) - xs = N1*(Xs+u1) - xm = N2*(Xm+u2) + xs = interpolate(N1, map(+,Xs,u1)) + xm = interpolate(N2, map(+,Xm,u2)) ge += w*vec((xm-xs)*Phi') end area += w diff --git a/src/solvers.jl b/src/solvers.jl index c0f6172..9670337 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -14,7 +14,7 @@ type Solver{S<:AbstractSolver} u :: Vector{Float64} la :: Vector{Float64} alpha :: Float64 # generalized alpha time integration coefficient - fields :: Dict{AbstractString, Field} + fields :: Dict{String, AbstractField} properties :: S end @@ -314,9 +314,9 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric end if !haskey(solver, "fint") - solver.fields["fint"] = Field(time => f) + solver.fields["fint"] = field(solver.time => f) else - update!(solver.fields["fint"], time => f) + update!(solver.fields["fint"], solver.time => f) end fint = solver.fields["fint"] @@ -327,7 +327,7 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric debug("Using generalized-α time integration, α=$alpha") K = (1-alpha)*K C1 = (1-alpha)*C1 - f = (1-alpha)*f + alpha*fint[end-1].data + f = (1-alpha)*f + alpha*fint.data[end-1].second end ndofs = solver.ndofs diff --git a/src/solvers_modal.jl b/src/solvers_modal.jl index cb11b5c..3b7b033 100644 --- a/src/solvers_modal.jl +++ b/src/solvers_modal.jl @@ -327,7 +327,7 @@ function update_xdmf!(solver::Solver{Modal}) elcon_arrays = Dict() @timeit "create topology arrays" for element_type in element_types - elements = filter_by_element_type(element_type, all_elements) + elements = collect(filter_by_element_type(element_type, all_elements)) nelements = length(elements) eldim = length(element_type) element_conn = zeros(Int, eldim, nelements) @@ -384,7 +384,7 @@ function update_xdmf!(solver::Solver{Modal}) for element_type in element_types timeit("save topology of element type $element_type") do - elements = filter_by_element_type(element_type, all_elements) + elements = collect(filter_by_element_type(element_type, all_elements)) nelements = length(elements) element_ids = map(get_element_id, elements) element_conn = elcon_arrays[element_type] diff --git a/test/test_add_elements.jl b/test/test_add_elements.jl deleted file mode 100644 index 10ee53d..0000000 --- a/test/test_add_elements.jl +++ /dev/null @@ -1,13 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using Base.Test - -@testset "add elements to problem" begin - problem = Problem(Elasticity, "test", 2) - element = Element(Quad4, [1, 2, 3, 4]) - elements = [element] - add_elements!(problem, elements) - @test problem.elements[1] == element -end diff --git a/test/test_common_failures.jl b/test/test_common_failures.jl deleted file mode 100644 index 774d044..0000000 --- a/test/test_common_failures.jl +++ /dev/null @@ -1,15 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -@testset "geometry missing" begin - el = Element(Quad4, [1, 2, 3, 4]) - pr = Problem(Elasticity, "problem", 2) - add_elements!(pr, [el]) - # this throws KeyError: geometry not found. - # it's descriptive enough to give hint to user - # what went wrong - @test_throws KeyError assemble!(pr) -end diff --git a/test/test_create_nodal_elements.jl b/test/test_create_nodal_elements.jl deleted file mode 100644 index f722699..0000000 --- a/test/test_create_nodal_elements.jl +++ /dev/null @@ -1,16 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using Base.Test - -using JuliaFEM -using JuliaFEM.Preprocess - -@testset "test create nodal elements" begin - m = Mesh() - add_node!(m, 1, [0.0, 0.0]) - add_node_to_node_set!(m, :test, 1) - els = create_nodal_elements(m, "test") - fel = first(els) - @test fel.connectivity == [1] -end diff --git a/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl b/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl index d3f6c5c..6fe3564 100644 --- a/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl +++ b/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl @@ -4,7 +4,8 @@ # http://ahojukka5.github.io/posts/finite-element-solution-for-one-element-problem/ using JuliaFEM -using JuliaFEM.Testing +using JuliaFEM: add_elements! +using Base.Test @testset "test 2d linear elasticity local matrices" begin element = Element(Quad4, [1, 2, 3, 4]) @@ -20,8 +21,9 @@ using JuliaFEM.Testing 4 => [0.0, 0.0]) update!(element, "geometry", X) update!(element, "displacement", u) - update!(element, "youngs modulus" => 288.0, "poissons ratio" => 1/3) - update!(element, "displacement load", DCTI([4.0, 8.0])) + update!(element, "youngs modulus", 288.0) + update!(element, "poissons ratio", 1/3) + update!(element, "displacement load", [4.0, 8.0]) problem = Problem(Elasticity, "[0x1] x [0x1] block", 2) update!(problem.properties, "formulation" => "plane_stress") diff --git a/test/test_elasticity_3d_nonlinear_with_surface_load.jl b/test/test_elasticity_3d_nonlinear_with_surface_load.jl index 11826d3..68aaa35 100644 --- a/test/test_elasticity_3d_nonlinear_with_surface_load.jl +++ b/test/test_elasticity_3d_nonlinear_with_surface_load.jl @@ -6,22 +6,22 @@ using JuliaFEM.Testing @testset "test continuum nonlinear elasticity with surface load" begin - nodes = Dict{Int64, Node}( - 1 => [0.0, 0.0, 0.0], - 2 => [1.0, 0.0, 0.0], - 3 => [1.0, 1.0, 0.0], - 4 => [0.0, 1.0, 0.0], - 5 => [0.0, 0.0, 1.0], - 6 => [1.0, 0.0, 1.0], - 7 => [1.0, 1.0, 1.0], - 8 => [0.0, 1.0, 1.0]) + X = Dict( + 1 => [0.0, 0.0, 0.0], + 2 => [1.0, 0.0, 0.0], + 3 => [1.0, 1.0, 0.0], + 4 => [0.0, 1.0, 0.0], + 5 => [0.0, 0.0, 1.0], + 6 => [1.0, 0.0, 1.0], + 7 => [1.0, 1.0, 1.0], + 8 => [0.0, 1.0, 1.0]) element1 = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8]) element2 = Element(Quad4, [5, 6, 7, 8]) - update!([element1, element2], "geometry", nodes) + update!([element1, element2], "geometry", X) update!([element1], "youngs modulus", 900.0) update!([element1], "poissons ratio", 0.25) - update!([element2], "displacement traction force", Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]) + update!([element2], "displacement traction force", [0.0, 0.0, -100.0]) elasticity_problem = Problem(Elasticity, "solve continuum block", 3) elasticity_problem.properties.finite_strain = true @@ -31,7 +31,7 @@ using JuliaFEM.Testing symxy = Element(Quad4, [1, 2, 3, 4]) symxz = Element(Quad4, [1, 2, 6, 5]) symyz = Element(Quad4, [1, 4, 8, 5]) - update!([symxy, symxz, symyz], "geometry", nodes) + update!([symxy, symxz, symyz], "geometry", X) symxy["displacement 3"] = 0.0 symxz["displacement 2"] = 0.0 symyz["displacement 1"] = 0.0 diff --git a/test/test_elasticity_hollow_sphere_with_surface_pressure.jl b/test/test_elasticity_hollow_sphere_with_surface_pressure.jl index 2d4fdd3..10e5ca1 100644 --- a/test/test_elasticity_hollow_sphere_with_surface_pressure.jl +++ b/test/test_elasticity_hollow_sphere_with_surface_pressure.jl @@ -108,8 +108,8 @@ function test_wedge_sphere(model, u_CA, S_CA) solver = LinearSolver(body, bc, lo) solver() - X = lo("geometry") - u = lo("displacement") + X = lo("geometry", 0.0) + u = lo("displacement", 0.0) nids = sort(collect(keys(X))) umag = Float64[norm(u[id]) for id in nids] um = mean(umag) diff --git a/test/test_elasticity_tet10_mass_matrix.jl b/test/test_elasticity_tet10_mass_matrix.jl index 90890eb..7f05a5c 100644 --- a/test/test_elasticity_tet10_mass_matrix.jl +++ b/test/test_elasticity_tet10_mass_matrix.jl @@ -2,8 +2,8 @@ # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md using JuliaFEM -using JuliaFEM: assemble_mass_matrix! -using JuliaFEM.Testing +using JuliaFEM: assemble_mass_matrix!, add_elements! +using Base.Test @testset "test tet10 mass matrix" begin X = Dict( diff --git a/test/test_elasticity_tet10_stiffness_matrix.jl b/test/test_elasticity_tet10_stiffness_matrix.jl index 52f97f5..e9b7851 100644 --- a/test/test_elasticity_tet10_stiffness_matrix.jl +++ b/test/test_elasticity_tet10_stiffness_matrix.jl @@ -2,7 +2,8 @@ # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md using JuliaFEM -using JuliaFEM.Testing +using JuliaFEM: add_elements! +using Base.Test @testset "test tet10 stiffness matrix" begin el = Element(Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) diff --git a/test/test_elasticity_tet4_stiffness_matrix.jl b/test/test_elasticity_tet4_stiffness_matrix.jl index cd3e89c..f59f85b 100644 --- a/test/test_elasticity_tet4_stiffness_matrix.jl +++ b/test/test_elasticity_tet4_stiffness_matrix.jl @@ -2,7 +2,8 @@ # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md using JuliaFEM -using JuliaFEM.Testing +using JuliaFEM: add_elements! +using Base.Test @testset "test tet4 stiffness matrix" begin el = Element(Tet4, [1, 2, 3, 4]) diff --git a/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl index 5eee961..5351699 100644 --- a/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl +++ b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl @@ -5,7 +5,8 @@ using JuliaFEM using JuliaFEM.Preprocess using JuliaFEM.Testing -# @testset "2d nonlinear elasticity: test nonhomogeneous boundary conditions and stress calculation" begin +#= +@testset "2d nonlinear elasticity: test nonhomogeneous boundary conditions and stress calculation" begin # field problem block = Problem(Elasticity, "BLOCK", 2) @@ -55,5 +56,7 @@ using JuliaFEM.Testing u3 = reshape(block.assembly.u, 2, 4)[:, 3] info("u3 = $u3") - #@test isapprox(u3, u3_expected, atol=1.0e-5) -# end + @test isapprox(u3, u3_expected, atol=1.0e-5) +end + +=# diff --git a/test/test_elasticplastic_3d_linear_with_surface_load.jl b/test/test_elasticplastic_3d_linear_with_surface_load.jl index 7bbe24b..ba3d358 100644 --- a/test/test_elasticplastic_3d_linear_with_surface_load.jl +++ b/test/test_elasticplastic_3d_linear_with_surface_load.jl @@ -5,7 +5,9 @@ using JuliaFEM using JuliaFEM.Preprocess using JuliaFEM.Testing -#@testset "test continuum 3d linear elasticity with surface load" begin +#= + +@testset "test continuum 3d linear elasticity with surface load" begin nodes = Dict{Int64, Node}( 1 => [0.0, 0.0, 0.0], 2 => [1.0, 0.0, 0.0], @@ -58,8 +60,10 @@ using JuliaFEM.Testing disp = element("displacement", [1.0, 1.0, 1.0], 1.0) info("displacement at tip: $disp") u_expected = 2.0 * [-1/3, -1/3, 1.0] - # @test isapprox(disp, u_expected) -#end + @test isapprox(disp, u_expected) +end + +=# # function solve_rod_model_elasticity(eltype) # fn = @__DIR__() * "/testdata/rod_short.med" diff --git a/test/test_elements.jl b/test/test_elements.jl deleted file mode 100644 index ee59dc4..0000000 --- a/test/test_elements.jl +++ /dev/null @@ -1,75 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -using JuliaFEM: group_by_element_type - -@testset "add time dependent field to element" begin - el = Element(Seg2, [1, 2]) - u1 = Vector{Float64}[[0.0, 0.0], [0.0, 0.0]] - u2 = Vector{Float64}[[1.0, 1.0], [1.0, 1.0]] - update!(el, "displacement", 0.0 => u1) - update!(el, "displacement", 1.0 => u2) - @test length(el["displacement"]) == 2 - @test isapprox(el("displacement", [0.0], 0.0), [0.0, 0.0]) - @test isapprox(el("displacement", [0.0], 0.5), [0.5, 0.5]) - @test isapprox(el("displacement", [0.0], 1.0), [1.0, 1.0]) - el2 = Element(Poi1, [1]) - update!(el2, "force 1", 0.0 => 1.0) -end - -@testset "add CVTV field to element" begin - el = Element(Seg2, [1, 2]) - f(xi, time) = xi[1]*time - update!(el, "my field", f) - v = el("my field", [1.0], 2.0) - @test isapprox(v, 2.0) -end - -@testset "add DCTI to element" begin - el = Element(Quad4, [1, 2, 3, 4]) - update!(el, "displacement load", DCTI([4.0, 8.0])) - @test isa(el["displacement load"], DCTI) - @test !isa(el["displacement load"].data, DCTI) - update!(el, "displacement load 2", [4.0, 8.0]) - @test isa(el["displacement load 2"], DCTI) - update!(el, "temperature", [1.0, 2.0, 3.0, 4.0]) - @test isa(el["temperature"], DVTI) - @test isapprox(el("displacement load", [0.0, 0.0], 0.0), [4.0, 8.0]) -end - -@testset "interpolate DCTI from element" begin - el = Element(Seg2, [1, 2]) - update!(el, "foobar", 1.0) - fb = el("foobar", [0.0], 0.0) - @test isa(fb, Float64) - @test isapprox(fb, 1.0) -end - -@testset "add elements to elements" begin - el1 = Element(Seg2, [1, 2]) - el2 = Element(Seg2, [3, 4]) - update!(el1, "master elements", [el2]) - lst = el1("master elements", 0.0) - @test isa(lst, Vector) -end - -@testset "extend basis" begin - el = Element(Quad4, [1, 2, 3, 4]) - expected = [ - 0.25 0.00 0.25 0.00 0.25 0.00 0.25 0.00 - 0.00 0.25 0.00 0.25 0.00 0.25 0.00 0.25] - @test isapprox(el([0.0, 0.0], 0.0, 2), expected) -end - -@testset "group elements" begin - e1 = Element(Seg2, [1, 2]) - e2 = Element(Quad4, [1, 2, 3, 4]) - elements = [e1, e2] - r = group_by_element_type(elements) - @test length(r) == 2 - @test first(r[Element{Seg2}]) == e1 - @test first(r[Element{Quad4}]) == e2 -end diff --git a/test/test_elements_2.jl b/test/test_elements_2.jl deleted file mode 100644 index e38f5b3..0000000 --- a/test/test_elements_2.jl +++ /dev/null @@ -1,53 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -@testset "inverse isoparametric mapping" begin - el = Element(Quad4, [1, 2, 3, 4]) - X = Dict{Int64, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [1.0, 1.0], - 4 => [0.0, 1.0]) - update!(el, "geometry", X) - time = 0.0 - X1 = el("geometry", [0.1, 0.2], time) - xi = get_local_coordinates(el, X1, time) - X2 = el("geometry", xi, time) - info("X1 = $X1, X2 = $X2") - @test isapprox(X1, X2) -end - -@testset "inside of linear element" begin - el = Element(Quad4, [1, 2, 3, 4]) - X = Dict{Int64, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [1.0, 1.0], - 4 => [0.0, 1.0]) - update!(el, "geometry", X) - time = 0.0 - @test inside(el, [0.5, 0.5], time) == true - @test inside(el, [1.0, 0.5], time) == true - @test inside(el, [1.0, 1.0], time) == true - @test inside(el, [1.01, 1.0], time) == false - @test inside(el, [1.0, 1.01], time) == false -end - -@testset "inside of quadratic element" begin - el = Element(Tri6, [1, 2, 3, 4, 5, 6]) - X = Dict{Int64, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [0.0, 1.0], - 4 => [0.5, 0.2], - 5 => [0.8, 0.6], - 6 => [-0.2, 0.5]) - update!(el, "geometry", X) - p = [0.94, 0.3] # visually checked to be inside - @test inside(el, p, 0.0) == true - p = [-0.2, 0.8] # visually checked to be outside - @test inside(el, p, 0.0) == false -end diff --git a/test/test_elements_add_fields.jl b/test/test_elements_add_fields.jl deleted file mode 100644 index 594ef2e..0000000 --- a/test/test_elements_add_fields.jl +++ /dev/null @@ -1,18 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -@testset "dict field" begin - el = Element(Seg2, [1, 2]) - X = Dict{Int64, Vector{Float64}}(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [0.5, 0.5]) - f = Field(X) - debug("field = $f") - #update!(el, "geometry", X) - el["geometry"] = f - @test isapprox(el("geometry")[1], [0.0, 0.0]) - @test isapprox(el("geometry", 0.0)[1], [0.0, 0.0]) - @test isapprox(el("geometry", 0.0)[3], [0.5, 0.5]) - @test isapprox(el("geometry", [0.0], 0.0), [0.5, 0.0]) -end diff --git a/test/test_extrapolate_to_nodes.jl b/test/test_extrapolate_to_nodes.jl deleted file mode 100644 index 677e71f..0000000 --- a/test/test_extrapolate_to_nodes.jl +++ /dev/null @@ -1,34 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Postprocess -using JuliaFEM.Testing - -@testset "extrapolate stress from gauss points to nodes" begin - X = Dict{Int, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [6.0, 0.0], - 3 => [6.0, 6.0], - 4 => [0.0, 6.0], - 5 => [12.0, 0.0], - 6 => [12.0, 6.0]) - el1 = Element(Quad4, [1, 2, 3, 4]) - el2 = Element(Quad4, [2, 5, 6, 3]) - el1.id = 1 - el2.id = 2 - elements = [el1, el2] - time = 0.0 - update!(elements, "geometry", X) - update!(get_integration_points(el1), "stress", time => [1.0, 2.0, 3.0]) - update!(get_integration_points(el2), "stress", time => [2.0, 3.0, 4.0]) - field_name = "stress" - field_dim = 3 - calc_nodal_values!(elements, field_name, field_dim, time) - s1 = el1("stress", [0.0, 0.0], time) - s2 = el2("stress", [0.0, 0.0], time) - # visually checked, see blog post "Postprocessing stress" - @test isapprox(s1, [1.125, 2.125, 3.125]) - @test isapprox(s2, [1.875, 2.875, 3.875]) -end - diff --git a/test/test_fields.jl b/test/test_fields.jl deleted file mode 100644 index fba93ac..0000000 --- a/test/test_fields.jl +++ /dev/null @@ -1,175 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -@testset "discrete, constant, time invariant field" begin - @test DCTI(0.0).data == 0.0 - @test isa(Field(0.0), DCTI) - f = DCTI(0.0) - update!(f, 1.0) - @test isapprox(f, DCTI(1.0)) - @test isapprox(f, 1.0) - @test 2*f == 2.0 # multiply by constant - @test f(1.0) == 1.0 # time interpolation - @test isapprox(reshape([2.0],1,1)*f, 2.0) # wanted behavior? -end - -@testset "discrete, variable, time invariant field" begin - @test DVTI([1.0, 2.0]).data == [1.0, 2.0] - @test isa(Field([1.0, 2.0]), DVTI) - - f = DVTI(zeros(2)) - update!(f, [2.0, 3.0]) - @test isapprox(f.data, [2.0, 3.0]) - @test length(f) == 2 - - # slicing - @test isapprox(f[1], 2.0) - @test isapprox(f[[1, 2]], [2.0, 3.0]) - - # boolean comparison and multiplying by a constant - @test f == DVTI([2.0, 3.0]) - @test isapprox(2*f, [4.0, 6.0]) - - f3 = 2*f - @test isa(f3, DVTI) - @test f3+f == 3*f - @test f3-f == f - - # spatial interpolation - N = [1.0, 2.0] - @test isapprox(N*f, 8.0) - - # time interpolation - @test isapprox(f(1.0), [2.0, 3.0]) - - # spatial interpolation of vector valued variable field - f2 = DVTI(Vector[[1.0, 2.0], [3.0, 4.0]]) - @test isapprox(f2[1], [1.0, 2.0]) - @test isapprox(f2[2], [3.0, 4.0]) - @test length(f2) == 2 - @test isapprox(N*f2, [1.0, 2.0] + [6.0, 8.0]) - - # iteration of DVTI field - s = zeros(2) - for j in f2 - s += j - end - @test isapprox(s, [4.0, 6.0]) - - @test vec(f2) == [1.0, 2.0, 3.0, 4.0] - @test isapprox([1.0 2.0]*f, [8.0]'') - - new_data = [2.0, 3.0, 4.0, 5.0] - f4 = similar(f2, new_data) - @test isa(f4, DVTI) - @test isapprox(f4.data[1], [2.0, 3.0]) - @test isapprox(f4.data[2], [4.0, 5.0]) -end - -@testset "discrete, constant, time-variant field" begin - f = Field(0.0 => 1.0) - @test isa(f, DCTV) - @test last(f).time == 0.0 - @test last(f).data == 1.0 - update!(f, 0.0 => 2.0) - @test last(f).time == 0.0 - @test last(f).data == 2.0 - @test length(f) == 1 - update!(f, 1.0 => 3.0) - @test last(f).time == 1.0 - @test last(f).data == 3.0 - @test length(f) == 2 - - @testset "interpolation in time direction" begin - @test isa(f(0.0), DCTI) # converts to time-invariant after time interpolation - @test isapprox(f(-1.0), 2.0) - @test isapprox(f(0.0), 2.0) - @test isapprox(f(0.5), 2.5) - @test isapprox(f(1.0), 3.0) - @test isapprox(f(2.0), 3.0) - end - - # create several time steps at once - f = DCTV(0.0 => 1.0, 1.0 => 2.0) - @test isapprox(f(0.5), 1.5) - -end - -@testset "discrete, variable, time-variant field" begin - f = Field(0.0 => [1.0, 2.0]) - @test isa(f, DVTV) - @test last(f).time == 0.0 - @test last(f).data == [1.0, 2.0] - update!(f, 0.0 => [2.0, 3.0]) - @test last(f).time == 0.0 - @test last(f).data == [2.0, 3.0] - @test length(f) == 1 - update!(f, 1.0 => [3.0, 4.0]) - @test last(f).time == 1.0 - @test last(f).data == [3.0, 4.0] - @test length(f) == 2 - - @testset "interpolation in time direction" begin - @test isa(f(0.0), DVTI) # converts to time-invariant after time interpolation - @test isapprox(f(-1.0), [2.0, 3.0]) - @test isapprox(f(0.0), [2.0, 3.0]) - @test isapprox(f(0.5), [2.5, 3.5]) - @test isapprox(f(1.0), [3.0, 4.0]) - @test isapprox(f(2.0), [3.0, 4.0]) - end - - # create several time steps at once - f = DVTV(0.0 => [1.0, 2.0], 1.0 => [2.0, 3.0]) - @test isapprox(f(0.5), [1.5, 2.5]) -end - -@testset "continuous, constant, time-invariant field" begin - f = Field(() -> 2.0) - @test isapprox(f([1.0], 2.0), 2.0) - -end - -@testset "continuous, constant, time variant field" begin - f = Field((time::Float64) -> 2.0*time) - @test isapprox(f([1.0], 2.0), 4.0) - -end - -@testset "continuous, variable, time invariant field" begin - f = Field((xi::Vector) -> sum(xi)) - @test isapprox(f([1.0, 2.0], 2.0), 3.0) -end - -@testset "continuous, variable, time variant field" begin - f = Field((xi::Vector, t::Float64) -> xi[1]*t) - @test isapprox(f([1.0], 2.0), 2.0) -end - -@testset "unknown function argument for continuous field" begin - @test_throws ErrorException Field((a, b, c) -> a*b*c) -end - -@testset "dictionary fields" begin - f1 = Dict{Int64, Vector{Float64}}(1 => [0.0, 0.0], 2 => [0.0, 0.0]) - f2 = Dict{Int64, Vector{Float64}}(1 => [1.0, 1.0], 2 => [1.0, 1.0]) - f = Field(0.0 => f1, 1.0 => f2) - @test isa(f, DVTV) - @test isapprox(f(0.0)[1], [0.0, 0.0]) - @test isapprox(f(1.0)[2], [1.0, 1.0]) - - f = Field(0.0 => f1) - update!(f, 1.0 => f2) - @test isa(f, DVTV) - @test isapprox(f(0.0)[1], [0.0, 0.0]) - @test isapprox(f(1.0)[2], [1.0, 1.0]) - - f = Field(f1) - @test isapprox(f(0.0)[1], [0.0, 0.0]) - @test isapprox(f[1], [0.0, 0.0]) - - f = Field(f1) - @test isa(f, DVTI) -end diff --git a/test/test_fields_time_interpolation.jl b/test/test_fields_time_interpolation.jl deleted file mode 100644 index ed76854..0000000 --- a/test/test_fields_time_interpolation.jl +++ /dev/null @@ -1,23 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -@testset "test interpolation of discrete constant time-variant field" begin - f = DCTV(0.0 => 0.0, 1.0 => 1.0) - @test isapprox(f(-1.0), DCTI(0.0)) - @test isapprox(f( 0.0), DCTI(0.0)) - @test isapprox(f( 0.3), DCTI(0.3)) - @test isapprox(f( 0.5), DCTI(0.5)) - @test isapprox(f( 0.9), DCTI(0.9)) - @test isapprox(f( 1.0), DCTI(1.0)) - @test isapprox(f( 1.5), DCTI(1.0)) - - f2 = DCTV(0.0 => 0.0, 0.25 => -0.1, 0.50 => -0.1) - @test isapprox(f2(0.0), DCTI(0.0)) - @test isapprox(f2(0.25), DCTI(-0.1)) - @test isapprox(f2(0.50), DCTI(-0.1)) - @test isapprox(f2(0.35), DCTI(-0.1)) -end - diff --git a/test/test_heat.jl b/test/test_heat.jl index d80e7e8..99c3aef 100644 --- a/test/test_heat.jl +++ b/test/test_heat.jl @@ -207,7 +207,7 @@ end 7 => [1.74360055518019E+04, -4.73227515822118E+02, -1.75280965396335E+02], 8 => [1.74447696000717E+04, -4.72904678032179E+02, -1.75280965396335E+02]) - T = p1("temperature") + T = p1("temperature", 0.0) for j in sort(collect(keys(T))) T1 = T[j][1] diff --git a/test/test_heat_3.jl b/test/test_heat_3.jl index bc64cf8..d247d79 100644 --- a/test/test_heat_3.jl +++ b/test/test_heat_3.jl @@ -35,7 +35,7 @@ using JuliaFEM.Testing T_fem = Float64[] T_acc = Float64[] - for (nid, X) in field("geometry") + for (nid, X) in field("geometry", 0.0) push!(T_fem, field("temperature", X)[1]) push!(T_acc, 1.0 + X[1]^2 + 2*X[2]^2) end diff --git a/test/test_integration_points.jl b/test/test_integration_points.jl deleted file mode 100644 index c7a8a5a..0000000 --- a/test/test_integration_points.jl +++ /dev/null @@ -1,15 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -@testset "test integration point" begin - a = sqrt(1.0/3.0) - ip = IP(1, 1.0, (a, -a)) - strain = [1.0 2.0; 3.0 4.0] - update!(ip, "strain", 0.0 => strain) - @test isapprox(ip("strain", 0.0), strain) - @test isapprox(ip("strain"), strain) -end - diff --git a/test/test_io.jl b/test/test_io.jl deleted file mode 100644 index 9f02d79..0000000 --- a/test/test_io.jl +++ /dev/null @@ -1,188 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing -using LightXML - -@testset "create new Xdmf object" begin - r = Xdmf() - expected = "" - @test string(r.xml) == expected -end - -@testset "put and get to Xdmf, low level" begin - io = Xdmf() - # h5 - write(io.hdf, "/Xdmf/Domain/Geometry", [1 2 3]) - @test isapprox(read(io.hdf, "/Xdmf/Domain/Geometry"), [1 2 3]) - # xml - obj = new_child(io.xml, "Domain") - set_attribute(obj, "Name", "Test Domain") - obj2 = find_element(io.xml, "Domain") - @test attribute(obj2, "Name") == "Test Domain" -end - -@testset "write data to HDF, automatically generate path" begin - xdmf = Xdmf() - di1 = new_dataitem(xdmf, [1 2 3]) - di2 = new_dataitem(xdmf, [4 5 6]) - @test contains(content(di1), "DataItem_1") - @test contains(content(di2), "DataItem_2") -end - -@testset "Xdmf filtering" begin - grid1 = new_element("Grid") - add_text(grid1, "I am first grid") - grid2 = new_element("Grid") - add_text(grid2, "I am second grid") - set_attribute(grid2, "Name", "Frame 2") - grid3 = new_element("Grid") - add_text(grid3, "I am third grid") - grids = [grid1, grid2, grid3] - - @test content(xdmf_filter(grids, "Grid")) == "I am first grid" - @test content(xdmf_filter(grids, "Grid[1]")) == "I am first grid" - @test content(xdmf_filter(grids, "Grid[2]")) == "I am second grid" - @test content(xdmf_filter(grids, "Grid[3]")) == "I am third grid" - @test content(xdmf_filter(grids, "Grid[end]")) == "I am third grid" - @test content(xdmf_filter(grids, "Grid[@Name=Frame 2]")) == "I am second grid" - @test xdmf_filter(grids, "Grid[0]") == nothing - @test xdmf_filter(grids, "Grid[4]") == nothing - @test xdmf_filter(grids, "Grid[@Name=Frame 3]") == nothing - @test xdmf_filter(grids, "Domain/Grid[@Name=Frame 3]") == nothing - @test xdmf_filter(grids, "Domain") == nothing -end - -@testset "XML traverse" begin - xdmf = Xdmf() - domain = new_child(xdmf.xml, "Domain") - grid = new_child(domain, "Grid") - set_attribute(grid, "CollectionType", "Temporal") - set_attribute(grid, "GridType", "Collection") - - frame1 = new_child(grid, "Grid") - time1 = new_child(frame1, "Time") - set_attribute(time1, "Value", 0.0) - X1 = new_child(frame1, "Geometry") - set_attribute(X1, "Type", "XY") - - frame2 = new_child(grid, "Grid") - set_attribute(frame2, "Name", "Frame 2") - time2 = new_child(frame2, "Time") - set_attribute(time2, "Value", 1.0) - X2 = new_child(frame2, "Geometry") - set_attribute(X2, "Type", "XY") - - add_child(grid, frame1) - add_child(grid, frame2) - - dataitem = new_dataitem(xdmf, "/Domain/Grid/Grid/2/Geometry", [1.0, 2.0]) - add_child(X2, dataitem) - - println(xdmf.xml) - @test read(xdmf, "/Domain/Grid/Grid/Time/Value") == "0.0" - @test read(xdmf, "/Domain/Grid/Grid[2]/Time/Value") == "1.0" - @test read(xdmf, "/Domain/Grid/Grid[end]/Time/Value") == "1.0" - @test read(xdmf, "/Domain/Grid/Grid[@Name=Frame 2]/Time/Value") == "1.0" - @test isapprox(read(xdmf, "/Domain/Grid/Grid[2]/Geometry/DataItem"), [1.0, 2.0]) -end - -@testset "write fields from different problems to Xdmf file" begin - - X = Dict( - 1 => [0.0, 0.0, 0.0], - 2 => [1.0, 0.0, 0.0], - 3 => [1.0, 1.0, 0.0], - 4 => [0.0, 1.0, 0.0], - 5 => [0.0, 0.0, 0.5], - 6 => [1.0, 0.0, 0.5], - 7 => [1.0, 1.0, 0.5], - 8 => [0.0, 1.0, 0.5], - 9 => [0.0, 0.0, 1.0], - 10 => [1.0, 0.0, 1.0], - 11 => [1.0, 1.0, 1.0], - 12 => [0.0, 1.0, 1.0]) - - u = Dict() - u[0] = Dict( - 1 => [0.0, 0.0, 0.0], - 2 => [0.0, 0.0, 0.0], - 3 => [0.0, 0.0, 0.0], - 4 => [0.0, 0.0, 0.0], - 5 => [0.0, 0.0, 0.0], - 6 => [0.0, 0.0, 0.0], - 7 => [0.0, 0.0, 0.0], - 8 => [0.0, 0.0, 0.0], - 9 => [0.0, 0.0, 0.0], - 10 => [0.0, 0.0, 0.0], - 11 => [0.0, 0.0, 0.0], - 12 => [0.0, 0.0, 0.0]) - u[1] = Dict( - 1 => [0.0, 0.0, 0.0], - 2 => [0.0, 0.0, 0.0], - 3 => [0.0, 0.0, 0.0], - 4 => [0.0, 0.0, 0.0], - 5 => [0.0, 0.0, -0.1], - 6 => [0.0, 0.0, -0.1], - 7 => [0.0, 0.0, -0.1], - 8 => [0.0, 0.0, -0.1], - 9 => [0.0, 0.0, -0.2], - 10 => [0.0, 0.0, -0.2], - 11 => [0.0, 0.0, -0.2], - 12 => [0.0, 0.0, -0.2]) - - T = Dict( - 1 => 10.0, - 2 => 10.0, - 3 => 10.0, - 4 => 10.0, - 5 => 20.0, - 6 => 20.0, - 7 => 20.0, - 8 => 20.0, - 9 => 30.0, - 10 => 30.0, - 11 => 30.0, - 12 => 30.0) - - rf = Dict() - rf[0] = Dict( - 1 => [0.0, 0.0, 0.0], - 2 => [0.0, 0.0, 0.0], - 3 => [0.0, 0.0, 0.0], - 4 => [0.0, 0.0, 0.0]) - rf[1] = Dict( - 1 => [0.0, 0.0, 1.0], - 2 => [0.0, 0.0, 1.0], - 3 => [0.0, 0.0, 1.0], - 4 => [0.0, 0.0, 1.0]) - - e1 = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8]) - e2 = Element(Hex8, [5, 6, 7, 8, 9, 10, 11, 12]) - e3 = Element(Quad4, [1, 2, 3, 4]) - update!([e1, e2, e3], "geometry", X) - update!([e1, e2, e3], "displacement", 0.0 => u[0]) - update!([e1, e2, e3], "displacement", 1.0 => u[1]) - update!([e1, e2, e3], "temperature", T) - update!(e3, "reaction force", 0.0 => rf[0]) - update!(e3, "reaction force", 1.0 => rf[1]) - - p1 = Problem(Elasticity, "lower", 3) - p1.elements = [e1] - p2 = Problem(Elasticity, "upper", 3) - p2.elements = [e2] - p3 = Problem(Dirichlet, "bc", 3, "displacement") - p3.elements = [e3] - - xdmf = Xdmf() - xdmf.format = "XML" - update_xdmf!(xdmf, p1, 0.0, ["displacement", "temperature"]) - update_xdmf!(xdmf, p2, 0.0, ["displacement"]) - update_xdmf!(xdmf, p3, 0.0, ["reaction force"]) - update_xdmf!(xdmf, p1, 1.0, ["displacement", "temperature"]) - update_xdmf!(xdmf, p2, 1.0, ["displacement"]) - update_xdmf!(xdmf, p3, 1.0, ["reaction force"]) - @test read(xdmf, "/Domain/Grid/Grid/Time/Value") == "0.0" - @test read(xdmf, "/Domain/Grid/Grid[2]/Time/Value") == "1.0" -end diff --git a/test/test_modal_analysis.jl b/test/test_modal_analysis.jl index 8d981b7..5e05957 100644 --- a/test/test_modal_analysis.jl +++ b/test/test_modal_analysis.jl @@ -19,12 +19,12 @@ function get_model() e2 = Element(Tri3, [1, 2, 3]) update!([e1, e2], "geometry", X) update!([e1, e2], "displacement", 0.0 => u) - update!(e1, "youngs modulus" => 96.0) - update!(e1, "poissons ratio" => 1.0/3.0) - update!(e1, "density" => 420.0) - update!(e2, "displacement 1" => 0.0) - update!(e2, "displacement 2" => 0.0) - update!(e2, "displacement 3" => 0.0) + update!(e1, "youngs modulus", 96.0) + update!(e1, "poissons ratio", 1.0/3.0) + update!(e1, "density", 420.0) + update!(e2, "displacement 1", 0.0) + update!(e2, "displacement 2", 0.0) + update!(e2, "displacement 3", 0.0) p1 = Problem(Elasticity, "test problem", 3) p1.properties.finite_strain = false p1.properties.geometric_stiffness = false diff --git a/test/test_mortar_2d_calculate_projection.jl b/test/test_mortar_2d_calculate_projection.jl index d957982..ac5aa33 100644 --- a/test/test_mortar_2d_calculate_projection.jl +++ b/test/test_mortar_2d_calculate_projection.jl @@ -3,22 +3,26 @@ using JuliaFEM using JuliaFEM.Testing +using JuliaFEM: calculate_normals function get_test_2d_model() - X = Dict{Int64, Vector{Float64}}( - 7 => [0.0, 1.0], - 8 => [5/4, 1.0], - 9 => [2.0, 1.0], - 10 => [0.0, 1.0], - 11 => [3/4, 1.0], - 12 => [2.0, 1.0]) + X = Dict( + 7 => [0.0, 1.0], + 8 => [5/4, 1.0], + 9 => [2.0, 1.0], + 10 => [0.0, 1.0], + 11 => [3/4, 1.0], + 12 => [2.0, 1.0]) mel1 = Element(Seg2, [7, 8]) mel2 = Element(Seg2, [8, 9]) sel1 = Element(Seg2, [10, 11]) sel2 = Element(Seg2, [11, 12]) update!([mel1, mel2, sel1, sel2], "geometry", X) update!([sel1, sel2], "master elements", [sel1, sel2]) - calculate_normals!([sel1, sel2], 0.0, Val{1}) + slave_elements = [sel1, sel2] + time = 0.0 + normals, tangents = calculate_normals(slave_elements, time, Val{1}) + update!(slave_elements, "normal", time => normals) return [sel1, sel2], [mel1, mel2] end @@ -28,6 +32,8 @@ end time = 0.0 X1 = sel1("geometry", [-1.0], time) n1 = sel1("normal", [-1.0], time) + println("X1 = ", X1) + println("n1 = ", n1) xi2 = project_from_slave_to_master(mel1, X1, n1, time) @test isapprox(xi2, -1.0) @@ -62,7 +68,10 @@ end mel1 = Element(Seg2, [3, 4]) update!([sel1, mel1], "geometry", X) time = 0.0 - calculate_normals!([sel1], time, Val{1}) + slave_elements = [sel1] + time = 0.0 + normals, tangents = calculate_normals(slave_elements, time, Val{1}) + update!(slave_elements, "normal", time => normals) X2 = mel1("geometry", [-1.0], time) xi = project_from_master_to_slave(sel1, X2, time) diff --git a/test/test_nodal_constraints.jl b/test/test_nodal_constraints.jl deleted file mode 100644 index 7e546b0..0000000 --- a/test/test_nodal_constraints.jl +++ /dev/null @@ -1,89 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -function get_model() - - X = Dict{Int, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [1.0, 1.0], - 4 => [0.0, 1.0]) - - body = Problem(Elasticity, "body", 2) - body.properties.formulation = :plane_stress - body.elements = [Element(Quad4, [1, 2, 3, 4])] - update!(body.elements, "geometry", X) - update!(body.elements, "youngs modulus", 288.0) - update!(body.elements, "poissons ratio", 1/3) - - # boundary conditions - bc_13 = Problem(Dirichlet, "symmetry 13", 2, "displacement") - bc_13.properties.dual_basis = true - bc_13.elements = [Element(Seg2, [1, 2])] - update!(bc_13.elements, "geometry", X) - update!(bc_13.elements, "displacement 2", 0.0) - - bc_23 = Problem(Dirichlet, "symmetry 23", 2, "displacement") - bc_23.properties.dual_basis = true - bc_23.elements = [Element(Seg2, [4, 1])] - update!(bc_23.elements, "geometry", X) - update!(bc_23, "displacement 1", 0.0) - - push!(bc_13.assembly.removed_dofs, 1, 2) - - solver = Solver(Nonlinear, "1x1 plane stress quad4 block") - push!(solver, body, bc_13, bc_23) - - return solver -end - -@testset "test dirichlet spc in point" begin - X = Dict{Int, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [1.0, 1.0], - 4 => [0.0, 1.0]) - solver = get_model() - update!(solver["symmetry 13"], "displacement 1", 0.0) - update!(solver["symmetry 23"], "displacement 2", 0.0) - nodal_bc = Problem(Dirichlet, "dx=0.5", 2, "displacement") - nodal_bc.elements = [Element(Poi1, [3])] - update!(nodal_bc, "geometry", X) - update!(nodal_bc, "displacement 1", 0.5) - update!(nodal_bc, "displacement 2", 0.0) - push!(solver, nodal_bc) - initialize!(solver["symmetry 13"]) - initialize!(solver["symmetry 23"]) - assemble!(solver["symmetry 13"]) - assemble!(solver["symmetry 23"]) - println(sparse(solver["symmetry 13"].assembly.C2)) - println(sparse(solver["symmetry 23"].assembly.C2)) - solver() - pel = nodal_bc.elements[1] - la = pel("lambda", [0.0], 0.0) - info("lambda: $la") - info(solver["body"].assembly.u) - @test isapprox(pel("displacement", [], 0.0), [0.5, 0.0]) -end - -@testset "test nodal point force" begin - X = Dict{Int, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [1.0, 1.0], - 4 => [0.0, 1.0]) - solver = get_model() - update!(solver["symmetry 13"], "displacement 1", 0.0) - update!(solver["symmetry 23"], "displacement 2", 0.0) - point_load = Element(Poi1, [3]) - update!(point_load, "geometry", X) - update!(point_load, "displacement traction force 1", 72.0) - update!(point_load, "displacement traction force 2", 27.0) - push!(solver["body"], point_load) - solver() - @test isapprox(point_load("displacement", [], 0.0), [0.5, 0.0]) -end - diff --git a/test/test_node_dof_mapping.jl b/test/test_node_dof_mapping.jl deleted file mode 100644 index 431256a..0000000 --- a/test/test_node_dof_mapping.jl +++ /dev/null @@ -1,26 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM.Testing - -#= -@testset "find dofs given a set of nodes" begin - nodes = [1, 3] - dim = 3 - dofs = find_dofs_by_nodes(dim, nodes) - @test dofs == [1, 2, 3, 7, 8, 9] -end - -@testset "find nodes given a set of dofs" begin - dofs = [2, 8, 9] - dim = 3 - nodes = find_nodes_by_dofs(dim, dofs) - @test nodes == [1, 3] - - dofs = [2, 12] - dim = 2 - nodes = find_nodes_by_dofs(dim,dofs) - @test nodes == [1, 6] - -end -=# diff --git a/test/test_preprocess.jl b/test/test_preprocess.jl deleted file mode 100644 index f270bd6..0000000 --- a/test/test_preprocess.jl +++ /dev/null @@ -1,89 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Preprocess -using JuliaFEM.Postprocess -using JuliaFEM.Testing - -datadir = first(splitext(basename(@__FILE__))) - -@testset "renumber element nodes" begin - mesh = Mesh() - add_element!(mesh, 1, :Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) - mapping = Dict{Symbol, Vector{Int}}( - :Tet10 => [1, 2, 4, 3, 5, 6, 7, 8, 9, 10]) - reorder_element_connectivity!(mesh, mapping) - @test mesh.elements[1] == [1, 2, 4, 3, 5, 6, 7, 8, 9, 10] - invmapping = Dict{Symbol, Vector{Int}}() - invmapping[:Tet10] = invperm(mapping[:Tet10]) - reorder_element_connectivity!(mesh, invmapping) - @test mesh.elements[1] == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] -end - -@testset "add_nodes! and add_elements!" begin - mesh = Mesh() - dic = Dict(1 => [1.,1.,1.], 2 => [2.,2.,2]) - add_nodes!(mesh, dic) - @test mesh.nodes == dic - - vec = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] - JuliaFEM.Preprocess.add_elements!(mesh,Dict(1=>(:Tet10,vec), - 11=>(:Tet10,vec))) - @test mesh.elements[1] == vec - @test mesh.elements[11] == vec -end - -@testset "find nearest nodes from mesh" begin - meshfile = joinpath(datadir, "block_2d.med") - mesh = aster_read_mesh(meshfile) - create_node_set_from_element_set!(mesh, "LOWER_LEFT", "UPPER_BOTTOM") - # nid 1 coords = (0.0, 0.5), nid 13 coords = (0.0, 0.5) - nid = find_nearest_node(mesh, [0.0, 0.5]; node_set="LOWER_LEFT") - @test first(nid) == 1 - nid = find_nearest_node(mesh, [0.0, 0.5]; node_set="UPPER_BOTTOM") - @test first(nid) == 13 -end - -@testset "test filter by element set" begin - mesh = aster_read_mesh(joinpath(datadir, "block_2d_1elem_quad4.med")) - mesh2 = filter_by_element_set(mesh, :BLOCK) - @test haskey(mesh2.element_sets, :BLOCK) - @test length(mesh2.elements) == 1 -end - -function calculate_volume(mesh_name, eltype) - mesh_file = joinpath(datadir, "primitives.med") - mesh = aster_read_mesh(mesh_file, mesh_name) - elements = create_elements(mesh; element_type=eltype) - V = 0.0 - time = 0.0 - for element in elements - for ip in get_integration_points(element) - detJ = element(ip, time, Val{:detJ}) - detJ > 0 || warn("negative determinant for element $eltype !") - V += ip.weight*detJ - end - end - info("volume of $eltype is $V") - return V -end - -@testset "calculate volume for 1 element models" begin - @test isapprox(calculate_volume("TRIANGLE_TRI3_1", :Tri3), 1/2) - @test isapprox(calculate_volume("TRIANGLE_TRI6_1", :Tri6), 1/2) - @test isapprox(calculate_volume("TRIANGLE_TRI7_1", :Tri7), 1/2) - @test isapprox(calculate_volume("SQUARE_QUAD4_1", :Quad4), 2^2) - @test isapprox(calculate_volume("SQUARE_QUAD8_1", :Quad8), 2^2) - @test isapprox(calculate_volume("SQUARE_QUAD9_1", :Quad9), 2^2) - @test isapprox(calculate_volume("TETRA_TET4_1", :Tet4), 1/6) - @test isapprox(calculate_volume("TETRA_TET10_1", :Tet10), 1/6) -# @test isapprox(calculate_volume("TETRA_TET14_1", :Tet14), 1/6) - @test isapprox(calculate_volume("CUBE_HEX8_1", :Hex8), 2^3) - @test isapprox(calculate_volume("CUBE_HEX20_1", :Hex20), 2^3) - @test isapprox(calculate_volume("CUBE_HEX27_1", :Hex27), 2^3) - @test isapprox(calculate_volume("WEDGE_WEDGE6_1", :Wedge6), 1) -# @test isapprox(calculate_volume("WEDGE_WEDGE15_1", :Wedge15, 1/2)) -# @test isapprox(calculate_volume("PYRAMID_PYRAMID5_1", :Pyramid5, ?)) -# @test isapprox(calculate_volume("PYRAMID_PYRAMID13_1", :Pyramid13, ?)) -end diff --git a/test/test_problem.jl b/test/test_problem.jl deleted file mode 100644 index a6cc7e8..0000000 --- a/test/test_problem.jl +++ /dev/null @@ -1,82 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -@testset "test initialize scalar field problem" begin - el = Element(Seg2, [1, 2]) - pr = Problem(Heat, "heat problem", 1) - push!(pr, el) - initialize!(pr) - @test haskey(el, "temperature") - # one timestep in field "temperature" - @test length(el["temperature"]) == 1 - # this way we access to field at default time t=0.0, it's different than ^! - @test length(el("temperature", 0.0)) == 2 - @test length(last(el, "temperature").data) == 2 -end - -@testset "test initialize vector field problem" begin - el = Element(Seg2, [1, 2]) - pr = Problem(Elasticity, "elasticity problem", 2) - push!(pr, el) - initialize!(pr) - @test haskey(el, "displacement") - @test length(el["displacement"]) == 1 - # this way we access to field at default time t=0.0, it's different than ^! - @test length(el("displacement", 0.0)) == 2 - @test length(last(el, "displacement").data) == 2 -end - -@testset "test initialize boundary problem" begin - el = Element(Seg2, [1, 2]) - pr = Problem(Dirichlet, "bc", 1, "temperature") - push!(pr, el) - initialize!(pr) - @test haskey(el, "lambda") - @test haskey(el, "temperature") -end - -#= -@testset "dict field depending from problems" begin - p1 = Problem(Elasticity, "Body 1", 2) - p2 = Problem(Elasticity, "Body 2", 2) - X = Dict{Int64, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [1.0, 1.0], - 4 => [0.0, 1.0]) - update!([p1, p2], "geometry", 0.0 => X) - @test isapprox(p1("geometry", 0.0)[1], [0.0, 0.0]) - @test isapprox(p2("geometry", 0.0)[1], [0.0, 0.0]) - p1("geometry", 0.0)[1] = [1.0, 2.0] - @test isapprox(p2("geometry", 0.0)[1], [1.0, 2.0]) -end - -@testset "dict field depending from problems" begin - X = Dict{Int64, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [1.0, 1.0], - 4 => [0.0, 1.0], - 5 => [0.0, 2.0], - 6 => [1.0, 2.0], - 7 => [1.0, 3.0], - 8 => [0.0, 3.0]) - p1 = Problem(Elasticity, "Body 1", 2) - p2 = Problem(Elasticity, "Body 2", 2) - e1 = Element(Quad4, [1, 2, 3, 4]) - e2 = Element(Quad4, [5, 6, 7, 8]) - push!(p1, e1) - push!(p2, e2) - update!(p1, "geometry", 0.0 => X) - update!(p2, "geometry", 0.0 => X) - @test isapprox(p1("geometry", 0.0)[1], [0.0, 0.0]) - @test isapprox(p2("geometry", 0.0)[1], [0.0, 0.0]) - p1("geometry", 0.0)[1] = [1.0, 2.0] - @test isapprox(p2("geometry", 0.0)[1], [1.0, 2.0]) - @test isapprox(e1("geometry", 0.0)[1], [1.0, 2.0]) -end - -=# diff --git a/test/test_solvers.jl b/test/test_solvers.jl deleted file mode 100644 index ab66fde..0000000 --- a/test/test_solvers.jl +++ /dev/null @@ -1,78 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -@testset "test linearsolver + xdmf writing" begin - el1 = Element(Quad4, [1, 2, 3, 4]) - el2 = Element(Seg2, [1, 2]) - el3 = Element(Seg2, [3, 4]) - X = Dict( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [1.0, 1.0], - 4 => [0.0, 1.0]) - update!([el1, el2, el3], "geometry", X) - update!(el1, "thermal conductivity", 6.0) - update!(el1, "density", 36.0) - - update!(el2, "heat flux", 0.0 => 0.0) - update!(el2, "heat flux", 1.0 => 600.0) - - problem = Problem(Heat, "test problem", 1) - problem.properties.formulation = "2D" - push!(problem.elements, el1, el2) - - update!(el3, "temperature 1", 0.0) - - bc = Problem(Dirichlet, "fixed", 1, "temperature") - - push!(bc.elements, el3) - - # Create a solver for a set of problems - solver = Solver(Linear, "solve heat problem") - push!(solver, problem, bc) - - # Solve problem at time t=1.0 and update fields - solver.time = 1.0 - solver.xdmf = Xdmf() - solver() - - # Postprocess. - # Interpolate temperature field along boundary of Γ₁ at time t=1.0 - xi = (0.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 - -@testset "problem not found from solver" begin - s = Solver(Linear, "demo solver") - @test_throws KeyError getindex(s, "not_found") -end - -@testset "automatic determination of problem dimension if not spesified" begin - s = Solver(Linear, "demo solver") - p = Problem(Elasticity, "demo problem", 2) - push!(s, p) - get_field_assembly(s) - @test s.ndofs == 0 - add!(p.assembly.K, [4], [4], reshape([4.0],1,1)) - get_field_assembly(s) - @test s.ndofs == 4 -end - -@testset "test for error when overdetermined system and requesting boundary assembly" begin - s = Solver(Linear, "demo solver") - @test_throws AssertionError get_boundary_assembly(s) # ndofs = 0 - p1 = Problem(Dirichlet, "bc1", 2, "displacement") - p2 = Problem(Dirichlet, "bc2", 2, "displacement") - # third dofs constrained - add!(p1.assembly.C2, [3], [3], reshape([1.0],1,1)) - add!(p2.assembly.C2, [3], [4], reshape([1.0],1,1)) - s.ndofs = 4 - push!(s, p1, p2) - @test_throws ErrorException get_boundary_assembly(s) -end diff --git a/test/test_solvers_geometry_not_defined.jl b/test/test_solvers_geometry_not_defined.jl deleted file mode 100644 index 0501074..0000000 --- a/test/test_solvers_geometry_not_defined.jl +++ /dev/null @@ -1,14 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -@testset "get nodal field from boundary condition if geometry is not defined" begin - bc = Problem(Dirichlet, "bc without geometry", 3, "displacement") - bc.elements = [Element(Poi1, [1])] - @test bc("geometry", 0.0) == nothing - s = Solver(Linear, "test solver") - s.problems = [bc] - @test length(s("geometry", 0.0)) == 0 -end diff --git a/test/test_sparse.jl b/test/test_sparse.jl deleted file mode 100644 index edd0ecf..0000000 --- a/test/test_sparse.jl +++ /dev/null @@ -1,47 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM -using JuliaFEM.Testing - -@testset "Add to SparseMatrixCOO" begin - A = SparseMatrixCOO() - A2 = reshape(collect(1:9), 3, 3) - add!(A, sparse(A2)) - @test isapprox(full(A), full(A2)) -end - -@testset "Add to SparseVectorCOO" begin - b = SparseVectorCOO() - b2 = collect(1:3) - add!(b, sparse(b2)) - @test isapprox(full(b), full(b2)) -end - -@testset "Failure to add data to sparse vector due dimensino mismatch" begin - b = SparseVectorCOO() - @test_throws ErrorException add!(b, [1, 2], [1.0, 2.0, 3.0]) -end - -@testset "Test combining of SparseMatrixCOO" begin - k = convert(Matrix{Float64}, reshape(collect(1:9), 3, 3)) - dofs1 = [1, 2, 3] - dofs2 = [2, 3, 4] - A = SparseMatrixCOO() - add!(A, dofs1, dofs1, k) - add!(A, dofs2, dofs2, k) - A1 = full(A) - optimize!(A) - A2 = full(A) - @test isapprox(A1, A2) -end - -@testset "resize of sparse matrix and sparse vector" begin - A = sparse(rand(3, 3)) - B = resize_sparse(A, 4, 4) - @test size(B) == (4, 4) - a = sparse(rand(3)) - b = resize_sparsevec(a, 4) - @test size(b) == (4, ) -end -