From ce3b6b9d19bf8716db2169ef657451ff2eb8bf3a Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 25 Feb 2016 10:44:55 +0200 Subject: [PATCH] 3d mortar autodiff not working --- src/elements.jl | 2 +- src/lagrange.jl | 7 +- src/mortar.jl | 12 ++ src/mortar_2d_autodiff.jl | 11 +- src/mortar_3d.jl | 89 ++++++---- src/mortar_3d_autodiff.jl | 312 +++++++++++++++++++++++++++++++++ src/preprocess_aster_reader.jl | 4 +- 7 files changed, 398 insertions(+), 39 deletions(-) create mode 100644 src/mortar_3d_autodiff.jl diff --git a/src/elements.jl b/src/elements.jl index 4d59834..28e03ea 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -124,7 +124,7 @@ function get_basis{E}(::Type{Element{E}}, xi::Vector{Float64}) return get_basis(E, xi) end -function get_basis{E}(element::Element{E}, xi::Vector{Float64}) +function get_basis{E}(element::Element{E}, xi::Vector) return get_basis(E, xi) end diff --git a/src/lagrange.jl b/src/lagrange.jl index 4a33f8c..af10ab0 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -141,10 +141,15 @@ end xi[2]^2, xi[3]^2, xi[1]*xi[2], xi[2]*xi[3], xi[3]*xi[1]]) -# some helpers to make accessing 1d basis functions more easily +# some helpers to make accessing 1d basis functions more easy function get_basis{T<:Real, E<:Union{Seg2,Seg3}}(::Type{E}, xi::T) get_basis(E, [xi]) end + function get_dbasis{T<:Real, E<:Union{Seg2,Seg3}}(::Type{E}, xi::T) get_dbasis(E, [xi]) end + +function get_reference_element_midpoint{E}(element::Element{E}) + get_reference_element_midpoint(E) +end diff --git a/src/mortar.jl b/src/mortar.jl index 9e3999b..fa3d06f 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -53,9 +53,21 @@ macro debug(msg) return msg end +function assemble!(problem::Problem{Mortar}, time::Real) + elements = get_elements(problem) + if length(elements) == 0 + info("$(typeof(problem)) : forget to add elements?") + return + end + # returns 3 if eldim 2 (tri3, quad4, ...) for 3d problems etc. + eldim = size(elements[1], 1)+1 + assemble!(problem, time, Val{eldim}) +end + include("mortar_2d.jl") include("mortar_2d_autodiff.jl") include("mortar_3d.jl") +include("mortar_3d_autodiff.jl") """ Remove inactive inequality constraints by using primal-dual active set strategy. """ function boundary_assembly_posthook!(solver::Solver, problem::Problem{Mortar}, C1, C2, D, g) diff --git a/src/mortar_2d_autodiff.jl b/src/mortar_2d_autodiff.jl index d8ea24b..a8e67df 100644 --- a/src/mortar_2d_autodiff.jl +++ b/src/mortar_2d_autodiff.jl @@ -73,7 +73,8 @@ function project_from_slave_to_master{E<:MortarElements2D}( end -function assemble!(problem::Problem{Mortar}, time::Real) +""" Assemble Mortar problem for two-dimensional problems, i.e. for Seg2 and Seg3 elements. """ +function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}) props = problem.properties field_dim = get_unknown_field_dimension(problem) @@ -99,7 +100,7 @@ function assemble!(problem::Problem{Mortar}, time::Real) push!(S, conn...) gdofs = get_gdofs(element, field_dim) X_el = element("geometry", time) - u_el = Field(Vector[u[:, i] for i in conn]) + u_el = Field(Vector[u[:,i] for i in conn]) x_el = X_el + u_el for ip in get_integration_points(element, Val{3}) dN = get_dbasis(element, ip) @@ -111,6 +112,7 @@ function assemble!(problem::Problem{Mortar}, time::Real) for i in 1:size(normals,2) normals[:,i] /= norm(normals[:,i]) end + # swap element normals in 2d if they point to inside of body if props.rotate_normals for i=1:size(normals,2) normals[:,i] = -normals[:,i] @@ -127,8 +129,8 @@ function assemble!(problem::Problem{Mortar}, time::Real) x1 = X1 + u1 la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) n1 = Field(Vector[normals[:,i] for i in slave_element_nodes]) - nnodes = size(slave_element, 2) + update!(slave_element, "normals", time => ForwardDiff.get_value(n1)) # 3. loop all master elements for master_element in slave_element["master elements"] @@ -202,7 +204,7 @@ function assemble!(problem::Problem{Mortar}, time::Real) gap[1,slave_element_nodes] += w*gn*Phi' #gap[1,slave_element_nodes] += w*gn*N1' - end + end # done integrating segment end # master elements done @@ -230,7 +232,6 @@ function assemble!(problem::Problem{Mortar}, time::Real) C[1,j] += gap[1, j] C[2,j] += lat else - #info("set node $j inactive") C[:,j] = la[:,j] end end diff --git a/src/mortar_3d.jl b/src/mortar_3d.jl index b4194e4..58591a9 100644 --- a/src/mortar_3d.jl +++ b/src/mortar_3d.jl @@ -45,21 +45,25 @@ Notes [1](http://stackoverflow.com/questions/8942950/how-do-i-find-the-orthogonal-projection-of-a-point-onto-a-plane) """ -function project_point_to_auxiliary_plane(p::Vector, x0::Vector, Q::Matrix) +function project_vertex_to_auxiliary_plane(p::Vector, x0::Vector, Q::Matrix) n = Q[:,1] ph = p - dot(p-x0, n)*n qproj = Q'*(ph-x0) - if !isapprox(qproj[1], 0.0; atol=1.0e-12) - info("project_point_to_auxiliary_plane(): point not projected correctly.") - info("p: $p") - info("x0: $x0") - info("Q: \n$Q") - info("qproj: $qproj") - error("Failed to project point to auxiliary plane.") + if abs(qproj[1]) > 1.0e-2 + # we should have something very little for normal direction if projected + # properly + info("project_point_to_auxiliary_plane(): vertex not projected correctly.") + info("p: $(ForwardDiff.get_value(p))") + info("x0: $(ForwardDiff.get_value(x0))") + info("Q: \n$(ForwardDiff.get_value(Q))") + info("qproj: $(ForwardDiff.get_value(qproj))") + error("Failed to project vertex to auxiliary plane.") end return qproj[2:3] end +project_point_to_auxiliary_plane = project_vertex_to_auxiliary_plane + """ Find edge intersections of two planar arbitrary shape polygons. @@ -114,7 +118,7 @@ function get_edge_intersections(S::Matrix, M::Matrix) for j=1:nm b = M[:,j]-S[:,i] A = [S[:,mod(i,ns)+1]-S[:,i] -M[:,mod(j,nm)+1]+M[:,j]] - if rank(A) == 2 + if rank(ForwardDiff.get_value(A)) == 2 r = A\b if (r[1]>=0) & (r[1]<=1) & (r[2]>=0) & (r[2]<=1) # intersection found k += 1 @@ -206,12 +210,12 @@ end """ function uniquetol(P, dim::Int; args...) @assert dim == 2 - items = Vector{Float64}[P[:,i] for i=1:size(P,dim)] - new_items = Vector{Float64}[] + items = Vector[P[:,i] for i=1:size(P,dim)] + new_items = Vector[] for item in items has_found = false for new_item in new_items - if isapprox(item, new_item; args...) + if isapprox(ForwardDiff.get_value(item), ForwardDiff.get_value(new_item); args...) has_found = true break end @@ -317,7 +321,7 @@ function calculate_polygon_centerpoint(P::Matrix) Cx += 1/(6*A)*(P[1,i] + P[1,inext])*(P[1,i]*P[2,inext] - P[1,inext]*P[2,i]) Cy += 1/(6*A)*(P[2,i] + P[2,inext])*(P[1,i]*P[2,inext] - P[1,inext]*P[2,i]) end - return Float64[Cx, Cy] + return [Cx, Cy] end """ @@ -401,33 +405,58 @@ julia> xquad*basis(theta[2:3]) 1.0 """ -function project_point_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix, element::Element{E}, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9) +function project_point_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix, + element::Element{E}, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9) + x = element("geometry", time) + return project_point_from_plane_to_surface(p, x0, Q, element, x, time; + max_iterations=max_iterations, iter_tol=iter_tol) +end +function project_vertex_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix, + element::Element{E}, x, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9) basis(xi) = get_basis(E, xi) dbasis(xi) = get_dbasis(E, xi) - x = element("geometry", time) ph = Q*[0; p] + x0 - theta = Float64[0.0, 0.0, 0.0] n = Q[:,1] + b(theta) = ph + theta[1]*n - basis(theta[2:3])*x + J(theta) = [n -dbasis(theta[2:3])*x] + theta = zeros(3) + dtheta = zeros(3) for i=1:max_iterations - b = ph + theta[1]*n - basis(theta[2:3])*x - J = [n -dbasis(theta[2:3])*x] - dtheta = J \ -b + # FIXME: gives NaN if partials in J + dtheta = ForwardDiff.get_value(J(theta)) \ -b(theta) + theta += dtheta + if norm(ForwardDiff.get_value(dtheta)) < iter_tol + return theta + end + end + info("failed to project vertex from auxiliary plane back to surface") + info("element type: $E") + info("element connectivity: $(get_connectivity(element))") + info("auxiliary plane: x0 = $(ForwardDiff.get_value(x0)), Q = $(ForwardDiff.get_value(Q))") + info("point coordinates on plane: $(ForwardDiff.get_value(p))") + info("element geometry: $(ForwardDiff.get_value(x.data))") + info("ph: $(ForwardDiff.get_value(ph))") + info("normal direction: $(ForwardDiff.get_value(n))") + info("parameter vector before giving up: $(ForwardDiff.get_value(theta))") + info("increment in parameter vector before giving up: $(ForwardDiff.get_value(dtheta))") + info("b([0.0, 0.0, 0.0]) = $(ForwardDiff.get_value(b([0.0, 0.0, 0.0])))") + info("J([0.0, 0.0, 0.0]) = $(ForwardDiff.get_value(J([0.0, 0.0, 0.0])))") + + info("iterations were") + theta = zeros(3) + dtheta = zeros(3) + for i=1:max_iterations + info("iter $i, theta = $(ForwardDiff.get_value(theta))") + info("b = $(ForwardDiff.get_value(b(theta)))") + info("J = $(ForwardDiff.get_value(J(theta)))") + dtheta = ForwardDiff.get_value(J(theta)) \ -b(theta) + info("dtheta = $(ForwardDiff.get_value(dtheta))") theta += dtheta if norm(dtheta) < iter_tol return theta end end - begin - info("projecting point from auxiliary plane back to surface didn't go very well.") - info("element type: $E") - info("element connectivity: $(get_connectivity(element))") - info("auxiliary plane: x0 = $x0, Q = $Q") - info("point coordinates on plane: $p") - info("element geometry: $x") - info("ph: $ph") - info("normal direction: $n") - info("parameter vector before giving up: $theta") - end + error("project_point_to_surface: did not converge in $max_iterations iterations!") end diff --git a/src/mortar_3d_autodiff.jl b/src/mortar_3d_autodiff.jl new file mode 100644 index 0000000..0275f08 --- /dev/null +++ b/src/mortar_3d_autodiff.jl @@ -0,0 +1,312 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +""" Assemble Mortar problem for three-dimensional problems, i.e. for Tri3, Tri6, Quad4, Quad8, Quad9 elements. """ +function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{3}}) + + props = problem.properties + field_dim = get_unknown_field_dimension(problem) + field_name = get_parent_field_name(problem) + + function calculate_interface(x::Vector) + + ndofs = round(Int, length(x)/2) # x = [u; la] + nnodes = round(Int, ndofs/field_dim) + u = reshape(x[1:ndofs], field_dim, nnodes) + la = reshape(x[ndofs+1:end], field_dim, nnodes) + fc = zeros(u) + gap = zeros(u) + gap_added = zeros(size(u)...) + C = zeros(la) + all_slave_nodes = Set{Int64}() + + # 1. calculate and average node normals for slave element nodes + normal = zeros(u) + tangent1 = zeros(u) + tangent2 = zeros(u) + for element in get_elements(problem) + haskey(element, "master elements") || continue + conn = get_connectivity(element) + push!(all_slave_nodes, conn...) + gdofs = get_gdofs(element, field_dim) + X_el = element("geometry", time) + u_el = Field(Vector[u[:,i] for i in conn]) + x_el = X_el + u_el + for ip in get_integration_points(element, Val{3}) + dN = get_dbasis(element, ip) + N = element(ip, time) + j = transpose(sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])) + #info("size of j = $(size(j))") + n = reshape(cross(j[:,1], j[:,2]), 3, 1) + normal[:, conn] += ip.weight*n*N + end + end + # calculate tangents + for i in 1:size(normal, 2) + i in all_slave_nodes || continue + normal[:,i] /= norm(normal[:,i]) + u1 = normal[:,i] + j = indmax(abs(u1)) + v2 = zeros(3) + v2[mod(j,3)+1] = 1.0 + u2 = v2 - dot(u1,v2)/dot(v2,v2)*v2 + u3 = cross(u1,u2) + tangent1[:,i] = u2/norm(u2) + tangent2[:,i] = u3/norm(u3) + end + if props.rotate_normals + for i=1:size(normal, 2) + normal[:,i] = -normal[:,i] + end + end + + # 2. loop slave elements and find contact segments + for slave_element in get_elements(problem) + haskey(slave_element, "master elements") || continue + + slave_element_nodes = get_connectivity(slave_element) + X1 = slave_element("geometry", time) + u1 = Field(Vector[u[:,i] for i in slave_element_nodes]) + x1 = X1 + u1 + la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) + n1 = Field(Vector[normal[:,i] for i in slave_element_nodes]) + t1 = Field(Vector[tangent1[:,i] for i in slave_element_nodes]) + t2 = Field(Vector[tangent2[:,i] for i in slave_element_nodes]) + nnodes = size(slave_element, 2) + update!(slave_element, "normals", time => ForwardDiff.get_value(n1.data)) + + # create auxiliary plane (x0, Q) + xi = get_reference_element_midpoint(slave_element) + N = vec(get_basis(slave_element, xi)) + x0 = N*x1 + Q = [N*n1 N*t1 N*t2] + + # project slave nodes to auxiliary plane + S = hcat([project_vertex_to_auxiliary_plane(p, x0, Q) for p in x1]...) + + # 3. loop all master elements + for master_element in slave_element["master elements"] + + master_element_nodes = get_connectivity(master_element) + X2 = master_element("geometry", time) + u2 = Field(Vector[u[:,i] for i in master_element_nodes]) + x2 = X2 + u2 + + x1_midpoint = mean(x1) + x2_midpoint = mean(x2) + distance = ForwardDiff.get_value(norm(x2_midpoint - x1_midpoint)) + distance > props.maximum_distance && continue + + # project master nodes to auxiliary plane + M = hcat([project_vertex_to_auxiliary_plane(p, x0, Q) for p in x2]...) + + # create polygon clipping on auxiliary plane +#= + P = nothing + neighbours = nothing + try + catch + info("polygon clipping failed") + info("S = ") + dump(ForwardDiff.get_value(S)) + info("M = ") + dump(ForwardDiff.get_value(M)) + error("cannot continue") + end +=# + P, neighbours = clip_polygon(S, M) + isa(P, Void) && continue # no clipping + info("polygon clip found: S = $(ForwardDiff.get_value(S)), M = $(ForwardDiff.get_value(M)), P = $(ForwardDiff.get_value(P))") + # special case, shared edge but no shared volume + size(P, 2) < 3 && continue + gap_added += 1 + + # clip polygon centerpoint + #C0 = calculate_polygon_centerpoint(P) + C0 = vec(mean(P, 2)) + npts = size(P, 2) # number of vertices in polygon + for pnt=1:npts # loop integration cells + x_cell = Field(Vector[C0, P[:,pnt], P[:,mod(pnt,npts)+1]]) +#= + try + catch + info("centerpoint: $(ForwardDiff.get_value(C0))") + info("P1 = $(ForwardDiff.get_value(P[:,pnt]))") + info("P2 = $(ForwardDiff.get_value(P[:,mod(pnt,npts)+1]))") + data = Vector[C0, P[:,pnt], P[:,mod(pnt,npts)+1]] + info("data = $(ForwardDiff.get_value(data))") + rethrow() + end +=# + # create dual basis + De = zeros(nnodes, nnodes) + Me = zeros(nnodes, nnodes) + for ip in get_integration_points(Tri3, Val{5}) + N = vec(get_basis(Tri3, ip.xi)) + x_g = N*x_cell + theta = zeros(3) + try + theta = project_vertex_from_plane_to_surface(x_g, x0, Q, slave_element, x1, time) + catch + info("creating dual basis did fail for finding projection back to slave surface.") + info("cell coords in auxiliary plane are:") + info(ForwardDiff.get_value(x_cell.data)) + info("interpolated value x_g is $(ForwardDiff.get_value(x_g))") + info("clip polygon centerpoint is $(ForwardDiff.get_value(C0))") + info("clip polygon is $(ForwardDiff.get_value(P))") + info("slave side nodes projected onto a plane are $(ForwardDiff.get_value(S))") + info("master side nodes projected onto a plane are $(ForwardDiff.get_value(M))") + rethrow() + end + xi_slave = theta[2:3] + N1 = slave_element(xi_slave, time) + + # jacobian determinant on integration cell + dNC = get_dbasis(Tri3, ip.xi) + JC = sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)]) + wC = ip.weight*det(JC) + + De += wC*diagm(vec(N1)) + Me += wC*N1'*N1 + end + Ae = De*inv(Me) + + # loop integration points of cell + for ip in get_integration_points(Tri3, Val{5}) + N = vec(get_basis(Tri3, ip.xi)) + x_g = N*x_cell + # project gauss point back to element surfaces + theta1 = project_vertex_from_plane_to_surface(x_g, x0, Q, slave_element, x1, time) + theta2 = project_vertex_from_plane_to_surface(x_g, x0, Q, master_element, x2, time) + xi_slave = theta1[2:3] + xi_master = theta2[2:3] + + # evaluate shape functions, calculate contact force and gap + N1 = vec(get_basis(slave_element, xi_slave)) + N2 = vec(get_basis(master_element, xi_master)) + Phi = Ae*N1 + + # jacobian determinant of integration cell + dNC = get_dbasis(Tri3, ip.xi) + JC = sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)]) + wC = ip.weight*det(JC) + + x_s = N1*x1 + n_s = N1*n1 + x_m = N2*x2 + la_s = Phi*la1 + gn = props.gap_sign*dot(n_s, x_s - x_m) + fc[:,slave_element_nodes] += wC*la_s*N1' + fc[:,master_element_nodes] -= wC*la_s*N2' + wg = wC*gn*Phi' + if any(isnan(wg)) + info("gap has NaNs!") + info("wC = $(ForwardDiff.get_value(wC))") + info("gn = $(ForwardDiff.get_value(gn))") + info("Phi = $(ForwardDiff.get_value(Phi))") + info("xi_slave = $(ForwardDiff.get_value(xi_slave))") + info("N1 = $(ForwardDiff.get_value(N1))") + info("Ae = $(ForwardDiff.get_value(Ae))") + info("De = $(ForwardDiff.get_value(De))") + info("Me = $(ForwardDiff.get_value(Me))") + info("x_cell(data) = $(ForwardDiff.get_value(x_cell.data))") + info("C0 = $(ForwardDiff.get_value(C0))") + info("P = $(ForwardDiff.get_value(P))") + error("fix this") + end + gap[1,slave_element_nodes] += wg + gap_added[1,slave_element_nodes] += 1 + end # done integrating cell + + end # done for all cells in this segment + + end # done all master elements for this slave element + + end # done all slave elements + + # like in 2d, check contact in nodes based on a complementarity condition + + nzgap = sort(nonzeros(sparse(ForwardDiff.get_value(gap)))) + all_slave_nodes = sort(collect(all_slave_nodes)) + info("gap: $nzgap") + info("size of normal = $(size(normal))") + info("size of la = $(size(la))") + info("size of C = $(size(C))") + info("S = $all_slave_nodes") + + for (i, j) in enumerate(all_slave_nodes) + if j in props.always_inactive + info("special node $j always inactive") + C[:,j] = la[:,j] + continue + end + n = normal[:,j] + t1 = tangent1[:,j] + t2 = tangent2[:,j] + lan = dot(n, la[:,j]) + +# if lan - gap[1, j] > 0 + info("set node $j active, normal direction = $(ForwardDiff.get_value(n)), tangent plane = $(ForwardDiff.get_value(t1)) x $(ForwardDiff.get_value(t2))") + C[1,j] += gap[1, j] + C[2,j] += dot(t1, la[:,j]) + C[3,j] += dot(t2, la[:,j]) + # else + # C[:,j] = la[:,j] + # end + end + + for (i, j) in enumerate(all_slave_nodes) + Ci = ForwardDiff.get_value(C[:,j]) + gapi = ForwardDiff.get_value(gap[:,j]) + fci = ForwardDiff.get_value(fc[:,j]) + lai = ForwardDiff.get_value(la[:,j]) + ui = ForwardDiff.get_value(u[:,j]) + ni = ForwardDiff.get_value(normal[:,j]) + gai = gap_added[:,j] + info("$i/$j: C = $Ci, f = $fci, gap = $gapi, la = $lai, u = $ui, n = $ni, gai = $gai") + end + + return vec([fc C]) + end + + # x doesn't mean deformed configuration here + x = [problem.assembly.u; problem.assembly.la] + ndofs = round(Int, length(x)/2) + A, allresults = ForwardDiff.jacobian(calculate_interface, x, + ForwardDiff.AllResults, cache=autodiffcache) + b = -ForwardDiff.value(allresults) + + A = sparse(A) + b = sparse(b) + SparseMatrix.droptol!(A, 1.0e-12) + SparseMatrix.droptol!(b, 1.0e-12) + + K = A[1:ndofs,1:ndofs] + C1 = transpose(A[1:ndofs,ndofs+1:end]) + C2 = A[ndofs+1:end,1:ndofs] + D = A[ndofs+1:end,ndofs+1:end] + f = b[1:ndofs] + g = b[ndofs+1:end] + + slaves = [101,108,111,112,113,120,123,124,125,126,129,130,149,150,151,152] + for j in slaves + dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3] + info("slave node $j, dofs $dofs") + info("Stiffness: $(K[dofs,:])") + info("force fc: $(C1[dofs,:])") + info("constraint: $(C2[dofs,:])") + info("lambdas: $(D[dofs,:])") + info("f = $(f[dofs]), g = $(g[dofs])") + end + + empty!(problem.assembly) + add!(problem.assembly.K, K) + add!(problem.assembly.C1, C1) + add!(problem.assembly.C2, C2) + add!(problem.assembly.D, D) + add!(problem.assembly.f, f) + add!(problem.assembly.g, g) + + return problem.assembly +end + diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl index 2a0b1a7..94c509a 100644 --- a/src/preprocess_aster_reader.jl +++ b/src/preprocess_aster_reader.jl @@ -3,13 +3,13 @@ using HDF5 using JuliaFEM -using JuliaFEM.Core: Element, Quad4, Tri3, Seg2, Hex8, update! +using JuliaFEM.Core: Element, Quad4, Tri3, Tet4, Seg2, Hex8, update! # TODO: this should be elsewhere function aster_create_elements(mesh, element_set, element_type=nothing) elements = Element[] - mapping = Dict(:QU4 => Quad4, :TR3 => Tri3, :SE2 => Seg2, :HE8 => Hex8) + mapping = Dict(:QU4 => Quad4, :TR3 => Tri3, :SE2 => Seg2, :HE8 => Hex8, :TE4 => Tet4) for (elid, (eltype, elset, elcon)) in mesh["connectivity"] if !haskey(mapping, eltype) error("aster_create_elements: unknown element mapping $eltype")