mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-19 09:54:55 +00:00
mortar assembly for 3d problems
This commit is contained in:
Jukka Aho
File diff suppressed because one or more lines are too long
+27
-6
@@ -54,24 +54,46 @@ function get_integration_points(::Type{Seg3})
|
||||
return get_integration_points(Seg3, Val{3})
|
||||
end
|
||||
|
||||
### 2d elements
|
||||
### 2d triangular elements
|
||||
|
||||
function get_integration_points(::Type{Tri3})
|
||||
# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
|
||||
|
||||
typealias TriangularElements Union{Type{Tri3}, Type{Tri6}}
|
||||
|
||||
function get_integration_points(::TriangularElements, ::Type{Val{1}})
|
||||
# http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
|
||||
[
|
||||
IntegrationPoint([1.0/3.0, 1.0/3.0], 0.5)
|
||||
]
|
||||
end
|
||||
|
||||
function get_integration_points(::Type{Tri6})
|
||||
function get_integration_points(::TriangularElements, ::Type{Val{2}})
|
||||
# http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
|
||||
[
|
||||
IntegrationPoint([2.0/3.0, 1.0/6.0], 1.0/6.0)
|
||||
IntegrationPoint([1.0/6.0, 2.0/3.0], 1.0/6.0)
|
||||
IntegrationPoint([2.0/3.0, 1.0/6.0], 1.0/6.0),
|
||||
IntegrationPoint([1.0/6.0, 2.0/3.0], 1.0/6.0),
|
||||
IntegrationPoint([1.0/6.0, 1.0/6.0], 1.0/6.0)
|
||||
]
|
||||
end
|
||||
|
||||
function get_integration_points(::TriangularElements, ::Type{Val{5}})
|
||||
# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
|
||||
[
|
||||
IntegrationPoint([0.33333333333333, 0.33333333333333], 0.22500000000000),
|
||||
IntegrationPoint([0.47014206410511, 0.47014206410511], 0.13239415278851),
|
||||
IntegrationPoint([0.47014206410511, 0.05971587178977], 0.13239415278851),
|
||||
IntegrationPoint([0.05971587178977, 0.47014206410511], 0.13239415278851),
|
||||
IntegrationPoint([0.10128650732346, 0.10128650732346], 0.12593918054483),
|
||||
IntegrationPoint([0.10128650732346, 0.79742698535309], 0.12593918054483),
|
||||
IntegrationPoint([0.79742698535309, 0.10128650732346], 0.12593918054483)
|
||||
]
|
||||
end
|
||||
|
||||
function get_integration_points(::Type{Tri3})
|
||||
return get_integration_points(Tri3, Val{1})
|
||||
end
|
||||
|
||||
|
||||
function get_integration_points(::Type{Quad4})
|
||||
[
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
|
||||
@@ -102,4 +124,3 @@ function get_integration_points(::Type{Tet10})
|
||||
IntegrationPoint([b, b, b], w)
|
||||
]
|
||||
end
|
||||
|
||||
|
||||
+430
-3
@@ -124,6 +124,434 @@ function project_from_master_to_slave{S,M}(slave::Element{S}, master::Element{M}
|
||||
error("find projection from master to slave: did not converge")
|
||||
end
|
||||
|
||||
### Mortar projection calculation for 3d cases
|
||||
|
||||
"""
|
||||
Construct auxiliary plane for surface.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x::Array{Float64, 2}
|
||||
Node coordinates
|
||||
ximp::Array{Float64, 1}
|
||||
Element mid-point in dimensionless mother element coordinates ξ
|
||||
normals::Array{Float64, 2}
|
||||
Normal directions in nodes
|
||||
|
||||
Returns
|
||||
-------
|
||||
x0, Q
|
||||
x0::Array{Float64, 1} - origo of auxiliary plane
|
||||
Q::Array{Float64, 2} - orthogonal basis, first vector is normal direction
|
||||
and two rest vectors create orthonormal right-handed basis.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Calculate auxiliary plane given nodal coordinates, midpoint of mother element,
|
||||
node normals and suitable function space:
|
||||
|
||||
julia> xquad = [
|
||||
... -2.5 2.5 2.0 -2.0
|
||||
... -2.0 -2.0 2.3 2.0
|
||||
... 1.0 0.7 0.0 1.0]
|
||||
julia> m_midpoint = [0.0, 0.0]
|
||||
julia> normals = [
|
||||
... 0.05989060 0.0590504 0.225612 0.2445800
|
||||
... -0.00748633 0.1670810 0.182034 -0.0305725
|
||||
... 0.99817700 0.9841730 0.957059 0.9691470]
|
||||
julia> basis(xi) = [
|
||||
... (1-xi[1])(1-xi[2])/4
|
||||
... (1+xi[1])(1-xi[2])/4
|
||||
... (1+xi[1])(1+xi[2])/4
|
||||
... (1-xi[1])(1+xi[2])/4]'
|
||||
julia> x0, Q = create_auxiliary_plane(xquad, mmidpoint, normals, basis)
|
||||
julia> x0
|
||||
3-element Array{Float64,1}:
|
||||
0.0
|
||||
0.075
|
||||
0.675
|
||||
julia> Q
|
||||
3x3 Array{Float64,2}:
|
||||
0.148586 0.988899 0.0
|
||||
0.0784519 -0.0117877 0.996848
|
||||
0.985783 -0.148118 -0.0793325
|
||||
|
||||
Notes
|
||||
-----
|
||||
- Midpoint in mother element typically (0, 0) for quadrangles and (1/3, 1/3)
|
||||
for triangles.
|
||||
- Uses Gram-Schmidt process to find orthogonal basis
|
||||
- [1](http://www.math.umn.edu/~olver/aims_/qr.pdf)
|
||||
- [2](http://www.ecs.umass.edu/ece/ece313/Online_help/gram.pdf)
|
||||
- [3](http://www.terathon.com/code/tangent.html)
|
||||
|
||||
"""
|
||||
# function create_auxiliary_plane(x, ximp, normals, basis)
|
||||
function create_auxiliary_plane(element::Element{Tri3}, time::Real)
|
||||
proj(u, v) = dot(v, u) / dot(u, u) * u
|
||||
xi = [1.0/3.0, 1.0/3.0]
|
||||
x0 = element("geometry", xi, time)
|
||||
n = element("nodal ntsys", xi, time)[:, 1]
|
||||
n /= norm(n)
|
||||
# gram-schmidt
|
||||
u1 = n
|
||||
j = indmax(abs(u1))
|
||||
v2 = zeros(3)
|
||||
v2[mod(j,3)+1] = 1.0
|
||||
u2 = v2 - proj(u1, v2)
|
||||
u3 = cross(u1, u2)
|
||||
t1 = u2/norm(u2)
|
||||
t2 = u3/norm(u3)
|
||||
new_basis = [n t1 t2]
|
||||
return x0, new_basis
|
||||
end
|
||||
|
||||
"""
|
||||
Project point q onto a plane given by a point p and normal n.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
q::Array{Float64, 2}
|
||||
point to project (row vector)
|
||||
x0::Array{Float64, 2}
|
||||
origo of plane
|
||||
n::Array{Float64, 2}
|
||||
normal vector of plane
|
||||
|
||||
Returns
|
||||
-------
|
||||
y::Array{Float64, 2}
|
||||
projected point
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> p = [-0.5 -1.0 4.0]'
|
||||
julia> x0 = [0.0 0.075 0.675]'
|
||||
julia> n = [0.1485860 0.0784519 0.9857830]'
|
||||
julia> project_node_to_auxiliary_plane(p, x0, n)
|
||||
3-element Array{Float64,1}:
|
||||
0.963455
|
||||
-1.2447
|
||||
0.925247
|
||||
|
||||
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)
|
||||
n = Q[:,1]
|
||||
ph = p - dot(p-x0, n)*n
|
||||
qproj = Q'*(ph-x0)
|
||||
@assert isapprox(qproj[1], 0.0)
|
||||
return qproj[2:3]
|
||||
end
|
||||
|
||||
"""
|
||||
Find edge intersections of two planar arbitrary shape polygons.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
S::Array{Float64,2}
|
||||
M::Array{Float64,2}
|
||||
|
||||
Matrices with size (2, n) where n is number of vertices of each polygon.
|
||||
|
||||
Returns
|
||||
-------
|
||||
P::Array{Float64,2}
|
||||
Intersection points of polygons
|
||||
n::Array{Float64,2}
|
||||
Neighbour info matrix with size (ns, mn). This keeps information which
|
||||
edges of polygons are intersecting. See further explanation in example
|
||||
below.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Find intersection points of two triangles:
|
||||
|
||||
julia> S = [0 0; 3 0; 0 3]'
|
||||
julia> M = [-1 1; 2 -1/2; 1 3/2]'
|
||||
julia> P, n = get_edge_intersections(S, M)
|
||||
julia> P
|
||||
2x4 Array{Float64,2}:
|
||||
1.0 1.75 0.0 0.0
|
||||
0.0 0.0 0.5 1.25
|
||||
julia> n
|
||||
3x3 Array{Int64,2}:
|
||||
1 1 0
|
||||
0 0 0
|
||||
1 0 1)
|
||||
|
||||
So intersection points are: (1.00, 0.00), (1.75, 0.00), (0.00, 0.50), (0.00, 1.25).
|
||||
"Neighbour matrix" can be interpreted as following:
|
||||
|
||||
1 1 0 <--> First edge of S intersects edges 1 and 2 of M
|
||||
0 0 0 <--> Second edge of S doesn't intersect at all
|
||||
1 0 1 <--> Third edge of S intersects with edges 1 and 3 of M
|
||||
|
||||
"""
|
||||
function get_edge_intersections(S::Matrix, M::Matrix)
|
||||
ns = size(S, 2)
|
||||
nm = size(M, 2)
|
||||
P = zeros(2, 0)
|
||||
n = zeros(Int64, ns, nm)
|
||||
k = 0
|
||||
for i=1:ns
|
||||
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
|
||||
r = A\b
|
||||
if (r[1]>=0) & (r[1]<=1) & (r[2]>=0) & (r[2]<=1) # intersection found
|
||||
k += 1
|
||||
f = S[:,i]+r[1]*(S[:,mod(i,ns)+1] - S[:,i])
|
||||
f = f''
|
||||
P = hcat(P, f)
|
||||
n[i, j] = 1
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
return P, n
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Find any points laying inside or border of triangle.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
Y::Array{Float64, 2}
|
||||
Triangle coordinates in 2×3 matrix
|
||||
X::Array{Float64, 2}
|
||||
List of points to test in 2×n matrix
|
||||
|
||||
Returns
|
||||
-------
|
||||
P::Array{Float64, 2}
|
||||
List of points in triangle in 2×m matrix, where m is number of points inside triangle
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> S = [0.0 0.0; 3.0 0.0; 0.0 3.0]' # triangle corner points
|
||||
julia> pts = [-1.0 1.0; 2.0 -0.5; 1.0 1.5; 0.5 1.5]' # points to tests
|
||||
julia> points_in_triangle(S, pts)
|
||||
2x2 Array{Float64,2}:
|
||||
1.0 0.5
|
||||
1.5 1.5
|
||||
"""
|
||||
function get_points_inside_triangle(Y::Matrix, X::Matrix)
|
||||
@assert size(Y, 2) == 3 # "Point in TRIANGLE..."
|
||||
P = zeros(2, 0)
|
||||
v0 = Y[:,2] - Y[:,1]
|
||||
v1 = Y[:,3] - Y[:,1] # find interior points of X in Y
|
||||
d00 = (v0'*v0)[1]
|
||||
d01 = (v0'*v1)[1]
|
||||
d11 = (v1'*v1)[1] # using baricentric coordinates
|
||||
id = 1/(d00*d11 - d01*d01)
|
||||
for i=1:size(X, 2)
|
||||
v2 = X[:,i] - Y[:,1]
|
||||
d02 = (v0'*v2)[1]
|
||||
d12 = (v1'*v2)[1]
|
||||
u = (d11*d02-d01*d12)*id
|
||||
v = (d00*d12-d01*d02)*id
|
||||
if (u>=0) & (v>=0) & (u+v<=1) # also include nodes on the boundary
|
||||
P = hcat(P, X[:,i]'')
|
||||
end
|
||||
end
|
||||
return P
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Make polygon clipping of shapes S and M.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
S::Array{Float64, 2}
|
||||
M::Array{Float64, 2}
|
||||
Shapes to clip. Needs to be triangles at the moment.
|
||||
|
||||
Returns
|
||||
-------
|
||||
Array{Float64, 2}, Array{Float64, 2}
|
||||
- Polygon vertices in 2×n matrix, sorted in counter-clockwise order.
|
||||
- 3×3 "neighbouring" matrix, see example.
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> S = [0 0; 3 0; 0 3]'
|
||||
julia> M = [-1 1; 2 -1/2; 2 2]'
|
||||
julia> P, n = clip_polygon(S, M)
|
||||
julia> P
|
||||
2x6 Array{Float64,2}:
|
||||
0.0 1.0 2.0 2.0 1.25 0.0
|
||||
0.5 0.0 0.0 1.0 1.75 1.33333,
|
||||
julia> n
|
||||
3x3 Array{Int64,2}:
|
||||
1 0 1 <- first edge of M ([-1 1; 2 -1/2]') intersects with edges 1 and 3 of S ([0 0; 3 0]' and [0 3; 0 0]')
|
||||
1 1 0 <- second edge of M ([2 -1/2; 2 2]') intersects with edges 1 and 2 of S
|
||||
0 1 1 <- third edge of M ([2 2; -1 1]') intersects with edgse 2 and 3 of S
|
||||
|
||||
"""
|
||||
function clip_polygon(S::Matrix, M::Matrix)
|
||||
P1, neighbours = get_edge_intersections(M, S)
|
||||
P2 = get_points_inside_triangle(M, S)
|
||||
P3 = get_points_inside_triangle(S, M)
|
||||
P = hcat(P1, P2, P3)
|
||||
meanval = mean(P, 2)
|
||||
tmp = P .- meanval
|
||||
angles = atan2(tmp[2,:], tmp[1,:])
|
||||
angles = reshape(angles, length(angles))
|
||||
order = sortperm(angles)
|
||||
P = copy(unique(P[:, order], 2))
|
||||
return P, neighbours
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Calculate polygon geometric center point
|
||||
|
||||
Parameters
|
||||
----------
|
||||
P::Array{Float64, 2}
|
||||
Polygon vertices in 2×n matrix
|
||||
|
||||
Returns
|
||||
-------
|
||||
Array{Float63, 2}
|
||||
Center point
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> P
|
||||
2x6 Array{Float64,2}:
|
||||
0.0 1.0 2.0 2.0 1.25 0.0
|
||||
0.5 0.0 0.0 1.0 1.75 1.33333,
|
||||
julia> C = get_polygon_cp(P)
|
||||
2x1 Array{Float64,2}:
|
||||
1.039740
|
||||
0.804701
|
||||
|
||||
"""
|
||||
function calculate_polygon_centerpoint(P::Matrix)
|
||||
n = size(P, 2)
|
||||
A = 0.0
|
||||
for i=1:n
|
||||
A += 1/2*(P[1,i]*P[2,mod(i,n)+1] - P[1,mod(i,n)+1]*P[2,i])
|
||||
end
|
||||
Cx = 0.0
|
||||
Cy = 0.0
|
||||
for i=1:n
|
||||
inext = mod(i, n)+1
|
||||
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]
|
||||
end
|
||||
|
||||
"""
|
||||
Project point from auxiliary plane to parametric surface given by (ξ₁, ξ₂)
|
||||
|
||||
Parameters
|
||||
----------
|
||||
p::Array{Float64,1}
|
||||
point in auxiliary plane, in (n,t1,t2) coordinate system
|
||||
x0::Array{Float64,1}
|
||||
origo of auxiliary plane cs
|
||||
Q::Array{Float64,2}
|
||||
basis of auxiliary plane cs
|
||||
x::Array{Float64,2}
|
||||
surface node coords
|
||||
basis::Array{Float64,2}
|
||||
surface basis functions
|
||||
dbasis::Array{Float64,2}
|
||||
partial derivatives of surface basis functions
|
||||
|
||||
Returns
|
||||
-------
|
||||
Array{Float64,2}
|
||||
solution vector (d, ξ₁, ξ₂) where d is distance to surface
|
||||
|
||||
Examples
|
||||
--------
|
||||
Define surface with node points, basis + dbasis
|
||||
|
||||
julia> xquad = [
|
||||
... -2.5 -2.0 1.0
|
||||
... 2.5 -2.0 0.7
|
||||
... 2.0 2.3 0.0
|
||||
... -2.0 2.0 1.0]'
|
||||
julia> basis(xi) = [
|
||||
... (1-xi[1])(1-xi[2])/4
|
||||
... (1+xi[1])(1-xi[2])/4
|
||||
... (1+xi[1])(1+xi[2])/4
|
||||
... (1-xi[1])(1+xi[2])/4]
|
||||
julia> dbasis(xi) = [
|
||||
... -(1-xi[2])/4 -(1-xi[1])/4
|
||||
... (1-xi[2])/4 -(1+xi[1])/4
|
||||
... (1+xi[2])/4 (1+xi[1])/4
|
||||
... -(1+xi[2])/4 (1-xi[1])/4]
|
||||
|
||||
We aim to find point p, which we first project to auxiliary plane defined as following
|
||||
julia> p = [-2.5 -2.0 1.0]'
|
||||
julia> x0 = [0.0 0.075 0.675]'
|
||||
julia> Q = [
|
||||
... 0.1485860 0.9888990 0.0000000
|
||||
... 0.0784519 -0.0117877 0.9968480
|
||||
... 0.9857830 -0.1481180 -0.0793325]
|
||||
|
||||
Our projected point is therefore
|
||||
julia> n = Q[:,1] # first component is normal direction
|
||||
julia> ph = project_node_to_auxiliary_plane(p, x0, n)
|
||||
julia> ph = Q'(ph-x0)
|
||||
julia> ph
|
||||
3x1 Array{Float64,2}:
|
||||
1.33264e-7
|
||||
-2.49593
|
||||
-2.09424
|
||||
|
||||
Our point ph is now in auxiliary plane in n,t1,t2 coordinate system. Next we
|
||||
project it back to surface defined by xquad*basis
|
||||
|
||||
julia> theta = project_point_from_plane_to_surface(ph, x0, Q, xquad, basis, dbasis)
|
||||
julia> theta
|
||||
3x1 Array{Float64,2}:
|
||||
-0.213874
|
||||
-0.999999
|
||||
-1.0
|
||||
|
||||
We see that our ξ₁ = ξ₂ = -1 so we found first point of xquad
|
||||
[-2.5 -2.0 1.0]' correctly.
|
||||
|
||||
julia> xquad*basis(theta[2:3])
|
||||
3-element Array{Float64,1}:
|
||||
-2.5
|
||||
-2.0
|
||||
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)
|
||||
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]
|
||||
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
|
||||
theta += dtheta
|
||||
if norm(dtheta) < iter_tol
|
||||
return theta
|
||||
end
|
||||
end
|
||||
error("project_point_to_auxiliary_plane: did not converge in $max_iterations iterations!")
|
||||
end
|
||||
|
||||
|
||||
### Mortar problem
|
||||
|
||||
"""
|
||||
@@ -142,7 +570,7 @@ end
|
||||
# Mortar assembly
|
||||
|
||||
function assemble!(assembly::Assembly, problem::BoundaryProblem{MortarProblem}, slave_element::Element, time::Number)
|
||||
|
||||
|
||||
# get dimension and name of PARENT field
|
||||
field_dim = problem.parent_field_dim
|
||||
field_name = problem.parent_field_name
|
||||
@@ -167,7 +595,7 @@ function assemble!(assembly::Assembly, problem::BoundaryProblem{MortarProblem},
|
||||
# projected integration point
|
||||
xi_projected = project_from_slave_to_master(slave_element, master_element, xi_gauss)
|
||||
|
||||
# add contribution to left hand side
|
||||
# add contribution to left hand side
|
||||
N1 = slave_element(xi_gauss, time)
|
||||
N2 = master_element(xi_projected, time)
|
||||
S = w*N1'*N1
|
||||
@@ -184,4 +612,3 @@ function assemble!(assembly::Assembly, problem::BoundaryProblem{MortarProblem},
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
+121
-4
@@ -5,10 +5,18 @@ module MortarTests
|
||||
|
||||
using JuliaFEM.Test
|
||||
|
||||
using JuliaFEM.Core: Element, Seg2, Quad4, MortarProblem, Assembly, assemble!
|
||||
using JuliaFEM.Core: Element, Seg2, Quad4, Tri3, MortarProblem, Assembly, assemble!
|
||||
using JuliaFEM.Core: PlaneStressElasticityProblem, DirichletProblem, DirectSolver
|
||||
|
||||
# 2d stuff
|
||||
using JuliaFEM.Core: project_from_slave_to_master, project_from_master_to_slave
|
||||
|
||||
# 3d stuff
|
||||
using JuliaFEM.Core: create_auxiliary_plane, project_point_to_auxiliary_plane,
|
||||
get_edge_intersections, get_points_inside_triangle,
|
||||
clip_polygon, calculate_polygon_centerpoint,
|
||||
project_point_from_plane_to_surface
|
||||
|
||||
function get_test_2d_model()
|
||||
# this is hand calculated and given as an example in my thesis
|
||||
N = Vector[
|
||||
@@ -118,7 +126,7 @@ function test_create_flat_2d_assembly()
|
||||
info("size of B = $(size(B))")
|
||||
info("B matrix in first slave element = \n$(B[10:11,:])")
|
||||
info("B matrix expected = \n$(B_expected[10:11,:])")
|
||||
@test isapprox(B, B_expected)
|
||||
@test isapprox(B, B_expected)
|
||||
|
||||
fill!(B_expected, 0.0)
|
||||
empty!(assembly)
|
||||
@@ -184,7 +192,7 @@ function test_2d_mortar_multiple_bodies_multiple_dirichlet_bc()
|
||||
dy1 = Seg2([1, 2])
|
||||
dy1["geometry"] = Vector[N[1], N[2]]
|
||||
dy1["displacement 2"] = 0.0
|
||||
|
||||
|
||||
boundary2 = DirichletProblem("displacement", 2)
|
||||
push!(boundary2, dy1)
|
||||
|
||||
@@ -282,7 +290,7 @@ function test_2d_mortar_three_bodies_shared_nodes()
|
||||
dy1 = Seg2([1, 2])
|
||||
dy1["geometry"] = Vector[N[1], N[2]]
|
||||
dy1["displacement 2"] = 0.0
|
||||
|
||||
|
||||
bc2 = DirichletProblem("displacement", 2)
|
||||
push!(bc2, dy1)
|
||||
|
||||
@@ -343,4 +351,113 @@ function test_2d_mortar_three_bodies_shared_nodes()
|
||||
end
|
||||
#test_2d_mortar_three_bodies_shared_nodes()
|
||||
|
||||
function test_auxiliary_plane_transforms()
|
||||
nodes = Vector{Float64}[
|
||||
[0.0, 0.0, 0.0],
|
||||
[1.0, 0.0, 0.0],
|
||||
[0.0, 1.0, 0.0]]
|
||||
e1 = Tri3([1, 2, 3])
|
||||
# local coordinate system N, T1, T2 in node
|
||||
R = [0.0 1.0 0.0
|
||||
0.0 0.0 1.0
|
||||
1.0 0.0 0.0]
|
||||
e1["geometry"] = Vector{Float64}[nodes[1], nodes[2], nodes[3]]
|
||||
e1["nodal ntsys"] = Matrix{Float64}[R, R, R]
|
||||
time::Real = 0.0
|
||||
x0, Q = create_auxiliary_plane(e1, time)
|
||||
info("x0 = $x0")
|
||||
info("Q = $Q")
|
||||
@test isapprox(x0, [1.0/3.0, 1.0/3.0, 0.0])
|
||||
@test isapprox(Q, R)
|
||||
p1 = Float64[1.0/3.0+0.1, 1.0/3.0+0.1, 1.0]
|
||||
p2 = project_point_to_auxiliary_plane(p1, x0, Q)
|
||||
info("point in auxiliary plane p2 = $p2")
|
||||
@test isapprox(p2, [0.1, 0.1])
|
||||
theta = project_point_from_plane_to_surface(p2, x0, Q, e1, time)
|
||||
info("theta = $theta")
|
||||
@test isapprox(theta[1], 0.0)
|
||||
X = e1("geometry", theta[2:3], time)
|
||||
info("projected point = $X")
|
||||
@test isapprox(X, Float64[1.0/3.0+0.1, 1.0/3.0+0.1, 0.0])
|
||||
end
|
||||
test_auxiliary_plane_transforms()
|
||||
|
||||
|
||||
function test_get_edge_intersections()
|
||||
# first case, two triangles
|
||||
S = [ 0.0 0.0; 3.0 0.0; 0.0 3.0]'
|
||||
M = [-1.0 1.0; 2.0 -0.5; 1.0 1.5]'
|
||||
P, n = get_edge_intersections(S, M)
|
||||
P_expected = [
|
||||
1.00 1.75 0.00 0.00
|
||||
0.00 0.00 0.50 1.25]
|
||||
n_expected = [
|
||||
1 1 0
|
||||
0 0 0
|
||||
1 0 1]
|
||||
@test isapprox(P, P_expected)
|
||||
@test isapprox(n, n_expected)
|
||||
|
||||
# slave 4 vertices non-convex, master triangle
|
||||
S = [ 0.0 0.0; 2.5 0.0; 1.0 1.0; 0.0 2.0]'
|
||||
M = [-1.0 1.0; 2.0 -0.5; 1.0 1.5]'
|
||||
P, n = get_edge_intersections(S, M)
|
||||
P_expected = [
|
||||
1.0 1.75 1.375 0.60 0.00 0.00
|
||||
0.0 0.00 0.750 1.40 0.50 1.25]
|
||||
n_expected = [
|
||||
1 1 0
|
||||
0 1 0
|
||||
0 0 1
|
||||
1 0 1]
|
||||
@test isapprox(P, P_expected)
|
||||
@test isapprox(n, n_expected)
|
||||
|
||||
# slave 3 triangle, master 4 vertices
|
||||
S = [ 0.0 0.0; 3.0 0.0; 0.0 3.0]'
|
||||
M = [-1.0 1.0; 2.0 -0.5; 1.0 1.5; -1.0 2.0]'
|
||||
P, n = get_edge_intersections(S, M)
|
||||
P_expected = [
|
||||
1.00 1.75 0.00 0.00
|
||||
0.00 0.00 0.50 1.75]
|
||||
n_expected = [
|
||||
1 1 0 0
|
||||
0 0 0 0
|
||||
1 0 1 0]
|
||||
@test isapprox(P, P_expected)
|
||||
@test isapprox(n, n_expected)
|
||||
end
|
||||
#test_get_edge_intersections()
|
||||
|
||||
|
||||
function test_get_points_inside_triangle()
|
||||
S = [0.0 0.0; 3.0 0.0; 0.0 3.0]'
|
||||
pts = [-1.0 1.0; 2.0 -0.5; 1.0 1.5; 0.5 1.5]'
|
||||
P = get_points_inside_triangle(S, pts)
|
||||
@test isapprox(P, [1.0 1.5; 0.5 1.5]')
|
||||
end
|
||||
#test_get_points_inside_triangle()
|
||||
|
||||
|
||||
function test_polygon_clipping()
|
||||
S = [0 0; 3 0; 0 3]'
|
||||
M = [-1 1; 2 -1/2; 2 2]'
|
||||
P, n = clip_polygon(S, M)
|
||||
@test isapprox(P, [0.0 0.5; 1.0 0.0; 2.0 0.0; 2.0 1.0; 1.25 1.75; 0.0 4/3]')
|
||||
@test isapprox(n, [1 0 1; 1 1 0; 0 1 1])
|
||||
end
|
||||
#test_polygon_clipping()
|
||||
|
||||
|
||||
function test_calculate_polygon_centerpoint()
|
||||
P = [
|
||||
0.0 1.0 2.0 2.0 1.25 0.0
|
||||
0.5 0.0 0.0 1.0 1.75 1.33333]
|
||||
C = calculate_polygon_centerpoint(P)
|
||||
info("Polygon centerpoint: $C")
|
||||
@test isapprox(C, [1.0397440690338993, 0.8047003412233396])
|
||||
end
|
||||
#test_calculate_polygon_centerpoint()
|
||||
|
||||
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user