diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index bfbfa35..386e2c4 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -30,9 +30,6 @@ module Testing end -include("io.jl") -export Xdmf, h5file, xmffile, xdmf_filter, new_dataitem, get_temporal_collection - include("fields.jl") export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment @@ -42,7 +39,8 @@ export AbstractPoint, Point, IntegrationPoint, IP, Node ### ELEMENTS ### include("elements.jl") # common element routines export Node, AbstractElement, Element, update!, get_connectivity, get_basis, - get_dbasis, inside, get_local_coordinates + get_dbasis, inside, get_local_coordinates, get_element_type, + filter_by_element_type, get_element_id include("elements_lagrange.jl") # Continuous Galerkin (Lagrange) elements export get_reference_coordinates, @@ -99,6 +97,8 @@ export calculate_normals, calculate_normals!, project_from_slave_to_master, project_from_master_to_slave, Mortar, get_slave_elements, get_polygon_clip +include("io.jl") +export Xdmf, h5file, xmffile, xdmf_filter, new_dataitem, update_xdmf!, save! ### ASSEMBLY + SOLVE ### include("assembly.jl") @@ -141,17 +141,11 @@ export aster_create_elements, parse_aster_med_file, is_aster_mail_keyword, end module Postprocess - include("postprocess_utils.jl") export calc_nodal_values!, get_nodal_vector, get_nodal_dict, copy_field!, calculate_area, calculate_center_of_mass, calculate_second_moment_of_mass, extract - -include("postprocess_xdmf.jl") -export XDMF, xdmf_new_result!, xdmf_save_field!, xdmf_save!, DataFrame - end -export Postprocessor module Abaqus include("abaqus.jl") diff --git a/src/abaqus.jl b/src/abaqus.jl index 2002a80..cebb5b5 100644 --- a/src/abaqus.jl +++ b/src/abaqus.jl @@ -564,8 +564,8 @@ function process_output_request(model::Model, solver::Solver, output_request::Ab haskey(code_mapping, code) || continue field_name = code_mapping[code] abbr = get(abbr_mapping, code, code) - table = solver(DataFrame, field_name, abbr, solver.time) - push!(tables, table) + #table = solver(DataFrame, field_name, abbr, solver.time) + #push!(tables, table) end length(tables) != 0 || continue results = join(tables..., on=:NODE, kind=:outer) diff --git a/src/api.jl b/src/api.jl deleted file mode 100644 index 8d72e42..0000000 --- a/src/api.jl +++ /dev/null @@ -1,5 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -include("api_types.jl") -include("api_functions.jl") diff --git a/src/io.jl b/src/io.jl index 5ba64fe..be6022c 100644 --- a/src/io.jl +++ b/src/io.jl @@ -8,6 +8,8 @@ type Xdmf name :: String xml :: XMLElement hdf :: HDF5File + hdf_counter :: Int + format :: String end function Xdmf() @@ -23,7 +25,7 @@ function xmffile(xdmf::Xdmf) end """ Initialize a new Xdmf object. """ -function Xdmf(name::String; overwrite=false) +function Xdmf(name::String; version="3.0", overwrite=false) xdmf = new_element("Xdmf") h5file = "$name.h5" xmlfile = "$name.xmf" @@ -47,16 +49,14 @@ function Xdmf(name::String; overwrite=false) end set_attribute(xdmf, "xmlns:xi", "http://www.w3.org/2001/XInclude") - set_attribute(xdmf, "Version", "2.1") + set_attribute(xdmf, "Version", version) flag = isfile(h5file) ? "r+" : "w" hdf = h5open(h5file, flag) - return Xdmf(name, xdmf, hdf) + return Xdmf(name, xdmf, hdf, 1, "HDF") end """ Return the basic structure of Xdmf document. Creates a new TemporalCollection if not found. - Basic structure for XML part of Xdmf file is - @@ -64,7 +64,6 @@ Basic structure for XML part of Xdmf file is - """ function get_temporal_collection(xdmf::Xdmf) domain = find_element(xdmf.xml, "Domain") @@ -194,6 +193,13 @@ function traverse(xdmf::Xdmf, x::XMLElement, attr_name::String) if '/' in attr_name items = split(attr_name, '/') new_item = xdmf_filter(childs, first(items)) + if new_item == nothing + info("traverse: childs:") + for child in childs + info(LightXML.name(child)) + end + error("traverse: failed, items = $items, xdmf_filter not find child") + end new_path = join(items[2:end], '/') return traverse(xdmf, new_item, new_path) end @@ -211,11 +217,14 @@ function read(xdmf::Xdmf, path::String) result = traverse(xdmf, xdmf.xml, path) if endswith(path, "DataItem") format = attribute(result, "Format"; required=true) - @assert format == "HDF" - h5file, path = map(String, split(content(result), ':')) - h5file = dirname(xdmf.name) * "/" * h5file - isfile(h5file) || throw("Xdmf: h5 file $h5file not found!") - return read(xdmf.hdf, path) + if format == "HDF" + h5file, path = map(String, split(content(result), ':')) + h5file = dirname(xdmf.name) * "/" * h5file + isfile(h5file) || throw("Xdmf: h5 file $h5file not found!") + return read(xdmf.hdf, path) + else + error("Read from Xdmf, reading from $format not implemented") + end else return result end @@ -228,14 +237,14 @@ function save!(xdmf::Xdmf) save_file(doc, xmffile(xdmf)) end -function new_dataitem{T,N}(xdmf::Xdmf, path::String, data::Array{T,N}; format="HDF") +function new_dataitem{T,N}(xdmf::Xdmf, path::String, data::Array{T,N}) dataitem = new_element("DataItem") datatype = replace("$T", "64", "") dimensions = join(reverse(size(data)), " ") set_attribute(dataitem, "DataType", datatype) set_attribute(dataitem, "Dimensions", dimensions) - set_attribute(dataitem, "Format", format) - if format == "HDF" + set_attribute(dataitem, "Format", xdmf.format) + if xdmf.format == "HDF" hdf = basename(h5file(xdmf)) if exists(xdmf.hdf, path) info("Xdmf: $path already existing in h5 file, not overwriting.") @@ -243,9 +252,194 @@ function new_dataitem{T,N}(xdmf::Xdmf, path::String, data::Array{T,N}; format="H write(xdmf.hdf, path, data) end add_text(dataitem, "$hdf:$path") - elseif format == "XML" - add_text(dataitem, strip(string(data), ['[', ']'])) + elseif xdmf.format == "XML" + text_data = string(data') + text_data = strip(text_data, ['[', ']']) + text_data = replace(text_data, ';', '\n') + text_data = "\n" * text_data * "\n" + add_text(dataitem, text_data) + else + error("Unsupported Xdmf big data format $(xdmf.format)") end return dataitem end +""" Create a new DataItem element, hdf path automatically determined. """ +function new_dataitem{T,N}(xdmf::Xdmf, data::Array{T,N}) + if xdmf.format == "XML" + # Path can be whatever as XML format does not store to HDF at all + return new_dataitem(xdmf, "/whatever", data) + else + debug("Determining path for HDF file automatically.") + path = "/DataItem_$(xdmf.hdf_counter)" + while exists(xdmf.hdf, path) + xdmf.hdf_counter += 1 + path = "/DataItem_$(xdmf.hdf_counter)" + end + debug("HDF path automatically determined to be $path") + return new_dataitem(xdmf, path, data) + end +end + +global const xdmf_element_mapping = Dict( + "Poi1" => "Polyvertex", + "Seg2" => "Polyline", + "Tri3" => "Triangle", + "Quad4" => "Quadrilateral", + "Tet4" => "Tetrahedron", + "Pyramid5" => "Pyramid", + "Wedge6" => "Wedge", + "Hex8" => "Hexahedron", + "Seg3" => "Edge_3", + "Tri6" => "Tri_6", + "Quad8" => "Quad_8", + "Tet10" => "Tet_10", + "Pyramid13" => "Pyramid_13", + "Wedge15" => "Wedge_15", + "Hex20" => "Hex_20") + +""" Write new fields to Xdmf file. + +Examples +-------- + +To write displacement and temperature fields from p1 at time t=0.0: + +julia> update_xdmf!(p1, 0.0, ["displacement", "temperature"]) +""" +function update_xdmf!(xdmf::Xdmf, problem::Problem, time::Float64, fields::Vector) + + info("Xdmf: storing fields $fields of problem $(problem.name) at time $time") + + # 1. find domain + xml = xdmf.xml + domain = find_element(xml, "Domain") + if domain == nothing + info("Xdmf: Domain not found, creating.") + domain = new_child(xml, "Domain") + else + debug("Xdmf: Domain already defined, skipping.") + end + + # 2. find for TemporalCollection + temporal_collection = find_element(domain, "Grid") + if temporal_collection == nothing + info("Xdmf: Temporal collection not found, creating.") + temporal_collection = new_child(domain, "Grid") + set_attribute(temporal_collection, "GridType", "Collection") + set_attribute(temporal_collection, "Name", "Time") + set_attribute(temporal_collection, "CollectionType", "Temporal") + else + debug("Xdmf: Temporal collection found, skipping.") + end + + # 2.1 make sure that Grid element we found really is TemporalCollection + collection_type = attribute(temporal_collection, "CollectionType"; required=true) + @assert collection_type == "Temporal" + + # 3. find for SpatialCollection at given time + spatial_collection = nothing + spatial_collection_exists = false + for spatial_collection in get_elements_by_tagname(temporal_collection, "Grid") + time_element = find_element(spatial_collection, "Time") + time_value = parse(attribute(time_element, "Value"; required=true)) + if isapprox(time_value, time) + info("Xdmf: SpatialCollection for time $time already exists.") + spatial_collection_exists = true + break + end + end + + if !spatial_collection_exists + info("Xdmf: SpatialCollection for time $time not found, creating.") + spatial_collection = new_child(temporal_collection, "Grid") + set_attribute(spatial_collection, "GridType", "Collection") + set_attribute(spatial_collection, "Name", "Problems") + set_attribute(spatial_collection, "CollectionType", "Spatial") + time_element = new_child(spatial_collection, "Time") + set_attribute(time_element, "Value", time) + end + + # 3.1 make sure that Grid element we found really is SpatialCollection + collection_type = attribute(spatial_collection, "CollectionType"; required=true) + @assert collection_type == "Spatial" + + for frame in get_elements_by_tagname(spatial_collection, "Grid") + frame_name = attribute(frame, "Name") + if frame_name == problem.name + warn("Xdmf: Already found Grid with name $frame_name for time $time, skipping.") + return + end + end + + frame_name = problem.name + info("Xdmf: Creating Grid for problem $frame_name") + frame = new_child(spatial_collection, "Grid") + set_attribute(frame, "Name", frame_name) + + # 4. save geometry + X_dict = problem("geometry", time) + node_ids = sort(collect(keys(X_dict))) + node_mapping = Dict(j => i for (i, j) in enumerate(node_ids)) + X_array = hcat([X_dict[nid] for nid in node_ids]...) + ndim, nnodes = size(X_array) + geom_type = (ndim == 2 ? "XY" : "XYZ") + info("Xdmf: Creating geometry, type = $geom_type, number of nodes = $nnodes") + X_dataitem = new_dataitem(xdmf, X_array) + geometry = new_child(frame, "Geometry") + set_attribute(geometry, "Type", geom_type) + add_child(geometry, X_dataitem) + + # 5. save topology + all_elements = get_elements(problem) + nelements = length(all_elements) + element_types = unique(map(get_element_type, all_elements)) + nelement_types = length(element_types) + 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) + nelements = length(elements) + info("Xdmf: $nelements elements of type $element_type") + sort!(elements, by=get_element_id) + element_ids = map(get_element_id, elements) + element_conn = map(element -> [node_mapping[j]-1 for j in get_connectivity(element)], elements) + element_conn = hcat(element_conn...) + element_code = split(string(element_type), ".")[end] + topology_dataitem = new_dataitem(xdmf, element_conn) + + topology = new_child(frame, "Topology") + set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code]) + set_attribute(topology, "NumberOfElements", length(elements)) + add_child(topology, topology_dataitem) + end + + # 6. save requested fields + for field_name in fields + field_dict = problem(field_name, time) + field_center = "Node" + field_node_ids = sort(collect(keys(field_dict))) + @assert node_ids == field_node_ids + field_dim = length(field_dict[first(field_node_ids)]) + if field_dim == 2 + info("Xdmf: Field dimension = 2, extending to 3") + for nid in field_node_ids + field_dict[nid] = [field_dict[nid]; 0.0] + end + field_dim == 3 + end + field_type = Dict(1 => "Scalar", 3 => "Vector", 6 => "Tensor6")[field_dim] + info("Xdmf: Saving field $field_name, type = $field_type, dimension = $field_dim, center = $field_center") + + field_array = hcat([field_dict[nid] for nid in field_node_ids]...) + field_dataitem = new_dataitem(xdmf, field_array) + attribute = new_child(frame, "Attribute") + set_attribute(attribute, "Name", ucfirst(field_name)) + set_attribute(attribute, "Center", field_center) + set_attribute(attribute, "AttributeType", field_type) + add_child(attribute, field_dataitem) + end + + save!(xdmf) + info("Xdmf: all done.") +end diff --git a/src/postprocess_xdmf.jl b/src/postprocess_xdmf.jl deleted file mode 100644 index 3ad67ab..0000000 --- a/src/postprocess_xdmf.jl +++ /dev/null @@ -1,194 +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 LightXML - -# element codes: http://www.paraview.org/pipermail/paraview/2013-July/028859.html -# > from ./VTK/ThirdParty/xdmf2/vtkxdmf2/libsrc/XdmfTopology.h -# > -# > // Topologies -# > #define XDMF_NOTOPOLOGY 0x0 -# > #define XDMF_POLYVERTEX 0x1 -# > #define XDMF_POLYLINE 0x2 -# > #define XDMF_POLYGON 0x3 -# > #define XDMF_TRI 0x4 -# > #define XDMF_QUAD 0x5 -# > #define XDMF_TET 0x6 -# > #define XDMF_PYRAMID 0x7 -# > #define XDMF_WEDGE 0x8 -# > #define XDMF_HEX 0x9 -# > #define XDMF_EDGE_3 0x0022 -# > #define XDMF_TRI_6 0x0024 -# > #define XDMF_QUAD_8 0x0025 -# > #define XDMF_QUAD_9 0x0023 -# > #define XDMF_TET_10 0x0026 -# > #define XDMF_PYRAMID_13 0x0027 -# > #define XDMF_WEDGE_15 0x0028 -# > #define XDMF_WEDGE_18 0x0029 -# > #define XDMF_HEX_20 0x0030 -# > #define XDMF_HEX_24 0x0031 -# > #define XDMF_HEX_27 0x0032 -# > #define XDMF_MIXED 0x0070 -# > #define XDMF_2DSMESH 0x0100 -# > #define XDMF_2DRECTMESH 0x0101 -# > #define XDMF_2DCORECTMESH 0x0102 -# > #define XDMF_3DSMESH 0x1100 -# > #define XDMF_3DRECTMESH 0x1101 -# > #define XDMF_3DCORECTMESH 0x1102 - -get_xdmf_element_code(element::Element{Poi1}) = 0x0001 -get_xdmf_element_code(element::Element{Seg2}) = 0x0002 -get_xdmf_element_code(element::Element{Seg3}) = 0x0003 -get_xdmf_element_code(element::Element{Tri3}) = 0x0004 -get_xdmf_element_code(element::Element{Quad4}) = 0x0005 -get_xdmf_element_code(element::Element{Tet4}) = 0x0006 -get_xdmf_element_code(element::Element{Hex8}) = 0x0009 -get_xdmf_element_code(element::Element{Tet10}) = 0x0026 - -type XDMF - dimension :: Int - use_hdf :: Bool - xdoc :: XMLDocument - domain :: XMLElement - temporal_collection :: XMLElement - current_grid - permutation :: Vector{Int} -end - -function XDMF() - xdoc = XMLDocument() - xroot = create_root(xdoc, "Xdmf") - set_attribute(xroot, "xmlns:xi", "http://www.w3.org/2001/XInclude") - set_attribute(xroot, "Version", "2.1") - domain = new_child(xroot, "Domain") - temporal_collection = new_child(domain, "Grid") - set_attribute(temporal_collection, "CollectionType", "Temporal") - set_attribute(temporal_collection, "GridType", "Collection") - set_attribute(temporal_collection, "Name", "Collection") - return XDMF(3, false, xdoc, domain, temporal_collection, Union{}, []) -end - -function xdmf_new_result!(xdmf::XDMF, elements::Vector, time) - grid = new_child(xdmf.temporal_collection, "Grid") - set_attribute(grid, "Name", "Grid") - time_ = new_child(grid, "Time") - set_attribute(time_, "Value", time) - xdmf.current_grid = grid - - # 1. calculate permutation - nids = Set() - X = Dict{Int64, Vector{Float64}}() - for element in elements - conn = get_connectivity(element) - push!(nids, conn...) - X_el = element["geometry"](time) - for (i, c) in enumerate(conn) - X[c] = X_el[i] - end - end - xdmf.permutation = sort(collect(nids)) - iperm = Dict{Int64, Int64}() - for (i, j) in enumerate(xdmf.permutation) - iperm[j] = i - end - - # 2. write nodes - geometry = new_child(grid, "Geometry") - set_attribute(geometry, "Type", xdmf.dimension == 3 ? "XYZ" : "XY") - dataitem = new_child(geometry, "DataItem") - set_attribute(dataitem, "DataType", "Float") - set_attribute(dataitem, "Format", "XML") - #set_attribute(dataitem, "Precision", 8) - s = [] - ndim = 0 - for i in xdmf.permutation - ndim += length(X[i]) - push!(s, join(round(X[i], 5), " ")) - end - set_attribute(dataitem, "Dimensions", ndim) - add_text(dataitem, "\n"*join(s, "\n")*"\n") - - # 3. write elements - topology = new_child(grid, "Topology") - set_attribute(topology, "TopologyType", "Mixed") - set_attribute(topology, "NumberOfElements", length(elements)) - dataitem = new_child(topology, "DataItem") - set_attribute(dataitem, "Format", "XML") - set_attribute(dataitem, "DataType", "Int") -# set_attribute(dataitem, "Precision", 8) - s = [] - eldim = 0 - for element in elements - eltype = get_xdmf_element_code(element) - # note: id numbers start from 0 in Xdmf - conn = [iperm[j] for j in get_connectivity(element)] - 1 - data = [eltype; conn] - eldim += length(data) - push!(s, join(data, " ")) - end - set_attribute(dataitem, "Dimensions", eldim) - add_text(dataitem, "\n"*join(s, "\n")*"\n") -end - -function xdmf_save_field!(xdmf, elements::Vector, time, field_name; field_type="Scalar", debug=false) - f = Dict() - field_dim = 0 - for element in elements - haskey(element, field_name) || continue - g = element[field_name](time) - conn = get_connectivity(element) - for (i, c) in enumerate(conn) - gi = g[i] - if (field_type == "Vector") && (length(gi) < 3) - # paraview goes crazy if 2d model with 2d displacement vector - gi = [gi; 0.0] - end - if field_dim == 0 - field_dim = length(gi) - end - field_dim == length(gi) || error("several dimensions in field, dim = $field_dim.") - f[c] = gi - end - end - - if length(f) == 0 - warn("xdmf_save_field!(): field $field_name was not found from set of elements") - return - end - - attribute = new_child(xdmf.current_grid, "Attribute") - set_attribute(attribute, "Center", "Node") - set_attribute(attribute, "Name", ucfirst(field_name)) - set_attribute(attribute, "Type", field_type) - dataitem = new_child(attribute, "DataItem") - set_attribute(dataitem, "DataType", "Float") - set_attribute(dataitem, "Format", "XML") - #set_attribute(dataitem, "Precision", 8) - debug && info("field dim = $field_dim") - debug && info(f) - s = [] - dim = 0 - for i in xdmf.permutation - gi = zeros(field_dim) - if haskey(f, i) - gi = f[i] - end - push!(s, join(round(gi, 5), " ")) - dim += length(gi) - end - set_attribute(dataitem, "Dimensions", dim) - add_text(dataitem, "\n"*join(s, "\n")*"\n") -end - -function xdmf_save_field!(xdmf, problem::Problem, time, field_name; field_type="Scalar") - xdmf_save_field!(xdmf, problem.elements, time, field_name; field_type=field_type) -end - -function xdmf_new_result!(xdmf, problem::Problem, time) - xdmf_new_result!(xdmf, problem.elements, time) -end - -function xdmf_save!(xdmf, filename) - save_file(xdmf.xdoc, filename) -end diff --git a/src/problems.jl b/src/problems.jl index 1e2a5f3..a9b04ee 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -89,6 +89,7 @@ type Problem{P<:AbstractProblem} dofmap :: Dict{Element, Vector{Int64}} # connects element local dofs to global dofs assembly :: Assembly fields :: Dict{AbstractString, Field} + postprocess_fields :: Vector{String} properties :: P end @@ -103,10 +104,10 @@ julia> prob2 = Problem(Elasticity, 3) """ function Problem{P<:FieldProblem}(::Type{P}, name::AbstractString, dimension::Int64) - return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), Dict(), P()) + return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), Dict(), Vector(), P()) end function Problem{P<:FieldProblem}(::Type{P}, dimension::Int64) - return Problem{P}("$P problem", dimension, "none", [], Dict(), Assembly(), Dict(), P()) + return Problem(P, "$P problem", dimension) end """ Construct a new boundary problem. @@ -119,28 +120,27 @@ julia> bc1 = Problem(Dirichlet, "support", 3, "displacement") solver. """ function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name) - return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), P()) + return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), Vector(), P()) end function Problem{P<:BoundaryProblem}(::Type{P}, main_problem::Problem) name = "$P problem" dimension = get_unknown_field_dimension(main_problem) parent_field_name = get_unknown_field_name(main_problem) - return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), P()) + return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), Vector(), P()) end -function get_formulation_type{P<:FieldProblem}(problem::Problem{P}) +function get_formulation_type(problem::Problem) return :incremental end -function get_formulation_type{P<:BoundaryProblem}(problem::Problem{P}) - return :incremental +function get_unknown_field_name{P<:BoundaryProblem}(::Type{P}) + return "lambda" end function get_assembly(problem) return problem.assembly end - """ Initialize element ready for calculation. """ function initialize!(problem::Problem, element::Element, time::Float64) field_name = get_unknown_field_name(problem) @@ -175,7 +175,7 @@ function initialize!(problem::Problem, time::Float64=0.0) end """ Update problem solution vector for assembly. """ -function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector; verbose=false) +function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector) # resize & fill with zeros vectors if length mismatch with current solution @@ -186,7 +186,7 @@ function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector; ve end if length(la) != length(assembly.la) - info("resizing lagrange multipliers vector u") + info("resizing lagrange multiplier vector la") resize!(assembly.la, length(la)) fill!(assembly.la, 0.0) end @@ -199,15 +199,12 @@ function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector; ve assembly.la_prev = copy(assembly.la) if get_formulation_type(problem) == :total - verbose && info("$(problem.name): total formulation, replacing solution vector with new values") assembly.u = u assembly.la = la elseif get_formulation_type(problem) == :incremental - verbose && info("$(problem.name): incremental formulation, adding increment to solution vector") assembly.u += u assembly.la = la elseif get_formulation_type(problem) == :forwarddiff - verbose && info("$(problem.name): forwarddiff formulation, adding increment to solution vector and reaction force vector") assembly.u += u assembly.la += la else @@ -259,13 +256,12 @@ end function update!{P<:BoundaryProblem}(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64) u, la = get_global_solution(problem, assembly) parent_field_name = get_parent_field_name(problem) # displacement - field_name = get_unknown_field_name(problem) # reaction force - # update solution u and reaction force λ for boundary elements + field_name = get_unknown_field_name(problem) # lambda + # update solution and lagrange multipliers for boundary elements for element in elements connectivity = get_connectivity(element) update!(element, parent_field_name, time => u[connectivity]) - # FIXME - update!(element, field_name, time => -la[connectivity]) + update!(element, field_name, time => la[connectivity]) end end diff --git a/src/problems_contact.jl b/src/problems_contact.jl index bac0010..46d5871 100644 --- a/src/problems_contact.jl +++ b/src/problems_contact.jl @@ -58,22 +58,6 @@ function Contact() default_fields) end -function get_unknown_field_name(problem::Problem{Contact}) - return "reaction force" -end - -function get_formulation_type(problem::Problem{Contact}) - #= - if problem.properties.use_forwarddiff - return :forwarddiff - else - return :incremental - end - =# - return :incremental - #return :forwarddiff -end - function assemble!(problem::Problem{Contact}, time::Real) if problem.properties.dimension == -1 problem.properties.dimension = dim = size(first(problem.elements), 1) diff --git a/src/problems_contact_2d.jl b/src/problems_contact_2d.jl index cac52d3..327e2f1 100644 --- a/src/problems_contact_2d.jl +++ b/src/problems_contact_2d.jl @@ -77,7 +77,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::T nsl = length(slave_element) X1 = slave_element("geometry", time) u1 = slave_element("displacement", time) - la1 = slave_element("reaction force", time) + la1 = slave_element("lambda", time) n1 = slave_element("normal", time) t1 = slave_element("tangent", time) x1 = X1 + u1 @@ -160,8 +160,8 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::T Te += w*reshape(kron(N2, n_s, Phi), 2, 4) He += w*reshape(kron(N1, t_s, Phi), 2, 4) ge += w*Phi*dot(n_s, x_m-x_s) - ce += w*N1*dot(n_s, -la_s) - Rn += w*dot(n_s, -la_s) + ce += w*N1*dot(n_s, la_s) + Rn += w*dot(n_s, la_s) contact_area += w contact_error += 1/2*w*dot(n_s, x_s-x_m)^2 diff --git a/src/problems_contact_3d.jl b/src/problems_contact_3d.jl index 191b790..7b67b25 100644 --- a/src/problems_contact_3d.jl +++ b/src/problems_contact_3d.jl @@ -117,7 +117,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time u1 = slave_element("displacement", time) x1 = X1 + u1 n1 = slave_element("normal", time) - la = slave_element("reaction force", time) + la = slave_element("lambda", time) Q3 = create_rotation_matrix(slave_element, time) @@ -277,7 +277,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time #u1 = sub_slave_element("displacement", time) #x1 = X1 + u1 n1 = sub_slave_element("normal", time) - #la = sub_slave_element("reaction force", time) + #la = sub_slave_element("lambda", time) # create auxiliary plane xi = mean(get_reference_coordinates(sub_slave_element)) @@ -599,16 +599,16 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T end complementarity_condition[j] = contact_pressure[j] - weighted_gap[j] - if complementarity_condition[j][1] < 0.0 - is_inactive[j] = 1 - is_active[j] = 0 - is_slip[j] = 0 - is_stick[j] = 0 - else + if complementarity_condition[j][1] > 0.0 is_inactive[j] = 0 is_active[j] = 1 is_slip[j] = 1 is_stick[j] = 0 + else + is_inactive[j] = 1 + is_active[j] = 0 + is_slip[j] = 0 + is_stick[j] = 0 end end @@ -670,3 +670,18 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T problem.assembly.g = g end + +function postprocess!(problem::Problem{Contact}, time::Float64, ::Type{Val{Symbol("contact pressure")}}) + n = problem("normal", time) + la = problem("lambda", time) + node_ids = keys(n) + cp = Dict(nid => dot(n[nid], la[nid]) for nid in node_ids) + # FIXME: have to define zero contact pressure & lambda to master elements + # elements because interface.elements = [slave_elements; master_elements] + for nid in keys(la) + if !haskey(cp, nid) + cp[nid] = 0.0 + end + end + update!(problem, "contact pressure", time => cp) +end diff --git a/src/problems_dirichlet.jl b/src/problems_dirichlet.jl index 6dda4a4..b9ff261 100644 --- a/src/problems_dirichlet.jl +++ b/src/problems_dirichlet.jl @@ -15,10 +15,6 @@ function Dirichlet() Dirichlet(:incremental, false, false, 1) end -function get_unknown_field_name(::Type{Dirichlet}) - return "reaction force" -end - function get_formulation_type(problem::Problem{Dirichlet}) return problem.properties.formulation end @@ -142,53 +138,8 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, end -#= - -function assemble!(assembly::Assembly, problem::Problem{DirichletProblem}, element::Element, time::Real) - - # get dimension and name of PARENT field - field_dim = problem.parent_field_dim - field_name = problem.parent_field_name - - gdofs = get_gdofs(element, field_dim) - for ip in get_integration_points(element, Val{2}) - w = ip.weight - J = get_jacobian(element, ip, time) - JT = transpose(J) - if size(JT, 2) == 1 # plane problem - w *= norm(JT) - else - w *= norm(cross(JT[:,1], JT[:,2])) - end - N = element(ip, time) - A = w*N'*N - - if haskey(element, field_name) - # add all dimensions at once if defined - # element["blaa"] = 0.0 - # or - # element["blaa"] = Vector{Float64}[[0.1, 0.2], [0.3, 0.4]] - g = element(field_name, ip, time) - if length(g) != length(N) - g = g*ones(length(N)) - end - for i=1:field_dim - ldofs = gdofs[i:field_dim:end] - add!(assembly.C1, ldofs, ldofs, A) - add!(assembly.C2, ldofs, ldofs, A) - end - add!(assembly.g, gdofs, w*g*N) - end - - for i=1:field_dim - if haskey(element, field_name*" $i") - g = element(field_name*" $i", ip, time) - ldofs = gdofs[i:field_dim:end] - add!(assembly.C1, ldofs, ldofs, A) - add!(assembly.C2, ldofs, ldofs, A) - add!(assembly.g, ldofs, w*g*N) - end - end - end +function postprocess!(problem::Problem{Dirichlet}, time::Float64, ::Type{Val{Symbol("reaction force")}}) + la = problem("lambda", time) + rf = Dict(nid => -lai for (nid, lai) in la) + update!(problem, "reaction force", time => rf) end -=# diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index b01184a..da28467 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -534,3 +534,80 @@ function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity}, end return Km, Kg, f end + +""" Return strain tensor. """ +function get_strain_tensor(problem, element, ip, time) + gradu = element("displacement", ip, time, Val{:Grad}) + eps = 0.5*(gradu' + gradu) + return eps +end + +""" Return stress tensor. """ +function get_stress_tensor(problem, element, ip, time) + eps = get_strain_tensor(problem, element, ip, time) + E = element("youngs modulus", ip, time) + nu = element("poissons ratio", ip, time) + mu = E/(2.0*(1.0+nu)) + la = E*nu/((1.0+nu)*(1.0-2.0*nu)) + S = la*trace(eps)*I + 2.0*mu*eps + return S +end + +""" Return stain vector in "ABAQUS" order 11, 22, 33, 12, 23, 13. """ +function get_strain_vector(problem, element, ip, time) + eps = get_strain_tensor(problem, element, ip, time) + return [eps[1,1], eps[2,2], eps[3,3], eps[1,2], eps[2,3], eps[1,3]] +end + +""" Return stress vector in "ABAQUS" order 11, 22, 33, 12, 23, 13. """ +function get_stress_vector(problem, element, ip, time) + S = get_stress_tensor(problem, element, ip, time) + return [S[1,1], S[2,2], S[3,3], S[1,2], S[2,3], S[1,3]] +end + +""" Make least squares fit for some field to nodes. """ +function lsq_fit(problem, elements, field, time) + A = SparseMatrixCOO() + b = SparseMatrixCOO() + volume = 0.0 + for element in elements + gdofs = get_connectivity(element) + for ip in get_integration_points(element) + detJ = element(ip, time, Val{:detJ}) + w = ip.weight*detJ + N = element(ip, time) + f = field(problem, element, ip, time) + add!(A, gdofs, gdofs, w*kron(N', N)) + for i=1:length(f) + add!(b, gdofs, w*f[i]*N, i) + end + volume += w + end + end + debug("Mass matrix for least-squares fit is assembled. Total volume to fit: $volume") + A = sparse(A) + b = sparse(b) + A = 1/2*(A + A') + + nz = get_nonzero_rows(A) + F = ldltfact(A[nz,nz]) + + x = F \ b[nz, :] + + nodal_values = Dict(node_id => vec(full(x[idx,:])) for (idx, node_id) in enumerate(nz)) + return nodal_values +end + + +""" Postprocessing, extrapolate strain to nodes using least-squares fit. """ +function postprocess!(problem::Problem{Elasticity}, time::Float64, ::Type{Val{:strain}}) + elements = get_elements(problem) + strain = lsq_fit(problem, elements, get_strain_vector, time) + update!(elements, "strain", time => strain) +end + +function postprocess!(problem::Problem{Elasticity}, time::Float64, ::Type{Val{:stress}}) + elements = get_elements(problem) + stress = lsq_fit(problem, elements, get_stress_vector, time) + update!(elements, "stress", time => stress) +end diff --git a/src/problems_mortar.jl b/src/problems_mortar.jl index 351c12a..73e8178 100644 --- a/src/problems_mortar.jl +++ b/src/problems_mortar.jl @@ -52,18 +52,6 @@ function Mortar() return Mortar(-1, false, false, false, false, Inf, true, true, true, 0.0, 1.0e-9, default_fields) end -function get_unknown_field_name(problem::Problem{Mortar}) - return "reaction force" -end - -function get_formulation_type(problem::Problem{Mortar}) - if problem.properties.use_forwarddiff - return :forwarddiff - else - return :incremental - end -end - function assemble!(problem::Problem{Mortar}, time::Float64) if length(problem.elements) == 0 warn("No elements defined in interface $(problem.name), this will result empty assembly!") diff --git a/src/solvers.jl b/src/solvers.jl index 94edf53..4f14ad2 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -77,7 +77,7 @@ If several field problems exists, they are simply summed together, so problems must have unique node ids. """ -function get_field_assembly(solver::Solver; show_info=true) +function get_field_assembly(solver::Solver) problems = get_field_problems(solver) M = SparseMatrixCOO() @@ -96,7 +96,7 @@ function get_field_assembly(solver::Solver; show_info=true) if solver.ndofs == 0 solver.ndofs = size(K, 1) - show_info && info("automatically determined problem dimension, ndofs = $(solver.ndofs)") + info("automatically determined problem dimension, ndofs = $(solver.ndofs)") end M = sparse(M, solver.ndofs, solver.ndofs) @@ -364,8 +364,8 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric end """ Default assembler for solver. """ -function assemble!(solver::Solver; show_info=true, timing=true, with_mass_matrix=false) - show_info && info("Assembling problems ...") +function assemble!(solver::Solver; timing=true, with_mass_matrix=false) + info("Assembling problems ...") function do_assemble(problem) t00 = Base.time() @@ -391,7 +391,7 @@ function assemble!(solver::Solver; show_info=true, timing=true, with_mass_matrix solver.ndofs = ndofs t1 = round(Base.time()-t0, 2) - show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.") + info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.") if timing info("Assembly times:") for (i, problem) in enumerate(solver.problems) @@ -423,12 +423,12 @@ function get_unknown_field_dimension(solver::Solver) end """ Default initializer for solver. """ -function initialize!(solver::Solver; show_info=true) +function initialize!(solver::Solver) if solver.initialized - show_info && info("initialize!(): solver already initialized") + warn("initialize!(): solver already initialized") return end - show_info && info("Initializing solver ...") + info("Initializing solver ...") problems = get_problems(solver) length(problems) != 0 || error("Empty solver, add problems to solver using push!") t0 = Base.time() @@ -459,7 +459,7 @@ function initialize!(solver::Solver; show_info=true) # initialize(problem, ....) end t1 = round(Base.time()-t0, 2) - show_info && info("Initialized solver in $t1 seconds.") + info("Initialized solver in $t1 seconds.") solver.initialized = true end @@ -468,7 +468,7 @@ function get_all_elements(solver::Solver) return [elements...;] end -function (solver::Solver)(field_name::AbstractString, time::Float64) +function (solver::Solver)(field_name::String, time::Float64) fields = [] for problem in get_problems(solver) field = problem(field_name, time) @@ -482,14 +482,14 @@ function (solver::Solver)(field_name::AbstractString, time::Float64) end """ Default update for solver. """ -function update!{S}(solver::Solver{S}; show_info=true) +function update!{S}(solver::Solver{S}) u = solver.u la = solver.la - show_info && info("Updating problems ...") + info("Updating problems ...") t0 = Base.time() - for problem in solver.problems + for problem in get_problems(solver) assembly = get_assembly(problem) elements = get_elements(problem) # update solution, first for assembly (u,la) ... @@ -498,123 +498,45 @@ function update!{S}(solver::Solver{S}; show_info=true) update!(problem, assembly, elements, solver.time) end - # if io is attached to solver, update hdf / xml also - if !isnull(solver.xdmf) - update_xdmf!(solver) - end - t1 = round(Base.time()-t0, 2) - show_info && info("Updated problems in $t1 seconds.") + info("Updated problems in $t1 seconds.") end -function update_xdmf!{S}(solver::Solver{S}; show_info=true) - xdmf = get(solver.xdmf) - temporal_collection = get_temporal_collection(xdmf) - - # 1. save geometry - X_ = solver("geometry", solver.time) - node_ids = sort(collect(keys(X_))) - X = hcat([X_[nid] for nid in node_ids]...) - ndim, nnodes = size(X) - geom_type = (ndim == 2 ? "XY" : "XYZ") - data_node_ids = new_dataitem(xdmf, "/Node IDs", node_ids) - data_geometry = new_dataitem(xdmf, "/Geometry", X) - geometry = new_element("Geometry") - set_attribute(geometry, "Type", geom_type) - add_child(geometry, data_geometry) - - # 2. save topology - nid_mapping = Dict(j=>i for (i, j) in enumerate(node_ids)) - all_elements = get_all_elements(solver) - nelements = length(all_elements) - debug("Saving topology: $nelements elements total.") - element_types = unique(map(get_element_type, all_elements)) - - xdmf_element_mapping = Dict( - "Poi1" => "Polyvertex", - "Seg2" => "Polyline", - "Tri3" => "Triangle", - "Quad4" => "Quadrilateral", - "Tet4" => "Tetrahedron", - "Pyramid5" => "Pyramid", - "Wedge6" => "Wedge", - "Hex8" => "Hexahedron", - "Seg3" => "Edge_3", - "Tri6" => "Tri_6", - "Quad8" => "Quad_8", - "Tet10" => "Tet_10", - "Pyramid13" => "Pyramid_13", - "Wedge15" => "Wedge_15", - "Hex20" => "Hex_20") - - topology = [] - for element_type in element_types - elements = filter_by_element_type(element_type, all_elements) - nelements = length(elements) - info("Xdmf save: $nelements elements of type $element_type") - sort!(elements, by=get_element_id) - element_ids = map(get_element_id, elements) - element_conn = map(element -> [nid_mapping[j]-1 for j in get_connectivity(element)], elements) - element_conn = hcat(element_conn...) - element_code = split(string(element_type), ".")[end] - dataitem = new_dataitem(xdmf, "/Topology/$element_code/Element IDs", element_ids) - dataitem = new_dataitem(xdmf, "/Topology/$element_code/Connectivity", element_conn) - topology_ = new_element("Topology") - set_attribute(topology_, "TopologyType", xdmf_element_mapping[element_code]) - set_attribute(topology_, "NumberOfElements", length(elements)) - add_child(topology_, dataitem) - push!(topology, topology_) - end - - # 3. save solved field - frame = new_element("Grid") - time = new_child(frame, "Time") - set_attribute(time, "Value", solver.time) - add_child(frame, geometry) - for topo in topology - add_child(frame, topo) - end - - unknown_field_name = get_unknown_field_name(solver) - U_ = solver(unknown_field_name, solver.time) - node_ids2 = sort(collect(keys(U_))) - @assert node_ids == node_ids2 - - ndim = length(U_[first(node_ids)]) - field_type = ndim == 1 ? "Scalar" : "Vector" - field_center = "Node" - if ndim == 2 - for nid in node_ids - U_[nid] = [U_[nid]; 0.0] +""" Default postprocess for solver. Loop all problems and run postprocess +functions to calculate secondary fields, i.e. contact pressure, stress, +heat flux, reaction force etc. quantities. +""" +function postprocess!(solver::Solver) + info("Running postprocess scripts for solver...") + for problem in get_problems(solver) + for field_name in problem.postprocess_fields + field = Val{Symbol(field_name)} + info("Running postprocess for problem $(problem.name), field $field_name") + postprocess!(problem, solver.time, field) end - ndim = 3 end - U = zeros(X) - for nid in node_ids - loc = nid_mapping[nid] - U[:,loc] = U_[nid] - end - unknown_field_name = ucfirst(unknown_field_name) - time = solver.time - path = "" - if S == Nonlinear - iteration = solver.properties.iteration - path = "/Results/Time $time/Iteration $iteration/Nodal Fields/$unknown_field_name" - elseif S == Linear - path = "/Results/Time $time/Nodal Fields/$unknown_field_name" - end - attribute = new_child(frame, "Attribute") - set_attribute(attribute, "Name", unknown_field_name) - set_attribute(attribute, "Center", field_center) - set_attribute(attribute, "AttributeType", field_type) - add_child(attribute, new_dataitem(xdmf, path, U)) - add_child(frame, attribute) - if (S == Linear) || ((S == Nonlinear) && has_converged(solver)) - add_child(temporal_collection, frame) - end - save!(xdmf) end +""" Default xdmf update for solver. Loop all problems and write them individually +to Xdmf file. By default write the main unknown field (displacement, temperature, +...) and any fields requested separately in `problem.postprocess_fields` vector +(stress, strain, ...) +""" +function update_xdmf!(solver::Solver) + if isnull(solver.xdmf) + info("update_xdmf: xdmf not attached to solver, not writing output to file.") + info("turn Xdmf writing on to solver by typing: solver.xdmf = Xdmf(\"results\")") + return + end + xdmf = get(solver.xdmf) + for problem in get_problems(solver) + fields = [get_unknown_field_name(problem); problem.postprocess_fields] + if is_boundary_problem(problem) + fields = [fields; get_parent_field_name(problem)] + end + update_xdmf!(xdmf, problem, solver.time, fields) + end +end ### Nonlinear quasistatic solver @@ -671,18 +593,23 @@ function (solver::Solver{Nonlinear})() info("Increment time t=$(round(solver.time, 3))") info(repeat("-", 80)) - # 2.1 update linearized assemblies + # 2.1 update assemblies assemble!(solver) + # 2.2 call solver for linearized system solve!(solver) + # 2.3 update solution back to elements update!(solver) # 2.4 check convergence - if has_converged(solver) + if properties.iteration >= properties.min_iterations && has_converged(solver) info("Converged in $(properties.iteration) iterations.") - properties.iteration >= properties.min_iterations && return true - info("Convergence criteria met, but iteration < min_iterations, continuing...") + # 2.4.1 run any postprocessing of problems + postprocess!(solver) + # 2.4.2 update Xdmf output + update_xdmf!(solver) + return true end end @@ -721,8 +648,8 @@ Main differences in this solver, compared to nonlinear solver are: type Linear <: AbstractSolver end -function assemble!(solver::Solver{Linear}; show_info=true) - show_info && info("Assembling problems ...") +function assemble!(solver::Solver{Linear}) + info("Assembling problems ...") tic() nproblems = 0 ndofs = 0 @@ -731,27 +658,27 @@ function assemble!(solver::Solver{Linear}; show_info=true) assemble!(problem, solver.time) nproblems += 1 else - show_info && info("$(problem.name) already assembled, skipping.") + info("$(problem.name) already assembled, skipping.") end ndofs = max(ndofs, size(problem.assembly.K, 2)) end solver.ndofs = ndofs t1 = round(toq(), 2) - show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.") + info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.") end -function (solver::Solver{Linear})(; show_info=true) +function (solver::Solver{Linear})() t0 = Base.time() - show_info && info(repeat("-", 80)) - show_info && info("Starting linear solver") - show_info && info("Increment time t=$(round(solver.time, 3))") - show_info && info(repeat("-", 80)) + info(repeat("-", 80)) + info("Starting linear solver") + info("Increment time t=$(round(solver.time, 3))") + info(repeat("-", 80)) initialize!(solver) assemble!(solver) solve!(solver) update!(solver) t1 = round(Base.time()-t0, 2) - show_info && info("Linear solver ready in $t1 seconds.") + info("Linear solver ready in $t1 seconds.") end """ Convenience function to call linear solver. """ @@ -769,60 +696,3 @@ function LinearSolver(name::AbstractString, problems::Problem...) end ### End of linear quasistatic solver - -### Postprocessor - -type Postprocessor <: AbstractSolver - assembly :: Assembly - F :: Union{Factorization, Void} -end - -function Postprocessor() - Postprocessor(Assembly(), nothing) -end - -function assemble!(solver::Solver{Postprocessor}; show_info=true) - show_info && info("Assembling problems ...") - tic() - nproblems = 0 - ndofs = 0 - assembly = solver.properties.assembly - empty!(assembly) - for problem in get_problems(solver) - for element in get_elements(problem) - postprocess!(assembly, problem, element, solver.time) - end - nproblems += 1 - ndofs = max(ndofs, size(problem.assembly.K, 2)) - end - solver.ndofs = ndofs - t1 = round(toq(), 2) - show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.") -end - -function (solver::Solver{Postprocessor})(; show_info=true) - t0 = Base.time() - show_info && info(repeat("-", 80)) - show_info && info("Starting postprocessor") - show_info && info("Increment time t=$(round(solver.time, 3))") - show_info && info(repeat("-", 80)) - initialize!(solver) - assemble!(solver) - assembly = solver.properties.assembly - M = sparse(assembly.M) - f = sparse(assembly.f) - F = cholfact(M) - q = F \ f - t1 = round(Base.time()-t0, 2) - show_info && info("Postprocess of results ready in $t1 seconds.") - return q -end - -""" Convenience function to call postprocessor. """ -function Postprocessor(problems::Problem...) - solver = Solver(Postprocessor, "default postprocessor") - if length(problems) != 0 - push!(solver, problems...) - end - return solver -end diff --git a/src/solvers_modal.jl b/src/solvers_modal.jl index db7923b..18c4b88 100644 --- a/src/solvers_modal.jl +++ b/src/solvers_modal.jl @@ -209,8 +209,7 @@ Parameters sigma Shift stiffness matrix by adding diagonal term, i.e. K_shifted = K + sigma*I """ -function (solver::Solver{Modal})(; show_info=true, debug=false, - bc_invertible=false, P=nothing, symmetric=true, +function (solver::Solver{Modal})(; bc_invertible=false, P=nothing, symmetric=true, empty_assemblies_before_solution=true, dense=false, info_matrices=false, sigma=0.0) info(repeat("-", 80)) @@ -273,13 +272,6 @@ function (solver::Solver{Modal})(; show_info=true, debug=false, info("Calculate $(props.nev) eigenvalues...") - if debug && length(nz) < 100 - info("Stiffness matrix:") - dump(round(full(K_red))) - info("Mass matrix:") - dump(round(full(M_red))) - end - tic() if symmetric @@ -358,14 +350,18 @@ function (solver::Solver{Modal})(; show_info=true, debug=false, end end - if !isnull(solver.xdmf) - update_xdmf!(solver) - end + update_xdmf!(solver) + return true end -function update_xdmf!(solver::Solver{Modal}; show_info=true) +function update_xdmf!(solver::Solver{Modal}) + + if isnull(solver.xdmf) + info("update_xdmf: xdmf not attached to solver, not writing file output.") + return + end if maximum(abs(imag(solver.properties.eigvals))) > 1.0e-9 error("Writing imaginary eigenvalues for Xdmf not supported.") diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl index d6995f7..4f47455 100644 --- a/test/test_elasticity_2d_linear_with_surface_load.jl +++ b/test/test_elasticity_2d_linear_with_surface_load.jl @@ -10,12 +10,6 @@ using JuliaFEM.Testing - solve 2d plane stress problem with known solution: surface traction force in 2d volume load in 2d - reaction force -- test postprocessing of nodal fields: - geometry - displacement - reaction force - concentrated force =# @testset "test 2d linear elasticity with surface + volume load" begin meshfile = "/geometry/2d_block/BLOCK_1elem.med" @@ -60,75 +54,4 @@ using JuliaFEM.Testing u3 = block("displacement", 0.0)[3] u3_expected = f/E*[-nu, 1] + g/(2*E)*[-nu, 1] @test isapprox(u3, u3_expected) -#= - # fetch nodal results X + u and join them into one table using DataFrames - X = solver(DataFrame, "geometry", :COOR) - u = solver(DataFrame, "displacement", :U) - la = solver(DataFrame, "reaction force", :RF) - f = solver(DataFrame, "concentrated force", :CF) - results = join(X, u, on=:NODE, kind=:outer) - results = join(results, la, on=:NODE, kind=:outer) - length(f) != 0 && (results = join(results, f, on=:NODE, kind=:outer)) - sort!(results, cols=[:NODE]) - println(results) - - u3 = extract(results, NODE=:N3, :U1, :U2) - @test isapprox(u3, u3_expected) - - # element details - el = first(block.elements) - X = el("geometry", 0.0) - debug("X = $X") - S1 = block(el, [0.0, 0.0], 0.0, Val{:S}) - S1 = S1[[1,4,2]] - E1 = block(el, [0.0, 0.0], 0.0, Val{:E}) - E1 = E1[[1,4,2]] - C1 = block(el, [0.0, 0.0], 0.0, Val{:COORD}) - info("strain = $E1, stress = $S1, at $C1") - @test isapprox(E1, [-2/3, 2.0, 0.0]) - @test isapprox(S1, [0.0, 576.0, 0.0]) - @test isapprox(C1, [0.5, 0.5]) - - S1 = block(DataFrame, 0.0, Val{:S}) - E1 = block(DataFrame, 0.0, Val{:E}) - C1 = block(DataFrame, 0.0, Val{:COORD}) - println(S1) - println(E1) - println(C1) - S = solver(DataFrame, 0.0, Val{:S}) - println(S) - - info(solver("displacement", 0.0)) - solver() - info(solver("displacement", 0.0)) - u = solver("displacement", 0.0)[3] - info("u3 = $u") - @test isapprox(u, u3_expected) - -# info("calling nonlinear solver") -# solver2 = NonlinearSolver(block, traction, bc_sym_23, bc_sym_13) -# solver2() -# u = solver2("displacement", 0.0)[3] -# info("nlsolver u3 = $u, expected = $u3_expected") -# @test isapprox(u, u3_expected; rtol=1.0e-5) -=# end - -#= TODO: to other file -@testset "test dump model to disk and read back before and after solution" begin - solver = get_model("test 2d linear elasticity with surface + volume load") - save("/tmp/model.jld", "linear_model", solver) - solver2 = load("/tmp/model.jld")["linear_model"] - solver2() - save("/tmp/model.jld", "results", solver2) - solver3 = load("/tmp/model.jld")["results"] - block = solver3["BLOCK"] - u3 = reshape(block.assembly.u, 2, 4)[:,3] - f = 288.0 - g = 576.0 - E = 288.0 - nu = 1/3 - u3_expected = f/E*[-nu, 1] + g/(2*E)*[-nu, 1] - @test isapprox(u3, u3_expected) -end -=# diff --git a/test/test_heat.jl b/test/test_heat.jl index c9bb627..874d291 100644 --- a/test/test_heat.jl +++ b/test/test_heat.jl @@ -91,24 +91,25 @@ end @test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [2.0, 2.0]) end -function T_acc(x) - # accurate solution - a = 0.01 - L = 0.20 - k = 50.0 - Tᵤ = 20.0 - h = 10.0 - P = 4*a - A = a^2 - α = h - β = sqrt((h*P)/(k*A)) - T̂ = 100.0 - C = [1.0 1.0; (α+k*β)*exp(β*L) (α-k*β)*exp(-β*L)] \ [T̂-Tᵤ, 0.0] - return dot(C, [exp(β*x), exp(-β*x)]) + Tᵤ -end - #= @testset "test 1d heat problem" begin + + function T_acc(x) + # accurate solution + a = 0.01 + L = 0.20 + k = 50.0 + Tᵤ = 20.0 + h = 10.0 + P = 4*a + A = a^2 + α = h + β = sqrt((h*P)/(k*A)) + T̂ = 100.0 + C = [1.0 1.0; (α+k*β)*exp(β*L) (α-k*β)*exp(-β*L)] \ [T̂-Tᵤ, 0.0] + return dot(C, [exp(β*x), exp(-β*x)]) + Tᵤ + end + X = Dict{Int, Vector{Float64}}( 1 => [0.0, 0.0, 0.0], 2 => [0.1, 0.0, 0.0], @@ -203,13 +204,6 @@ end 7 => [1.74360055518019E+04, -4.73227515822118E+02, -1.75280965396335E+02], 8 => [1.74447696000717E+04, -4.72904678032179E+02, -1.75280965396335E+02]) - postprocessor = Postprocessor(p1) - flux = full(postprocessor()) - fluxd = Dict{Int64, Vector{Float64}}() - for j=1:8 - fluxd[j] = vec(flux[j,:]) - end - T = p1("temperature") for j in sort(collect(keys(T))) @@ -220,23 +214,6 @@ end @test rtol < 1.0e-9 end - for j=1:8 - q1 = get_integration_points(first(rod))[j]("heat flux", 0.0) - q2 = FLUX_ELGA[j] - rtol = norm(q1-q2)/max(norm(q1),norm(q2))*100.0 - @printf "ip %i flux, JF: (% e,% e,% e), CA: (% e,% e,% e), rtol: %10.6f %%\n" j q1... q2... rtol - # @test rtol < 0.05 - # testing in integration points makes no sense because they are in different order in CA - end - - for j in sort(collect(keys(fluxd))) - q1 = fluxd[j] - q2 = FLUX_NOEU[j] - rtol = norm(q1-q2)/max(norm(q1),norm(q2))*100.0 - @printf "node %i flux, JF: (% e,% e,% e), CA: (% e,% e,% e), rtol: %10.6f %%\n" j q1... q2... rtol - @test rtol < 1.0e-9 - end - end @testset "compare simple 3d heat problem to analytical solution" begin diff --git a/test/test_io.jl b/test/test_io.jl index 13537fb..9f02d79 100644 --- a/test/test_io.jl +++ b/test/test_io.jl @@ -7,7 +7,7 @@ using LightXML @testset "create new Xdmf object" begin r = Xdmf() - expected = "" + expected = "" @test string(r.xml) == expected end @@ -23,6 +23,14 @@ end @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") @@ -79,3 +87,102 @@ end @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 8f297de..7bac02f 100644 --- a/test/test_modal_analysis.jl +++ b/test/test_modal_analysis.jl @@ -39,7 +39,7 @@ end @testset "test eigenvalues for single tet4 element" begin solver = get_model() - solver(; debug=true) + solver() @test isapprox(solver.properties.eigvals, [4/3, 1/3]) end @@ -49,7 +49,7 @@ end # problem.properties.finite_strain = true problem.properties.geometric_stiffness = true solver.properties.geometric_stiffness = true - solver(; debug=true) + solver() @test isapprox(solver.properties.eigvals, [5/3, 2/3]) end @@ -119,4 +119,3 @@ end solver() @test isapprox(solver.properties.eigvals[1], 1.0) end - diff --git a/test/test_mortar_2d_contact.jl b/test/test_mortar_2d_contact.jl index 83f1459..964ceef 100644 --- a/test/test_mortar_2d_contact.jl +++ b/test/test_mortar_2d_contact.jl @@ -73,8 +73,8 @@ end master = first(contact.elements) slave = last(contact.elements) u = master("displacement", [0.0], 0.0) - la = slave("reaction force", [0.0], 0.0) + la = slave("lambda", [0.0], 0.0) info("u = $u, la = $la") @test isapprox(u, [-0.2, -0.15]) - @test isapprox(la, [0.0, -30.375]) + @test isapprox(la, [0.0, 30.375]) end diff --git a/test/test_mortar_2d_mesh_tie.jl b/test/test_mortar_2d_mesh_tie.jl index 66ee551..54884bd 100644 --- a/test/test_mortar_2d_mesh_tie.jl +++ b/test/test_mortar_2d_mesh_tie.jl @@ -174,8 +174,8 @@ end solver() slave_elements = get_slave_elements(interface) - node_ids, la = get_nodal_vector(slave_elements, "reaction force", 0.0) + node_ids, la = get_nodal_vector(slave_elements, "lambda", 0.0) for lai in la - @test isapprox(lai, [0.0, -10.0]) + @test isapprox(lai, [0.0, 10.0]) end end diff --git a/test/test_mortar_3d_mesh_tie_modal.jl b/test/test_mortar_3d_mesh_tie_modal.jl index cfe2804..8e331eb 100644 --- a/test/test_mortar_3d_mesh_tie_modal.jl +++ b/test/test_mortar_3d_mesh_tie_modal.jl @@ -151,25 +151,10 @@ FIN() freqs_jf = sqrt(solver.properties.eigvals)/(2*pi) freqs_ca = [1.12946E+00, 1.13141E+00, 2.93779E+00, 2.94143E+00, 4.51684E+00] - freq_jf = freqs_jf[1] - freq_ca = freqs_ca[1] - rtol = norm(freq_jf - freq_ca)/max(freq_jf, freq_ca) - info("rtol = $rtol") for (i, freq) in enumerate(freqs_jf) @printf "mode %i | freq JuliaFEM %8.3f | freq Code Aster %8.3f\n" i freqs_jf[i] freqs_ca[i] end - if rtol > 1.0e-3 - outfile = tempname() * ".xmf" - info("Something went wrong, results are saved to $outfile") - result = XDMF() - elems = [body1.elements; body2.elements] - for (i, freq) in enumerate(freqs_jf) - xdmf_new_result!(result, elems, freq) - xdmf_save_field!(result, elems, freq, "displacement"; field_type="Vector") - end - xdmf_save!(result, outfile) - end - @test rtol < 0.05 + @test isapprox(freqs_ca, freqs_jf; rtol=0.04) end diff --git a/test/test_mortar_autodiff.jl b/test/test_mortar_autodiff.jl deleted file mode 100644 index d42d4c9..0000000 --- a/test/test_mortar_autodiff.jl +++ /dev/null @@ -1,50 +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_testproblems(u, la) - nodes = Dict{Int64, Node}( - 1 => [0.0, 0.0], - 2 => [2.0, 0.0], - 3 => [1.0, 2.0], - 4 => [0.0, 2.0], - 5 => [2.5, 0.0], - 6 => [4.5, 0.0], - 7 => [4.5, 1.0], - 8 => [2.5, 1.0]) - displacement = Dict{Int64, Vector{Float64}}() - reaction_force = Dict{Int64, Vector{Float64}}() - for i=1:8 - displacement[i] = u[:,i] - reaction_force[i] = la[:,i] - end - bc5 = Seg2([3, 2]) - bc6 = Seg2([8, 5]) - update!([bc5, bc6], "geometry", nodes) - update!([bc5, bc6], "displacement", displacement) - update!([bc5, bc6], "reaction force", reaction_force) - bc5["master elements"] = [bc6] - contact1 = Problem(Mortar, "contact between bodies", 2, "displacement") - contact2 = Problem(Mortar, "contact between bodies", 2, "displacement") - contact2.properties.formulation = :forwarddiff - contact2.assembly.u = vec(u) - contact2.assembly.la = vec(la) - push!(contact1, bc5, bc6) - push!(contact2, bc5, bc6) - return contact1, contact2 -end - -#= TODO: Fix test -@testset "test linearization of contact force in undeformed state" begin - u = zeros(2, 8) - la = zeros(2, 8) - contact1, contact2 = get_testproblems(u, la) - assemble!(contact1, 0.0) - assemble!(contact2, 0.0) - @test isapprox(full(contact1.assembly.C1), full(contact2.assembly.C1)) - @test isapprox(full(contact1.assembly.K), full(contact2.assembly.K)) -end -=# - diff --git a/test/test_nodal_constraints.jl b/test/test_nodal_constraints.jl index f0aa399..7e546b0 100644 --- a/test/test_nodal_constraints.jl +++ b/test/test_nodal_constraints.jl @@ -63,8 +63,8 @@ end println(sparse(solver["symmetry 23"].assembly.C2)) solver() pel = nodal_bc.elements[1] - la = pel("reaction force", [0.0], 0.0) - info("reaction force: $la") + 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 diff --git a/test/test_postprocess_xdmf.jl b/test/test_postprocess_xdmf.jl deleted file mode 100644 index 05d6c5e..0000000 --- a/test/test_postprocess_xdmf.jl +++ /dev/null @@ -1,102 +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 - -testdata = """\ - - - - - - - - - -""" - -@testset "write simple xmf file" 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]) - u = Dict{Int64, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [0.0, 0.0], - 3 => [0.5, 1.0], - 4 => [0.0, 0.0]) - n = Dict{Int64, Vector{Float64}}( - 2 => [1.0, 0.0], - 3 => [1.0, 0.0]) - el1 = Element(Quad4, [1, 2, 3, 4]) - el2 = Element(Seg2, [2, 3]) - update!([el1, el2], "geometry", X) - update!([el1, el2], "displacement", u) - update!(el2, "normal", n) - xdmf = XDMF() - xdmf.dimension = 2 - xdmf_new_result!(xdmf, [el1, el2], 0.0) - xdmf_save_field!(xdmf, [el1, el2], 0.0, "displacement"; field_type="Vector") - xdmf_save_field!(xdmf, [el1, el2], 0.0, "normal"; field_type="Vector") - xdmf_save!(xdmf, "/tmp/test.xmf") - # TODO: how to test? -end - diff --git a/test/test_problem.jl b/test/test_problem.jl index c3d8dd9..76c45cc 100644 --- a/test/test_problem.jl +++ b/test/test_problem.jl @@ -34,7 +34,7 @@ end pr = Problem(Dirichlet, "bc", 1, "temperature") push!(pr, el) initialize!(pr) - @test haskey(el, "reaction force") + @test haskey(el, "lambda") @test haskey(el, "temperature") end diff --git a/test/test_problems_contact_2d.jl b/test/test_problems_contact_2d.jl index ccd973d..9752a74 100644 --- a/test/test_problems_contact_2d.jl +++ b/test/test_problems_contact_2d.jl @@ -60,9 +60,9 @@ datadir = first(splitext(basename(@__FILE__))) push!(solver, upper, lower, bc_fixed, bc_sym_23, load, contact) solver() slaves = get_slave_elements(contact) - node_ids, la = get_nodal_vector(slaves, "reaction force", 0.0) + node_ids, la = get_nodal_vector(slaves, "lambda", 0.0) node_ids, n = get_nodal_vector(slaves, "normal", 0.0) - pres = [dot(ni, -lai) for (ni, lai) in zip(n, la)] + pres = [dot(ni, lai) for (ni, lai) in zip(n, la)] #@test isapprox(maximum(pres), 4060.010799583303) # 12 % error in maximum pressure # integrate pressure in normal and tangential direction @@ -75,15 +75,15 @@ datadir = first(splitext(basename(@__FILE__))) w = ip.weight*sel(ip, time, Val{:detJ}) n = sel("normal", ip, time) t = Q'*n - la = sel("reaction force", ip, time) - Rn += w*dot(n, -la) + la = sel("lambda", ip, time) + Rn += w*dot(n, la) Rt += w*dot(t, la) end end info("2d hertz: Rn = $Rn, Rt = $Rt") info("2d hertz: maximum pressure pmax = ", maximum(pres)) @test isapprox(maximum(pres), 3585.0; rtol = 0.13) - # under 0.15 % error in reaction force + # under 0.15 % error in resultant force @test isapprox(Rn, 35.0e3; rtol=0.020) @test isapprox(Rt, 0.0; atol=200.0) end @@ -136,9 +136,9 @@ end push!(solver, upper, lower, bc_fixed, bc_sym_23, load, contact) solver() slaves = get_slave_elements(contact) - node_ids, la = get_nodal_vector(slaves, "reaction force", 0.0) + node_ids, la = get_nodal_vector(slaves, "lambda", 0.0) node_ids, n = get_nodal_vector(slaves, "normal", 0.0) - pres = [dot(ni, -lai) for (ni, lai) in zip(n, la)] + pres = [dot(ni, lai) for (ni, lai) in zip(n, la)] #@test isapprox(maximum(pres), 4060.010799583303) # 12 % error in maximum pressure # integrate pressure in normal and tangential direction @@ -151,15 +151,15 @@ end w = ip.weight*sel(ip, time, Val{:detJ}) n = sel("normal", ip, time) t = Q'*n - la = sel("reaction force", ip, time) - Rn += w*dot(n, -la) + la = sel("lambda", ip, time) + Rn += w*dot(n, la) Rt += w*dot(t, la) end end info("2d hertz: Rn = $Rn, Rt = $Rt") info("2d hertz: maximum pressure pmax = ", maximum(pres)) @test isapprox(maximum(pres), 3585.0; rtol = 0.13) - # under 0.15 % error in reaction force + # under 0.15 % error in resultant force @test isapprox(Rn, 35.0e3; rtol=0.020) @test isapprox(Rt, 0.0; atol=200.0) end @@ -219,9 +219,9 @@ end node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0) node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0) - node_ids, reaction_force = get_nodal_vector(get_slave_elements(interface), "reaction force", 0.0) + node_ids, lambda = get_nodal_vector(get_slave_elements(interface), "lambda", 0.0) u2 = [u[2] for u in displacement] - f2 = [f[2] for f in reaction_force] + f2 = [f[2] for f in lambda] maxabsu2 = maximum(abs(u2)) stdabsu2 = std(abs(u2)) info("max(abs(u2)) = $maxabsu2, std(abs(u2)) = $stdabsu2") @@ -242,9 +242,9 @@ end node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0) node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0) - node_ids, reaction_force = get_nodal_vector(get_slave_elements(interface), "reaction force", 0.0) + node_ids, lambda = get_nodal_vector(get_slave_elements(interface), "lambda", 0.0) u2 = [u[2] for u in displacement] - f2 = [f[2] for f in reaction_force] + f2 = [f[2] for f in lambda] maxabsu2 = maximum(abs(u2)) stdabsu2 = std(abs(u2)) info("max(abs(u2)) = $maxabsu2, std(abs(u2)) = $stdabsu2") diff --git a/test/test_problems_contact_2d_autodiff.jl b/test/test_problems_contact_2d_autodiff.jl index 6f26896..9e88e6e 100644 --- a/test/test_problems_contact_2d_autodiff.jl +++ b/test/test_problems_contact_2d_autodiff.jl @@ -76,9 +76,9 @@ end node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 1.0) node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 1.0) - node_ids, reaction_force = get_nodal_vector(get_slave_elements(interface), "reaction force", 1.0) + node_ids, la = get_nodal_vector(get_slave_elements(interface), "lambda", 1.0) u2 = [u[2] for u in displacement] - f2 = [f[2] for f in reaction_force] + f2 = [f[2] for f in la] maxabsu2 = maximum(abs(u2)) stdabsu2 = std(abs(u2)) info("max(abs(u2)) = $maxabsu2, std(abs(u2)) = $stdabsu2") diff --git a/test/test_problems_contact_3d.jl b/test/test_problems_contact_3d.jl index 2c27f9d..3dd3c7b 100644 --- a/test/test_problems_contact_3d.jl +++ b/test/test_problems_contact_3d.jl @@ -76,12 +76,18 @@ end solver() node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0) node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0) - # node_ids, pressure = get_nodal_vector(interface.elements, "contact pressure", 0.0) + # test postprocess of fields + postprocess!(interface, 0.0, Val{Symbol("contact pressure")}) + node_ids, contact_pressure = get_nodal_vector(interface.elements, "contact pressure", 0.0) u3 = [u[3] for u in displacement] maxabsu3 = maximum(abs(u3)) stdabsu3 = std(abs(u3)) + maxpres = maximum(contact_pressure) + stdpres = std(contact_pressure) info("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3") + info("max(contact_pressure) = $maxpres, std(contact_pressure) = $stdpres") @test isapprox(stdabsu3, 0.0; atol=1.0e-12) + @test isapprox(maxpres, 172.8; atol=1.0e-6) end @testset "small sliding contact patch test, tet4 + dual basis" begin @@ -116,7 +122,7 @@ end @testset "small sliding contact patch test, tet10 + dual basis, alpha=0.2" begin solver = get_model(tet10_meshfile) - solver.xdmf = Xdmf("contact_sl_quad_disp_results"; overwrite=true) + solver.xdmf = Xdmf("contact_dl_quad_disp_results"; overwrite=true) interface = solver["LOWER_TO_UPPER"] interface.properties.dual_basis = true interface.properties.alpha = 0.2 diff --git a/test/test_problems_elasticity.jl b/test/test_problems_elasticity.jl new file mode 100644 index 0000000..781f056 --- /dev/null +++ b/test/test_problems_elasticity.jl @@ -0,0 +1,54 @@ +# 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.Testing + +@testset "test postprocessing of strain and stress" 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, 1.0], + 6 => [1.0, 0.0, 1.0], + 7 => [1.0, 1.0, 1.0], + 8 => [0.0, 1.0, 1.0]) + + u = Dict( + 1 => [0.0, 0.0, 0.0], + 2 => [-1/3, 0.0, 0.0], + 3 => [-1/3, -1/3, 0.0], + 4 => [0.0, -1/3, 0.0], + 5 => [0.0, 0.0, 1.0], + 6 => [-1/3, 0.0, 1.0], + 7 => [-1/3, -1/3, 1.0], + 8 => [0.0, -1/3, 1.0]) + + element = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8]) + update!(element, "geometry", 0.0 => X) + update!(element, "displacement", 0.0 => u) + update!(element, "youngs modulus", 288.0) + update!(element, "poissons ratio", 1/3) + + body = Problem(Elasticity, "[0,1]³ elastic unit block", 3) + body.elements= [element] + postprocess!(body, 0.0, Val{:strain}) + postprocess!(body, 0.0, Val{:stress}) + + geom = element("geometry", 0.0) + strain = element("strain", 0.0) + stress = element("stress", 0.0) + + for i in 1:length(geom) + a = geom[i] + b = strain[i] + c = stress[i] + debug("$a -> $b -> $c") + @test isapprox(b, [-1/3, -1/3, 1.0, 0.0, 0.0, 0.0]) + @test isapprox(c, [0.0, 0.0, 288.0, 0.0, 0.0, 0.0]) + end + +end diff --git a/test/test_solvers_postprocess.jl b/test/test_solvers_postprocess.jl new file mode 100644 index 0000000..8b277f6 --- /dev/null +++ b/test/test_solvers_postprocess.jl @@ -0,0 +1,94 @@ +# 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 +using JuliaFEM.Abaqus: create_surface_elements + +datadir = first(splitext(basename(@__FILE__))) + +@testset "test postprocessing of secondary fields" begin + mesh = abaqus_read_mesh(joinpath(datadir, "tet4.inp")) + + upper = Problem(Elasticity, "UPPER", 3) + upper.elements = create_elements(mesh, "UPPER") + update!(upper, "youngs modulus", 3*288.0) + update!(upper, "poissons ratio", 1/3) + push!(upper.postprocess_fields, "strain", "stress") + info("upper postprocess: $(upper.postprocess_fields)") + + lower = Problem(Elasticity, "LOWER", 3) + lower.elements = create_elements(mesh, "LOWER") + update!(lower, "youngs modulus", 288.0) + update!(lower, "poissons ratio", 1/3) + push!(lower.postprocess_fields, "strain", "stress") + info("lower postprocess: $(lower.postprocess_fields)") + + bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement") + bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP") + update!(bc_upper, "displacement 3", -0.4) + + bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement") + bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM") + update!(bc_lower, "displacement 3", 0.0) + + # point-wise boundary conditions to prevent free body move + nid1 = find_nearest_nodes(mesh, [0.0, 0.0, 0.0])[1] + nid2 = find_nearest_nodes(mesh, [1.0, 0.0, 0.0])[1] + nid3 = find_nearest_nodes(mesh, [0.0, 1.0, 0.0])[1] + nid4 = find_nearest_nodes(mesh, [0.0, 0.0, 1.0])[1] + nid5 = find_nearest_nodes(mesh, [1.0, 0.0, 1.0])[1] + nid6 = find_nearest_nodes(mesh, [0.0, 1.0, 1.0])[1] + + bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement") + # nodes in X2=0 plane + bc_sym13.elements = [Element(Poi1, [j]) for j in [nid1, nid2, nid4, nid5]] + update!(bc_sym13, "geometry", mesh.nodes) + update!(bc_sym13, "displacement 2", 0.0) + + bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement") + # nodes in X1=0 plane + bc_sym23.elements = [Element(Poi1, [j]) for j in [nid1, nid3, nid4, nid6]] + update!(bc_sym23, "geometry", mesh.nodes) + update!(bc_sym23, "displacement 1", 0.0) + + for bc in [bc_upper, bc_lower, bc_sym13, bc_sym23] + push!(bc.postprocess_fields, "reaction force") + end + + interface = Problem(Contact, "LOWER_TO_UPPER", 3, "displacement") + interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER") + interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER") + update!(interface_slave_elements, "master elements", interface_master_elements) + interface.elements = [interface_slave_elements; interface_master_elements] + interface.properties.dual_basis = true + interface.properties.contact_state_in_first_iteration = :AUTO + push!(interface.postprocess_fields, "contact pressure") + + solver = NonlinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface) + + solver.xdmf = Xdmf("contact_two_blocks_postprocess"; overwrite=true) + solver() + + node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0) + node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0) + # node_ids, pressure = get_nodal_vector(interface.elements, "contact pressure", 0.0) + u3 = [u[3] for u in displacement] + maxabsu3 = maximum(abs(u3)) + stdabsu3 = std(abs(u3)) + info("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3") + @test isapprox(stdabsu3, 0.0; atol=1.0e-12) + upper_X = [0.5, 0.5, 0.75] + lower_X = [0.5, 0.5, 0.25] + strain_upper = upper("strain", upper_X) + stress_upper = upper("stress", upper_X) + strain_lower = lower("strain", lower_X) + stress_lower = lower("stress", lower_X) + info("strain at $upper_X = $strain_upper") + info("stress at $upper_X = $stress_upper") + info("strain at $lower_X = $strain_lower") + info("stress at $lower_X = $stress_lower") +end + diff --git a/test/test_solvers_postprocess/tet4.inp b/test/test_solvers_postprocess/tet4.inp new file mode 100644 index 0000000..fabfbb8 --- /dev/null +++ b/test/test_solvers_postprocess/tet4.inp @@ -0,0 +1,589 @@ +**NSET COUNT = 116 +*NODE +127, 1.00000, 0.00000, 0.75000 +128, 1.00000, 0.00000, 1.00000 +129, 0.75000, 0.00000, 1.00000 +133, 0.75000, 0.00000, 0.75000 +136, 1.00000, 0.00000, 0.50000 +139, 0.75000, 0.00000, 0.50000 +142, 0.50000, 0.00000, 1.00000 +145, 0.50000, 0.00000, 0.75000 +149, 0.50000, 0.00000, 0.50000 +152, 0.25000, 0.00000, 1.00000 +155, 0.25000, 0.00000, 0.75000 +159, 0.25000, 0.00000, 0.50000 +162, 0.00000, 0.00000, 1.00000 +165, 0.00000, 0.00000, 0.75000 +169, 0.00000, 0.00000, 0.50000 +172, 0.75000, 1.00000, 1.00000 +173, 1.00000, 1.00000, 1.00000 +174, 1.00000, 1.00000, 0.75000 +178, 0.75000, 1.00000, 0.75000 +181, 1.00000, 1.00000, 0.50000 +184, 0.75000, 1.00000, 0.50000 +187, 0.50000, 1.00000, 1.00000 +190, 0.50000, 1.00000, 0.75000 +194, 0.50000, 1.00000, 0.50000 +197, 0.25000, 1.00000, 1.00000 +200, 0.25000, 1.00000, 0.75000 +204, 0.25000, 1.00000, 0.50000 +207, 0.00000, 1.00000, 1.00000 +210, 0.00000, 1.00000, 0.75000 +214, 0.00000, 1.00000, 0.50000 +217, 1.00000, 0.75000, 1.00000 +220, 1.00000, 0.75000, 0.75000 +224, 1.00000, 0.75000, 0.50000 +227, 1.00000, 0.50000, 1.00000 +230, 1.00000, 0.50000, 0.75000 +234, 1.00000, 0.50000, 0.50000 +237, 1.00000, 0.25000, 1.00000 +240, 1.00000, 0.25000, 0.75000 +244, 1.00000, 0.25000, 0.50000 +252, 0.00000, 0.75000, 1.00000 +255, 0.00000, 0.75000, 0.75000 +259, 0.00000, 0.75000, 0.50000 +262, 0.00000, 0.50000, 1.00000 +265, 0.00000, 0.50000, 0.75000 +269, 0.00000, 0.50000, 0.50000 +272, 0.00000, 0.25000, 1.00000 +275, 0.00000, 0.25000, 0.75000 +279, 0.00000, 0.25000, 0.50000 +288, 0.75000, 0.75000, 1.00000 +292, 0.50000, 0.75000, 1.00000 +296, 0.25000, 0.75000, 1.00000 +302, 0.75000, 0.50000, 1.00000 +306, 0.50000, 0.50000, 1.00000 +310, 0.25000, 0.50000, 1.00000 +316, 0.75000, 0.25000, 1.00000 +320, 0.50000, 0.25000, 1.00000 +324, 0.25000, 0.25000, 1.00000 +337, 0.75000, 0.75000, 0.50000 +341, 0.50000, 0.75000, 0.50000 +345, 0.25000, 0.75000, 0.50000 +351, 0.75000, 0.50000, 0.50000 +355, 0.50000, 0.50000, 0.50000 +359, 0.25000, 0.50000, 0.50000 +365, 0.75000, 0.25000, 0.50000 +369, 0.50000, 0.25000, 0.50000 +373, 0.25000, 0.25000, 0.50000 +438, 1.00000, 0.00000, 0.25000 +439, 1.00000, 0.00000, 0.50000 +440, 0.66667, 0.00000, 0.50000 +444, 0.66667, 0.00000, 0.25000 +447, 1.00000, 0.00000, 0.00000 +450, 0.66667, 0.00000, 0.00000 +453, 0.33333, 0.00000, 0.50000 +456, 0.33333, 0.00000, 0.25000 +460, 0.33333, 0.00000, 0.00000 +463, 0.00000, 0.00000, 0.50000 +466, 0.00000, 0.00000, 0.25000 +470, 0.00000, 0.00000, 0.00000 +473, 0.66667, 1.00000, 0.50000 +474, 1.00000, 1.00000, 0.50000 +475, 1.00000, 1.00000, 0.25000 +479, 0.66667, 1.00000, 0.25000 +482, 1.00000, 1.00000, 0.00000 +485, 0.66667, 1.00000, 0.00000 +488, 0.33333, 1.00000, 0.50000 +491, 0.33333, 1.00000, 0.25000 +495, 0.33333, 1.00000, 0.00000 +498, 0.00000, 1.00000, 0.50000 +501, 0.00000, 1.00000, 0.25000 +505, 0.00000, 1.00000, 0.00000 +508, 1.00000, 0.66667, 0.50000 +511, 1.00000, 0.66667, 0.25000 +515, 1.00000, 0.66667, 0.00000 +518, 1.00000, 0.33333, 0.50000 +521, 1.00000, 0.33333, 0.25000 +525, 1.00000, 0.33333, 0.00000 +533, 0.00000, 0.66667, 0.50000 +536, 0.00000, 0.66667, 0.25000 +540, 0.00000, 0.66667, 0.00000 +543, 0.00000, 0.33333, 0.50000 +546, 0.00000, 0.33333, 0.25000 +550, 0.00000, 0.33333, 0.00000 +559, 0.66667, 0.66667, 0.50000 +563, 0.33333, 0.66667, 0.50000 +569, 0.66667, 0.33333, 0.50000 +573, 0.33333, 0.33333, 0.50000 +584, 0.66667, 0.66667, 0.00000 +588, 0.33333, 0.66667, 0.00000 +594, 0.66667, 0.33333, 0.00000 +598, 0.33333, 0.33333, 0.00000 +608, 0.66006, 0.47128, 0.70078 +609, 0.19762, 0.64074, 0.81889 +610, 0.65771, 0.82813, 0.74829 +611, 0.59993, 0.15953, 0.76106 +612, 0.40915, 0.16311, 0.74937 +613, 0.50000, 0.50000, 0.25000 +** +**ELSET COUNT = 172 +**HWCOLOR COMP 54 0 +*ELEMENT, TYPE=C3D4, ELSET=UPPER + 570, 252, 262, 255, 609 + 571, 252, 310, 262, 609 + 572, 252, 296, 310, 609 + 573, 200, 252, 255, 609 + 574, 296, 306, 310, 609 + 575, 262, 310, 265, 609 + 576, 259, 265, 269, 359 + 577, 184, 337, 220, 610 + 578, 259, 345, 265, 359 + 579, 288, 608, 292, 610 + 580, 187, 288, 292, 610 + 581, 172, 288, 187, 610 + 582, 178, 172, 187, 610 + 583, 178, 220, 172, 610 + 584, 172, 220, 288, 610 + 585, 230, 240, 237, 316 + 586, 178, 184, 220, 610 + 587, 178, 190, 184, 610 + 588, 178, 187, 190, 610 + 589, 187, 292, 190, 610 + 590, 190, 292, 609, 610 + 591, 190, 609, 341, 610 + 592, 190, 341, 194, 610 + 593, 184, 190, 194, 610 + 594, 184, 194, 337, 610 + 595, 194, 341, 337, 610 + 596, 337, 341, 608, 610 + 597, 230, 337, 608, 610 + 598, 230, 608, 288, 610 + 599, 220, 230, 288, 610 + 600, 220, 337, 230, 610 + 601, 341, 609, 608, 610 + 602, 292, 608, 609, 610 + 603, 306, 316, 320, 611 + 604, 306, 608, 316, 611 + 605, 365, 608, 369, 611 + 606, 240, 608, 365, 611 + 607, 240, 316, 608, 611 + 608, 240, 133, 316, 611 + 609, 240, 365, 133, 611 + 610, 365, 139, 133, 611 + 611, 365, 369, 139, 611 + 612, 369, 149, 139, 611 + 613, 149, 145, 611, 612 + 614, 240, 127, 128, 316 + 615, 320, 129, 142, 611 + 616, 316, 129, 320, 611 + 617, 316, 133, 129, 611 + 618, 129, 133, 142, 611 + 619, 142, 133, 145, 611 + 620, 139, 145, 133, 611 + 621, 139, 149, 145, 611 + 622, 310, 320, 324, 612 + 623, 306, 320, 310, 612 + 624, 306, 611, 320, 612 + 625, 306, 608, 611, 612 + 626, 306, 609, 608, 612 + 627, 306, 310, 609, 612 + 628, 265, 609, 310, 612 + 629, 265, 310, 324, 612 + 630, 265, 324, 275, 612 + 631, 265, 275, 373, 612 + 632, 265, 373, 359, 612 + 633, 265, 359, 609, 612 + 634, 355, 609, 359, 612 + 635, 355, 608, 609, 612 + 636, 355, 369, 608, 612 + 637, 369, 611, 608, 612 + 638, 214, 255, 259, 345 + 639, 142, 611, 145, 612 + 640, 320, 611, 142, 612 + 641, 320, 142, 324, 612 + 642, 324, 142, 152, 612 + 643, 324, 152, 155, 612 + 644, 275, 324, 155, 612 + 645, 275, 155, 373, 612 + 646, 373, 155, 159, 612 + 647, 373, 159, 149, 612 + 648, 369, 373, 149, 612 + 649, 359, 373, 369, 612 + 650, 355, 359, 369, 612 + 651, 369, 149, 611, 612 + 652, 149, 155, 145, 612 + 653, 149, 159, 155, 612 + 654, 152, 145, 155, 612 + 655, 142, 145, 152, 612 + 656, 237, 240, 128, 316 + 657, 244, 136, 127, 365 + 658, 240, 365, 127, 133 + 659, 240, 351, 244, 365 + 660, 220, 224, 230, 337 + 661, 227, 288, 230, 302 + 662, 200, 207, 252, 296 + 663, 240, 127, 316, 133 + 664, 230, 337, 234, 351 + 665, 275, 162, 165, 155 + 666, 172, 220, 217, 288 + 667, 200, 252, 210, 255 + 668, 190, 292, 197, 296 + 669, 230, 234, 240, 351 + 670, 204, 210, 214, 345 + 671, 210, 255, 214, 345 + 672, 275, 165, 373, 155 + 673, 172, 174, 217, 220 + 674, 200, 204, 341, 345 + 675, 240, 244, 127, 365 + 676, 220, 230, 227, 288 + 677, 272, 275, 324, 155 + 678, 272, 152, 162, 155 + 679, 200, 255, 210, 345 + 680, 190, 200, 194, 341 + 681, 197, 207, 200, 296 + 682, 190, 197, 200, 296 + 683, 279, 165, 169, 159 + 684, 200, 207, 210, 252 + 685, 224, 234, 230, 337 + 686, 187, 197, 190, 292 + 687, 275, 165, 279, 373 + 688, 265, 275, 269, 373 + 689, 272, 324, 152, 155 + 690, 262, 265, 310, 324 + 691, 265, 272, 275, 324 + 692, 184, 224, 220, 337 + 693, 272, 162, 275, 155 + 694, 262, 272, 265, 324 + 695, 178, 174, 172, 220 + 696, 136, 365, 139, 133 + 697, 227, 230, 237, 302 + 698, 178, 184, 181, 220 + 699, 172, 174, 173, 217 + 700, 194, 200, 204, 341 + 701, 234, 244, 240, 351 + 702, 200, 210, 204, 345 + 703, 230, 237, 302, 316 + 704, 128, 316, 127, 129 + 705, 217, 220, 227, 288 + 706, 127, 365, 136, 133 + 707, 265, 269, 359, 373 + 708, 279, 373, 165, 159 + 709, 165, 159, 373, 155 + 710, 269, 275, 279, 373 + 711, 127, 129, 316, 133 + 712, 181, 220, 184, 224 + 713, 178, 181, 174, 220 + 714, 292, 302, 306, 608 + 715, 288, 302, 292, 608 + 716, 230, 302, 288, 608 + 717, 230, 316, 302, 608 + 718, 230, 240, 316, 608 + 719, 230, 351, 240, 608 + 720, 230, 337, 351, 608 + 721, 302, 316, 306, 608 + 722, 240, 351, 365, 608 + 723, 351, 355, 365, 608 + 724, 341, 355, 351, 608 + 725, 337, 341, 351, 608 + 726, 355, 369, 365, 608 + 727, 341, 355, 608, 609 + 728, 292, 608, 306, 609 + 729, 292, 306, 296, 609 + 730, 190, 292, 296, 609 + 731, 190, 200, 341, 609 + 732, 190, 296, 200, 609 + 733, 200, 296, 252, 609 + 734, 200, 255, 345, 609 + 735, 200, 345, 341, 609 + 736, 341, 345, 355, 609 + 737, 345, 359, 355, 609 + 738, 265, 359, 345, 609 + 739, 255, 265, 345, 609 + 740, 255, 265, 259, 345 + 741, 255, 262, 265, 609 +** +**ELSET COUNT = 92 +**HWCOLOR COMP 1 0 +*ELEMENT, TYPE=C3D4, ELSET=LOWER + 742, 456, 460, 466, 598 + 743, 444, 450, 456, 598 + 744, 444, 525, 447, 594 + 745, 536, 540, 588, 598 + 746, 453, 456, 463, 543 + 747, 533, 543, 536, 573 + 748, 444, 456, 453, 573 + 749, 466, 550, 546, 598 + 750, 536, 543, 546, 573 + 751, 491, 536, 495, 588 + 752, 438, 447, 444, 525 + 753, 444, 594, 450, 598 + 754, 444, 521, 525, 594 + 755, 438, 444, 440, 521 + 756, 540, 546, 550, 598 + 757, 453, 456, 543, 573 + 758, 488, 491, 559, 563 + 759, 505, 536, 540, 588 + 760, 440, 444, 453, 573 + 761, 491, 495, 584, 588 + 762, 444, 447, 450, 594 + 763, 485, 491, 495, 584 + 764, 440, 569, 444, 573 + 765, 495, 501, 505, 536 + 766, 473, 479, 511, 559 + 767, 438, 440, 439, 521 + 768, 491, 533, 536, 563 + 769, 444, 569, 521, 594 + 770, 473, 479, 475, 511 + 771, 491, 498, 501, 533 + 772, 495, 536, 505, 588 + 773, 438, 444, 521, 525 + 774, 473, 511, 508, 559 + 775, 440, 521, 444, 569 + 776, 511, 515, 521, 584 + 777, 491, 498, 533, 563 + 778, 515, 525, 521, 584 + 779, 460, 550, 466, 598 + 780, 439, 521, 440, 569 + 781, 456, 463, 543, 546 + 782, 479, 482, 511, 515 + 783, 533, 536, 563, 573 + 784, 479, 515, 511, 584 + 785, 450, 460, 456, 598 + 786, 521, 569, 559, 594 + 787, 559, 594, 569, 613 + 788, 559, 584, 594, 613 + 789, 521, 559, 584, 594 + 790, 536, 546, 540, 598 + 791, 511, 521, 559, 584 + 792, 479, 559, 491, 584 + 793, 491, 559, 563, 584 + 794, 511, 521, 518, 559 + 795, 491, 584, 563, 588 + 796, 563, 588, 584, 613 + 797, 559, 563, 584, 613 + 798, 559, 569, 563, 613 + 799, 563, 569, 573, 613 + 800, 536, 563, 573, 588 + 801, 563, 573, 588, 613 + 802, 573, 598, 588, 613 + 803, 491, 501, 495, 536 + 804, 536, 573, 546, 598 + 805, 536, 588, 573, 598 + 806, 588, 598, 594, 613 + 807, 584, 588, 594, 613 + 808, 444, 569, 594, 598 + 809, 521, 584, 525, 594 + 810, 444, 573, 569, 598 + 811, 569, 598, 573, 613 + 812, 569, 594, 598, 613 + 813, 444, 456, 573, 598 + 814, 439, 518, 521, 569 + 815, 456, 546, 543, 573 + 816, 479, 485, 482, 515 + 817, 508, 511, 518, 559 + 818, 456, 546, 573, 598 + 819, 488, 498, 491, 563 + 820, 479, 485, 515, 584 + 821, 491, 563, 536, 588 + 822, 473, 488, 479, 559 + 823, 479, 491, 485, 584 + 824, 491, 533, 501, 536 + 825, 479, 511, 559, 584 + 826, 479, 488, 491, 559 + 827, 475, 479, 482, 511 + 828, 518, 559, 521, 569 + 829, 473, 475, 508, 511 + 830, 473, 475, 474, 508 + 831, 460, 470, 466, 550 + 832, 456, 466, 463, 546 + 833, 456, 466, 546, 598 +** +**ELSET COUNT = 18 +*SURFACE,TYPE=ELEMENT,NAME=LOWER_TO_UPPER +799,S1 +798,S1 +764,S2 +814,S2 +747,S2 +828,S2 +819,S2 +822,S2 +777,S3 +746,S4 +757,S4 +760,S4 +780,S4 +783,S4 +817,S4 +758,S4 +774,S4 +830,S4 +** +**ELSET COUNT = 18 +*SURFACE,TYPE=ELEMENT,NAME=LOWER_BOTTOM +806,S1 +807,S1 +831,S2 +779,S2 +785,S2 +778,S2 +753,S3 +762,S3 +744,S3 +745,S3 +809,S3 +761,S3 +820,S3 +816,S3 +756,S4 +759,S4 +772,S4 +763,S4 +** +**ELSET COUNT = 12 +*SURFACE,TYPE=ELEMENT,NAME=LOWER_SYM13 +831,S1 +742,S1 +832,S1 +746,S1 +785,S1 +743,S1 +748,S1 +760,S1 +762,S1 +752,S1 +755,S1 +767,S1 +** +**ELSET COUNT = 12 +*SURFACE,TYPE=ELEMENT,NAME=LOWER_SYM23 +749,S1 +756,S1 +790,S1 +750,S1 +747,S1 +759,S1 +831,S3 +832,S3 +781,S3 +765,S3 +824,S3 +771,S3 +** +**ELSET COUNT = 32 +*SURFACE,TYPE=ELEMENT,NAME=UPPER_TOP +678,S1 +689,S1 +642,S1 +641,S1 +615,S1 +616,S1 +622,S1 +623,S1 +603,S1 +721,S1 +571,S1 +572,S1 +574,S1 +729,S1 +714,S1 +715,S1 +580,S1 +581,S1 +704,S2 +694,S2 +661,S2 +681,S2 +686,S2 +703,S3 +662,S3 +668,S3 +656,S4 +690,S4 +697,S4 +705,S4 +666,S4 +699,S4 +** +**ELSET COUNT = 32 +*SURFACE,TYPE=ELEMENT,NAME=UPPER_TO_LOWER +647,S1 +648,S1 +612,S1 +611,S1 +696,S1 +649,S1 +650,S1 +726,S1 +723,S1 +737,S1 +736,S1 +724,S1 +725,S1 +595,S1 +594,S1 +708,S2 +657,S2 +701,S2 +578,S2 +685,S2 +692,S2 +707,S3 +659,S3 +664,S3 +674,S3 +683,S4 +710,S4 +576,S4 +638,S4 +670,S4 +700,S4 +712,S4 +** +**ELSET COUNT = 16 +*SURFACE,TYPE=ELEMENT,NAME=UPPER_SYM13 +653,S1 +652,S1 +654,S1 +655,S1 +621,S1 +620,S1 +619,S1 +618,S1 +709,S2 +711,S2 +683,S3 +665,S3 +678,S3 +696,S4 +706,S4 +704,S4 +** +**ELSET COUNT = 16 +*SURFACE,TYPE=ELEMENT,NAME=UPPER_SYM23 +683,S1 +687,S1 +665,S1 +693,S1 +710,S1 +688,S1 +691,S1 +694,S1 +576,S1 +740,S1 +741,S1 +570,S1 +638,S1 +671,S1 +667,S3 +684,S3 +** +**Property Definitions +** +*SOLID SECTION, ELSET=UPPER, MATERIAL=Def_Material +*SOLID SECTION, ELSET=LOWER, MATERIAL=Def_Material +** +**Material Definitions +** +**Material:Def_Material +*MATERIAL,NAME=Def_Material +*ELASTIC,TYPE=ISO +2.08000e+005,3.00000e-001 +*DENSITY +7.80000e-009, +*SPECIFIC HEAT +5.00000e-001 +*CONDUCTIVITY +4.98100e-002 +** \ No newline at end of file diff --git a/test/test_types.jl b/test/test_types.jl deleted file mode 100644 index d48f606..0000000 --- a/test/test_types.jl +++ /dev/null @@ -1,6 +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 -