mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-20 18:18:31 +00:00
2d tie contact working.
This commit is contained in:
Jukka Aho
@@ -0,0 +1,401 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# 2d tie contact\n",
|
||||
"\n",
|
||||
"Author: Jukka Aho\n",
|
||||
"\n",
|
||||
"Abstract: 2d tie contact.\n",
|
||||
"\n",
|
||||
"Model:\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"Each element is modelled as own \"body\" and they are connected using tie contacts. Segments 5-6 and 9-10 and 6-7 are slave surfaces, so node 6 or 9 is on at least two tie contacts as slave node. Moreover this model has dirichlet boundary $y=0$ at bottom of body 1 and $x=0$ on left. To get the accurate solution one needs to minimize \n",
|
||||
"\\begin{equation}\n",
|
||||
"\\frac{15}{2}u_{1}^{4} + 60 u_{1}^{3} + \\frac{15}{4}u_{1}^{2} u_{2}^{2} + 15 u_{1}^{2} u_{2} + 120 u_{1}^{2} + 15 u_{1} u_{2}^{2} + 60 u_{1} u_{2} + \\frac{15}{2}u_{2}^{4} + 60 u_{2}^{3} + 120 u_{2}^{2} + 50 u_{2}\n",
|
||||
",\n",
|
||||
"\\end{equation}\n",
|
||||
"which gives approximate $u_1 = 0.0634862$ and $u_2 = -0.277183$ for the displacement of upper right corner.\n",
|
||||
"[Wolfram](http://www.wolframalpha.com/input/?i=local+minimum+15*x^4%2F2+%2B+60*x^3+%2B+15*x^2*y^2%2F4+%2B+15*x^2*y+%2B+120*x^2+%2B+15*x*y^2+%2B+60*x*y+%2B+15*y^4%2F2+%2B+60*y^3+%2B+120*y^2+%2B+50*y)."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using JuliaFEM\n",
|
||||
"using JuliaFEM: Element, Seg2, Quad4, PlaneStressElasticityProblem, DirichletProblem, MortarProblem, DirectSolver"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"nodes = Dict{Int64, Vector{Float64}}(\n",
|
||||
" 1 => [0.0, 0.0],\n",
|
||||
" 2 => [2.0, 0.0],\n",
|
||||
" 3 => [2.0, 1.0],\n",
|
||||
" 4 => [0.0, 1.0],\n",
|
||||
" 5 => [0.0, 1.0],\n",
|
||||
" 6 => [1.0, 1.0],\n",
|
||||
" 7 => [1.0, 2.0],\n",
|
||||
" 8 => [0.0, 2.0],\n",
|
||||
" 9 => [1.0, 1.0],\n",
|
||||
" 10 => [2.0, 1.0],\n",
|
||||
" 11 => [2.0, 2.0],\n",
|
||||
" 12 => [1.0, 2.0]);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"connectivity = Dict{Int64, Vector{Int64}}(\n",
|
||||
" 1 => [1, 2, 3, 4],\n",
|
||||
" 2 => [5, 6, 7, 8],\n",
|
||||
" 3 => [9, 10, 11, 12]);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"3"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"elements = Element[]\n",
|
||||
"for c in values(connectivity)\n",
|
||||
" element = Quad4(c)\n",
|
||||
" element[\"geometry\"] = Vector{Float64}[nodes[i] for i in c]\n",
|
||||
" element[\"youngs modulus\"] = 900.0\n",
|
||||
" element[\"poissons ratio\"] = 0.25\n",
|
||||
" push!(elements, element)\n",
|
||||
"end\n",
|
||||
"length(elements)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Create three bodies, each containing one element."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"body1 = PlaneStressElasticityProblem()\n",
|
||||
"body2 = PlaneStressElasticityProblem()\n",
|
||||
"body3 = PlaneStressElasticityProblem()\n",
|
||||
"push!(body1, elements[1])\n",
|
||||
"push!(body2, elements[2])\n",
|
||||
"push!(body3, elements[3]);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Surface traction to the top of bodies 2 and 3:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"t2 = Seg2([8, 7])\n",
|
||||
"t2[\"geometry\"] = Vector{Float64}[nodes[8], nodes[7]]\n",
|
||||
"t2[\"displacement traction force\"] = Vector{Float64}[[0.0, -100.0], [0.0, -100.0]]\n",
|
||||
"t3 = Seg2([12, 11])\n",
|
||||
"t3[\"geometry\"] = Vector{Float64}[nodes[12], nodes[11]]\n",
|
||||
"t3[\"displacement traction force\"] = Vector{Float64}[[0.0, -100.0], [0.0, -100.0]]\n",
|
||||
"push!(body2, t2)\n",
|
||||
"push!(body3, t3);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Boundary conditions: $x=0$ for left boundary."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"dx1 = Seg2([1, 4])\n",
|
||||
"dx1[\"geometry\"] = Vector[nodes[1], nodes[4]]\n",
|
||||
"dx1[\"displacement 1\"] = 0.0\n",
|
||||
"dx2 = Seg2([5, 8])\n",
|
||||
"dx2[\"geometry\"] = Vector[nodes[5], nodes[8]]\n",
|
||||
"dx2[\"displacement 1\"] = 0.0\n",
|
||||
"bc1 = DirichletProblem(\"displacement\", 2)\n",
|
||||
"push!(bc1, dx1)\n",
|
||||
"push!(bc1, dx2);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"$y=0$ for bottom of model"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"dy1 = Seg2([1, 2])\n",
|
||||
"dy1[\"geometry\"] = Vector[nodes[1], nodes[2]]\n",
|
||||
"dy1[\"displacement 2\"] = 0.0\n",
|
||||
"bc2 = DirichletProblem(\"displacement\", 2)\n",
|
||||
"push!(bc2, dy1);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Mortar boundary conditions: tie contact between body 1 and body 2"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)]\n",
|
||||
"\n",
|
||||
"master1 = Seg2([4, 3])\n",
|
||||
"master1[\"geometry\"] = Vector[nodes[4], nodes[3]]\n",
|
||||
"slave1 = Seg2([5, 6])\n",
|
||||
"slave1[\"geometry\"] = Vector[nodes[5], nodes[6]]\n",
|
||||
"slave1[\"master elements\"] = Element[master1]\n",
|
||||
"slave1[\"nodal ntsys\"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]\n",
|
||||
"contact1 = MortarProblem(\"displacement\", 2)\n",
|
||||
"push!(contact1, slave1);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Tie contact between body 1 and body 3"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"slave2 = Seg2([9, 10])\n",
|
||||
"slave2[\"geometry\"] = Vector[nodes[9], nodes[10]]\n",
|
||||
"slave2[\"nodal ntsys\"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]\n",
|
||||
"slave2[\"master elements\"] = Element[master1]\n",
|
||||
"contact2 = MortarProblem(\"displacement\", 2)\n",
|
||||
"push!(contact2, slave2);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Tie contact between body 2 and body 3"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"master2 = Seg2([6, 7])\n",
|
||||
"master2[\"geometry\"] = Vector[nodes[6], nodes[7]]\n",
|
||||
"slave3 = Seg2([9, 12])\n",
|
||||
"slave3[\"geometry\"] = Vector[nodes[9], nodes[12]]\n",
|
||||
"slave3[\"nodal ntsys\"] = Matrix[rotation_matrix(0.0), rotation_matrix(0.0)]\n",
|
||||
"slave3[\"master elements\"] = Element[master2]\n",
|
||||
"contact3 = MortarProblem(\"displacement\", 2)\n",
|
||||
"push!(contact3, slave3);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"All defined. Solve it."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"solver = DirectSolver()\n",
|
||||
"push!(solver, body1)\n",
|
||||
"push!(solver, body2)\n",
|
||||
"push!(solver, body3)\n",
|
||||
"push!(solver, bc1)\n",
|
||||
"push!(solver, bc2)\n",
|
||||
"push!(solver, contact1)\n",
|
||||
"push!(solver, contact2)\n",
|
||||
"push!(solver, contact3);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"INFO: # of field problems: 3\n",
|
||||
"INFO: # of boundary problems: 5\n",
|
||||
"INFO: Starting iteration 1\n",
|
||||
"INFO: # of dofs: 24, # of interface dofs: 15\n",
|
||||
"INFO: solved. length of solution vector = 48\n",
|
||||
"INFO: Iteration took 9.311098465 seconds\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(5,true)"
|
||||
]
|
||||
},
|
||||
"execution_count": 13,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"INFO: Starting iteration 2\n",
|
||||
"INFO: # of dofs: 24, # of interface dofs: 15\n",
|
||||
"INFO: solved. length of solution vector = 48\n",
|
||||
"INFO: Iteration took 0.003787437 seconds\n",
|
||||
"INFO: Starting iteration 3\n",
|
||||
"INFO: # of dofs: 24, # of interface dofs: 15\n",
|
||||
"INFO: solved. length of solution vector = 48\n",
|
||||
"INFO: Iteration took 0.020931551 seconds\n",
|
||||
"INFO: Starting iteration 4\n",
|
||||
"INFO: # of dofs: 24, # of interface dofs: 15\n",
|
||||
"INFO: solved. length of solution vector = 48\n",
|
||||
"INFO: Iteration took 0.003763852 seconds\n",
|
||||
"INFO: Starting iteration 5\n",
|
||||
"INFO: # of dofs: 24, # of interface dofs: 15\n",
|
||||
"INFO: solved. length of solution vector = 48\n",
|
||||
"INFO: Iteration took 0.003679408 seconds\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"iterations, converged = call(solver, 0.0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"INFO: displacement at [2.0,2.0] = [0.06348623177789343,-0.27718303785565257]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using JuliaFEM.Test\n",
|
||||
"\n",
|
||||
"@test converged\n",
|
||||
"\n",
|
||||
"X = elements[2](\"geometry\", [1.0, 1.0], 0.0)\n",
|
||||
"u = elements[2](\"displacement\", [1.0, 1.0], 0.0)\n",
|
||||
"info(\"displacement at $X = $u\")\n",
|
||||
"@test isapprox(u, [0.0634862, -0.277183], atol=1.0e-5)"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Julia 0.4.0",
|
||||
"language": "julia",
|
||||
"name": "julia-0.4"
|
||||
},
|
||||
"language_info": {
|
||||
"file_extension": ".jl",
|
||||
"mimetype": "application/julia",
|
||||
"name": "julia",
|
||||
"version": "0.4.1"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 0
|
||||
}
|
||||
+19
-2
@@ -3,9 +3,26 @@
|
||||
|
||||
# Functions to handle global assembly of problem
|
||||
|
||||
function assemble!(assembly::Assembly, problem::Problem, time::Number=0.0)
|
||||
empty!(assembly)
|
||||
function assemble!(assembly::Assembly, problem::Problem, time::Number=0.0, empty_assembly::Bool=true)
|
||||
if empty_assembly
|
||||
empty!(assembly)
|
||||
end
|
||||
for equation in get_equations(problem)
|
||||
assemble!(assembly, equation, time, problem)
|
||||
end
|
||||
end
|
||||
|
||||
function assemble(problem::Problem, time::Number=0.0)
|
||||
assembly = Assembly()
|
||||
for equation in get_equations(problem)
|
||||
assemble!(assembly, equation, time, problem)
|
||||
end
|
||||
return assembly
|
||||
end
|
||||
|
||||
function Base.(:+)(ass1::Assembly, ass2::Assembly)
|
||||
mass_matrix = ass1.mass_matrix + ass2.mass_matrix
|
||||
stiffness_matrix = ass1.stiffness_matrix + ass2.stiffness_matrix
|
||||
force_vector = ass1.force_vector + ass2.force_vector
|
||||
return Assembly(mass_matrix, stiffness_matrix, force_vector)
|
||||
end
|
||||
|
||||
+91
-60
@@ -6,6 +6,7 @@
|
||||
type DirectSolver <: Solver
|
||||
field_problems :: Vector{FieldProblem}
|
||||
boundary_problems :: Vector{BoundaryProblem}
|
||||
parallel :: Bool
|
||||
nonlinear_problem :: Bool
|
||||
max_iterations :: Int64
|
||||
tol :: Float64
|
||||
@@ -21,90 +22,120 @@ end
|
||||
|
||||
""" Default initializer. """
|
||||
function DirectSolver()
|
||||
DirectSolver([], [], true, 10, 1.0e-6)
|
||||
DirectSolver([], [], false, true, 10, 1.0e-6)
|
||||
end
|
||||
|
||||
""" Call solver to solve a set of problems. """
|
||||
function call(solver::DirectSolver, time::Number=0.0)
|
||||
@assert length(solver.field_problems) == 1
|
||||
@assert length(solver.boundary_problems) == 1
|
||||
#@assert length(solver.field_problems) == 1
|
||||
info("# of field problems: $(length(solver.field_problems))")
|
||||
info("# of boundary problems: $(length(solver.boundary_problems))")
|
||||
@assert solver.nonlinear_problem == true
|
||||
|
||||
problem1 = solver.field_problems[1]
|
||||
problem2 = solver.boundary_problems[1]
|
||||
# check that all problems are "same kind"
|
||||
field_name = get_unknown_field_name(solver.field_problems[1])
|
||||
field_dim = get_unknown_field_dimension(solver.field_problems[1])
|
||||
for field_problem in solver.field_problems
|
||||
get_unknown_field_name(field_problem) == field_name || error("several different fields not supported yet")
|
||||
get_unknown_field_dimension(field_problem) == field_dim || error("several different field dimensions not supported yet")
|
||||
end
|
||||
|
||||
x = zeros(3)
|
||||
dx = zeros(3)
|
||||
dims = nothing
|
||||
# create initial fields for this increment
|
||||
# i.e., copy last known values as initial guess
|
||||
# for this increment
|
||||
|
||||
for field_problem in solver.field_problems
|
||||
for equation in get_equations(field_problem)
|
||||
element = get_element(equation)
|
||||
gdofs = get_gdofs(field_problem, equation)
|
||||
if !isapprox(last(element[field_name]).time, time)
|
||||
last_data = copy(last(element[field_name]).data)
|
||||
push!(element[field_name], time => last_data)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
for boundary_problem in solver.boundary_problems
|
||||
for equation in get_equations(boundary_problem)
|
||||
element = get_element(equation)
|
||||
gdofs = get_gdofs(boundary_problem, equation)
|
||||
eqdim = size(equation)[2]
|
||||
data = Vector{Float64}[zeros(field_dim) for i in 1:eqdim]
|
||||
if !isapprox(last(element["reaction force"]).time, time)
|
||||
push!(element["reaction force"], time => data)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
dim = 0
|
||||
|
||||
for iter=1:solver.max_iterations
|
||||
tic()
|
||||
info("Starting iteration $iter")
|
||||
assembly1 = Assembly()
|
||||
assemble!(assembly1, problem1, time)
|
||||
assembly2 = Assembly()
|
||||
assemble!(assembly2, problem2, time)
|
||||
|
||||
A1 = sparse(assembly1.stiffness_matrix)
|
||||
dims = size(A1)
|
||||
b1 = sparse(assembly1.force_vector, dims[1], 1)
|
||||
A2 = sparse(assembly2.stiffness_matrix, dims[1], dims[2])
|
||||
b2 = sparse(assembly2.force_vector, dims[1], 1)
|
||||
|
||||
# create a saddle point problem
|
||||
A = [A1 A2; A2' zeros(A2)]
|
||||
b = [b1; b2]
|
||||
mapper = solver.parallel ? pmap : map
|
||||
|
||||
if length(b) != length(x)
|
||||
info("iter $iter: resizing solution vector")
|
||||
resize!(x, length(b))
|
||||
resize!(dx, length(b))
|
||||
fill!(x, 0.0)
|
||||
fill!(dx, 0.0)
|
||||
end
|
||||
# assemble boundary problems
|
||||
boundary_assembly = sum(mapper((p)->assemble(p, time), solver.boundary_problems))
|
||||
boundary_dofs = unique(boundary_assembly.stiffness_matrix.I)
|
||||
|
||||
# solve problem, update solution vector
|
||||
# assemble field problems
|
||||
# in principle if we want to static condensation we need to pass boundary dofs
|
||||
# to field problems in order to know which dofs are interior dofs and can be
|
||||
# condensated.
|
||||
field_assembly = sum(mapper((p)->assemble(p, time), solver.field_problems))
|
||||
field_dofs = unique(field_assembly.stiffness_matrix.I)
|
||||
info("# of dofs: $(length(field_dofs)), # of interface dofs: $(length(boundary_dofs))")
|
||||
|
||||
# create sparse matrices and saddle point problem
|
||||
K = sparse(field_assembly.stiffness_matrix)
|
||||
dim = size(K, 1)
|
||||
r = sparse(field_assembly.force_vector, dim, 1)
|
||||
C = sparse(boundary_assembly.stiffness_matrix, dim, dim)
|
||||
g = sparse(boundary_assembly.force_vector, dim, 1)
|
||||
A = [K C'; C spzeros(dim, dim)]
|
||||
b = [r; g]
|
||||
|
||||
# solve increment for linearized problem
|
||||
nz = unique(rowvals(A)) # take only non-zero rows
|
||||
dx[nz] = lufact(A[nz,nz]) \ full(b[nz])
|
||||
x += dx
|
||||
sol = zeros(b)
|
||||
sol[nz] = lufact(A[nz,nz]) \ full(b[nz])
|
||||
info("solved. length of solution vector = $(length(sol))")
|
||||
#info(full(sol[nz]))
|
||||
|
||||
# get "problem-wise" solution vectors
|
||||
x1 = x[1:dims[1]]
|
||||
x2 = x[dims[1]+1:end]
|
||||
|
||||
# update field for elements in problem 1
|
||||
for equation in get_equations(problem1)
|
||||
element = get_element(equation)
|
||||
field_name = get_unknown_field_name(problem1)
|
||||
gdofs = get_gdofs(problem1, equation)
|
||||
local_sol = vec(full(x1[gdofs]))
|
||||
eqsize = size(equation)
|
||||
if eqsize[1] != 1
|
||||
# update elements in field problems
|
||||
for field_problem in solver.field_problems
|
||||
for equation in get_equations(field_problem)
|
||||
element = get_element(equation)
|
||||
gdofs = get_gdofs(field_problem, equation)
|
||||
eqsize = size(equation)
|
||||
local_sol = vec(full(sol[gdofs])) # incremental data for element
|
||||
local_sol = reshape(local_sol, eqsize)
|
||||
local_sol = Vector{Float64}[local_sol[:,i] for i=1:size(local_sol,2)]
|
||||
last(element[field_name]).data += local_sol # <-- added
|
||||
end
|
||||
#info("problem1: pushing to $field_name")
|
||||
push!(element[field_name], time => local_sol)
|
||||
end
|
||||
|
||||
# update field for elements in problem 2 (Dirichlet boundary)
|
||||
for equation in get_equations(problem2)
|
||||
element = get_element(equation)
|
||||
field_name = "reaction force" #get_unknown_field_name(problem2)
|
||||
gdofs = get_gdofs(problem2, equation)
|
||||
local_sol = vec(full(x1[gdofs]))
|
||||
eqsize = size(equation)
|
||||
if eqsize[1] != 1
|
||||
local_sol = reshape(local_sol, eqsize)
|
||||
# update elements in boundary problems
|
||||
for boundary_problem in solver.boundary_problems
|
||||
for equation in get_equations(boundary_problem)
|
||||
element = get_element(equation)
|
||||
gdofs = get_gdofs(boundary_problem, equation) + dim
|
||||
eqsize = size(equation)
|
||||
local_sol = vec(full(sol[gdofs]))
|
||||
#info("local sol = $local_sol")
|
||||
local_sol = reshape(local_sol, field_dim, eqsize[2])
|
||||
local_sol = Vector{Float64}[local_sol[:,i] for i=1:size(local_sol,2)]
|
||||
last(element["reaction force"]).data = local_sol # <-- replaced
|
||||
end
|
||||
#info("problem2: pushing to $field_name")
|
||||
push!(element[field_name], time => local_sol)
|
||||
end
|
||||
|
||||
if norm(dx[1:dims[1]]) < solver.tol
|
||||
return (iter, true)
|
||||
end
|
||||
|
||||
info("Iteration took $(toq()) seconds")
|
||||
|
||||
if norm(sol[1:dim]) < solver.tol
|
||||
return (iter, true)
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
info("Warning: did not coverge in $(solver.max_iterations) iterations!")
|
||||
|
||||
+4
-2
@@ -39,8 +39,10 @@ function Base.size(equation::DBC2D2)
|
||||
end
|
||||
|
||||
function Base.convert(::Type{DirichletEquation}, element::Seg2)
|
||||
integration_points = line3()
|
||||
haskey(element, "reaction force") || (element["reaction force"] = 0.0 => zeros(2))
|
||||
integration_points = get_integration_points(element, Val{3})
|
||||
if !haskey(element, "reaction force")
|
||||
element["reaction force"] = (0.0 => Vector{Float64}[])
|
||||
end
|
||||
DBC2D2(element, integration_points)
|
||||
end
|
||||
|
||||
|
||||
+2
-2
@@ -122,7 +122,7 @@ function Base.size(equation::CPS4)
|
||||
end
|
||||
|
||||
function Base.convert(::Type{PlaneStressElasticityEquation}, element::Quad4)
|
||||
integration_points = get_default_integration_points(element)
|
||||
integration_points = get_integration_points(element)
|
||||
if !haskey(element, "displacement")
|
||||
element["displacement"] = 0.0 => [zeros(2) for i=1:4]
|
||||
end
|
||||
@@ -140,7 +140,7 @@ function Base.size(equation::CPS2)
|
||||
end
|
||||
|
||||
function Base.convert(::Type{PlaneStressElasticityEquation}, element::Seg2)
|
||||
integration_points = get_default_integration_points(element)
|
||||
integration_points = get_integration_points(element)
|
||||
if !haskey(element, "displacement")
|
||||
element["displacement"] = 0.0 => [zeros(2) for i=1:2]
|
||||
end
|
||||
|
||||
+3
-3
@@ -67,9 +67,9 @@ function Base.getindex(element::Element, field_name)
|
||||
return element.fields[field_name]
|
||||
end
|
||||
|
||||
function get_integration_points(element)
|
||||
return get_default_integration_points(element)
|
||||
end
|
||||
#function get_integration_points(element)
|
||||
# return get_default_integration_points(element)
|
||||
#end
|
||||
|
||||
"""Add new Field to element.
|
||||
|
||||
|
||||
@@ -11,14 +11,10 @@ type Assembly
|
||||
mass_matrix :: SparseMatrixIJV
|
||||
stiffness_matrix :: SparseMatrixIJV
|
||||
force_vector :: SparseMatrixIJV
|
||||
lhs :: SparseMatrixIJV
|
||||
rhs :: SparseMatrixIJV
|
||||
end
|
||||
|
||||
function Assembly()
|
||||
return Assembly(
|
||||
SparseMatrixIJV(),
|
||||
SparseMatrixIJV(),
|
||||
SparseMatrixIJV(),
|
||||
SparseMatrixIJV(),
|
||||
SparseMatrixIJV())
|
||||
@@ -28,8 +24,6 @@ function Base.empty!(assembly::Assembly)
|
||||
empty!(assembly.mass_matrix)
|
||||
empty!(assembly.stiffness_matrix)
|
||||
empty!(assembly.force_vector)
|
||||
empty!(assembly.lhs)
|
||||
empty!(assembly.rhs)
|
||||
end
|
||||
|
||||
function get_mass_matrix
|
||||
@@ -184,8 +178,6 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro
|
||||
return R
|
||||
end
|
||||
|
||||
#info("field = $field")
|
||||
#info("vec(field) = $(vec(field))")
|
||||
jacobian, allresults = ForwardDiff.jacobian(calc_R, vec(field), AllResults, cache=autodiffcache)
|
||||
add!(assembly.stiffness_matrix, gdofs, gdofs, jacobian)
|
||||
add!(assembly.force_vector, gdofs, -ForwardDiff.value(allresults))
|
||||
|
||||
@@ -208,6 +208,18 @@ function Base.similar{T}(field::DVTI, data::Vector{T})
|
||||
return typeof(field)(newdata)
|
||||
end
|
||||
|
||||
function Base.start(::DVTI)
|
||||
return 1
|
||||
end
|
||||
|
||||
function Base.next(f::DVTI, state)
|
||||
return f.data[state], state+1
|
||||
end
|
||||
|
||||
function Base.done(f::DVTI, s)
|
||||
return s > length(f.data)
|
||||
end
|
||||
|
||||
### Accessing continuous fields
|
||||
|
||||
function Base.call(field::CVTI, xi::Vector)
|
||||
|
||||
+2
-2
@@ -93,13 +93,13 @@ end
|
||||
# Conversions element -> equation
|
||||
|
||||
function Base.convert(::Type{HeatEquation}, element::Quad4)
|
||||
integration_points = get_default_integration_points(element)
|
||||
integration_points = get_integration_points(element)
|
||||
haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(4))
|
||||
DC2D4(element, integration_points)
|
||||
end
|
||||
|
||||
function Base.convert(::Type{HeatEquation}, element::Seg2)
|
||||
integration_points = get_default_integration_points(element)
|
||||
integration_points = get_integration_points(element)
|
||||
haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(2))
|
||||
DC2D2(element, integration_points)
|
||||
end
|
||||
|
||||
+13
-10
@@ -1,8 +1,9 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# Let's drop here all integration schemes and some defaults for different element types
|
||||
|
||||
function get_default_integration_points(element::Quad4)
|
||||
function get_integration_points(Quad4::Element)
|
||||
[
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
|
||||
@@ -11,21 +12,22 @@ function get_default_integration_points(element::Quad4)
|
||||
]
|
||||
end
|
||||
|
||||
typealias LineElement Union{Seg2, Seg3}
|
||||
|
||||
function line1()
|
||||
function get_integration_points(element::LineElement, ::Type{Val{1}})
|
||||
[
|
||||
IntegrationPoint([0.0], 2.0)
|
||||
]
|
||||
end
|
||||
|
||||
function line2()
|
||||
function get_integration_points(element::LineElement, ::Type{Val{2}})
|
||||
[
|
||||
IntegrationPoint([-sqrt(1/3)], 1)
|
||||
IntegrationPoint([+sqrt(1/3)], 1)
|
||||
]
|
||||
end
|
||||
|
||||
function line3()
|
||||
function get_integration_points(element::LineElement, ::Type{Val{3}})
|
||||
[
|
||||
IntegrationPoint([0.0], 8/9),
|
||||
IntegrationPoint([-sqrt(3/5)], 5/9),
|
||||
@@ -33,7 +35,7 @@ function line3()
|
||||
]
|
||||
end
|
||||
|
||||
function line4()
|
||||
function get_integration_points(element::LineElement, ::Type{Val{4}})
|
||||
[
|
||||
IntegrationPoint([+sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
|
||||
IntegrationPoint([-sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
|
||||
@@ -42,7 +44,7 @@ function line4()
|
||||
]
|
||||
end
|
||||
|
||||
function line5()
|
||||
function get_integration_points(element::LineElement, ::Type{Val{5}})
|
||||
[
|
||||
IntegrationPoint([-1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900),
|
||||
IntegrationPoint([-1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),
|
||||
@@ -52,10 +54,11 @@ function line5()
|
||||
]
|
||||
end
|
||||
|
||||
function get_default_integration_points(element::Seg2)
|
||||
return line1()
|
||||
function get_integration_points(element::Seg2)
|
||||
return get_integration_points(element, Val{1})
|
||||
end
|
||||
|
||||
function get_default_integration_points(element::MSeg2)
|
||||
return line3()
|
||||
function get_integration_points(element::Seg3)
|
||||
return get_integration_points(element, Val{2})
|
||||
end
|
||||
|
||||
|
||||
+148
-24
@@ -1,17 +1,132 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# Mortar equations
|
||||
# Mortar projection calculation for 2d
|
||||
|
||||
""" Find projection from slave nodes to master element, i.e. find xi2 from
|
||||
master element corresponding to the xi1.
|
||||
"""
|
||||
function project_from_slave_to_master(slave::Element, master::Element, xi1::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9)
|
||||
# slave_basis = get_basis(slave)
|
||||
|
||||
# slave side geometry and normal direction at xi1
|
||||
X1 = slave("geometry", xi1, time)
|
||||
N1 = slave("nodal ntsys", xi1, time)[:,1]
|
||||
|
||||
# master side geometry at xi2
|
||||
master_basis = master.basis.data.basis
|
||||
master_dbasis = master.basis.data.dbasis
|
||||
master_geometry = master("geometry")(time)
|
||||
|
||||
function X2(xi2)
|
||||
N = master_basis([xi2])
|
||||
return sum([N[i]*master_geometry[i] for i=1:length(N)])
|
||||
end
|
||||
|
||||
function dX2(xi2)
|
||||
dN = master_dbasis([xi2])
|
||||
return sum([dN[i]*master_geometry[i] for i=1:length(dN)])
|
||||
end
|
||||
|
||||
# master_basis = get_basis(master)
|
||||
# X2(xi2) = master_basis("geometry", [xi2], time)
|
||||
# dX2(xi2) = dmaster_basis("geometry", xi2, time)
|
||||
|
||||
# equation to solve
|
||||
R(xi2) = det([X2(xi2)-X1 N1]')
|
||||
dR(xi2) = det([dX2(xi2) N1]')
|
||||
# dR = ForwardDiff.derivative(R)
|
||||
|
||||
# go!
|
||||
xi2 = 0.0
|
||||
for i=1:max_iterations
|
||||
dxi2 = -R(xi2) / dR(xi2)
|
||||
xi2 += dxi2
|
||||
if norm(dxi2) < tol
|
||||
return Float64[xi2]
|
||||
end
|
||||
end
|
||||
error("find projection from slave to master: did not converge")
|
||||
end
|
||||
|
||||
""" Find projection from master surface to slave point, i.e. find xi1 from slave
|
||||
element corresponding to the xi2. """
|
||||
function project_from_master_to_slave(slave::Element, master::Element, xi2::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9)
|
||||
# slave_basis = get_basis(slave)
|
||||
|
||||
# slave side geometry and normal direction at xi1
|
||||
|
||||
slave_geometry = slave("geometry")(time)
|
||||
slave_normals = slave("nodal ntsys")(time)
|
||||
slave_basis = slave.basis.data.basis
|
||||
slave_dbasis = slave.basis.data.dbasis
|
||||
|
||||
function X1(xi1)
|
||||
N = slave_basis([xi1])
|
||||
return sum([N[i]*slave_geometry[i] for i=1:length(N)])
|
||||
end
|
||||
|
||||
function dX1(xi1)
|
||||
dN = slave_dbasis([xi1])
|
||||
return sum([dN[i]*slave_geometry[i] for i=1:length(dN)])
|
||||
end
|
||||
|
||||
function N1(xi1)
|
||||
N = slave_basis([xi1])
|
||||
return sum([N[i]*slave_normals[i] for i=1:length(N)])[:,1]
|
||||
end
|
||||
|
||||
function dN1(xi1)
|
||||
dN = slave_dbasis([xi1])
|
||||
return sum([dN[i]*slave_normals[i] for i=1:length(dN)])[:,1]
|
||||
end
|
||||
|
||||
#X1(xi1) = slave_basis("geometry", [xi1], time)
|
||||
#N1(xi1) = slave_basis("nodal ntsys", [xi1], time)[:,1]
|
||||
|
||||
#master_basis = get_basis(master)
|
||||
|
||||
# master side geometry at xi2
|
||||
#X2 = master_basis("geometry", xi2, time)
|
||||
X2 = master("geometry", xi2, time)
|
||||
|
||||
# equation to solve
|
||||
R(xi1) = det([X1(xi1)-X2 N1(xi1)]')
|
||||
dR(xi1) = det([dX1(xi1) N1(xi1)]') + det([X1(xi1)-X2 dN1(xi1)]')
|
||||
|
||||
#=
|
||||
info("R(-1.0) = $(R(-1.0))")
|
||||
info("R( 0.0) = $(R(0.0))")
|
||||
info("R( 1.0) = $(R(1.0))")
|
||||
info("R( 1.5) = $(R(1.5))")
|
||||
info("dR(-1.0) = $(dR(-1.0))")
|
||||
info("dR( 0.0) = $(dR(0.0))")
|
||||
info("dR( 1.0) = $(dR(1.0))")
|
||||
info("dR( 1.5) = $(dR(1.5))")
|
||||
=#
|
||||
|
||||
#dR = ForwardDiff.derivative(R)
|
||||
|
||||
# go!
|
||||
xi1 = 0.0
|
||||
for i=1:max_iterations
|
||||
dxi1 = -R(xi1) / dR(xi1)
|
||||
xi1 += dxi1
|
||||
if norm(dxi1) < tol
|
||||
return Float64[xi1]
|
||||
end
|
||||
end
|
||||
error("find projection from master to slave: did not converge")
|
||||
end
|
||||
|
||||
|
||||
### Mortar equations
|
||||
|
||||
abstract MortarEquation <: Equation
|
||||
|
||||
function get_unknown_field_name(equation::MortarEquation)
|
||||
return "reaction force"
|
||||
end
|
||||
|
||||
""" Mortar boundary condition element for 2-dimensional problem, 2 node line segment. """
|
||||
type MBC2D2 <: MortarEquation
|
||||
element :: MSeg2
|
||||
element :: Seg2
|
||||
integration_points :: Vector{IntegrationPoint}
|
||||
end
|
||||
|
||||
@@ -19,11 +134,16 @@ function Base.size(equation::MBC2D2)
|
||||
return (1, 2)
|
||||
end
|
||||
|
||||
function Base.convert(::Type{MortarEquation}, element::MSeg2)
|
||||
return MBC2D2(element, get_default_integration_points(element))
|
||||
function Base.convert(::Type{MortarEquation}, element::Seg2)
|
||||
integration_points = get_integration_points(element, Val{3})
|
||||
if !haskey(element, "reaction force")
|
||||
element["reaction force"] = (0.0 => Vector{Float64}[])
|
||||
end
|
||||
MBC2D2(element, integration_points)
|
||||
end
|
||||
|
||||
# Mortar problem
|
||||
|
||||
### Mortar problem
|
||||
|
||||
"""
|
||||
Parameters
|
||||
@@ -38,26 +158,24 @@ type MortarProblem <: BoundaryProblem
|
||||
equations :: Vector{MortarEquation}
|
||||
end
|
||||
|
||||
function MortarProblem(dimension::Int=1, equations=[])
|
||||
MortarProblem("reaction force", dimension, equations)
|
||||
function MortarProblem(unknown_field_name, unknown_field_dimension::Int=1)
|
||||
MortarProblem(unknown_field_name, unknown_field_dimension, [])
|
||||
end
|
||||
|
||||
# Mortar projection calculation
|
||||
|
||||
""" Find master or "mortar" elements for this slave element. """
|
||||
function get_master_elements(element::MortarElement)
|
||||
return element.master_elements
|
||||
end
|
||||
# Mortar assembly
|
||||
|
||||
function assemble!(assembly::Assembly, equation::MortarEquation, time::Number=0.0, problem=nothing)
|
||||
isa(problem, Void) && error("Mortar boundary problem needs problem to be defined")
|
||||
field_dim = problem.unknown_field_dimension
|
||||
field_name = problem.unknown_field_name
|
||||
|
||||
slave_element = get_element(equation)
|
||||
master_elements = get_master_elements(slave_element)
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
slave_basis = get_basis(slave_element)
|
||||
detJ = det(slave_basis)
|
||||
dim = size(equation, 1) # number of nodes
|
||||
slave_dofs = get_gdofs(slave_element, dim)
|
||||
for master_element in master_elements
|
||||
master_dofs = get_gdofs(master_element, dim)
|
||||
|
||||
for master_element in slave_element["master elements"]
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
xi1a = project_from_master_to_slave(slave_element, master_element, [-1.0])
|
||||
xi1b = project_from_master_to_slave(slave_element, master_element, [ 1.0])
|
||||
xi1 = clamp([xi1a xi1b], -1.0, 1.0)
|
||||
@@ -78,8 +196,14 @@ function assemble!(assembly::Assembly, equation::MortarEquation, time::Number=0.
|
||||
# add contribution to left hand side
|
||||
N1 = slave_basis(xi_gauss, time)
|
||||
N2 = master_basis(xi_projected, time)
|
||||
add!(assembly.lhs, slave_dofs, slave_dofs, w*N1'*N1)
|
||||
add!(assembly.lhs, slave_dofs, master_dofs, -w*N1'*N2)
|
||||
S = w*N1'*N1
|
||||
M = w*N1'*N2
|
||||
for i=1:field_dim
|
||||
sd = slave_dofs[i:field_dim:end]
|
||||
md = master_dofs[i:field_dim:end]
|
||||
add!(assembly.stiffness_matrix, sd, sd, S)
|
||||
add!(assembly.stiffness_matrix, sd, md, -M)
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
@@ -21,64 +21,4 @@ function MSeg2(connectivity, master_elements=[], biorthogonal=false)
|
||||
return MSeg2(connectivity, Basis(basis, dbasisdxi), FieldSet(), master_elements)
|
||||
end
|
||||
|
||||
""" Find projection from slave nodes to master element, i.e. find xi2 from
|
||||
master element corresponding to the xi1.
|
||||
"""
|
||||
function project_from_slave_to_master(slave::MortarElement, master::MortarElement, xi1::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9)
|
||||
slave_basis = get_basis(slave)
|
||||
master_basis = get_basis(master)
|
||||
|
||||
# slave side geometry and normal direction at xi1
|
||||
X1 = slave_basis("geometry", xi1, time)
|
||||
N1 = slave_basis("nodal ntsys", xi1, time)[:,1]
|
||||
|
||||
# master side geometry at xi2
|
||||
X2(xi2) = master_basis("geometry", [xi2], time)
|
||||
# dX2(xi2) = dmaster_basis("geometry", xi2, time)
|
||||
|
||||
# equation to solve
|
||||
R(xi2) = det([X2(xi2)-X1 N1]')
|
||||
# dR(xi2) = det([dX2(xi2) N1]')
|
||||
dR = ForwardDiff.derivative(R)
|
||||
|
||||
# go!
|
||||
xi2 = 0.0
|
||||
for i=1:max_iterations
|
||||
dxi2 = -R(xi2) / dR(xi2)
|
||||
xi2 += dxi2
|
||||
if norm(dxi2) < tol
|
||||
return Float64[xi2]
|
||||
end
|
||||
end
|
||||
error("find projection from slave to master: did not converge")
|
||||
end
|
||||
|
||||
""" Find projection from master surface to slave point, i.e. find xi1 from slave element corresponding to the xi2. """
|
||||
function project_from_master_to_slave(slave::MortarElement, master::MortarElement, xi2::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9)
|
||||
slave_basis = get_basis(slave)
|
||||
master_basis = get_basis(master)
|
||||
|
||||
# slave side geometry and normal direction at xi1
|
||||
X1(xi1) = slave_basis("geometry", [xi1], time)
|
||||
N1(xi1) = slave_basis("nodal ntsys", [xi1], time)[:,1]
|
||||
|
||||
# master side geometry at xi2
|
||||
X2 = master_basis("geometry", xi2, time)
|
||||
|
||||
# equation to solve
|
||||
R(xi1) = det([X1(xi1)-X2 N1(xi1)]')
|
||||
# dR(xi1) = det([dX1(xi1) N1(xi1)]') + det([X1(xi1)-X2 dN1(xi1)]')
|
||||
dR = ForwardDiff.derivative(R)
|
||||
|
||||
# go!
|
||||
xi1 = 0.0
|
||||
for i=1:max_iterations
|
||||
dxi1 = -R(xi1) / dR(xi1)
|
||||
xi1 += dxi1
|
||||
if norm(dxi1) < tol
|
||||
return Float64[xi1]
|
||||
end
|
||||
end
|
||||
error("find projection from master to slave: did not converge")
|
||||
end
|
||||
|
||||
|
||||
+1
-1
@@ -170,7 +170,7 @@ function call(solver::SimpleSolver, time::Number=0.0)
|
||||
local_sol = reshape(local_sol, eqsize)
|
||||
end
|
||||
#info("problem2: pushing to $field_name")
|
||||
push!(element[field_name], time => local_sol)
|
||||
#push!(element[field_name], time => local_sol)
|
||||
end
|
||||
|
||||
return norm(x1)
|
||||
|
||||
@@ -36,6 +36,21 @@ function Base.append!(A::SparseMatrixIJV, I::Vector{Int}, J::Vector{Int}, V::Vec
|
||||
append!(A.V, V)
|
||||
end
|
||||
|
||||
function Base.isempty(A::SparseMatrixIJV)
|
||||
return isempty(A.I) && isempty(A.J) && isempty(A.V)
|
||||
end
|
||||
|
||||
function Base.(:+)(A::SparseMatrixIJV, B::SparseMatrixIJV)
|
||||
if isempty(A)
|
||||
return B
|
||||
end
|
||||
if isempty(B)
|
||||
return A
|
||||
end
|
||||
C = SparseMatrixIJV([A.I;B.I], [A.J;B.J], [A.V;B.V])
|
||||
return C
|
||||
end
|
||||
|
||||
function Base.full(A::SparseMatrixIJV, args...)
|
||||
return full(sparse(A.I, A.J, A.V, args...))
|
||||
end
|
||||
|
||||
+233
-21
@@ -6,8 +6,9 @@ module MortarTests
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Test
|
||||
|
||||
using JuliaFEM: MSeg2, Seg2, MortarProblem, MortarEquation, MortarElement, Assembly, assemble!
|
||||
using JuliaFEM: get_basis, grad, project_from_slave_to_master, project_from_master_to_slave
|
||||
using JuliaFEM: Seg2, MortarProblem, MortarEquation, MortarElement, Assembly, assemble!, Element
|
||||
using JuliaFEM: get_basis, grad, project_from_slave_to_master, project_from_master_to_slave, Quad4
|
||||
using JuliaFEM: PlaneStressElasticityProblem, DirichletProblem, DirectSolver
|
||||
|
||||
function get_test_2d_model()
|
||||
# this is hand calculated and given as an example in my thesis
|
||||
@@ -17,22 +18,24 @@ function get_test_2d_model()
|
||||
[0.0, 1.0], [5/4, 1.0], [2.0, 1.0],
|
||||
[0.0, 1.0], [3/4, 1.0], [2.0, 1.0]]
|
||||
rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)]
|
||||
slave1 = MSeg2([10, 11])
|
||||
|
||||
master1 = Seg2([7, 8])
|
||||
master1["geometry"] = Vector[N[7], N[8]]
|
||||
master2 = Seg2([8, 9])
|
||||
master2["geometry"] = Vector[N[8], N[9]]
|
||||
|
||||
slave1 = Seg2([10, 11])
|
||||
slave1["geometry"] = Vector[N[10], N[11]]
|
||||
# should be n = [0 -1]' and t = [1 0]'
|
||||
slave1["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
|
||||
slave2 = MSeg2([11, 12])
|
||||
slave1["master elements"] = Element[master1, master2]
|
||||
|
||||
slave2 = Seg2([11, 12])
|
||||
slave2["geometry"] = Vector[N[11], N[12]]
|
||||
# should be n = [0 -1]' and t = [1 0]'
|
||||
slave2["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
|
||||
master1 = MSeg2([7, 8])
|
||||
master1["geometry"] = Vector[N[7], N[8]]
|
||||
master2 = MSeg2([8, 9])
|
||||
master2["geometry"] = Vector[N[8], N[9]]
|
||||
push!(slave1.master_elements, master1)
|
||||
push!(slave1.master_elements, master2)
|
||||
push!(slave2.master_elements, master1)
|
||||
push!(slave2.master_elements, master2)
|
||||
slave2["master elements"] = Element[master1, master2]
|
||||
|
||||
return [slave1, slave2], [master1, master2]
|
||||
end
|
||||
|
||||
@@ -57,13 +60,36 @@ function test_calc_flat_2d_projection()
|
||||
@test X1 == [5/4, 1.0]
|
||||
end
|
||||
|
||||
function test_calc_flat_2d_projection_rotated()
|
||||
master1 = Seg2([3, 4])
|
||||
master1["geometry"] = Vector{Float64}[[0.0, 1.0], [0.0, 0.0]]
|
||||
slave1 = Seg2([1, 2])
|
||||
slave1["geometry"] = Vector{Float64}[[0.0, 0.0], [0.0, 1.0]]
|
||||
slave1["nodal ntsys"] = Matrix{Float64}[[1.0 0.0; 0.0 1.0], [1.0 0.0; 0.0 1.0]]
|
||||
xi = project_from_master_to_slave(slave1, master1, [-1.0])
|
||||
info("xi = $xi")
|
||||
@test xi == [ 1.0]
|
||||
xi = project_from_master_to_slave(slave1, master1, [1.0])
|
||||
info("xi = $xi")
|
||||
@test xi == [-1.0]
|
||||
|
||||
xi = project_from_slave_to_master(slave1, master1, [-1.0])
|
||||
info("xi = $xi")
|
||||
@test xi == [ 1.0]
|
||||
xi = project_from_slave_to_master(slave1, master1, [1.0])
|
||||
info("xi = $xi")
|
||||
@test xi == [-1.0]
|
||||
|
||||
|
||||
end
|
||||
|
||||
function test_create_flat_2d_assembly()
|
||||
slaves, masters = get_test_2d_model()
|
||||
slave1, slave2 = slaves
|
||||
master1, master2 = masters
|
||||
|
||||
info("creating problem")
|
||||
problem = MortarProblem()
|
||||
problem = MortarProblem("temperature", 1)
|
||||
info("pushing slave elements to problem")
|
||||
push!(problem, slave1)
|
||||
push!(problem, slave2)
|
||||
@@ -77,7 +103,7 @@ function test_create_flat_2d_assembly()
|
||||
info("creating assembly")
|
||||
assembly = Assembly()
|
||||
assemble!(assembly, problem.equations[1], 0.0, problem)
|
||||
B = round(full(assembly.lhs, 12, 12), 6)
|
||||
B = round(full(assembly.stiffness_matrix, 12, 12), 6)
|
||||
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,:])")
|
||||
@@ -95,7 +121,7 @@ function test_create_flat_2d_assembly()
|
||||
B_expected[S3,S3] += [9/100 27/200; 27/200 39/100]
|
||||
B_expected[S3,M3] -= [3/20 3/40; 9/40 3/10]
|
||||
assemble!(assembly, problem.equations[2], 0.0, problem)
|
||||
B = full(assembly.lhs)
|
||||
B = full(assembly.stiffness_matrix)
|
||||
info("size of B = $(size(B))")
|
||||
info("B matrix in second slave element = \n$(B[11:12,:])")
|
||||
info("B matrix expected = \n$(B_expected[11:12,:])")
|
||||
@@ -103,16 +129,202 @@ function test_create_flat_2d_assembly()
|
||||
@test isapprox(B, B_expected)
|
||||
end
|
||||
|
||||
function test_patch_test_heat_2d()
|
||||
function test_2d_mortar_multiple_bodies_multiple_dirichlet_bc()
|
||||
N = Vector[
|
||||
[0.0, 0.0], [1.0, 0.0],
|
||||
[0.0, 1.0], [1.0, 1.0],
|
||||
[0.0, 1.0], [1.0, 1.0],
|
||||
[0.0, 2.0], [1.0, 2.0]]
|
||||
|
||||
e1 = Quad4([1, 2, 4, 3])
|
||||
e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
|
||||
e2 = Quad4([5, 6, 8, 7])
|
||||
e2["geometry"] = Vector[N[5], N[6], N[8], N[7]]
|
||||
for el in [e1, e2]
|
||||
el["youngs modulus"] = 900.0
|
||||
el["poissons ratio"] = 0.25
|
||||
end
|
||||
b1 = Seg2([7, 8])
|
||||
b1["geometry"] = Vector[N[7], N[8]]
|
||||
b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
|
||||
|
||||
|
||||
slaves, masters = get_test_2d_model()
|
||||
problem2 = MortarProblem()
|
||||
for slave in slaves:
|
||||
push!(problem2, slave)
|
||||
body1 = PlaneStressElasticityProblem()
|
||||
push!(body1, e1)
|
||||
|
||||
body2 = PlaneStressElasticityProblem()
|
||||
push!(body2, e2)
|
||||
push!(body2, b1)
|
||||
|
||||
# boundary elements for dirichlet dx=0
|
||||
dx1 = Seg2([1, 3])
|
||||
dx1["geometry"] = Vector[N[1], N[3]]
|
||||
dx2 = Seg2([5, 7])
|
||||
dx2["geometry"] = Vector[N[5], N[7]]
|
||||
for dx in [dx1, dx2]
|
||||
dx["displacement 1"] = 0.0
|
||||
end
|
||||
|
||||
boundary1 = DirichletProblem("displacement", 2)
|
||||
push!(boundary1, dx1)
|
||||
push!(boundary1, dx2)
|
||||
|
||||
# boundary elements for dirichlet dy=0
|
||||
dy1 = Seg2([1, 2])
|
||||
dy1["geometry"] = Vector[N[1], N[2]]
|
||||
dy1["displacement 2"] = 0.0
|
||||
|
||||
boundary2 = DirichletProblem("displacement", 2)
|
||||
push!(boundary2, dy1)
|
||||
|
||||
# mortar boundary between two bodies
|
||||
rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)]
|
||||
|
||||
master1 = Seg2([3, 4])
|
||||
master1["geometry"] = Vector[N[3], N[4]]
|
||||
|
||||
slave1 = Seg2([5, 6])
|
||||
slave1["geometry"] = Vector[N[5], N[6]]
|
||||
slave1["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
|
||||
slave1["master elements"] = Element[master1]
|
||||
|
||||
boundary3 = MortarProblem("displacement", 2)
|
||||
push!(boundary3, slave1)
|
||||
|
||||
solver = DirectSolver()
|
||||
push!(solver, body1)
|
||||
push!(solver, body2)
|
||||
push!(solver, boundary1)
|
||||
push!(solver, boundary2)
|
||||
push!(solver, boundary3)
|
||||
|
||||
# launch solver
|
||||
norm = solver(0.0)
|
||||
|
||||
disp = e2("displacement", [1.0, 1.0], 0.0)
|
||||
info("displacement at tip: $disp")
|
||||
# code aster verification, two_elements.comm
|
||||
@test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01])
|
||||
end
|
||||
|
||||
|
||||
function test_2d_mortar_three_bodies_shared_nodes()
|
||||
N = Vector[
|
||||
[0.0, 0.0], [2.0, 0.0],
|
||||
[0.0, 1.0], [2.0, 1.0],
|
||||
[0.0, 1.0], [1.0, 1.0],
|
||||
[0.0, 2.0], [1.0, 2.0],
|
||||
[1.0, 1.0], [2.0, 1.0],
|
||||
[1.0, 2.0], [2.0, 2.0]]
|
||||
|
||||
e1 = Quad4([1, 2, 4, 3])
|
||||
e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
|
||||
|
||||
e2 = Quad4([5, 6, 8, 7])
|
||||
e2["geometry"] = Vector[N[5], N[6], N[8], N[7]]
|
||||
|
||||
e3 = Quad4([9, 10, 12, 11])
|
||||
e3["geometry"] = Vector[N[9], N[10], N[12], N[11]]
|
||||
|
||||
for el in [e1, e2, e3]
|
||||
el["youngs modulus"] = 900.0
|
||||
el["poissons ratio"] = 0.25
|
||||
end
|
||||
|
||||
b1 = Seg2([7, 8])
|
||||
b1["geometry"] = Vector[N[7], N[8]]
|
||||
b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
|
||||
|
||||
b2 = Seg2([11, 12])
|
||||
b2["geometry"] = Vector[N[11], N[12]]
|
||||
b2["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
|
||||
|
||||
body1 = PlaneStressElasticityProblem()
|
||||
push!(body1, e1)
|
||||
|
||||
body2 = PlaneStressElasticityProblem()
|
||||
push!(body2, e2)
|
||||
push!(body2, b1)
|
||||
|
||||
body3 = PlaneStressElasticityProblem()
|
||||
push!(body3, e3)
|
||||
push!(body3, b2)
|
||||
|
||||
# boundary elements for dirichlet dx=0
|
||||
dx1 = Seg2([1, 3])
|
||||
dx1["geometry"] = Vector[N[1], N[3]]
|
||||
dx2 = Seg2([5, 7])
|
||||
dx2["geometry"] = Vector[N[5], N[7]]
|
||||
for dx in [dx1, dx2]
|
||||
dx["displacement 1"] = 0.0
|
||||
end
|
||||
|
||||
bc1 = DirichletProblem("displacement", 2)
|
||||
push!(bc1, dx1)
|
||||
push!(bc1, dx2)
|
||||
|
||||
# boundary elements for dirichlet dy=0
|
||||
dy1 = Seg2([1, 2])
|
||||
dy1["geometry"] = Vector[N[1], N[2]]
|
||||
dy1["displacement 2"] = 0.0
|
||||
|
||||
bc2 = DirichletProblem("displacement", 2)
|
||||
push!(bc2, dy1)
|
||||
|
||||
# mortar boundary between body 1 and body 2
|
||||
rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)]
|
||||
|
||||
master1 = Seg2([3, 4])
|
||||
master1["geometry"] = Vector[N[3], N[4]]
|
||||
|
||||
slave1 = Seg2([5, 6])
|
||||
slave1["geometry"] = Vector[N[5], N[6]]
|
||||
slave1["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
|
||||
slave1["master elements"] = Element[master1]
|
||||
bc3 = MortarProblem("displacement", 2)
|
||||
push!(bc3, slave1)
|
||||
|
||||
# mortar boundary between body 1 and body 3
|
||||
slave2 = Seg2([9, 10])
|
||||
slave2["geometry"] = Vector[N[9], N[10]]
|
||||
slave2["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
|
||||
slave2["master elements"] = Element[master1]
|
||||
bc4 = MortarProblem("displacement", 2)
|
||||
push!(bc4, slave2)
|
||||
|
||||
# mortar boundary between body 2 and body 3
|
||||
master2 = Seg2([9, 11])
|
||||
master2["geometry"] = Vector[N[9], N[11]]
|
||||
|
||||
slave3 = Seg2([6, 8])
|
||||
slave3["geometry"] = Vector[N[6], N[8]]
|
||||
#slave3["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
|
||||
slave3["nodal ntsys"] = Matrix[rotation_matrix(0.0), rotation_matrix(0.0)]
|
||||
slave3["master elements"] = Element[master2]
|
||||
bc5 = MortarProblem("displacement", 2)
|
||||
push!(bc5, slave3)
|
||||
|
||||
solver = DirectSolver()
|
||||
push!(solver, body1)
|
||||
push!(solver, body2)
|
||||
push!(solver, body3)
|
||||
|
||||
push!(solver, bc1)
|
||||
push!(solver, bc2)
|
||||
|
||||
push!(solver, bc3)
|
||||
push!(solver, bc4)
|
||||
push!(solver, bc5)
|
||||
|
||||
# launch solver
|
||||
norm = solver(0.0)
|
||||
|
||||
disp = e2("displacement", [1.0, 1.0], 0.0)
|
||||
info("displacement at tip: $disp")
|
||||
# code aster verification, two_elements.comm
|
||||
@test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01])
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||
end
|
||||
|
||||
+73
-4
@@ -67,6 +67,7 @@ function test_simplesolver()
|
||||
info("Temperature at point X = $X is T = $T")
|
||||
@test isapprox(T, 100.0)
|
||||
end
|
||||
#test_simplesolver()
|
||||
|
||||
function atest_direct_solver()
|
||||
|
||||
@@ -190,13 +191,13 @@ function test_solver_multiple_dirichlet_bc()
|
||||
b1["geometry"] = Vector[N[3], N[4]]
|
||||
b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
|
||||
|
||||
#free_dofs = [3, 5, 6, 8]
|
||||
free_dofs = [3, 6, 7, 8]
|
||||
problem = PlaneStressElasticityProblem()
|
||||
push!(problem, e1)
|
||||
push!(problem, b1)
|
||||
|
||||
# manually solve problem 1
|
||||
# free_dofs = [3, 5, 6, 8]
|
||||
# free_dofs = [3, 6, 7, 8]
|
||||
#solve!(problem, free_dofs, 0.0; max_iterations=10)
|
||||
#disp = e1("displacement", [1.0, 1.0], 0.0)
|
||||
#info("displacement at tip: $disp")
|
||||
@@ -214,21 +215,89 @@ function test_solver_multiple_dirichlet_bc()
|
||||
|
||||
problem2 = DirichletProblem("displacement", 2)
|
||||
push!(problem2, dx)
|
||||
push!(problem2, dy)
|
||||
|
||||
problem3 = DirichletProblem("displacement", 2)
|
||||
push!(problem3, dy)
|
||||
|
||||
solver = DirectSolver()
|
||||
push!(solver, problem)
|
||||
push!(solver, problem2)
|
||||
push!(solver, problem3)
|
||||
|
||||
# launch solver
|
||||
norm = solver(0.0)
|
||||
|
||||
# info(e1("displacement"))
|
||||
# info(last(e1["displacement"]))
|
||||
disp = e1("displacement", [1.0, 1.0], 0.0)
|
||||
info("displacement at tip: $disp")
|
||||
@test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
|
||||
|
||||
end
|
||||
|
||||
# test_solver_multiple_dirichlet_bc()
|
||||
function test_solver_multiple_bodies_multiple_dirichlet_bc()
|
||||
N = Vector[
|
||||
[0.0, 0.0], [1.0, 0.0],
|
||||
[0.0, 1.0], [1.0, 1.0],
|
||||
[0.0, 2.0], [1.0, 2.0]]
|
||||
|
||||
e1 = Quad4([1, 2, 4, 3])
|
||||
e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
|
||||
e2 = Quad4([3, 4, 6, 5])
|
||||
e2["geometry"] = Vector[N[3], N[4], N[6], N[5]]
|
||||
for el in [e1, e2]
|
||||
el["youngs modulus"] = 900.0
|
||||
el["poissons ratio"] = 0.25
|
||||
end
|
||||
b1 = Seg2([5, 6])
|
||||
b1["geometry"] = Vector[N[5], N[6]]
|
||||
b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
|
||||
|
||||
body1 = PlaneStressElasticityProblem()
|
||||
push!(body1, e1)
|
||||
|
||||
body2 = PlaneStressElasticityProblem()
|
||||
push!(body2, e2)
|
||||
push!(body2, b1)
|
||||
|
||||
# boundary elements for dirichlet dx=0
|
||||
dx1 = Seg2([1, 3])
|
||||
dx1["geometry"] = Vector[N[1], N[3]]
|
||||
dx2 = Seg2([3, 5])
|
||||
dx2["geometry"] = Vector[N[3], N[5]]
|
||||
for dx in [dx1, dx2]
|
||||
dx["displacement 1"] = 0.0
|
||||
end
|
||||
|
||||
boundary1 = DirichletProblem("displacement", 2)
|
||||
push!(boundary1, dx1)
|
||||
push!(boundary1, dx2)
|
||||
|
||||
# boundary elements for dirichlet dy=0
|
||||
dy1 = Seg2([1, 2])
|
||||
dy1["geometry"] = Vector[N[1], N[2]]
|
||||
dy1["displacement 2"] = 0.0
|
||||
|
||||
boundary2 = DirichletProblem("displacement", 2)
|
||||
push!(boundary2, dy1)
|
||||
|
||||
|
||||
solver = DirectSolver()
|
||||
push!(solver, body1)
|
||||
push!(solver, body2)
|
||||
push!(solver, boundary1)
|
||||
push!(solver, boundary2)
|
||||
|
||||
# launch solver
|
||||
norm = solver(0.0)
|
||||
|
||||
disp = e2("displacement", [1.0, 1.0], 0.0)
|
||||
info("displacement at tip: $disp")
|
||||
# code aster verification, two_elements.comm
|
||||
@test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01])
|
||||
|
||||
end
|
||||
|
||||
# test_solver_multiple_bodies_multiple_dirichlet_bc()
|
||||
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user