From 80fe76e282facf37c37bf82341e5ed9c0144a7e4 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 2 Aug 2017 19:50:53 +0300 Subject: [PATCH 1/2] Change calculate_mortar_assembly return arguments - function returns now s, m, D, M and user can find P however wants - function is also now exported to global namespace --- docs/src/theory.md | 4 ++-- src/Mortar2D.jl | 1 + src/calculate_mortar_assembly.jl | 26 +++++++++++++------------- test/test_calculate_mortar_assembly.jl | 11 +++++++---- 4 files changed, 23 insertions(+), 19 deletions(-) diff --git a/docs/src/theory.md b/docs/src/theory.md index 5df1d8e..3b73a53 100644 --- a/docs/src/theory.md +++ b/docs/src/theory.md @@ -151,8 +151,8 @@ and ``\boldsymbol{M}`` or ``\boldsymbol{P}=\boldsymbol{D}^{-1}\boldsymbol{M}``. Calculation projection matrix ``\boldsymbol{P}`` is implemented as function `calculate_mortar_assembly`: ```@example 0 -P = calculate_mortar_assembly(elements, element_types, coords, - slave_element_ids, master_element_ids) +s, m, D, M = calculate_mortar_assembly(elements, element_types, coords, + slave_element_ids, master_element_ids) ``` This last command combines everything above to single command to calculate diff --git a/src/Mortar2D.jl b/src/Mortar2D.jl index 98c5ffa..8138d39 100644 --- a/src/Mortar2D.jl +++ b/src/Mortar2D.jl @@ -7,4 +7,5 @@ include("calculate_projections.jl") include("calculate_segments.jl") include("calculate_mortar_matrices.jl") include("calculate_mortar_assembly.jl") +export calculate_mortar_assembly end diff --git a/src/calculate_mortar_assembly.jl b/src/calculate_mortar_assembly.jl index 029b1dd..69d52a2 100644 --- a/src/calculate_mortar_assembly.jl +++ b/src/calculate_mortar_assembly.jl @@ -8,14 +8,17 @@ slave_element_ids::Vector{Int}, master_element_ids::Vector{Int}) -Given data, calculate projection matrix `P`. This is the main function of -package. +Given data, calculate projection matrices `D` and `M`. This is the main +function of package. Relation between matrices is ``D u_s = M u_m``, where +``u_s`` is slave nodes and ``u_m`` master nodes. -Matrix ``P`` is defined as a projection between slave and master surface, -i.e. -```math - D u_s = M u_m \\Rightarrow u_s = D^{-1} M u_m = P u_m. +# Example + +```julia +s, m, D, M = calculate_mortar_assembly(elements, element_types, coords, + slave_element_ids, master_element_ids) ``` + """ function calculate_mortar_assembly(elements::Dict{Int, Vector{Int}}, element_types::Dict{Int, Symbol}, @@ -65,11 +68,8 @@ function calculate_mortar_assembly(elements::Dict{Int, Vector{Int}}, end end - D = sparse(D_I, D_J, D_V) - Df = cholfact(1/2*(D+D')[S,S]) - - M_ = sparse(M_I, M_J, M_V) - P = Df \ M_[S,M] - SparseArrays.droptol!(P, 1.0e-12) - return P + # return global matrices + slave and master dofs + Dg = sparse(D_I, D_J, D_V) + Mg = sparse(M_I, M_J, M_V) + return S, M, Dg, Mg end diff --git a/test/test_calculate_mortar_assembly.jl b/test/test_calculate_mortar_assembly.jl index dd2ef97..ef2fedc 100644 --- a/test/test_calculate_mortar_assembly.jl +++ b/test/test_calculate_mortar_assembly.jl @@ -19,18 +19,21 @@ using Mortar2D: calculate_mortar_assembly elements = merge(Es, Em) coords = merge(Xs, Xm) - P = calculate_mortar_assembly(elements, element_types, coords, - slave_element_ids, master_element_ids) + s, m, D, M = calculate_mortar_assembly(elements, + element_types, + coords, + slave_element_ids, + master_element_ids) # test case 1, constant pressure load um = ones(nm) - us = P*um + us = D[s,s] \ (M[s,m]*um) us_expected = ones(ns) @test isapprox(us, us_expected) # test case 2, linear load um = linspace(0, 1, nm) - us = P*um + us = D[s,s] \ (M[s,m]*um) us_expected = linspace(0, 1, ns) @test isapprox(us, us_expected) end From cb54d40ae40cf359e10da8b5b78ced4956d73b47 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 2 Aug 2017 19:51:55 +0300 Subject: [PATCH 2/2] Update README * add usage example --- README.md | 49 +++++++++++++++++- docs/src/figs/poisson_problem_discretized.png | Bin 0 -> 4596 bytes 2 files changed, 48 insertions(+), 1 deletion(-) create mode 100644 docs/src/figs/poisson_problem_discretized.png diff --git a/README.md b/README.md index a33fda8..99f340d 100644 --- a/README.md +++ b/README.md @@ -25,4 +25,51 @@ Testing package can be done using `Pkg.test`, i.e. julia> Pkg.test("Mortar2D") ``` -Probably the easiest way to test the functionality of package is to use [JuliaBox](https://juliabox.com/). +Probably the easiest way to test the functionality of package is to +use [JuliaBox](https://juliabox.com/). + +## Usage example + +Let us calculate projection matrices D and M for the following problem: + +![](docs/src/figs/poisson_problem_discretized.png) + +Problem setup: + +```julia +Xs = Dict(1 => [0.0, 1.0], 2 => [5/4, 1.0], 3 => [2.0, 1.0]) +Xm = Dict(4 => [0.0, 1.0], 5 => [1.0, 1.0], 6 => [2.0, 1.0]) +coords = merge(Xm , Xs) +Es = Dict(1 => [1, 2], 2 => [2, 3]) +Em = Dict(3 => [4, 5], 4 => [5, 6]) +elements = merge(Es, Em) +element_types = Dict(1 => :Seg2, 2 => :Seg2, 3 => :Seg2, 4 => :Seg2) +slave_element_ids = [1, 2] +master_element_ids = [3, 4] +``` + +Calculate projection matrices D and M + +```julia +s, m, D, M = calculate_mortar_assembly( + elements, element_types, coords, + slave_element_ids, master_element_ids) +``` + +According to theory, the interface should transfer constant without any +error. Let's test that: + +```julia +u_m = ones(3) +u_s = D[s,s] \ (M[s,m]*um) + +# output + +3-element Array{Float64,1}: + 1.0 + 1.0 + 1.0 +``` + +The rest of the story can be read from the [documentation](https://juliafem.github.io/Mortar2D.jl/latest/). +There's also brief review to the theory behind non-conforming finite element meshes. diff --git a/docs/src/figs/poisson_problem_discretized.png b/docs/src/figs/poisson_problem_discretized.png new file mode 100644 index 0000000000000000000000000000000000000000..8d3d528b1c24d93057d259f4e349417fde09300a GIT binary patch literal 4596 zcmeHLS5y`(^H{bg9tXZ?>TWj|0JxSMW%#ZVm^8x^X z4m-{4+^ZN1ctm4wn?vcL+rlm|>(S)f=|y*$<#lGyasaK)(g z%LxlqsaFwxe9CjVg?;UfvfoOi`mgW{@tym}%Xqd@n&o;q2k3^pwyio~+bb~vMFQll zjLgt=-1h9Zp`?y0GO+G-)6BgFIhC;#Qr0{?9C6E9 z#!3;C{yxUhT80Y-Ns@doL_L6 zbKRYB1hRN!D@UsS`OlE%3E{FNF_vb*hzY7l(nO^PH^+A0^O- zS5hx)wl_v_R8&M?3h_R_ClkB+b?t6d4R6KF_hoI+M_T#&8SFcQZsxyJnwI^;^+hhO;;F8OOZSy$DofJRGZzpEZC!q2s4+M^x+ zez-sFrBh<6_M(a}Xoz>NO`gNCN#_Oi{zW_5t%rrH(?=&tUA_o4s(B*ONH1EsK z@Z@A{@QUy*cfp?CYwuxM_3dQ;0>jbASeNm;x8baPy$dsY$lQ9p=YH|k7?Q^n>lGnApTX~e%_d1kd|m0-I>SMWz=0O{6w$u)t@nB#O*$XlHrtS z=JzY?O!NM;;F4;m&`-M_WrlafLON{@Wz|qVTC)foX2Cdf>>7a+-*|E6xXFq}5!Jfg z%%lGIdc;|@9A0#a`0a9_Yi{iGTGO3TrgCiBQhzJBghxQIygOED?(~+9mnwCAjq+Obs{xf+gaCUNCWaoA5crG>+xGSD*zkV%* z{!*#qKB$w{fxR--s>JLsP<>b!{y4R9(_(rezOP~Z3G}vhIAb3)g~bpwsZ^zsqT~6W zJ>n**K0h-=VUyEz6SoN436Z}Do(RpITex-5Y5e!o9}GuhQna)NM0T^&f3vI7GzS5;lqi`*HQDr7O@nKRSA zls9kQBun+kw!VfhJ2fa&Z} zm!2~Z@_Riesw4=H?@hOfu+s~P&Fp^(+m5V^eXTeKKnm`LIo@vo&ij z>eoRqhi=2Uwh5*Pk0!`LVtFeXKK+Vb7M)y8dROk-y2o?)>_#vis<#@~TiOXz@w_?D zOPN@oahP;YfsIr{_IyQVk2o%XPzB;EW&y)OIH{(JA4Vvu*8UJD4}4}Zz?@OSyft^2 zfl1)H!$r7nbC49zaHP#wt>zyLM2xzAPVR{T{2vVJfz2s?MA=Qh)Si>_-qHh`Of;mc z7$hA?ktQS_?T#&|<8D|p_&D>zW|WmKm5irsObA?CsxWf+nR~$DFxOfL1GrJ#EU*Z< zdrTjvH`JvizRI@j>u!2@#cn&5yEw{p+ z2_q`FFqT%CwWq*x3SJ}Ogv{6R|C{)4YYAi#Gd3_FIG!{y;fjNzI^C}%Eo?L4VwiH% zh3Asf(MS3$h7#;)VDYAav_Z{pC7Ld83uDYui8+Fpo%&IT6*z4v9a@b`66H-Gcy~YJ zf!?}xYusg+r98ij;rL$eHq%&Sn(F?mOZ61CYsxCg{^m_h1>ZFNa~UMI$cu}gJ)$Rn z^WNXAIs>`_{w;a!3(?kHzTk_ZV-B!{b}z6?T_AYT6a%7u&|xWD1vK-Qt>FQ%3oGrb zkZ$3vbP?o)p(#L3pysIDUgo(@jHOBXzir=dADEsMF))9LxDPD(D3ZC*%6A_RuP^OYlHju)&MpR2O!lL5!l z|8|6dlDAvZH*g}3)_ekA-Z1Q3Mjt7g^slt467Yx48v_L-PB&jG9LL2ECk#DLegpSW2lQ7F`0sd8K9m*yYF?nsn|BkA?>slcL zaZkGU+4cf}CQ_87?aXvf7q(@k-b=jGTAeHzCV~F{x z^+(sYc8_K#WI5=w>$J3p^%hm6c?3E)m!J&J4Dt@mr#5uIoP2l-UneNLN*Xgv!= zZA>)-DDDQEXQ!}!3LkR!pQxRDyj^U+!+f^?a3)ssadygDTWu9@%;EghtFr}}e~YDL zgKa`)Hzb}g9XH&RJ*c+D$1$7iln@ed_|A9FM7%L|ul(iHp&fJdP|gXSufbco+QF-d z4I3z$9ozT7NP*QSaei}Sy3lfi`*oi@BM$jP6_xBLOlbmnLHBUMw$RBizd}YJOTsX| z{$pHWFEhA?Y~Ir}-5s``h(;t(&)z%zpB@waFW~OuBQCE*@gD9CSJsaW0AOiqV^V4C G8T$_lp6`PI literal 0 HcmV?d00001