From fe55cb0faa27f1b99ab7ededdcc4f5fae25fa249 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Tue, 24 Nov 2015 03:06:56 +0200 Subject: [PATCH] 2d tie contact working. --- notebooks/2015-11-23-2d-tie-contact.ipynb | 401 ++++++++++++++++ src/assembly.jl | 21 +- src/directsolver.jl | 151 +++--- src/dirichlet.jl | 6 +- src/elasticity.jl | 4 +- src/elements.jl | 6 +- src/equations.jl | 8 - src/fields2.jl | 12 + src/heat.jl | 4 +- src/integrate.jl | 23 +- src/mortar.jl | 172 ++++++- src/mortar_elements.jl | 60 --- src/solvers.jl | 2 +- src/sparse.jl | 15 + test/test_mortar.jl | 254 +++++++++- test/test_solver.jl | 77 ++- .../2015-10-22-plane-stress/two_elements.comm | 57 +++ .../2015-10-22-plane-stress/two_elements.mail | 20 + .../2015-10-22-plane-stress/two_elements.resu | 437 ++++++++++++++++++ 19 files changed, 1531 insertions(+), 199 deletions(-) create mode 100644 notebooks/2015-11-23-2d-tie-contact.ipynb create mode 100644 verification/2015-10-22-plane-stress/two_elements.comm create mode 100644 verification/2015-10-22-plane-stress/two_elements.mail create mode 100644 verification/2015-10-22-plane-stress/two_elements.resu diff --git a/notebooks/2015-11-23-2d-tie-contact.ipynb b/notebooks/2015-11-23-2d-tie-contact.ipynb new file mode 100644 index 0000000..14b0c63 --- /dev/null +++ b/notebooks/2015-11-23-2d-tie-contact.ipynb @@ -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", + "![model](http://s4.postimg.org/u9yqeryul/Screenshot_from_2015_11_24_02_00_00.png)\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 +} diff --git a/src/assembly.jl b/src/assembly.jl index 0bbbe3b..793a082 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -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 diff --git a/src/directsolver.jl b/src/directsolver.jl index 8fea44b..91aa2f9 100644 --- a/src/directsolver.jl +++ b/src/directsolver.jl @@ -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!") diff --git a/src/dirichlet.jl b/src/dirichlet.jl index cb18a78..2ebada9 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -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 diff --git a/src/elasticity.jl b/src/elasticity.jl index 8452dc8..af85b30 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -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 diff --git a/src/elements.jl b/src/elements.jl index 63ffaaa..f273812 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -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. diff --git a/src/equations.jl b/src/equations.jl index 8c610a5..34a4f2e 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -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)) diff --git a/src/fields2.jl b/src/fields2.jl index 98c95eb..58b53f7 100644 --- a/src/fields2.jl +++ b/src/fields2.jl @@ -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) diff --git a/src/heat.jl b/src/heat.jl index e7d63cf..70caa93 100644 --- a/src/heat.jl +++ b/src/heat.jl @@ -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 diff --git a/src/integrate.jl b/src/integrate.jl index 2cb807c..bd2cf24 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -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 + diff --git a/src/mortar.jl b/src/mortar.jl index 900ff30..15cd9e4 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -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 diff --git a/src/mortar_elements.jl b/src/mortar_elements.jl index 52263c2..d10643c 100644 --- a/src/mortar_elements.jl +++ b/src/mortar_elements.jl @@ -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 diff --git a/src/solvers.jl b/src/solvers.jl index 3f1cbd7..d1992ef 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -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) diff --git a/src/sparse.jl b/src/sparse.jl index 1e622f1..95a0390 100644 --- a/src/sparse.jl +++ b/src/sparse.jl @@ -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 diff --git a/test/test_mortar.jl b/test/test_mortar.jl index 9a788ca..82a1cf1 100644 --- a/test/test_mortar.jl +++ b/test/test_mortar.jl @@ -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 diff --git a/test/test_solver.jl b/test/test_solver.jl index c6e6e04..beeba9d 100644 --- a/test/test_solver.jl +++ b/test/test_solver.jl @@ -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 diff --git a/verification/2015-10-22-plane-stress/two_elements.comm b/verification/2015-10-22-plane-stress/two_elements.comm new file mode 100644 index 0000000..528f622 --- /dev/null +++ b/verification/2015-10-22-plane-stress/two_elements.comm @@ -0,0 +1,57 @@ +DEBUT() + +MAIL = LIRE_MAILLAGE() + +MO = AFFE_MODELE(MAILLAGE = MAIL, + AFFE = _F(MAILLE=('E1', 'E2', 'E3'), PHENOMENE='MECANIQUE', MODELISATION='C_PLAN')) + +MAT = DEFI_MATERIAU(ELAS = _F(E=900.0, NU=0.25)) + +CHMAT = AFFE_MATERIAU( + MAILLAGE = MAIL, + AFFE = _F(MAILLE=('E1', 'E2'), MATER=MAT)) + +BC = AFFE_CHAR_MECA( + MODELE = MO, + DDL_IMPO = ( + _F(NOEUD = ('N1','N2'), DY=0), + _F(NOEUD = ('N1','N3','N5'), DX=0))) + +LO = AFFE_CHAR_MECA( + MODELE = MO, + FORCE_CONTOUR = _F(MAILLE='E3', FY=-100.0)) + +LIST = DEFI_LIST_REEL( + DEBUT = 0, + INTERVALLE = _F(JUSQU_A=1.0, NOMBRE=1)) + +STEP = DEFI_FONCTION( + NOM_PARA='INST', + VALE=(0,0,1,1)) + +RESU = STAT_NON_LINE( + MODELE=MO, + CHAM_MATER=CHMAT, + EXCIT=( + _F(CHARGE=BC), + _F(CHARGE=LO, FONC_MULT=STEP)), + NEWTON=_F(REAC_INCR=1, MATRICE='TANGENTE', REAC_ITER=1), + COMP_ELAS=_F(DEFORMATION='GROT_GDEP'), + INCREMENT=_F(LIST_INST=LIST), + CONVERGENCE=_F(RESI_GLOB_RELA=1.0E-12)) + +RESU = CALC_CHAMP( + reuse=RESU, + RESULTAT=RESU, + MODELE=MO, + DEFORMATION=('EPSI_ELNO', 'EPSG_ELNO', 'EPSI_ELGA', 'EPSG_ELGA'), + ENERGIE=('ENEL_ELEM', 'ENEL_NOEU', 'ETOT_ELEM', 'ETOT_NOEU'), + FORCE=('FORC_NODA'), + CONTRAINTE='SIGM_ELNO') + +IMPR_RESU( + MODELE = MO, + FORMAT = 'RESULTAT', + RESU = _F(RESULTAT = RESU)) + +FIN() diff --git a/verification/2015-10-22-plane-stress/two_elements.mail b/verification/2015-10-22-plane-stress/two_elements.mail new file mode 100644 index 0000000..9e520a9 --- /dev/null +++ b/verification/2015-10-22-plane-stress/two_elements.mail @@ -0,0 +1,20 @@ + + COOR_2D + N1 0.0 0.0 + N2 1.0 0.0 + N3 0.0 1.0 + N4 1.0 1.0 + N5 0.0 2.0 + N6 1.0 2.0 + FINSF + + QUAD4 + E1 N1 N2 N4 N3 + E2 N3 N4 N6 N5 + FINSF + + SEG2 + E3 N5 N6 + FINSF + + FIN diff --git a/verification/2015-10-22-plane-stress/two_elements.resu b/verification/2015-10-22-plane-stress/two_elements.resu new file mode 100644 index 0000000..9b92423 --- /dev/null +++ b/verification/2015-10-22-plane-stress/two_elements.resu @@ -0,0 +1,437 @@ + + + -- CODE_ASTER -- VERSION : EXPLOITATION (stable) -- + + Version 11.4.0 du 05/06/2013 + Copyright EDF R&D 1991 - 2015 + + Exécution du : Mon Nov 23 16:32:44 2015 + Nom de la machine : jukka-desktop + Architecture : 64bit + Type de processeur : x86_64 + Système d'exploitation : Linux 3.13.0-68-generic + Langue des messages : en (UTF-8) + + + !------------------------------------------------------------------------------------! + ! ! + ! ! + ! Vous utilisez une vieille version de Code_Aster. ! + ! ! + ! En mettant à jour votre version, vous bénéficierez des dernières améliorations ! + ! apportées au code depuis 15 mois. ! + ! Si vous avez des développements privés, vous risquez d'avoir un travail ! + ! important de portage si vous ne suivez pas les mises à jour. ! + ! ! + ! ! + ! Ceci est une alarme. Si vous ne comprenez pas le sens de cette ! + ! alarme, vous pouvez obtenir des résultats inattendus ! ! + !------------------------------------------------------------------------------------! + + Parallélisme MPI : inactif + Parallélisme OpenMP : actif + Nombre de processus utilisés : 1 + Version de la librairie HDF5 : 1.8.8 + Version de la librairie MED : 3.0.6 + Librairie MUMPS : installée + Version de la librairie SCOTCH : 5.1.10 + Mémoire limite pour l'exécution : 4096.00 Mo + consommée par l'initialisation : 197.46 Mo + par les objets du jeu de commandes : 1.54 Mo + reste pour l'allocation dynamique : 3897.01 Mo + Taille limite des fichiers d'échange : 48.00 Go + + + -------------------------------------------------------------------------------- + ASTER 11.04.00 CONCEPT RESU CALCULE LE 23/11/2015 A 16:32:44 DE TYPE EVOL_NOLI + + + ======> + + ------> + CHAMP AUX NOEUDS DE NOM SYMBOLIQUE DEPL + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + NOEUD DX DY + N1 0.00000000000000E+00 0.00000000000000E+00 + N2 0.00000000000000E+00 0.00000000000000E+00 + N3 0.00000000000000E+00 0.00000000000000E+00 + N4 0.00000000000000E+00 0.00000000000000E+00 + N5 0.00000000000000E+00 0.00000000000000E+00 + N6 0.00000000000000E+00 0.00000000000000E+00 + + + ------> + CHAMP PAR ELEMENT AUX POINTS DE GAUSS DE NOM SYMBOLIQUE EPSI_ELGA + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + E1 EPXX EPYY EPZZ EPXY + 1 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 2 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 3 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 4 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + E2 EPXX EPYY EPZZ EPXY + 1 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 2 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 3 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 4 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + + + ------> + CHAMP PAR ELEMENT AUX NOEUDS DE NOM SYMBOLIQUE EPSI_ELNO + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + E1 EPXX EPYY EPZZ EPXY + N1 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N2 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N4 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N3 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + E2 EPXX EPYY EPZZ EPXY + N3 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N4 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N6 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N5 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + + + ------> + CHAMP PAR ELEMENT AUX POINTS DE GAUSS DE NOM SYMBOLIQUE SIEF_ELGA + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + E1 SIXX SIYY SIZZ SIXY + 1 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 2 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 3 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 4 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + E2 SIXX SIYY SIZZ SIXY + 1 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 2 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 3 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 4 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + + + ------> + CHAMP PAR ELEMENT AUX NOEUDS DE NOM SYMBOLIQUE SIGM_ELNO + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + E1 SIXX SIYY SIZZ SIXY + N1 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N2 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N4 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N3 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + E2 SIXX SIYY SIZZ SIXY + N3 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N4 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N6 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N5 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + + + ------> + CHAMP AUX NOEUDS DE NOM SYMBOLIQUE FORC_NODA + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + NOEUD DX DY + N1 0.00000000000000E+00 0.00000000000000E+00 + N2 0.00000000000000E+00 0.00000000000000E+00 + N3 0.00000000000000E+00 0.00000000000000E+00 + N4 0.00000000000000E+00 0.00000000000000E+00 + N5 0.00000000000000E+00 0.00000000000000E+00 + N6 0.00000000000000E+00 0.00000000000000E+00 + + + ------> + CHAMP PAR ELEMENT AUX POINTS DE GAUSS DE NOM SYMBOLIQUE EPSG_ELGA + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + E1 EPXX EPYY EPZZ EPXY + 1 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 2 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 3 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 4 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + E2 EPXX EPYY EPZZ EPXY + 1 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 2 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 3 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + 4 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + + + ------> + CHAMP PAR ELEMENT AUX NOEUDS DE NOM SYMBOLIQUE EPSG_ELNO + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + E1 EPXX EPYY EPZZ EPXY + N1 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N2 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N4 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N3 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + E2 EPXX EPYY EPZZ EPXY + N3 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N4 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N6 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + N5 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 + + + ------> + CHAMP PAR ELEMENT AUX POINTS DE GAUSS DE NOM SYMBOLIQUE VARI_ELGA + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + E1 VARI + 1 0.00000000000000E+00 + 2 0.00000000000000E+00 + 3 0.00000000000000E+00 + 4 0.00000000000000E+00 + E2 VARI + 1 0.00000000000000E+00 + 2 0.00000000000000E+00 + 3 0.00000000000000E+00 + 4 0.00000000000000E+00 + + + ------> + CARTE DE NOM SYMBOLIQUE COMPORTEMENT + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + + ------> + CHAMP AUX NOEUDS DE NOM SYMBOLIQUE ENEL_NOEU + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + NOEUD TOTALE + N1 0.00000000000000E+00 + N2 0.00000000000000E+00 + N3 0.00000000000000E+00 + N4 0.00000000000000E+00 + N5 0.00000000000000E+00 + N6 0.00000000000000E+00 + + + ------> + CHAMP PAR ELEMENT CONSTANT SUR L'ELEMENT DE NOM SYMBOLIQUE ENEL_ELEM + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + E1 TOTALE + 0.00000000000000E+00 + E2 TOTALE + 0.00000000000000E+00 + + + ------> + CHAMP PAR ELEMENT CONSTANT SUR L'ELEMENT DE NOM SYMBOLIQUE ETOT_ELEM + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + E1 TOTALE + 0.00000000000000E+00 + E2 TOTALE + 0.00000000000000E+00 + + + ------> + CHAMP AUX NOEUDS DE NOM SYMBOLIQUE ETOT_NOEU + NUMERO D'ORDRE: 0 INST: 0.00000000000000E+00 + NOEUD TOTALE + N1 0.00000000000000E+00 + N2 0.00000000000000E+00 + N3 0.00000000000000E+00 + N4 0.00000000000000E+00 + N5 0.00000000000000E+00 + N6 0.00000000000000E+00 + + + ======> + + ------> + CHAMP AUX NOEUDS DE NOM SYMBOLIQUE DEPL + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + NOEUD DX DY + N1 2.58861162337177E-27 3.23117426778526E-27 + N2 3.17431158889468E-02 -3.23117426778526E-27 + N3 0.00000000000000E+00 -1.38591518927826E-01 + N4 3.17431158889468E-02 -1.38591518927826E-01 + N5 -2.42338070083895E-27 -2.77183037855653E-01 + N6 3.17431158889468E-02 -2.77183037855653E-01 + + + ------> + CHAMP PAR ELEMENT AUX POINTS DE GAUSS DE NOM SYMBOLIQUE EPSI_ELGA + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + E1 EPXX EPYY EPZZ EPXY + 1 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 1.73472347597681E-18 + 2 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 0.00000000000000E+00 + 3 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 0.00000000000000E+00 + 4 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 6.93889390390723E-18 + E2 EPXX EPYY EPZZ EPXY + 1 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 3.46944695195361E-18 + 2 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 3.46944695195361E-18 + 3 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 0.00000000000000E+00 + 4 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 0.00000000000000E+00 + + + ------> + CHAMP PAR ELEMENT AUX NOEUDS DE NOM SYMBOLIQUE EPSI_ELNO + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + E1 EPXX EPYY EPZZ EPXY + N1 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 -2.32408877239647E-19 + N2 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 6.22737709701874E-20 + N4 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 -3.23703807471396E-18 + N3 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 1.20807905608675E-17 + E2 EPXX EPYY EPZZ EPXY + N3 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 4.73935267345112E-18 + N4 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 4.73935267345113E-18 + N6 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 -1.26990572149751E-18 + N5 3.17431158889468E-02 -1.38591518927826E-01 3.56161343462932E-02 -1.26990572149751E-18 + + + ------> + CHAMP PAR ELEMENT AUX POINTS DE GAUSS DE NOM SYMBOLIQUE SIEF_ELGA + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + E1 SIXX SIYY SIZZ SIXY + 1 -6.26575285800744E-33 -9.39413593841838E+01 -8.24861551021306E-15 1.65027734475050E-15 + 2 8.24861551021306E-15 -9.39413593841838E+01 -8.24861551021306E-15 3.44319560568744E-16 + 3 -3.30517155972965E-32 -9.39413593841838E+01 0.00000000000000E+00 2.42991616671316E-15 + 4 -8.24861551021306E-15 -9.39413593841838E+01 -8.24861551021306E-15 3.73587395089492E-15 + E2 SIXX SIYY SIZZ SIXY + 1 8.24861551021306E-15 -9.39413593841838E+01 -8.24861551021306E-15 3.95353358111226E-15 + 2 -3.30517155972155E-32 -9.39413593841838E+01 0.00000000000000E+00 2.42991616658465E-15 + 3 -8.24861551021306E-15 -9.39413593841838E+01 0.00000000000000E+00 -1.74127704570419E-15 + 4 -8.24861551021306E-15 -9.39413593841838E+01 -8.24861551021306E-15 -2.17659631176580E-16 + + + ------> + CHAMP PAR ELEMENT AUX NOEUDS DE NOM SYMBOLIQUE SIGM_ELNO + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + E1 SIXX SIYY SIZZ SIXY + N1 2.36658271566304E-30 -9.39413593841838E+01 -7.14351057789485E-15 1.36490973013559E-15 + N2 1.42870211557897E-14 -9.39413593841838E+01 -1.23729232653196E-14 -8.97075504607280E-16 + N4 -7.88860905221012E-31 -9.39413593841838E+01 7.14351057789484E-15 2.71528378132807E-15 + N3 -1.42870211557897E-14 -9.39413593841838E+01 -1.23729232653196E-14 4.97726901607094E-15 + E2 SIXX SIYY SIZZ SIXY + N3 1.84113289108962E-14 -9.39413593841838E+01 -1.12678183330014E-14 6.03797894026866E-15 + N4 -1.10510493231821E-15 -9.39413593841838E+01 3.01920282278832E-15 3.39899616701011E-15 + N6 -1.01627134006832E-14 -9.39413593841838E+01 3.01920282278832E-15 -3.82572240486059E-15 + N5 -1.53921260881079E-14 -9.39413593841837E+01 -1.12678183330014E-14 -1.18673963160204E-15 + + + ------> + CHAMP AUX NOEUDS DE NOM SYMBOLIQUE FORC_NODA + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + NOEUD DX DY + N1 -2.27249029899045E-15 4.84616754209406E+01 + N2 1.67634515516756E-16 4.84616754209406E+01 + N3 8.24782431065233E-16 -2.13162820728030E-14 + N4 1.38833126914655E-16 0.00000000000000E+00 + N5 3.22406470732547E-15 -4.84616754209406E+01 + N6 -2.08282448183166E-15 -4.84616754209406E+01 + + + ------> + CHAMP PAR ELEMENT AUX POINTS DE GAUSS DE NOM SYMBOLIQUE EPSG_ELGA + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + E1 EPXX EPYY EPZZ EPXY + 1 3.22469285921164E-02 -1.28987714368465E-01 3.22469285921163E-02 1.22635478889943E-18 + 2 3.22469285921164E-02 -1.28987714368465E-01 3.22469285921163E-02 -5.90966979577186E-19 + 3 3.22469285921163E-02 -1.28987714368465E-01 3.22469285921163E-02 -1.07180290248795E-18 + 4 3.22469285921163E-02 -1.28987714368465E-01 3.22469285921163E-02 5.94968929391909E-18 + E2 EPXX EPYY EPZZ EPXY + 1 3.22469285921163E-02 -1.28987714368465E-01 3.22469285921163E-02 2.49400872401844E-18 + 2 3.22469285921163E-02 -1.28987714368465E-01 3.22469285921163E-02 2.39764404945753E-18 + 3 3.22469285921163E-02 -1.28987714368465E-01 3.22469285921163E-02 -1.10131056674545E-19 + 4 3.22469285921163E-02 -1.28987714368465E-01 3.22469285921163E-02 -1.37663821136368E-20 + + + ------> + CHAMP PAR ELEMENT AUX NOEUDS DE NOM SYMBOLIQUE EPSG_ELNO + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + E1 EPXX EPYY EPZZ EPXY + N1 3.22469285921164E-02 -1.28987714368465E-01 3.22469285921163E-02 -5.34546328115394E-19 + N2 3.22469285921164E-02 -1.28987714368465E-01 3.22469285921163E-02 -3.82928119133670E-19 + N4 3.22469285921163E-02 -1.28987714368465E-01 3.22469285921163E-02 -4.51507221340354E-18 + N3 3.22469285921163E-02 -1.28987714368465E-01 3.22469285921164E-02 1.09458208614060E-17 + E2 EPXX EPYY EPZZ EPXY + N3 3.22469285921163E-02 -1.28987714368465E-01 3.22469285921163E-02 3.44719003875771E-18 + N4 3.22469285921163E-02 -1.28987714368465E-01 3.22469285921163E-02 3.28028152636338E-18 + N6 3.22469285921163E-02 -1.28987714368465E-01 3.22469285921163E-02 -1.06331237141381E-18 + N5 3.22469285921163E-02 -1.28987714368465E-01 3.22469285921163E-02 -8.96403859019481E-19 + + + ------> + CHAMP PAR ELEMENT AUX POINTS DE GAUSS DE NOM SYMBOLIQUE VARI_ELGA + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + E1 VARI + 1 0.00000000000000E+00 + 2 0.00000000000000E+00 + 3 0.00000000000000E+00 + 4 0.00000000000000E+00 + E2 VARI + 1 0.00000000000000E+00 + 2 0.00000000000000E+00 + 3 0.00000000000000E+00 + 4 0.00000000000000E+00 + + + ------> + CARTE DE NOM SYMBOLIQUE COMPORTEMENT + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + + ------> + CHAMP AUX NOEUDS DE NOM SYMBOLIQUE ENEL_NOEU + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + NOEUD TOTALE + N1 4.90276611274910E+00 + N2 4.90276611274910E+00 + N3 4.90276611274910E+00 + N4 4.90276611274910E+00 + N5 4.90276611274909E+00 + N6 4.90276611274910E+00 + + + ------> + CHAMP PAR ELEMENT CONSTANT SUR L'ELEMENT DE NOM SYMBOLIQUE ENEL_ELEM + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + E1 TOTALE + 4.90276611274910E+00 + E2 TOTALE + 4.90276611274910E+00 + + + ------> + CHAMP PAR ELEMENT CONSTANT SUR L'ELEMENT DE NOM SYMBOLIQUE ETOT_ELEM + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + E1 TOTALE + 6.50973784359943E+00 + E2 TOTALE + 6.50973784359942E+00 + + + ------> + CHAMP AUX NOEUDS DE NOM SYMBOLIQUE ETOT_NOEU + NUMERO D'ORDRE: 1 INST: 1.00000000000000E+00 + NOEUD TOTALE + N1 6.50973784359942E+00 + N2 6.50973784359942E+00 + N3 6.50973784359942E+00 + N4 6.50973784359942E+00 + N5 6.50973784359942E+00 + N6 6.50973784359943E+00 + + FERMETURE DE LA BASE "GLOBALE" EFFECTUEE. + + Arrêt normal dans "FIN". + ARRET NORMAL DANS "FIN" PAR APPEL A "JEFINI". + + MEMOIRE JEVEUX MINIMALE REQUISE POUR L'EXECUTION : 21.14 Mo + MEMOIRE JEVEUX OPTIMALE REQUISE POUR L'EXECUTION : 27.32 Mo + MAXIMUM DE MEMOIRE UTILISEE PAR LE PROCESSUS LORS DE L'EXECUTION : 226.66 Mo + + ******************************************************************************** + * COMMAND : USER : SYSTEM : USER+SYS : ELAPSED * + ******************************************************************************** + * init (jdc) : 0.16 : 0.02 : 0.18 : 0.17 * + * . compile : 0.00 : 0.00 : 0.00 : 0.00 * + * . exec_compile : 0.05 : 0.01 : 0.06 : 0.06 * + * . report : 0.01 : 0.00 : 0.01 : 0.00 * + * . build : 0.00 : 0.00 : 0.00 : 0.00 * + * DEBUT : 0.00 : 0.02 : 0.02 : 0.03 * + * LIRE_MAILLAGE : 0.01 : 0.00 : 0.01 : 0.00 * + * AFFE_MODELE : 0.00 : 0.00 : 0.00 : 0.00 * + * DEFI_MATERIAU : 0.00 : 0.00 : 0.00 : 0.00 * + * AFFE_MATERIAU : 0.01 : 0.00 : 0.01 : 0.00 * + * AFFE_CHAR_MECA : 0.00 : 0.00 : 0.00 : 0.01 * + * AFFE_CHAR_MECA : 0.00 : 0.00 : 0.00 : 0.00 * + * DEFI_LIST_REEL : 0.00 : 0.00 : 0.00 : 0.00 * + * DEFI_FONCTION : 0.00 : 0.00 : 0.00 : 0.00 * + * STAT_NON_LINE : 0.07 : 0.00 : 0.07 : 0.06 * + * CALC_CHAMP : 0.04 : 0.00 : 0.04 : 0.04 * + * IMPR_RESU : 0.01 : 0.00 : 0.01 : 0.02 * + * FIN : 0.01 : 0.01 : 0.02 : 0.02 * + * . part Superviseur : 0.18 : 0.04 : 0.22 : 0.22 * + * . part Fortran : 0.14 : 0.01 : 0.15 : 0.15 * + ******************************************************************************** + * TOTAL_JOB : 0.32 : 0.05 : 0.37 : 0.37 * + ******************************************************************************** +