From e22d6ac36b548a1550e799bb50d203f212b933a2 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Fri, 29 Sep 2023 21:56:28 +0200 Subject: [PATCH] Starting to add Poisson example --- examples/.poisson.jl.swp | Bin 0 -> 28672 bytes examples/poisson.jl | 352 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 352 insertions(+) create mode 100644 examples/.poisson.jl.swp create mode 100644 examples/poisson.jl diff --git a/examples/.poisson.jl.swp b/examples/.poisson.jl.swp new file mode 100644 index 0000000000000000000000000000000000000000..b84e5208f1304e8c8b4b31be4fe07ce26a0706f1 GIT binary patch literal 28672 zcmeI44~!gFdB6t}+)zkNp-xmlVeT9?vzNWI&xs*e*yNh4jV~emvEv9f_G<2Sc6WTc zvv>2yw?~>kRTCPKsG+T(g$6{5K!i$^wx~);Nz~MU0<9{QMnF;9B(wp9N+?C7sDKa%R@xc#&b;J9UwB%Zby* z)4yOCP2cvGoM1_^Hai+lr!&2G@2hrCpNckp-Cc+5j^i)8-qBvKv+Ra5 zCp){GzTI8vIKi%!o*M)`kF?FIEb~eWlolAr0#VqUzVdwI@;$S=MSSPvg!z&eUq4RC zQl_*(X@SxLr3Fe0lolv0P+Fk0Kxu*hCl&}-o@d-jU!EKHCydX}$(_&5tYO4ZjUF zxE!{@cg{DAZ^6I9=iv`v72XW5gqOn8FEotL!H3}AZiQQ*4HNM11U?>w55cd(E8uSlc02-~g!|xbxC;VsU>{rrj}cV45AK8p zTnn#&ZSXw;B7YC}!`;w_7$P^fiCnd?1ZDffssp#{c5I(|3x#)@ zbDd6a$qpUUj{0uLwf%K-F={p)KQNR6i`1p{^P0ZnXswRbmwLTH#=%a{4qtgi)ss{$ zu#gcr?1h)_NoOD4JUgY8v~$YtL{6-Z<2968Vw?2uqH9T0I+>U4H)jpGOdt8;N#0V= z3tcaAj3gmByG_&UHJp0btBWb>a#8c#J<50fUS)N&O{B9K5q1*FDSE1C=3dF^{s9SH^?1r+sPSCD4 zbR(Ib#F8SVHOYwTd%xe@t?J%sUfqkl(3BD7c%knmLn?U?J9ZFO?M2;;Jb~0M52)(< zAR~S6EPJ-PII9bo9E-EDLZQ7FbySDzU7M!*6IRf+SDbn?@|MI7LB$fwsF7XOk5{+q zRjRAQ%Iq_#>|oDdv;D@r+icFXe7oV&!?<3j6Rxi~)@r4iZ^26$S`7@#cbauE@KB2` zGVvHl>6YtETu&~wptjXml6&AcT#soduv{;6P+6=}_oGfc!`L$KsOd469aIw~gNoVc zF>tfR2ifsWKO%NF2X@z!n&=^w%RNZxQq5NuPr_>{@ioV7wZmYPXnJ@^?W>X9pmUM$ zsK#1Ouj_>V`p`Dj2E9+_)eY^(4QIQYT*W;+JhXr*2Gm_lo)h0~P>HCbo)XGBzQc1? zcN>|9C8<@?}| zOdWTL1O)V4qPoV-(NDEh98b|Z&G}^N&b(GVC+gasSC?9;h0fXCptdiatX~-MX3@Hb zf+@jtvK9<_ol|KEb*;1Xh)||ViG9#sL64^Iq zk_+w+GJ7&InV~}N5ng*h4mrdLns*oxrs{^+;OQ?&@q6DphNi)9(0@(T!4M-dzfBI~*Sh-ulhlqru%6c^3!wL`qwy_|mrQ>SA2AyyI%-^8N$5jbtR%mXDkAVYVz^W%J$Z zwPEa>uyu6O6P90iaJoUp;ifVJ$b2g#6324Xdx=z6g-R6seo)qdchYg%aYlF6_$ zIj-U}#4{(~C=-Z0Fy@*YwoU3X_gzAN0yFZWz-j1VH|XT`-FbDFIFHaBH8W}_8l~>R z89j?e6fiu^8+2%I_Jx7UL>=~~xwY=5rPau$6%5tEPErC7}EuJruS zOnreFaZ+xty=ezB> z5AZ2?5BwUm;X1eiOn4!D7k~dB;WXS2?}N9)GQ0uy!ujxZeEmn^ez*&6f;wCcCTxQ* z;TBkgd5{>v z3*d2V^Duk_Zi7B_pbae$+d0sHI-Gzv!By~DH~176Ub!a3Cb3nb_Ba)eBO@jrSI!j+r+!8byEken`b;K&!XsA6j6S&8VzF?9$PRq=Ef|8?bG*8OC`1nH5q2%weuL z@+wn)DD8%M%0X6{>Q^?`>Fkj#PPz2bQp$#m-(*sa8M9Rh`neCoLeQQy=ex@TnD5`*=gy0$2s@tc@8( zW=T3p5!XAl5KfU$XrIxYYaFAny%<;xCtL_CQ)^=0tkqQ2x(N%BjBQjY@;GnEE7Dl{ zC3iAd?N4P#hFEMt7asQB|V8e z$MN*FIKn4Qu1TnDt;RGs>}nusK{gH3^(9%(8*)?uF>W@;>kQjeY7WOpBz3X2mMgHu zqLC=9!^){3^Eiu2dZtXGIeLziD?@6;vYQ8RWzJ0nl9NJPEKTWBgQ-7~P+z>BnKFi_ zev>8d9SLiz>%oUL+NmW!1XPI3nj?mz6qu7F6*aKbqGnaCd{)&0X;m$v>Ax$wlQOKt z|NjQ>=f`AyIJN%&C4BoY!h`TT@NT#be0U4I8g{_5`1ap`FTexvez+Z0;6|8-3HTwt z{a4`AAbSCBhc<}s|614q&xh~e-+vJvg!^C}j>2Bp0ngyye+@nle+aSgh5dZ%_@$LT(PQwH6`*0`Na0q7LVmJq$!GC`OJ_{d!J3wLp zH^MK$CGdUx`LDo3@KLx2x^N9#0zVJufdS9r%RdF5hSTsHAb$R3a1J~dzJ_oAF$kdz zd*O%p^gn>F!n;BC4OnmiJcWOM8|;FY!FKpI_f_ugr{F!X3OB%hn1J)(<-nL99GQ7x zkjjrtMry(k89iy#o7tyVa`es{nuW`+lC?S*3Nq!Umf))C)dpG5 zWY7%XSZkSLBNt=SMe+kn+S_bdRytV8EX#_%1zN=eT~m%d zofM|B_`6_`H0<>tDZcSs$E&!i;dmi;I(9UqW=S7x_=roXwBd`b+9)cZAAGTq5+AlS zSTB}RhP1)%lvw5T;!nk#7;d?^POmGn>B4KM)zY)v#-rY+HCMI`$qRjs&*19G?2IKE zsh96itaZ|m)bycV(kThoNaxCL)rfoq|9Ex`wcNt(LDCXg!;f9H+$=Xob`%+soj%w} zlze06a@F^1{P~)7jPm1N>KnDo$ZVY~RFmXF6?3-t60x9`$`RGno}`S=mRD=B)?y14 z_HJ@hvpvs?Zq|blb{^Bl#r7?>ZCmOVonc6>JY20 z(P190fEU3(q1Ru)d!Y~IQCgt1Kxu)}0;L5?3zQZpEl^tEr_}=DI>?a^ALOsckr<<* ze6_jj++~qosu7)${Jbh?_twlKMsmHe7Rkw{=C4I&ldxM+6_O@LVE|Dw2nLS`!Ab6z zVK@;j4#i1MacoR-iaJ=7rfQMBB`QwVO2^5DHWeW`M#st|w=jB)#@NSYw}_9jvlvR^tunxg4ASXKI=Kqq0*mqX3QiP9s`!tjOcP tAgapNW&Ty$^a1`7CG|gAt8sO_n%U+!+?|<`q(0d^nfpg?sU4cee*+rK2?YQE literal 0 HcmV?d00001 diff --git a/examples/poisson.jl b/examples/poisson.jl new file mode 100644 index 00000000..87880d81 --- /dev/null +++ b/examples/poisson.jl @@ -0,0 +1,352 @@ +import GalerkinToolkit as glk +import ForwardDiff +using StaticArrays +using LinearAlgebra +using SparseArrays +using WriteVTK + +function add_default_params(params_in,params_default) + UmD = setdiff(keys(params_in),keys(params_default)) + for k in UmD + @warn "Parameter with key :$k is unused" + end + UiD = intersect(keys(params_default),keys(params_in)) + DmU = setdiff(keys(params_default),keys(params_in)) + a = [ k=>params_in[k] for k in UiD] + b = [ k=>params_default[k] for k in DmU] + Dict(vcat(a,b)) +end + +function default_params_poisson() + params = Dict{Symbol,Any}() + params[:mesh] = glk.cartesian_mesh((0,1,0,1),(10,10)) + params[:u] = (x) -> sum(x) + params[:f] = (x) -> 0.0 + params[:g] = (x) -> 0.0 + params[:dirichlet_tags] = ["boundary"] + params[:neumann_tags] = String[] + params[:integration_degree] = 2 + params[:solve] = \ + params[:example_name] = "poisson" + params +end + +function poisson(params_in) + + # Process params + params_default = default_params_poisson() + params = add_default_params(params_in,params_default) + + # Get the parameters + mesh = params[:mesh] + u = params[:u] + f = params[:f] + g = params[:g] + dirichlet_tags = params[:dirichlet_tags] + degree = params[:integration_degree] + solve = params[:solve] + example_name = params[:example_name] + neumann_tags = params[:neumann_tags] + + # Dirichlet values + node_to_tag = zeros(glk.num_nodes(mesh)) + tag_to_name = dirichlet_tags + glk.classify_mesh_nodes!(node_to_tag,mesh,tag_to_name) + free_and_dirichlet_nodes = glk.partition_from_mask(i->i==0,node_to_tag) + node_to_x = glk.node_coordinates(mesh) + dirichlet_nodes = last(free_and_dirichlet_nodes) + x_dirichlet = view(node_to_x,dirichlet_nodes) + u_dirichlet = u.(x_dirichlet) + + # Reference cell + d = glk.num_dims(mesh) + ref_cells = glk.reference_faces(mesh,d) + cell_to_rid = glk.face_reference_id(mesh,d) + + # Integration rule + integration_rules = map(ref_cells) do ref_cell + glk.quadrature(glk.geometry(ref_cell),degree) + end + w = map(glk.weights,integration_rules) + q = map(glk.coordinates,integration_rules) + + # Shape functions + ab = map(integration_rules,ref_cells) do integration_rule,ref_cell + q = glk.coordinates(integration_rule) + shape_functions = glk.shape_functions(ref_cell) + m = length(q) + n = glk.num_nodes(ref_cell) + a1 = glk.tabulation_matrix!(shape_functions)(ForwardDiff.gradient,zeros(eltype(q),m,n),q) + b1 = glk.tabulation_matrix!(shape_functions)(glk.value,zeros(m,n),q) + (a1,b1) + end + ∇s = map(first,ab) + s = map(last,ab) + + # Count coo entries + n_coo = 0 + cell_to_nodes = glk.face_nodes(mesh,d) + ncells = glk.num_faces(mesh,d) + node_to_free_node = glk.permutation(free_and_dirichlet_nodes) + nfree = length(first(free_and_dirichlet_nodes)) + for cell in 1:ncells + nodes = cell_to_nodes[cell] + for nj in nodes + if node_to_free_node[nj] > nfree + continue + end + for ni in nodes + if node_to_free_node[ni] > nfree + continue + end + n_coo += 1 + end + end + end + + # Allocate coo values + I_coo = zeros(Int32,n_coo) + J_coo = zeros(Int32,n_coo) + V_coo = zeros(Float64,n_coo) + b = zeros(nfree) + + # Allocate auxiliary buffers + ∇x = map(i->similar(i,size(i,2)),∇s) + Aes = map(i->zeros(size(i,2),size(i,2)),∇s) + fes = map(i->zeros(size(i,2)),∇s) + ues = map(i->zeros(size(i,2)),∇s) + + Tx = SVector{d,Float64} + TJ = typeof(zero(Tx)*zero(Tx)') + + # Fill coo values + i_coo = 0 + for cell in 1:ncells + rid = cell_to_rid[cell] + nodes = cell_to_nodes[cell] + Ae = Aes[rid] + ue = ues[rid] + fe = fes[rid] + nl = length(nodes) + ∇se = ∇s[rid] + ∇xe = ∇x[rid] + se = s[rid] + we = w[rid] + nq = length(we) + fill!(Ae,zero(eltype(Ae))) + fill!(fe,zero(eltype(Ae))) + fill!(ue,zero(eltype(Ae))) + for k in 1:nl + nk = nodes[k] + gk = node_to_free_node[nk] + if gk <= nfree + continue + end + ue[k] = u_dirichlet[gk-nfree] + end + for iq in 1:nq + Jt = zero(TJ) + xint = zero(Tx) + for k in 1:nl + x = node_to_x[nodes[k]] + ∇sqx = ∇se[iq,k] + sqx = se[iq,k] + Jt += ∇sqx*x' + xint += sqx*x + end + detJt = det(Jt) + invJt = inv(Jt) + dV = abs(detJt)*we[iq] + for k in 1:nl + ∇xe[k] = invJt*∇se[iq,k] + end + for j in 1:nl + for i in 1:nl + Ae[i,j] += ∇xe[i]⋅∇xe[j]*dV + end + end + fx = f(xint) + for k in 1:nl + fe[k] += fx*se[iq,k]*dV + end + end + for i in 1:nl + for j in 1:nl + fe[i] -= Ae[i,j]*ue[j] + end + end + for i in 1:nl + ni = nodes[i] + gi = node_to_free_node[ni] + if gi > nfree + continue + end + b[gi] += fe[i] + end + for j in 1:nl + nj = nodes[j] + gj = node_to_free_node[nj] + if gj > nfree + continue + end + for i in 1:nl + ni = nodes[i] + gi = node_to_free_node[ni] + if gi > nfree + continue + end + i_coo += 1 + I_coo[i_coo] = gi + J_coo[i_coo] = gj + V_coo[i_coo] = Ae[i,j] + end + end + end + + # Neumann conditions + if length(neumann_tags) != 0 + + face_to_rid = glk.face_reference_id(mesh,d-1) + face_to_nodes = glk.face_nodes(mesh,d-1) + ref_faces = glk.reference_faces(mesh,d-1) + n_ref_faces = length(ref_faces) + + # Integration rule + integration_rules_f = map(ref_faces) do ref_face + glk.quadrature(glk.geometry(ref_face),degree) + end + w_f = map(glk.weights,integration_rules_f) + q_f = map(glk.coordinates,integration_rules_f) + + # Shape functions + ab = map(integration_rules_f,ref_faces) do integration_rule,ref_cell + q = glk.coordinates(integration_rule) + shape_functions = glk.shape_functions(ref_cell) + m = length(q) + n = glk.num_nodes(ref_cell) + a1 = glk.tabulation_matrix!(shape_functions)(ForwardDiff.gradient,zeros(eltype(q),m,n),q) + b1 = glk.tabulation_matrix!(shape_functions)(glk.value,zeros(m,n),q) + (a1,b1) + end + ∇s_f = map(first,ab) + s_f = map(last,ab) + + tag_to_groups = glk.physical_groups(mesh,d-1) + neum_face_to_face = reduce(union,map(tag->tag_to_groups[tag],neumann_tags)) + fes = map(i->zeros(size(i,2)),s_f) + + Ts = typeof(zero(SVector{d-1,Float64})*zero(Tx)') + + for face in neum_face_to_face + rid = face_to_rid[face] + nodes = face_to_nodes[face] + ∇se = ∇s_f[rid] + se = s_f[rid] + we = w_f[rid] + nq = length(we) + fe = fes[rid] + fill!(fe,zero(eltype(fe))) + nl = length(nodes) + for iq in 1:nq + Jt = zero(Ts) + xint = zero(Tx) + for k in 1:nl + x = node_to_x[nodes[k]] + ∇sqx = ∇se[iq,k] + sqx = se[iq,k] + Jt += ∇sqx*x' + xint += sqx*x + end + J = transpose(Jt) + dS = sqrt(det(Jt*J))*we[iq] + ∇ux = ForwardDiff.gradient(u,xint) + gx = g(xint) + for k in 1:nl + fe[k] += gx*se[iq,k]*dS + end + end + for i in 1:nl + ni = nodes[i] + gi = node_to_free_node[ni] + if gi > nfree + continue + end + b[gi] += fe[i] + end + end + end + + A = sparse(I_coo,J_coo,V_coo,nfree,nfree) + u_free = solve(A,b) + nnodes = glk.num_nodes(mesh) + uh = zeros(nnodes) + uh[first(free_and_dirichlet_nodes)] = u_free + uh[last(free_and_dirichlet_nodes)] = u_dirichlet + + eh1 = 0.0 + el2 = 0.0 + for cell in 1:ncells + rid = cell_to_rid[cell] + nodes = cell_to_nodes[cell] + ue = ues[rid] + nl = length(nodes) + ∇se = ∇s[rid] + ∇xe = ∇x[rid] + se = s[rid] + we = w[rid] + nq = length(we) + for k in 1:nl + nk = nodes[k] + ue[k] = uh[nk] + end + for iq in 1:nq + Jt = zero(TJ) + xint = zero(Tx) + for k in 1:nl + x = node_to_x[nodes[k]] + ∇sqx = ∇se[iq,k] + sqx = se[iq,k] + Jt += ∇sqx*x' + xint += sqx*x + end + detJt = det(Jt) + invJt = inv(Jt) + dV = abs(detJt)*we[iq] + for k in 1:nl + ∇xe[k] = invJt*∇se[iq,k] + end + ux = u(xint) + ∇ux = ForwardDiff.gradient(u,xint) + ∇uhx = zero(∇ux) + uhx = zero(ux) + for k in 1:nl + uek = ue[k] + ∇uhx += uek*∇xe[k] + uhx += uek*se[iq,k] + end + ∇ex = ∇ux - ∇uhx + ex = ux - uhx + eh1 += (∇ex⋅∇ex + ex*ex)*dV + el2 += ex*ex*dV + end + end + + results = Dict{Symbol,Any}() + eh1 = sqrt(eh1) + el2 = sqrt(el2) + results[:eh1] = eh1 + results[:el2] = el2 + results[:nnodes] = nnodes + results[:nfree] = nfree + results[:ncells] = ncells + + vtk_grid(example_name,glk.vtk_args(mesh)...) do vtk + glk.vtk_physical_groups!(vtk,mesh) + vtk["tag"] = node_to_tag + vtk["uh"] = uh + vtk["dim"] = glk.face_dim(mesh) + end + results +end + +