From b3f053174688f8bbebec275037b66529f33dd2f6 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Tue, 21 Mar 2017 08:36:18 +0200 Subject: [PATCH] Add postprocessing features (#100) * refactored code for solvers. * Added elementary tests for least-squares fitting of strain and stress fields * A more realistic postprocess + Xdmf writing test * removed debug keyword argument from test * Rewrite update_xdmf! New function to update Xdmf file no longer takes Solver object but xdmf, problem, time and fields to write, for example julia> update_xdmf!(xdmf, problem, 0.0, ["displacement", "temperature"]) All problems are written separately and put together into one SpatialCollection, allowing to have more structured Xdmf and making it easier to write complicated field configurations. Support for Xdmf API 3.0 added. * Support for Tensor6 field writing * moved update_xdmf! to io.jl * Removed some empty files * Not use old Postprocessor, obsolete code. * Not use old XDMF (obsolete code). Fixed test. * removed some postprocessing to pass test, maybe we should drop abaqus.jl from code as obsolete * add function get_temporal_collection back, it's used by update_xdmf of modal solver * postprocess of boundary problems also * added test for contact pressure. dl+quad test output was written in wrong file, fixed. * postprocess for contact pressure * contact pressure postprocess * with boundary problems always store also the primary unknown field * Change "reaction force" -> "lambda" * testing postprocess of reaction force also * sign convention --- src/JuliaFEM.jl | 14 +- src/abaqus.jl | 4 +- src/api.jl | 5 - src/io.jl | 226 ++++++- src/postprocess_xdmf.jl | 194 ------ src/problems.jl | 30 +- src/problems_contact.jl | 16 - src/problems_contact_2d.jl | 6 +- src/problems_contact_3d.jl | 31 +- src/problems_dirichlet.jl | 57 +- src/problems_elasticity.jl | 77 +++ src/problems_mortar.jl | 12 - src/solvers.jl | 258 ++------ src/solvers_modal.jl | 22 +- ..._elasticity_2d_linear_with_surface_load.jl | 77 --- test/test_heat.jl | 57 +- test/test_io.jl | 109 +++- test/test_modal_analysis.jl | 5 +- test/test_mortar_2d_contact.jl | 4 +- test/test_mortar_2d_mesh_tie.jl | 4 +- test/test_mortar_3d_mesh_tie_modal.jl | 17 +- test/test_mortar_autodiff.jl | 50 -- test/test_nodal_constraints.jl | 4 +- test/test_postprocess_xdmf.jl | 102 --- test/test_problem.jl | 2 +- test/test_problems_contact_2d.jl | 28 +- test/test_problems_contact_2d_autodiff.jl | 4 +- test/test_problems_contact_3d.jl | 10 +- test/test_problems_elasticity.jl | 54 ++ test/test_solvers_postprocess.jl | 94 +++ test/test_solvers_postprocess/tet4.inp | 589 ++++++++++++++++++ test/test_types.jl | 6 - 32 files changed, 1305 insertions(+), 863 deletions(-) delete mode 100644 src/api.jl delete mode 100644 src/postprocess_xdmf.jl delete mode 100644 test/test_mortar_autodiff.jl delete mode 100644 test/test_postprocess_xdmf.jl create mode 100644 test/test_problems_elasticity.jl create mode 100644 test/test_solvers_postprocess.jl create mode 100644 test/test_solvers_postprocess/tet4.inp delete mode 100644 test/test_types.jl 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 -