diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst index e0bc570..84281cc 100644 --- a/CONTRIBUTING.rst +++ b/CONTRIBUTING.rst @@ -78,15 +78,40 @@ Testing is made easy by using our `Makefile`. From there one founds convenient functions `make test`, `make test_file` and `make test_function` to make testing more rapid. -Git issues ----------- +One possible workflow: See +`this post _`: + +.. code:: bash + + while sleep_until_modified.sh src/equations.jl; do clear; make test_file FILE=test/test_equations.jl; done + +or something similar. Every time file `src/equations.jl` is changed, tests from +file `test/test_equations.jl` are run. + +Git spesific things +------------------- Have done local changes, want to get latest updates and get "Cannot pull with rebase: You have unstaged changes. Please commit or stash them." http://stackoverflow.com/questions/23517464/error-cannot-pull-with-rebase-you-have-unstaged-changes -Solution: git stash + git pull + git pop +.. code-block:: bash + + git stash + git pull + git stash pop + +Local modifictions done, want to get clean file from repo + +.. code-block:: bash + + git checkout -- + +Having untracked files not wanting to commit at all? Put them to `.gitignore` if +they are somehow generally unnecessary things. Or if they are "personal" (have +made e.g. notebooks not wanting to commit), add them to your private ignore rules +`.git/info/exclude` so that they don't show as untracked for you. Use of UTF-8 characters in program code --------------------------------------- @@ -140,5 +165,6 @@ Documentation ------------- We use restructured text to document this project. Information how to write rst format is described `here `_. See issue -`#49 `_. +`#49 `_. Keep line width max +80 characters. diff --git a/docs/build_notebooks.jl b/docs/build_notebooks.jl index 65812fe..6d7c427 100644 --- a/docs/build_notebooks.jl +++ b/docs/build_notebooks.jl @@ -63,24 +63,34 @@ function run_notebooks() end runtime = toc() bn = "tutorials/$(ipynb[1:end-6])" + #try + # run(`ipython nbconvert tutorials/$ipynb --to rst --output=$bn`) + #catch error + # warn("unable to convert notebook to rst format") + # Base.showerror(Base.STDOUT, error) + #end + try - run(`ipython nbconvert tutorials/$ipynb --to rst --output=$bn`) + run(`ipython nbconvert tutorials/$ipynb --to html --output=$bn`) catch error - warn("unable to convert notebook to rst format") + warn("unable to convert notebook to html format") Base.showerror(Base.STDOUT, error) end + try run(`ipython nbconvert tutorials/$ipynb --to latex --output=$bn`) catch error warn("unable to convert notebook to tex format") Base.showerror(Base.STDOUT, error) end + try - run(`lualatex $bn.tex`) + run(`lualatex --output-directory=tutorials $bn.tex`) catch error warn("unable to convert notebook from tex to pdf") Base.showerror(Base.STDOUT, error) end + data = Dict("author" => "unknown", "status" => status, "runtime" => runtime, "filename" => ipynb, "last_run" => time(), "description"=>"") res = parse_rst("$bn.rst") diff --git a/src/dirichlet.jl b/src/dirichlet.jl index ba9ae5a..68f50df 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -10,12 +10,29 @@ abstract DirichletEquation <: Equation type DirichletProblem <: BoundaryProblem unknown_field_name :: ASCIIString unknown_field_dimension :: Int - equations :: Array{DirichletEquation, 1} - element_mapping :: Dict{DataType, DataType} + equations :: Vector{DirichletEquation} + element_mapping :: Dict{Element, Equation} field_value :: Function end -function DirichletProblem(dimension::Int, field_value::Function=(X)->[0.0,0.0,0.0]) +""" Initialize new Dirichlet boundary condition. + +Parameters +---------- +dimension + dimension of unknown field +field_value + boundary function + +Examples +-------- + +Create u(X) = 0.0 boundary condition for three-dimensional elasticity problem: + +>>> u(X) = [0.0, 0.0, 0.0] +>>> bc = DirichletProblem(3, u) +""" +function DirichletProblem(dimension::Int=1, field_value::Function=(X)->[0.0,0.0,0.0]) element_mapping = nothing if dimension == 1 element_mapping = Dict( @@ -28,12 +45,10 @@ end """ Dirichlet boundary condition element for 2 node line segment """ type DBC2D2 <: DirichletEquation element :: Seg2 - integration_points :: Array{IntegrationPoint, 1} + integration_points :: Vector{IntegrationPoint} end function DBC2D2(element::Seg2) - integration_points = [ - IntegrationPoint([-sqrt(1/3)], 1.0), - IntegrationPoint([+sqrt(1/3)], 1.0)] + integration_points = default_integration_points(element) if !haskey(element, "reaction force") element["reaction force"] = zeros(1, 2) end @@ -41,9 +56,7 @@ function DBC2D2(element::Seg2) end Base.size(equation::DBC2D2) = (1, 2) -function calculate_local_assembly!(assembly::LocalAssembly, equation::DirichletEquation, - unknown_field_name::ASCIIString, time::Number=Inf, - problem=nothing) +function calculate_local_assembly!(assembly::LocalAssembly, equation::DirichletEquation, unknown_field_name::ASCIIString, time::Number=0.0, problem=nothing) initialize_local_assembly!(assembly, equation) element = get_element(equation) basis = get_basis(element) diff --git a/src/field_utils.jl b/src/field_utils.jl deleted file mode 100644 index 301dcdf..0000000 --- a/src/field_utils.jl +++ /dev/null @@ -1,140 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -# There are here for now until I figure a better place for them. - - -""" -Fit field s.t. || ∫ (Nᵢ(ξ)αᵢ - f(el, ξ)) dS || -> min! - -Parameters ----------- -f::Function - Needs to take (el::Element, xi::Vector) as argument -fixed_coeffs::Int[] - These coefficients are not changed during fitting -> constrained optimizatio -""" -function fit_field!(el::Element, field, f, fixed_coeffs=Int[]) - w = [ - 128/225, - (332+13*sqrt(70))/900, - (332+13*sqrt(70))/900, - (332-13*sqrt(70))/900, - (332-13*sqrt(70))/900] - xi = Vector[ - [0.0], - [ 1/3*sqrt(5 - 2*sqrt(10/7))], - [-1/3*sqrt(5 - 2*sqrt(10/7))], - [ 1/3*sqrt(5 + 2*sqrt(10/7))], - [-1/3*sqrt(5 + 2*sqrt(10/7))]] - n = get_number_of_basis_functions(el) - fld = get_field(el, field) - nfld = length(fld[1]) - #Logging.debug("dim of field $field: $nfld") - - M = zeros(n, n) - b = zeros(n, nfld) - for i=1:length(w) - detJ = get_detJ(el, xi[i]) - N = get_basis(el, xi[i]) - M += w[i]*N*N'*detJ - fi = f(el, xi[i]) - for j=1:nfld - b[:, j] += w[i]*N*fi[j]*detJ - end - end - - coeffs = zeros(n) - for j=1:nfld - for k=1:n - coeffs[k] = fld[k][j] - end - if length(fixed_coeffs) != 0 - # constrained problem, some coefficients are fixed - N = Int[] # rest of coeffs - S = Int[] # fixed coeffs - for i = 1:n - if i in fixed_coeffs - push!(S, i) - else - push!(N, i) - end - end - lhs = M[N,N] - rhs = b[N,j] - M[N,S]*coeffs[S] - coeffs[N] = lhs \ rhs - else - coeffs[:] = M \ b[:,j] - end - for k=1:n - fld[k][j] = coeffs[k] - end - end - set_field(el, field, fld) - return -end - - -""" -Fit field s.t. || ∫ ∂/∂ξ(∑Nᵢ(ξ)αᵢ)f(el, ξ) dS || -> min! -""" -function fit_derivative_field!(el::Element, field, f, fixed_coeffs=Int[]) - w = [ - 128/225, - (332+13*sqrt(70))/900, - (332+13*sqrt(70))/900, - (332-13*sqrt(70))/900, - (332-13*sqrt(70))/900] - xi = Vector[ - [0.0], - [ 1/3*sqrt(5 - 2*sqrt(10/7))], - [-1/3*sqrt(5 - 2*sqrt(10/7))], - [ 1/3*sqrt(5 + 2*sqrt(10/7))], - [-1/3*sqrt(5 + 2*sqrt(10/7))]] - n = get_number_of_basis_functions(el) - fld = get_field(el, field) - nfld = length(fld[1]) - #Logging.debug("dim of field $field: $nfld") - - M = zeros(n, n) - b = zeros(n, nfld) - for i=1:length(w) - detJ = get_detJ(el, xi[i]) - dNdxi = get_dbasisdxi(el, xi[i]) - dNdX = dNdxi / detJ - M += w[i]*dNdX*dNdX'*detJ - fi = f(el, xi[i]) - for j=1:nfld - b[:, j] += w[i]*dNdX*fi[j]*detJ - end - end - - coeffs = zeros(n) - for j=1:nfld - for k=1:n - coeffs[k] = fld[k][j] - end - if length(fixed_coeffs) != 0 - #Logging.info("constrained problem, some coefficients are fixed") - N = Int[] # rest of coeffs - S = Int[] # fixed coeffs - for i = 1:n - if i in fixed_coeffs - push!(S, i) - else - push!(N, i) - end - end - lhs = M[N,N] - rhs = b[N,j] - M[N,S]*coeffs[S] - coeffs[N] = lhs \ rhs - else - coeffs[:] = M \ b[:,j] - end - for k=1:n - fld[k][j] = coeffs[k] - end - end - set_field(el, field, fld) - return -end diff --git a/src/mortar.jl b/src/mortar.jl index e7e48f8..871080a 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -1,101 +1,80 @@ -""" -calculate "local" normals in elements, in a way that -n = Nᵢnᵢ gives some reasonable results for ξ ∈ [-1, 1] -""" -function calculate_normals!(el::Element, t, field_name=symbol("normals")) - new_field!(el, field_name, Vector) - for xi in Vector[[-1.0], [1.0]] - t = dinterpolate(el, :Geometry, xi) - n = [0 -1; 1 0]*t - n /= norm(n) - push_field!(el, field_name, n) - end -end +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +# Mortar projection integration + +abstract MortarEquation <: Equation """ -Alter normal field such that normals of adjacent elements are averaged. +Parameters +---------- +node_csys + coordinate system in node, normal + tangent + "binormal" +element_pairs + m x s matrix of boolean values, indicating elements sharing + common surface. s is number of slave elements and m is number + of master elements. """ -function average_normals!(elements, normal_field=symbol("normals")) - d = Dict() - for el in elements - c = get_connectivity(el) - n = get_field(el, normal_field) - for (ci, ni) in zip(c, n) - d[ci] = haskey(d, ci) ? d[ci] + ni : ni +type MortarProblem <: BoundaryProblem + unknown_field_name :: ASCIIString + unknown_field_dimension :: Int + equations :: Vector{MortarEquation} + element_mapping :: Dict{Element, MortarEquation} + master_elements :: Vector{Element} # mortar surface + node_csys :: Dict{Int, Matrix{Float64}} + element_pairs :: Matrix{Bool} +end + +function MortarProblem(dimension::Int=1, equations=[], master_elements=[]) + element_mapping = Dict( + Seg2 => MBC2D2, + ) + MortarProblem("reaction force", dimension, equations, element_mapping, master_elements, Dict(), zeros(0,0)) +end + +""" Mortar boundary condition element for 2-dimensional problem, 2 node line segment. """ +type MBC2D2 <: MortarEquation + element :: Seg2 # == non-mortar surface element + integration_points :: Vector{IntegrationPoint} +end +function MBC2D2(element::Seg2) + integration_points = default_integration_points(element) + if !haskey(element, "reaction force") + element["reaction force"] = zeros(1, 2) + end + MBC2D2(element, integration_points) +end +Base.size(equation::MBC2D2) = (1, 2) + +function find_master_elements(slave_element, problem) + # find slave element "position" in element pairs matrix + all_elements = map((equation) -> get_element(equation), problem.equations) + seid = findfirst(slave_element, all_elements) + info("slave element id = $seid") + # find master element "positions" in element pairs matrix + meids = find(problem.element_pairs[:, seid]) + info("master element ids = $meids") + # master elements + master_elements = problem.master_elements[meids] + return master_elements +end + +function calculate_local_assembly!(assembly::LocalAssembly, equation::MortarEquation, unknown_field_name::ASCIIString, time::Number=0.0, problem=nothing) + # slave element = non-mortar element where integration happens + # master element = mortar element projected to non-mortar side + isa(problem, Void) && error("Cannot create projection without problem") + initialize_local_assembly!(assembly, equation) + slave_element = get_element(equation) + basis = get_basis(slave_element) + detJ = det(basis) + master_elements = find_master_elements(equation, problem) + for master_element in master_elements + for ip in get_integration_points(slave_element) + mortar_basis = 0 # ... + assembly.stiffness_matrix += w*basis'*basis + assembly.force_vector += w*N'*gn end end - for (ci, ni) in d - d[ci] /= norm(d[ci]) - end - for el in elements - c = get_connectivity(el) - new_normals = [d[ci] for ci in c] - set_field(el, normal_field, new_normals) - end end -""" Find projection from slave nodes to master element. """ -function calc_projection_slave_nodes_to_master_element(sel, mel) - X1 = get_field(sel, :Geometry) - N1 = get_field(sel, :Normals) - X2(xi) = interpolate(mel, :Geometry, xi) - dX2(xi) = dinterpolate(mel, :Geometry, xi) - R(xi, k) = det([X2(xi) - X1[k] N1[k]]') - dR(xi, k) = det([dX2(xi) N1[k]]') - xi2 = Vector[[0.0], [0.0]] - for k=1:2 - xi = xi2[k] - for i=1:3 - dxi = -R(xi, k)/dR(xi, k) - xi += dxi - if abs(dxi) < 1.0e-9 - break - end - end - xi2[k] = xi - end - clamp!(xi2, -1, 1) - return xi2 -end - -""" Find projection from master nodes to slave element. """ -function calc_projection_master_nodes_to_slave_element(sel, mel) - X1(xi) = interpolate(sel, :Geometry, xi) - dX1(xi) = dinterpolate(sel, :Geometry, xi) - N1(xi) = interpolate(sel, :Normals, xi) - dN1(xi) = dinterpolate(sel, :Normals, xi) - X2 = get_field(mel, :Geometry) - R(xi, k) = det([X1(xi) - X2[k] N1(xi)]') - dR(xi, k) = det([dX1(xi) N1(xi)]') + det([X1(xi) - X2[k] dN1(xi)]') - xi1 = Vector[[0.0], [0.0]] - for k=1:2 - xi = xi1[k] - for i=1:3 - dxi = -R(xi, k)/dR(xi, k) - xi += dxi - if abs(dxi) < 1.0e-9 - break - end - end - xi1[k] = xi - end - clamp!(xi1, -1, 1) - return xi1 -end - -function has_projection(sel, mel) - xi1 = calc_projection_master_nodes_to_slave_element(sel, mel) - l = abs(xi1[2]-xi1[1])[1] - return l > 1.0e-9 -end - -""" -Calculate projection between 1d boundary elements -""" -function calc_projection(sel, mel) - xi1 = calc_projection_master_nodes_to_slave_element(sel, mel) - xi2 = calc_projection_slave_nodes_to_master_element(sel, mel) - return xi1, xi2 -end - diff --git a/test/test_mortar.jl b/test/test_mortar.jl new file mode 100644 index 0000000..dab6820 --- /dev/null +++ b/test/test_mortar.jl @@ -0,0 +1,83 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +module MortarTests + +using JuliaFEM +using JuliaFEM.Test + +function test_calc_flat_2d_assembly() + # this is hand calculated and given example in my thesis + N = Vector[ + [0.0, 2.0], [1.0, 2.0], [2.0, 2.0], + [0.0, 0.0], [1.0, 0.0], [2.0, 0.0], + [0.0, 1.0], [5/4, 1.0], [2.0, 1.0], + [0.0, 1.0], [3/4, 1.0], [2.0, 1.0]] + + slave1 = Seg2([10, 11]) + slave1["geometry"] = Vector[N10, N11] + + slave2 = Seg2([11, 12]) + slave2["geometry"] = Vector[N11, N12] + + master1 = Seg2([7, 8]) + master1["geometry"] = Vector[N7, N8] + + master2 = Seg2([8, 9]) + master2["geometry"] = Vector[N8, N9] + + problem = MortarProblem() + push!(problem, slave1) + push!(problem, slave2) + push!(problem.master_elements, master1) + push!(problem.master_elements, master2) + + rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)] + # should be n = [0 -1]' and t = [1 0]' + @test isapprox(rotation_matrix(-phi/2), [[0 -1]' [1 0]']) + + problem.node_csys = Dict( + 10 => rotation_matrix(-phi/2), + 11 => rotation_matrix(-phi/2), + 12 => rotation_matrix(-phi/2)) + + # first index = master element id + # second index = slave element id + problem.element_pairs = zeros(2, 2) + # first slave element connects to master element 1 + problem.element_pairs[1, 1] = true + # second slave element connects to master element 1 + problem.element_pairs[1, 2] = true + # second slave element connects to master element 2 + problem.element_pairs[2, 2] = true + + B_expected = zeros(12, 9) + + S1 = [10, 11] + M1 = [7, 8] + B_expected[S1,S1] += [1/4 1/8; 1/8 1/4] + B_expected[S1,M1] += [3/10 3/40; 9/40 3/20] + + la = initialize_local_assembly(problem) + calculate_local_assembly!(la, problem.equations[1], "reaction force", 0.0, problem=problem) + B = full(la.lhs) + @test isapprox(B, B_expected) + + fill!(B_expected, 0.0) + + S2 = [11, 12] + M2 = [7, 8] + B_expected[S2,S2] += [49/150 11/150; 11/150 2/75] + B_expected[S2,M2] += [13/150 47/150; 1/75 13/150] + S3 = [11, 12] + M3 = [8, 9] + B_expected[S3,S3] += [9/100 27/200; 27/200 39/100] + B_expected[S3,M3] += [3/20 3/40; 9/40 3/10] + + la = initialize_local_assembly(problem) + calculate_local_assembly!(la, problem.equations[1], "reaction force", 0.0, problem=problem) + B = full(la.lhs) + @test isapprox(B, B_expected) +end + +end