mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
mortar tests, bug finding
This commit is contained in:
Jukka Aho
@@ -1,9 +1,7 @@
|
||||
julia 0.4
|
||||
BaseTestNext
|
||||
ForwardDiff
|
||||
Logging
|
||||
LightXML
|
||||
Docile
|
||||
Lexicon
|
||||
FactCheck
|
||||
HDF5
|
||||
Lint
|
||||
JLD
|
||||
|
||||
|
||||
File diff suppressed because one or more lines are too long
@@ -6,6 +6,8 @@ This is JuliaFEM -- Finite Element Package
|
||||
"""
|
||||
module JuliaFEM
|
||||
|
||||
# include("common.jl")
|
||||
|
||||
""" JuliaFEM Core module. """
|
||||
module Core
|
||||
include("core.jl")
|
||||
@@ -19,6 +21,10 @@ include("api.jl")
|
||||
end
|
||||
|
||||
module Preprocess
|
||||
macro debug(msg)
|
||||
haskey(ENV, "DEBUG") || return
|
||||
return msg
|
||||
end
|
||||
include("abaqus_reader.jl")
|
||||
include("preprocess_aster_reader.jl")
|
||||
end
|
||||
|
||||
@@ -0,0 +1,52 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# common routines
|
||||
|
||||
"""
|
||||
A very simple debugging macro. It executes commands if environment variable DEBUG is set.
|
||||
|
||||
Usage
|
||||
-----
|
||||
Instead of starting session `julia file.jl`, do `DEBUG=1 julia file.jl`.
|
||||
Or set `export DEBUG=1` for your `.bashrc`.
|
||||
|
||||
Running inside code
|
||||
-------------------
|
||||
|
||||
julia> @debug info("moimoi heihei")
|
||||
|
||||
will get executed iff environment variable DEBUG is set.
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> @debug info("moimoi")
|
||||
(empty)
|
||||
julia> ENV["DEBUG"] = 1
|
||||
julia> @debug info("moimoi")
|
||||
INFO: moimoi
|
||||
julia> @debug begin
|
||||
... info("matrix is")
|
||||
... dump([1 2; 3 4])
|
||||
... end
|
||||
INFO: matrix is
|
||||
Array(Int64(2,2)) 2x2 Array{Int64,2}:
|
||||
1 2
|
||||
3 4
|
||||
|
||||
"""
|
||||
macro debug(msg)
|
||||
haskey(ENV, "DEBUG") || return
|
||||
return msg
|
||||
end
|
||||
|
||||
function set_debug_on!()
|
||||
ENV["DEBUG"] = 1;
|
||||
end
|
||||
|
||||
function set_debug_off!()
|
||||
pop!(ENV, "DEBUG");
|
||||
end
|
||||
|
||||
export @debug, set_debug_on!, set_debug_off!
|
||||
|
||||
+2
-23
@@ -1,31 +1,10 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM
|
||||
|
||||
import Base: +, -, /, *, push!, convert, getindex, setindex!, length, similar, call, vec, endof, append!
|
||||
|
||||
"""
|
||||
A very simple debugging macro. It prints debug message if environment variable
|
||||
JULIAFEM_DEBUG is found.
|
||||
|
||||
Usage: instead of starting session `julia file.jl` do `DEBUG=1 julia file.jl`.
|
||||
Or set `export DEBUG=1` for your `.bashrc`.
|
||||
"""
|
||||
macro debug(msg)
|
||||
haskey(ENV, "DEBUG") || return
|
||||
# return :( println("DEBUG: ", $msg) )
|
||||
return msg
|
||||
end
|
||||
|
||||
function set_debug_on!()
|
||||
ENV["DEBUG"] = 1;
|
||||
end
|
||||
|
||||
function set_debug_off!()
|
||||
pop!(ENV, "DEBUG");
|
||||
end
|
||||
|
||||
export @debug, set_debug_on!, set_debug_off!
|
||||
|
||||
using ForwardDiff
|
||||
autodiffcache = ForwardDiffCache()
|
||||
# export derivative, jacobian, hessian
|
||||
|
||||
@@ -143,6 +143,7 @@ end
|
||||
function call(solver::DirectSolver, time::Number=0.0)
|
||||
info("# of field problems: $(length(solver.field_problems))")
|
||||
info("# of boundary problems: $(length(solver.boundary_problems))")
|
||||
(length(solver.field_problems) != 0) || error("no field problems defined.")
|
||||
|
||||
timing = Dict{ASCIIString, Float64}()
|
||||
tic(timing, "solver")
|
||||
|
||||
+5
-2
@@ -3,8 +3,11 @@
|
||||
|
||||
abstract DirichletProblem <: AbstractProblem
|
||||
|
||||
function DirichletProblem(parent_field_name, parent_field_dim, dim=1, elements=[])
|
||||
return BoundaryProblem{DirichletProblem}(parent_field_name, parent_field_dim, dim, elements)
|
||||
function DirichletProblem(parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=Element[])
|
||||
return BoundaryProblem{DirichletProblem}("dirichlet boundary", parent_field_name, parent_field_dim, dim, elements)
|
||||
end
|
||||
function DirichletProblem(problem_name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=Element[])
|
||||
return BoundaryProblem{DirichletProblem}(problem_name, parent_field_name, parent_field_dim, dim, elements)
|
||||
end
|
||||
|
||||
function assemble!(assembly::Assembly, problem::BoundaryProblem{DirichletProblem}, element::Element, time::Number)
|
||||
|
||||
@@ -264,3 +264,13 @@ function calculate_normal_tangential_coordinates!{E}(element::Element{E}, time::
|
||||
element["normal-tangential coordinates"] = ntcoords
|
||||
end
|
||||
|
||||
""" Pick values from nodes and set to element according to connectivity. """
|
||||
function update(element::Element, field_name::ASCIIString, data::Union{Vector, Dict})
|
||||
element[field_name] = [data[i] for i in get_connectivity(element)]
|
||||
end
|
||||
function update(elements::Vector{Element}, field_name::ASCIIString, data::Union{Vector, Dict})
|
||||
# info("update $field_name for $(length(elements)) elements.")
|
||||
for element in elements
|
||||
update(element, field_name, data)
|
||||
end
|
||||
end
|
||||
|
||||
@@ -78,6 +78,19 @@ function get_integration_points(::TriangularElement, ::Type{Val{2}})
|
||||
]
|
||||
end
|
||||
|
||||
function get_integration_points(::TriangularElement, ::Type{Val{4}})
|
||||
# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
|
||||
# FIXME: something wrong here with weights ..?
|
||||
[
|
||||
IntegrationPoint([0.44594849091597, 0.44594849091597], 0.5*0.22338158967801),
|
||||
IntegrationPoint([0.44594849091597, 0.10810301816807], 0.5*0.22338158967801),
|
||||
IntegrationPoint([0.10810301816807, 0.44594849091597], 0.5*0.22338158967801),
|
||||
IntegrationPoint([0.09157621350977, 0.09157621350977], 0.5*0.10995174365532),
|
||||
IntegrationPoint([0.09157621350977, 0.81684757298046], 0.5*0.10995174365532),
|
||||
IntegrationPoint([0.81684757298046, 0.09157621350977], 0.5*0.10995174365532)
|
||||
]
|
||||
end
|
||||
|
||||
function get_integration_points(::TriangularElement, ::Type{Val{5}})
|
||||
# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
|
||||
# FIXME: something wrong here with weights ..?
|
||||
|
||||
@@ -5,8 +5,8 @@
|
||||
|
||||
abstract LinearElasticityProblem <: ElasticityProblem
|
||||
|
||||
function LinearElasticityProblem(dim::Int=3, elements=[])
|
||||
return Problem{LinearElasticityProblem}(dim, elements)
|
||||
function LinearElasticityProblem(name="linear elasticity", dim::Int=3, elements=[])
|
||||
return Problem{LinearElasticityProblem}(name, dim, elements)
|
||||
end
|
||||
|
||||
""" Elasticity equations, general 3D case. """
|
||||
@@ -57,13 +57,22 @@ function assemble!{E<:CG, P<:LinearElasticityProblem}(assembly::Assembly, proble
|
||||
L = w*T*N*norm(cross(JT[:,1], JT[:,2]))
|
||||
add!(assembly.force_vector, gdofs, vec(L))
|
||||
end
|
||||
for dim in 1:problem.dim
|
||||
if haskey(element, "displacement traction force $dim")
|
||||
T = element("displacement traction force $dim", ip, time)
|
||||
ldofs = gdofs[dim:problem.dim:end]
|
||||
JT = transpose(J)
|
||||
L = w*T*N*norm(cross(JT[:,1], JT[:,2]))
|
||||
add!(assembly.force_vector, ldofs, vec(L))
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
abstract PlaneStressLinearElasticityProblem <: LinearElasticityProblem
|
||||
|
||||
function PlaneStressLinearElasticityProblem(dim::Int=2, elements=[])
|
||||
return Problem{PlaneStressLinearElasticityProblem}(dim, elements)
|
||||
function PlaneStressLinearElasticityProblem(name="plane stress linear elasticity", dim::Int=2, elements=[])
|
||||
return Problem{PlaneStressLinearElasticityProblem}(name, dim, elements)
|
||||
end
|
||||
|
||||
""" Elasticity equations, plane stress. """
|
||||
|
||||
+151
-35
@@ -3,6 +3,13 @@
|
||||
|
||||
# Mortar projection calculation for 2d
|
||||
|
||||
|
||||
macro debug(msg)
|
||||
haskey(ENV, "DEBUG") || return
|
||||
return msg
|
||||
end
|
||||
|
||||
|
||||
""" Find projection from slave nodes to master element, i.e. find xi2 from
|
||||
master element corresponding to the xi1.
|
||||
"""
|
||||
@@ -617,7 +624,18 @@ function project_point_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix
|
||||
return theta
|
||||
end
|
||||
end
|
||||
error("project_point_to_auxiliary_plane: did not converge in $max_iterations iterations!")
|
||||
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
|
||||
|
||||
|
||||
@@ -632,8 +650,12 @@ node_csys
|
||||
"""
|
||||
abstract MortarProblem <: AbstractProblem
|
||||
|
||||
function MortarProblem(parent_field_name, parent_field_dim, dim=1, elements=[])
|
||||
return BoundaryProblem{MortarProblem}(parent_field_name, parent_field_dim, dim, elements)
|
||||
function MortarProblem(parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=[])
|
||||
return BoundaryProblem{MortarProblem}("mortar problem", parent_field_name, parent_field_dim, dim, elements)
|
||||
end
|
||||
|
||||
function MortarProblem(problem_name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=[])
|
||||
return BoundaryProblem{MortarProblem}(problem_name, parent_field_name, parent_field_dim, dim, elements)
|
||||
end
|
||||
|
||||
# Mortar assembly
|
||||
@@ -688,6 +710,31 @@ end
|
||||
|
||||
typealias MortarElements3D Union{Tri3, Quad4}
|
||||
|
||||
""" Find master elements from list of potential master elements. """
|
||||
function find_master_elements(slave_element::Element, time::Real)
|
||||
x0, Q = create_auxiliary_plane(slave_element, time)
|
||||
Sl = Vector{Float64}[]
|
||||
for p in slave_element("geometry", time)
|
||||
push!(Sl, project_point_to_auxiliary_plane(p, x0, Q))
|
||||
end
|
||||
S = hcat(Sl...)
|
||||
master_elements = Element[]
|
||||
|
||||
for master_element in slave_element["master elements"]
|
||||
M = Vector{Float64}[]
|
||||
for p in master_element("geometry", time)
|
||||
push!(M, project_point_to_auxiliary_plane(p, x0, Q))
|
||||
end
|
||||
M = hcat(M...)
|
||||
P, neighbours = clip_polygon(S, M)
|
||||
isa(P, Void) && continue # no clipping
|
||||
size(P, 2) < 3 && continue # shared edge, no contribution
|
||||
push!(master_elements, master_element)
|
||||
end
|
||||
|
||||
return master_elements
|
||||
end
|
||||
|
||||
function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryProblem{MortarProblem}, slave_element::Element{E}, time::Real)
|
||||
field_dim = problem.parent_field_dim
|
||||
field_name = problem.parent_field_name
|
||||
@@ -702,11 +749,14 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro
|
||||
for p in slave_element("geometry", time)
|
||||
push!(Sl, project_point_to_auxiliary_plane(p, x0, Q))
|
||||
end
|
||||
@debug info("auxiliary plane coords and basis: origo = $x0")
|
||||
@debug info("basis:")
|
||||
@debug dump(round(Q, 3))
|
||||
@debug begin
|
||||
info("auxiliary plane coords and basis: origo = $x0")
|
||||
info("basis:")
|
||||
dump(round(Q, 3))
|
||||
end
|
||||
#S = reshape([S...;], 2, size(slave_element)[2])
|
||||
S = hcat(Sl...)
|
||||
slave_geom = Field(Vector{Float64}[S[:,j] for j=1:size(S,2)])
|
||||
|
||||
for master_element in slave_element["master elements"]
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
@@ -717,11 +767,16 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro
|
||||
end
|
||||
#M = reshape([M...;], 2, size(master_element)[2])
|
||||
M = hcat(M...)
|
||||
master_geom = Field(Vector{Float64}[M[:,j] for j=1:size(M,2)])
|
||||
|
||||
P = nothing
|
||||
neighbours = nothing
|
||||
@debug info("applying polygon clip algorithm, S & M = ")
|
||||
@debug dump(round(S, 3))
|
||||
@debug dump(round(M, 3))
|
||||
@debug begin
|
||||
info("applying polygon clip algorithm, S & M = ")
|
||||
dump(round(S, 3))
|
||||
dump(round(M, 3))
|
||||
end
|
||||
|
||||
try
|
||||
P, neighbours = clip_polygon(S, M)
|
||||
catch
|
||||
@@ -735,53 +790,114 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro
|
||||
error("cannot continue")
|
||||
end
|
||||
isa(P, Void) && continue # no clipping
|
||||
@debug info("polygon on auxilyary plane: ")
|
||||
@debug dump(round(P, 3))
|
||||
C = calculate_polygon_centerpoint(P)
|
||||
@debug info("center point = $C")
|
||||
@debug begin
|
||||
info("polygon coords on auxilyary plane: ")
|
||||
dump(round(P, 3))
|
||||
end
|
||||
|
||||
if size(P, 2) < 3
|
||||
# shared edge but no shared volume. skipping
|
||||
continue
|
||||
info("this is not polygon at all.")
|
||||
info("clipping S")
|
||||
dump(S)
|
||||
info("clipping M")
|
||||
dump(M)
|
||||
error("size(P, 2) < 3")
|
||||
end
|
||||
C = calculate_polygon_centerpoint(P)
|
||||
npts = size(P, 2) # number of vertices in polygon
|
||||
@debug info("number of vectices in polygon: $npts")
|
||||
# S = zeros(3, 3)
|
||||
# M = zeros(3, 3)
|
||||
|
||||
@debug begin
|
||||
info("clip polygon info")
|
||||
theta = project_point_from_plane_to_surface(C, x0, Q, slave_element, time)
|
||||
CC = slave_element("geometry", theta[2:3], time)
|
||||
info("center point on slave: $CC")
|
||||
info("number of vectices in polygon: $npts")
|
||||
on_slave = zeros(3, 0)
|
||||
on_master = zeros(3, 0)
|
||||
for i=1:size(P, 2)
|
||||
theta = project_point_from_plane_to_surface(P[:,i], x0, Q, slave_element, time)
|
||||
on_slave = [on_slave slave_element("geometry", theta[2:3], time)]
|
||||
theta = project_point_from_plane_to_surface(P[:,i], x0, Q, master_element, time)
|
||||
on_master = [on_master master_element("geometry", theta[2:3], time)]
|
||||
end
|
||||
info("polygon coords projected to slave element")
|
||||
dump(round(on_slave, 3))
|
||||
info("polygon coords projected to master element")
|
||||
dump(round(on_master, 3))
|
||||
end
|
||||
|
||||
for i=1:npts # loop vertices and create temporary integrate cells
|
||||
xvec = [C[1], P[1, i], P[1, mod(i, npts)+1]]
|
||||
yvec = [C[2], P[2, i], P[2, mod(i, npts)+1]]
|
||||
X = hcat(xvec, yvec)'
|
||||
@debug info("cell $i, coords = ")
|
||||
@debug dump(round(X, 3))
|
||||
geom = Field(Vector{Float64}[X[:,j] for j=1:size(X,2)])
|
||||
|
||||
@debug begin
|
||||
on_slave = zeros(3, 0)
|
||||
on_master = zeros(3, 0)
|
||||
for j=1:size(X, 2)
|
||||
theta = project_point_from_plane_to_surface(X[:,j], x0, Q, slave_element, time)
|
||||
on_slave = [on_slave slave_element("geometry", theta[2:3], time)]
|
||||
theta = project_point_from_plane_to_surface(X[:,j], x0, Q, master_element, time)
|
||||
on_master = [on_master master_element("geometry", theta[2:3], time)]
|
||||
end
|
||||
info("cell $i coords projected to slave element")
|
||||
dump(round(on_slave, 3))
|
||||
info("cell $i coords projected to master element")
|
||||
dump(round(on_master, 3))
|
||||
end
|
||||
|
||||
# integration cell geometry, i.e., Tri3
|
||||
cell = Field(Vector{Float64}[X[:,j] for j=1:size(X,2)])
|
||||
# info("geom = $geom")
|
||||
for ip in get_integration_points(Tri3, Val{5})
|
||||
# calculate determiant of jacobian
|
||||
#dN = get_dbasis(E, ip.xi)
|
||||
dN = get_dbasis(Tri3, ip.xi)
|
||||
J = sum([kron(dN[:,j], geom[j]') for j=1:length(geom)])
|
||||
w = ip.weight*det(J)
|
||||
# gauss point in auxiliary plane
|
||||
#N = get_basis(E, ip.xi)
|
||||
N = get_basis(Tri3, ip.xi)
|
||||
x = vec(N*geom)
|
||||
xi = vec(N*cell) # xi defined in auxilary plane
|
||||
#xi = ip.xi
|
||||
# info("x = $x")
|
||||
# find projection of gauss point to master and slave elements
|
||||
theta1 = project_point_from_plane_to_surface(x, x0, Q, slave_element, time)
|
||||
theta2 = project_point_from_plane_to_surface(x, x0, Q, master_element, time)
|
||||
theta1 = project_point_from_plane_to_surface(xi, x0, Q, slave_element, time)
|
||||
theta2 = project_point_from_plane_to_surface(xi, x0, Q, master_element, time)
|
||||
xi_slave = theta1[2:3]
|
||||
xi_master = theta2[2:3]
|
||||
@debug begin
|
||||
X_slave = slave_element("geometry", xi_slave, time)
|
||||
X_master = master_element("geometry", xi_master, time)
|
||||
info("integration point on slave: $xi_slave => $X_slave")
|
||||
info("integration point on master: $xi_master => $X_master")
|
||||
end
|
||||
# evaluate shape functions values in gauss point and add contribution to matrices
|
||||
N1 = slave_element(theta1[2:3], time)
|
||||
N2 = master_element(theta2[2:3], time)
|
||||
Sm = w*N1'*N1
|
||||
N1 = slave_element(xi_slave, time)
|
||||
N2 = master_element(xi_master, time)
|
||||
|
||||
# calculate determiant of jacobian
|
||||
dNC = get_dbasis(Tri3, ip.xi)
|
||||
dNS = get_dbasis(Quad4, xi_slave)
|
||||
dNM = get_dbasis(Quad4, xi_master)
|
||||
JC = sum([kron(dNC[:,j], cell[j]') for j=1:length(cell)])
|
||||
JN = sum([kron(dNS[:,j], slave_geom[j]') for j=1:length(slave_geom)])
|
||||
JM = sum([kron(dNM[:,j], master_geom[j]') for j=1:length(master_geom)])
|
||||
wS = det(JN)
|
||||
wM = det(JM)
|
||||
wC = det(JC)
|
||||
@debug info("weight S = $wS, weight M = $wM, weight C = $wC")
|
||||
|
||||
Sm = ip.weight*N1'*N1*wC
|
||||
# FIXME: master side transpose -- why?
|
||||
Mm = w*(N1'*N2)'
|
||||
Mm = ip.weight*(N1'*N2)'*wC
|
||||
for k=1:field_dim
|
||||
sd = slave_dofs[k:field_dim:end]
|
||||
md = master_dofs[k:field_dim:end]
|
||||
add!(assembly.stiffness_matrix, sd, sd, Sm)
|
||||
add!(assembly.stiffness_matrix, sd, md, -Mm)
|
||||
# info("sd = $sd")
|
||||
# info("md = $md")
|
||||
end
|
||||
end
|
||||
# info("breaking on first")
|
||||
# break
|
||||
end
|
||||
# info("S = \n$S")
|
||||
# info("M = \n$M")
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
@@ -1,6 +1,8 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM
|
||||
|
||||
using HDF5
|
||||
|
||||
function aster_parse_nodes(section::ASCIIString; strip_characters=true)
|
||||
@@ -53,6 +55,144 @@ function parse(mesh::ASCIIString, ::Type{Val{:CODE_ASTER_MAIL}})
|
||||
return model
|
||||
end
|
||||
|
||||
|
||||
function aster_renumber_nodes_!(mesh, node_numbering)
|
||||
old_nodes = mesh["nodes"]
|
||||
new_nodes = typeof(old_nodes)()
|
||||
for (node_id, node_coords) in old_nodes
|
||||
new_node_id = node_numbering[node_id]
|
||||
new_nodes[new_node_id] = node_coords
|
||||
end
|
||||
mesh["nodes"] = new_nodes
|
||||
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
|
||||
new_elcon = [node_numbering[node_id] for node_id in elcon]
|
||||
mesh["connectivity"][elid] = (eltype, elset, new_elcon)
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
function aster_renumber_nodes!(mesh1, mesh2)
|
||||
|
||||
reserved_node_ids = Set(collect(keys(mesh1["nodes"])))
|
||||
@debug info("already reserved node ids: $reserved_node_ids")
|
||||
mesh2_node_numbering = Dict{Int64, Int64}()
|
||||
|
||||
# find new node ids assigned for mesh 2
|
||||
k = 1
|
||||
for node_id in sort(collect(keys(mesh2["nodes"])))
|
||||
@debug info("mesh2: processing node $node_id")
|
||||
# if node id is reserved in mesh 1, find new number
|
||||
if node_id in reserved_node_ids
|
||||
@debug info("node id conflict, $node_id already defined in mesh 1, renumbering")
|
||||
while k in reserved_node_ids
|
||||
k += 1
|
||||
end
|
||||
@debug info("mesh2: node $node_id -> $k")
|
||||
mesh2_node_numbering[node_id] = k
|
||||
push!(reserved_node_ids, k)
|
||||
else
|
||||
mesh2_node_numbering[node_id] = node_id
|
||||
end
|
||||
end
|
||||
@debug info("new node numering:")
|
||||
@debug println(mesh2_node_numbering)
|
||||
aster_renumber_nodes_!(mesh2, mesh2_node_numbering)
|
||||
|
||||
#=
|
||||
# create new nodes
|
||||
mesh2_old_nodes = mesh2["nodes"]
|
||||
mesh2_new_nodes = typeof(mesh2_old_nodes)()
|
||||
for (node_id, node_coords) in mesh2_old_nodes
|
||||
new_node_id = mesh2_node_numbering[node_id]
|
||||
mesh2_new_nodes[new_node_id] = node_coords
|
||||
end
|
||||
mesh2["nodes"] = mesh2_new_nodes
|
||||
|
||||
# update connectivity
|
||||
for (elid, (eltype, elset, elcon)) in mesh2["connectivity"]
|
||||
new_elcon = [mesh2_node_numbering[node_id] for node_id in elcon]
|
||||
mesh2["connectivity"][elid] = (eltype, elset, new_elcon)
|
||||
end
|
||||
=#
|
||||
|
||||
end
|
||||
|
||||
|
||||
function aster_renumber_elements!(mesh1, mesh2)
|
||||
|
||||
reserved_element_ids = Set(collect(keys(mesh1["connectivity"])))
|
||||
@debug info("already reserved element ids: $reserved_element_ids")
|
||||
mesh2_element_numbering = Dict{Int64, Int64}()
|
||||
|
||||
# find new element ids assigned for mesh 2
|
||||
k = 1
|
||||
for element_id in sort(collect(keys(mesh2["connectivity"])))
|
||||
@debug info("mesh2: processing element $element_id")
|
||||
# if node id is reserved in mesh 1, find new number
|
||||
if element_id in reserved_element_ids
|
||||
@debug info("element id conflict, $element_id already defined in mesh 1, renumbering")
|
||||
while k in reserved_element_ids
|
||||
k += 1
|
||||
end
|
||||
@debug info("mesh2: element $element_id -> $k")
|
||||
mesh2_element_numbering[element_id] = k
|
||||
push!(reserved_element_ids, k)
|
||||
else
|
||||
mesh2_element_numbering[element_id] = element_id
|
||||
end
|
||||
end
|
||||
@debug info("element numbering for mesh 2:")
|
||||
@debug info(mesh2_element_numbering)
|
||||
|
||||
# create new elements
|
||||
mesh2_old_elements = mesh2["connectivity"]
|
||||
mesh2_new_elements = typeof(mesh2_old_elements)()
|
||||
for (element_id, element_data) in mesh2_old_elements
|
||||
new_element_id = mesh2_element_numbering[element_id]
|
||||
mesh2_new_elements[new_element_id] = element_data
|
||||
end
|
||||
mesh2["connectivity"] = mesh2_new_elements
|
||||
|
||||
end
|
||||
|
||||
|
||||
function aster_combine_meshes(mesh1, mesh2)
|
||||
|
||||
# check that meshes are ready to be combined
|
||||
node_ids_mesh_1 = collect(keys(mesh1["nodes"]))
|
||||
node_ids_mesh_2 = collect(keys(mesh2["nodes"]))
|
||||
if length(intersect(node_ids_mesh_1, node_ids_mesh_2)) != 0
|
||||
error("nodes with same id number in both meshes, failed.")
|
||||
end
|
||||
element_ids_mesh_1 = collect(keys(mesh1["connectivity"]))
|
||||
element_ids_mesh_2 = collect(keys(mesh2["connectivity"]))
|
||||
if length(intersect(element_ids_mesh_1, element_ids_mesh_2)) != 0
|
||||
error("elements with same id number in both meshes, failed.")
|
||||
end
|
||||
@assert similar(mesh1) == similar(mesh2)
|
||||
@assert similar(mesh1["nodes"]) == similar(mesh2["nodes"])
|
||||
@assert similar(mesh1["connectivity"]) == similar(mesh2["connectivity"])
|
||||
|
||||
new_mesh = similar(mesh1)
|
||||
new_mesh["nodes"] = similar(mesh1["nodes"])
|
||||
new_mesh["connectivity"] = similar(mesh1["connectivity"])
|
||||
|
||||
for (node_id, node_coords) in mesh1["nodes"]
|
||||
new_mesh["nodes"][node_id] = node_coords
|
||||
end
|
||||
for (node_id, node_coords) in mesh2["nodes"]
|
||||
new_mesh["nodes"][node_id] = node_coords
|
||||
end
|
||||
for (element_id, element_data) in mesh1["connectivity"]
|
||||
new_mesh["connectivity"][element_id] = element_data
|
||||
end
|
||||
for (element_id, element_data) in mesh2["connectivity"]
|
||||
new_mesh["connectivity"][element_id] = element_data
|
||||
end
|
||||
return new_mesh
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Code Aster binary file (.med), which is exported from SALOME.
|
||||
"""
|
||||
|
||||
@@ -4,11 +4,13 @@
|
||||
abstract AbstractProblem
|
||||
|
||||
type Problem{T<:AbstractProblem}
|
||||
name :: ASCIIString
|
||||
dim :: Int
|
||||
elements :: Vector{Element}
|
||||
end
|
||||
|
||||
type BoundaryProblem{T<:AbstractProblem}
|
||||
name :: ASCIIString
|
||||
parent_field_name :: ASCIIString
|
||||
parent_field_dim :: Int
|
||||
dim :: Int
|
||||
@@ -36,3 +38,4 @@ end
|
||||
function Base.push!(problem::AllProblems, element::Element)
|
||||
push!(problem.elements, element)
|
||||
end
|
||||
|
||||
|
||||
+2
-3
@@ -1,7 +1,7 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using Base.Test
|
||||
using BaseTestNext
|
||||
|
||||
abstract TestResult
|
||||
|
||||
@@ -159,5 +159,4 @@ function print_test_statistics()
|
||||
return passed, failed, errors, critical
|
||||
end
|
||||
|
||||
export @test, run_test, print_test_statistics
|
||||
|
||||
export @test, @testset, @test_throws, run_test, print_test_statistics
|
||||
|
||||
+238
-2
@@ -5,7 +5,8 @@ module MortarTests
|
||||
|
||||
using JuliaFEM.Test
|
||||
|
||||
using JuliaFEM.Core: Element, Seg2, Quad4, Tri3, MortarProblem, Assembly, assemble!
|
||||
using JuliaFEM.Core: Element, Seg2, Quad4, Tri3, Hex8, MortarProblem, Assembly, assemble!,
|
||||
get_connectivity, update
|
||||
using JuliaFEM.Core: PlaneStressElasticityProblem, DirichletProblem, DirectSolver
|
||||
|
||||
# 2d stuff
|
||||
@@ -19,6 +20,7 @@ using JuliaFEM.Core: create_auxiliary_plane, project_point_to_auxiliary_plane,
|
||||
calculate_normal_tangential_coordinates!,
|
||||
is_point_inside_convex_polygon
|
||||
|
||||
using JuliaFEM.Core: LinearElasticityProblem
|
||||
|
||||
function get_test_2d_model()
|
||||
# this is hand calculated and given as an example in my thesis
|
||||
@@ -691,6 +693,240 @@ function test_assemble_3d_problem_quad4_3()
|
||||
@test isapprox(stiffness_matrix, B)
|
||||
|
||||
end
|
||||
test_assemble_3d_problem_quad4_3()
|
||||
#test_assemble_3d_problem_quad4_3()
|
||||
|
||||
|
||||
function test_3d_problem()
|
||||
nodes = Vector{Float64}[
|
||||
[0.0, 0.0, 0.0],
|
||||
[1.0, 0.0, 0.0],
|
||||
[1.0, 1.0, 0.0],
|
||||
[0.0, 1.0, 0.0],
|
||||
[0.0, 0.0, 0.5],
|
||||
[1.0, 0.0, 0.5],
|
||||
[1.0, 1.0, 0.5],
|
||||
[0.0, 1.0, 0.5],
|
||||
|
||||
[0.0, 0.0, 0.5],
|
||||
[1.0, 0.0, 0.5],
|
||||
[1.0, 1.0, 0.5],
|
||||
[0.0, 1.0, 0.5],
|
||||
[0.0, 0.0, 1.0],
|
||||
[1.0, 0.0, 1.0],
|
||||
[1.0, 1.0, 1.0],
|
||||
[0.0, 1.0, 1.0],
|
||||
]
|
||||
el1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8])
|
||||
el2 = Hex8([9, 10, 11, 12, 13, 14, 15, 16])
|
||||
sym121 = Quad4([1, 2, 3, 4])
|
||||
sym131 = Quad4([1, 2, 6, 5])
|
||||
sym132 = Quad4([9, 10, 14, 13])
|
||||
sym231 = Quad4([4, 1, 5, 8])
|
||||
sym232 = Quad4([12, 9, 13, 16])
|
||||
force = Quad4([14, 15, 16, 13])
|
||||
l2u = Quad4([5, 6, 7, 8])
|
||||
u2l = Quad4([9, 10, 11, 12])
|
||||
elements = Element[el1, el2, sym121, sym131, sym132, sym231, sym232, force, l2u, u2l]
|
||||
update(elements, "geometry", nodes)
|
||||
el1["youngs modulus"] = el2["youngs modulus"] = 900.0
|
||||
el1["poissons ratio"] = el2["poissons ratio"] = 0.25
|
||||
sym121["displacement 3"] = 0.0
|
||||
sym131["displacement 2"] = sym132["displacement 2"] = 0.0
|
||||
sym231["displacement 1"] = sym232["displacement 1"] = 0.0
|
||||
force["displacement traction force 3"] = -100.0
|
||||
l2u["master elements"] = Element[u2l]
|
||||
calculate_normal_tangential_coordinates!(l2u, 0.0)
|
||||
fb = LinearElasticityProblem("two elastic blocks")
|
||||
push!(fb, el1, el2, force)
|
||||
bc = DirichletProblem("symmetry boundaries", "displacement", 3)
|
||||
push!(bc, sym121, sym131, sym132, sym231, sym232)
|
||||
tie = MortarProblem("tie contact between bodies", "displacement", 3)
|
||||
push!(tie, l2u)
|
||||
solver = DirectSolver("solution of elasticity problem")
|
||||
push!(solver, fb)
|
||||
push!(solver, bc)
|
||||
push!(solver, tie)
|
||||
solver.nonlinear_problem = false
|
||||
solver.method = :UMFPACK
|
||||
call(solver, 0.0)
|
||||
X = el2("geometry", [1.0, 1.0, 1.0], 0.0)
|
||||
u = el2("displacement", [1.0, 1.0, 1.0], 0.0)
|
||||
info("displacement at $X = $u")
|
||||
@test isapprox(u, 1/36*[1, 1, -4])
|
||||
end
|
||||
#test_3d_problem()
|
||||
|
||||
#=
|
||||
@testset "plane quad4 projector tests" begin
|
||||
a = 1/2
|
||||
b = 1/3
|
||||
nodes = Dict{Int64, Vector{Float64}}(
|
||||
1 => [0.0, 0.0, 0.0],
|
||||
2 => [1/2, 0.0, 0.0],
|
||||
3 => [1.0, 0.0, 0.0],
|
||||
4 => [0.0, 1.0, 0.0],
|
||||
5 => [1/2, 1.0, 0.0],
|
||||
6 => [1.0, 1.0, 0.0],
|
||||
|
||||
7 => [0.0, 0.0, 0.0],
|
||||
8 => [1/3, 0.0, 0.0],
|
||||
9 => [2/3, 0.0, 0.0],
|
||||
10 => [1.0, 0.0, 0.0],
|
||||
11 => [0.0, 1/2, 0.0],
|
||||
12 => [1/3, 1/2, 0.0],
|
||||
13 => [2/3, 1/2, 0.0],
|
||||
14 => [1.0, 1/2, 0.0],
|
||||
15 => [0.0, 1.0, 0.0],
|
||||
16 => [1/3, 1.0, 0.0],
|
||||
17 => [2/3, 1.0, 0.0],
|
||||
18 => [1.0, 1.0, 0.0],
|
||||
)
|
||||
sel1 = Quad4([1, 2, 5, 4])
|
||||
sel2 = Quad4([2, 3, 6, 5])
|
||||
mel1 = Quad4([7, 8, 12, 11])
|
||||
mel2 = Quad4([8, 9, 13, 12])
|
||||
mel3 = Quad4([9, 10, 14, 13])
|
||||
mel4 = Quad4([11, 12, 16, 15])
|
||||
mel5 = Quad4([12, 13, 17, 16])
|
||||
mel6 = Quad4([13, 14, 18, 17])
|
||||
update(Element[sel1, sel2, mel1, mel2, mel3, mel4, mel5, mel6], "geometry", nodes)
|
||||
calculate_normal_tangential_coordinates!(sel1, 0.0)
|
||||
calculate_normal_tangential_coordinates!(sel2, 0.0)
|
||||
prob = MortarProblem("temperature", 1)
|
||||
push!(prob, sel1)
|
||||
push!(prob, sel2)
|
||||
|
||||
sel1["master elements"] = [mel1, mel2, mel4, mel5]
|
||||
sel2["master elements"] = [mel2, mel3, mel5, mel6]
|
||||
stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)*2592*6
|
||||
info("interface matrix:")
|
||||
dump(round(stiffness_matrix, 3))
|
||||
B = [
|
||||
864 432 0 432 216 0 -420 -375 -15 0 -504 -450 -18 0 -84 -75 -3 0
|
||||
432 1728 432 216 864 216 -120 -690 -690 -120 -144 -828 -828 -144 -24 -138 -138 -24
|
||||
0 432 864 0 216 432 0 -15 -375 -420 0 -18 -450 -504 0 -3 -75 -84
|
||||
432 216 0 864 432 0 -84 -75 -3 0 -504 -450 -18 0 -420 -375 -15 0
|
||||
216 864 216 432 1728 432 -24 -138 -138 -24 -144 -828 -828 -144 -120 -690 -690 -120
|
||||
0 216 432 0 432 864 0 -3 -75 -84 0 -18 -450 -504 0 -15 -375 -420
|
||||
]
|
||||
info("expected interface matrix:")
|
||||
dump(round(B, 3))
|
||||
@test isapprox(stiffness_matrix, B)
|
||||
end
|
||||
=#
|
||||
|
||||
@testset "plane quad4 projector master 3x3 slave 2x2" begin
|
||||
a = 1/2
|
||||
b = 1/3
|
||||
nodes = Dict{Int64, Vector{Float64}}(
|
||||
1 => [0*a, 0*a, 0.0],
|
||||
2 => [1*a, 0*a, 0.0],
|
||||
3 => [2*a, 0*a, 0.0],
|
||||
4 => [0*a, 1*a, 0.0],
|
||||
5 => [1*a, 1*a, 0.0],
|
||||
6 => [2*a, 1*a, 0.0],
|
||||
7 => [0*a, 2*a, 0.0],
|
||||
8 => [1*a, 2*a, 0.0],
|
||||
9 => [2*a, 2*a, 0.0],
|
||||
|
||||
10 => [0*b, 0*b, 0.0],
|
||||
11 => [1*b, 0*b, 0.0],
|
||||
12 => [2*b, 0*b, 0.0],
|
||||
13 => [3*b, 0*b, 0.0],
|
||||
14 => [0*b, 1*b, 0.0],
|
||||
15 => [1*b, 1*b, 0.0],
|
||||
16 => [2*b, 1*b, 0.0],
|
||||
17 => [3*b, 1*b, 0.0],
|
||||
18 => [0*b, 2*b, 0.0],
|
||||
19 => [1*b, 2*b, 0.0],
|
||||
20 => [2*b, 2*b, 0.0],
|
||||
21 => [3*b, 2*b, 0.0],
|
||||
22 => [0*b, 3*b, 0.0],
|
||||
23 => [1*b, 3*b, 0.0],
|
||||
24 => [2*b, 3*b, 0.0],
|
||||
25 => [3*b, 3*b, 0.0],
|
||||
)
|
||||
sel1 = Quad4([1, 2, 5, 4])
|
||||
sel2 = Quad4([2, 3, 6, 5])
|
||||
sel3 = Quad4([4, 5, 8, 7])
|
||||
sel4 = Quad4([5, 6, 9, 8])
|
||||
mel1 = Quad4([10, 11, 15, 14])
|
||||
mel2 = Quad4([11, 12, 16, 15])
|
||||
mel3 = Quad4([12, 13, 17, 16])
|
||||
mel4 = Quad4([14, 15, 19, 18])
|
||||
mel5 = Quad4([15, 16, 20, 19])
|
||||
mel6 = Quad4([16, 17, 21, 20])
|
||||
mel7 = Quad4([18, 19, 23, 22])
|
||||
mel8 = Quad4([19, 20, 24, 23])
|
||||
mel9 = Quad4([20, 21, 25, 24])
|
||||
update(Element[sel1, sel2, sel3, sel4, mel1, mel2, mel3,
|
||||
mel4, mel5, mel6, mel7, mel8, mel9], "geometry", nodes)
|
||||
calculate_normal_tangential_coordinates!(sel1, 0.0)
|
||||
calculate_normal_tangential_coordinates!(sel2, 0.0)
|
||||
calculate_normal_tangential_coordinates!(sel3, 0.0)
|
||||
calculate_normal_tangential_coordinates!(sel4, 0.0)
|
||||
prob = MortarProblem("temperature", 1)
|
||||
push!(prob, sel1)
|
||||
push!(prob, sel2)
|
||||
push!(prob, sel3)
|
||||
push!(prob, sel4)
|
||||
|
||||
master_elements = [mel1, mel2, mel3, mel4, mel5, mel6, mel7, mel8, mel9]
|
||||
sel1["master elements"] = master_elements
|
||||
sel2["master elements"] = master_elements
|
||||
sel3["master elements"] = master_elements
|
||||
sel4["master elements"] = master_elements
|
||||
B = sparse(assemble(prob, 0.0).stiffness_matrix, 25, 25)*46656
|
||||
B = full(B)
|
||||
D = B[1:9,1:9]
|
||||
M = B[1:9,10:end]
|
||||
info("interface matrix D:")
|
||||
dump(round(D, 3))
|
||||
info("interface matrix M:")
|
||||
dump(round(M, 3))
|
||||
D_expected = [
|
||||
1296 648 0 648 324 0 0 0 0
|
||||
648 2592 648 324 1296 324 0 0 0
|
||||
0 648 1296 0 324 648 0 0 0
|
||||
648 324 0 2592 1296 0 648 324 0
|
||||
324 1296 324 1296 5184 1296 324 1296 324
|
||||
0 324 648 0 1296 2592 0 324 648
|
||||
0 0 0 648 324 0 1296 648 0
|
||||
0 0 0 324 1296 324 648 2592 648
|
||||
0 0 0 0 324 648 0 648 1296]
|
||||
|
||||
M_expected = [
|
||||
-784 -700 -28 0 -700 -625 -25 0 -28 -25 -1 0 0 0 0 0
|
||||
-224 -1288 -1288 -224 -200 -1150 -1150 -200 -8 -46 -46 -8 0 0 0 0
|
||||
0 -28 -700 -784 0 -25 -625 -700 0 -1 -25 -28 0 0 0 0
|
||||
-224 -200 -8 0 -1288 -1150 -46 0 -1288 -1150 -46 0 -224 -200 -8 0
|
||||
-64 -368 -368 -64 -368 -2116 -2116 -368 -368 -2116 -2116 -368 -64 -368 -368 -64
|
||||
0 -8 -200 -224 0 -46 -1150 -1288 0 -46 -1150 -1288 0 -8 -200 -224
|
||||
0 0 0 0 -28 -25 -1 0 -700 -625 -25 0 -784 -700 -28 0
|
||||
0 0 0 0 -8 -46 -46 -8 -200 -1150 -1150 -200 -224 -1288 -1288 -224
|
||||
0 0 0 0 0 -1 -25 -28 0 -25 -625 -700 0 -28 -700 -784]
|
||||
|
||||
|
||||
info("D - D_expected")
|
||||
dump(D - D_expected)
|
||||
info("M - M_expected")
|
||||
dump(M - M_expected)
|
||||
|
||||
@test isapprox(D, D_expected)
|
||||
@test isapprox(M, M_expected)
|
||||
#=
|
||||
B = [
|
||||
864 432 0 432 216 0 -420 -375 -15 0 -504 -450 -18 0 -84 -75 -3 0
|
||||
432 1728 432 216 864 216 -120 -690 -690 -120 -144 -828 -828 -144 -24 -138 -138 -24
|
||||
0 432 864 0 216 432 0 -15 -375 -420 0 -18 -450 -504 0 -3 -75 -84
|
||||
432 216 0 864 432 0 -84 -75 -3 0 -504 -450 -18 0 -420 -375 -15 0
|
||||
216 864 216 432 1728 432 -24 -138 -138 -24 -144 -828 -828 -144 -120 -690 -690 -120
|
||||
0 216 432 0 432 864 0 -3 -75 -84 0 -18 -450 -504 0 -15 -375 -420
|
||||
]
|
||||
info("expected interface matrix:")
|
||||
dump(round(B, 3))
|
||||
@test isapprox(stiffness_matrix, B)
|
||||
=#
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@@ -7,7 +7,8 @@ using JuliaFEM
|
||||
using JuliaFEM.Test
|
||||
|
||||
#using JuliaFEM: parse
|
||||
using JuliaFEM.Preprocess: aster_parse_nodes
|
||||
using JuliaFEM.Preprocess: aster_parse_nodes, aster_renumber_nodes!,
|
||||
aster_renumber_elements!, aster_combine_meshes
|
||||
|
||||
function test_read_mesh()
|
||||
mesh = """
|
||||
@@ -69,6 +70,48 @@ function test_parse_nodes()
|
||||
@test nodes[8] == Float64[0.0, 1.0, 1.0]
|
||||
@test length(nodes) == 8
|
||||
end
|
||||
test_parse_nodes()
|
||||
#test_parse_nodes()
|
||||
|
||||
@testset "test combining meshes" begin
|
||||
|
||||
mesh1 = Dict{ASCIIString, Any}(
|
||||
"nodes" => Dict{Int64, Vector{Float64}}(
|
||||
1 => [1.0, 2.0],
|
||||
2 => [2.0, 3.0]
|
||||
),
|
||||
"connectivity" => Dict{Int64, Tuple{Symbol, Symbol, Vector{Int64}}}(
|
||||
1 => (:SE2, :GRP1, [1, 2])
|
||||
),
|
||||
)
|
||||
|
||||
mesh2 = Dict{ASCIIString, Any}(
|
||||
"nodes" => Dict{Int64, Vector{Float64}}(
|
||||
1 => [3.0, 4.0],
|
||||
2 => [4.0, 5.0]
|
||||
),
|
||||
"connectivity" => Dict{Int64, Tuple{Symbol, Symbol, Vector{Int64}}}(
|
||||
1 => (:SE2, :GRP1, [1, 2])
|
||||
),
|
||||
)
|
||||
|
||||
aster_renumber_nodes!(mesh1, mesh2)
|
||||
@test length(mesh2["nodes"]) == 2
|
||||
@test 3 in keys(mesh2["nodes"])
|
||||
@test 4 in keys(mesh2["nodes"])
|
||||
@test mesh2["connectivity"][1] == (:SE2, :GRP1, [3, 4])
|
||||
|
||||
aster_renumber_elements!(mesh1, mesh2)
|
||||
@test length(mesh2["connectivity"]) == 1
|
||||
@test mesh2["connectivity"][2] == (:SE2, :GRP1, [3, 4])
|
||||
|
||||
mesh = aster_combine_meshes(mesh1, mesh2)
|
||||
@test length(mesh["nodes"]) == 4
|
||||
@test mesh["nodes"][1] == [1.0, 2.0]
|
||||
@test mesh["nodes"][2] == [2.0, 3.0]
|
||||
@test mesh["nodes"][3] == [3.0, 4.0]
|
||||
@test mesh["nodes"][4] == [4.0, 5.0]
|
||||
@test mesh["connectivity"][1] == (:SE2, :GRP1, [1, 2])
|
||||
@test mesh["connectivity"][2] == (:SE2, :GRP1, [3, 4])
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user