From d5d7ae8ca73dc5aa018ae00e475e7749de40c873 Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Mon, 24 Jun 2024 16:08:31 +0200 Subject: [PATCH] updated version requirements, several small improvements, pngs instead of svgs in examples, version bump --- Project.toml | 6 +- assets/example280.png | Bin 0 -> 7169 bytes examples/Example200_LowLevelPoisson.jl | 31 +++--- .../Example205_LowLevelSpaceTimePoisson.jl | 20 ++-- examples/Example210_LowLevelNavierStokes.jl | 19 +--- examples/Example280_BasisPlotter.jl | 6 +- examples/Example281_DiscontinuousPlot.jl | 10 +- .../Example290_InterpolationBetweenMeshes.jl | 4 +- src/ExtendableFEMBase.jl | 3 +- src/dofmaps.jl | 98 +++++++++++++++++- src/fematrix.jl | 3 + src/finiteelements.jl | 25 ++--- test/Project.toml | 1 + test/runtests.jl | 3 + test/test_fematrix_and_vector.jl | 39 +++++++ 15 files changed, 192 insertions(+), 76 deletions(-) create mode 100644 assets/example280.png create mode 100644 test/test_fematrix_and_vector.jl diff --git a/Project.toml b/Project.toml index f67107c..7828a98 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ExtendableFEMBase" uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c" authors = ["Christian Merdon "] -version = "0.4" +version = "0.5" [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" @@ -20,8 +20,8 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] DiffResults = "1" DocStringExtensions = "0.8,0.9" -ExtendableGrids = "1.3" -ExtendableSparse = "^1.2" +ExtendableGrids = "1.8" +ExtendableSparse = "^1.4" ForwardDiff = "^0.10.35" Term = "2" UnicodePlots = "^3.6" diff --git a/assets/example280.png b/assets/example280.png new file mode 100644 index 0000000000000000000000000000000000000000..2296fbbb619a22e32331855fdd65f38a53782921 GIT binary patch literal 7169 zcmZ{J1yq#XxBkljiVWS|f^-SeB?u~lATZJ)-67oy(k0R*-3;9f(jX<00>glSGy~E( zcl7st|Np&f-TSUJ@64RF_c{Blv+F!h#7i}0d|Vn_008ixswip#00sp8yBPYy3)n9&=h+3Z{cIJCuxuJ+&zQ30)?PBa|6}-uc&BV}s_(gjsBWEIA020< zS41hBQQsy*M7AS{1ViKwnFR<(sN2sa5O4L1B;&Dapt@{l=h{|daa*ys7h?Zx>rub; zLy&?>x4AjM%qEzMP4uz#wGh*tq@<+M(o%D>Tx-(?J@qUJr6;~U!Gd;HNX24L;4ATXeOM7^N;-^e0o{IHHu3US zk4jx^6Iv<%^F*T{fY{|#0VIvW29BAdFn?!K>nthwC_!m(?^y;<$H?nakH7=$PuesC zd`AAAv@DU$gZNQdt9@bgWeb9hDUQB~=Jmp9ivKNx41Bn_ciPoI9Immoux z!p~$e%oHG`)Lf1B-(M?>q$4FH-dO}EE$AV|FJ?*KS5>qsKZ8jwPFQf;+fMrh92879 z5l|;sO*inNs#dLMvCOCIlkpg~HYXAis_r3BIndl(3}yV>2NYf8`*RBqk37Qi0#B7^ z8*}CMg-c5f=Bv29o&xHiF8eZuRxE~7z57o!K>z81llVk>n6z0V4%u~N>g7I z_^hR>t~#zC70pJQ=J)%?8&K-Z%~Ky6tkhmuk{=+SGQ;Kh!S2Pl*WxUH7Q3 zo}CTRnM-{Rg>ts-JSlN}|13Y9fe{M;QbsKy_c4QnbY2-WRh}fw6vbR*Hw=c4O5_TW zvtt6OwW?f&LUJMLX%^mNG2-T*d?SzXsf+OSLUl4=;;g(JQ=DK@UwsN0T}zA0%V#lq zm8AwbMrlLKHFZrE$)lf^cHZI*X^yBVd0y8QzY$H2xAR(TC0(WZASOgd#H}27xg0zR z0P)71`3>DZ2ldnZ->)t%WZpEc+GBt?u6FVF?kH33r=PT282LJ-~ z5z0LIVI92SoYgrbObk8b+i1j|Jkhh1tozF-Hz|ucQX|deck+lJTU*1t55XZ79Ot6+^HRW0IA=uA7!VXt$eHG zlQdK>U%9+7`Ga+sN4N^As!{%mYnFDn*|E@9HY?igaaxw1W`^xdyo;=$zvtrm6Lp_V zU9IJ3)I~ptf(GY7U(Gu-9(8D8ba((kM;x%yN(cYmqk?KX5ZSNoBHOy}*g#{~SdY~J zND;u}&L~SMHOpE}k&*TAzym!gS#DBx&W(%H@u@#0K6fa*lF1@mrH4b`#&^k?=LAr- z<(X~|VKse~I$5inlf5Nud2&wSrTIWLI48d(IyS72G^&mxgsKZL`4I^F36=_P{q^TWp0Q+%n~>(5J*+Y zCX0-CwquZjG)R`~|M6J`DwI|>g)1~h>BJZu9fcvghgxo4%n@^Q*LzwhlMt4bcjrn75{+2xd0-1>BIBiib>GZ`-a3CgBvKgJMyk#!_p?yEkOt+y zfh2Wg_WgYDe(blx{!cRuE43A5Zd`kTU|YwFa*ow*W_}Z2LUM6+b9hrsv@2g^t(-xm z9b8if3AUcLp;c|A!^jc-qwy!5>GsbiigTH#S@eIcC8Ma5tMHD2g(qBrs6Xs6R`~74 zmGh2A;?V2wq5$Aki6IY{htm`a{?$R~)k4 z$L(~Z_gW-Tx)_i{Ei<3A z_FZ=-$r*i1?Pv=fz1lM{*?zFbjfdV}Y>6^cBI#@{d_9x#Kywy39e2IZcZ5g;Umpax z3&dSqvIxeFkUr}p?re}yHQj#jvrgoGLn-kIPEi&7B`kpw1$XXMHTI=`HQOFfAOAylH8z8OKX5R~%( zX^)LMnU2xswwOHM36*#LPST@66K9P%)`MJ#G25rchp)~zXPtG^dayB$ao(LTNf&;b z>eM;^wp+p+=rWwcKX!l) z@eHNN#xLYVU*d#ymr-AW(ZWnN{52Eh<3b{77%P65oNd^(k~>9s1xzKu$U;x6&~l|o zk!r47OtTnWnTLfe%u2}e20RN^OsYH{1pP|diA{#yPL?g6)BwOJHaAMs=r-l|Z%pid zdukx4>-9nLbUC72@2o7Lz3~tRD7+)>B428EZ8$aC2tZC&@3JYO93+7$=<8P^Qx}V) zKXR2UF6{LBC30xeAANr?V$ZIowH;L;c0~A!TzSbmV6d6&rMSoatod9Z9n51` zSZ8XP*v~WirRzM`-A@Q1HKkaItul8%%_8kNhy48-8DBMdDZgK`JYxYz#<13VzzH_@ zY01D@95$yY07TSdmrPOKE8PA-_%<%9l!jIiJ(rKMs_wFl=yY>88p!y_>XWu72k2N*+|ur%sVN zD%3>3e0rXQFg`>IjIc)VER&YNRIKap4x@bfH0rf+CsqSR-;GKtB7EV-9VzA)iQay` zXYE75{*HuF)8cY>-X|)?eHrNdj}CWcWWIGd+B;D80^&#_Y>(m!uE6Qwq?RQ4fMR@` zi@Q3NKAL;km5;J9mGPO`#K6WYxeFlEZq?X$95}4LBO_XgM`efENBZL9YSY=c`4ZA^ z>%NuJcl))is~S5=^%C1OYbaS{bP0DUTMD+$bZZ%8w1A?;)(-AH1f-3LG*oS~jDU_p z%90~apXCl77(JbC{<)k{cy&sSG!o|Uv`RM;f8E%I?uU+TI)9wy`$JLZ%S#Vk5eNS| zX0usNYv;}A56*A^_(`IF$~P4Ay!m1*O*ZrB&U05IC%Qln9@YB}J;!@wu>@M=duFRI8zChy8-1KDMYR4c`@EE%%NafpZQ6H`i074%SapNf?;-Ujk$q-%P~Wj`t|;vQgM0=;AN?2 z6K9;XuQwIwK29Ks19)VQxPH%F%J9AC-S?k0dEW!a^U3fZG1&~HAY@M4qU9U1U-k*0rH9<~MUXATauCvUhfqk)Rcne9S?|rIX zgvHQgo}~7FKwQgp_&UM#81+vzyG^W+GMD*bY_k;B3$$L*L`j}RDq&Fn9Y)iZK5n6| z4f+$=u#waj*WI6Wd(y=plnEPi0=Kdjg0~})uKw+wV{f`m&>m8@V<=J%LPrylDV_i| zPc^33&O+*!%#3my|KJt{ZPL-~$6Rl|v&p3;=+(#$ zlJL+Hdv;&fhD`0gs+W~WYw!B*NQDhcepiT_jyH1zQlQ;COiOR0-%q*OE%Qm3GnZ+% zWckGXF@r3f+U#vf>>ly*_AluB-@Beea=wU0IdeCp%0?;1jfUfBOH14~n{%f?t?mab z*?FN61fENz{XYGcl*nhvue_OODN*|7b$Isz>hGW85RuTuLJuA|fgf)C+MnoCo?$)E zpr&u*+MPI?Vk-J9==+G(X*gz|LEd0&Y$V;m{FB5k@xhQJ@};>nh!I5{llR!(GR?(b z7k95*c-|N5tra!;Qc5f4$j{s>h@QB_5zG|wApAGa7gc!F*UdR?3x4&;SgWvI?G~!H z$&{{?lmj3E_g@#}?V1mScC)w z5JGn`F4-$}Iwp11T~3C3m(ogxWsX%&#aRs2onOr>=UHGP)_Xh^6h~}t5J34>Nr#h+ zS_T$HqVGs99`I-+VM~$oxY_A7bJSafb(fkD z>-ZB-CMm2WP(HM568;@P`9x8J;y|bR+*a?L2`!?47MYKJq6)zcRZiBmtT11`;luma z-^7rCtf_1jxD)~0u^|$S1f9Y(H4yZEf~P1Cj~%Mzt;u<=b?EN7%g_Z_Jv4)&Ed`bG z|A7_$xD|tv14eGF+L5$~NuXb@^AlXnF8EoWMd3>^tl$1h8Wcp^n-3AMFi!_4aV-j% z$PEIF4RmlrA2T5+|M4s5i6-pxQu@5Giv*>SQK4SxU!C$7_E$MM&vKnj;)1BJrwAr` zt_N9)EUABD9wAKp8jA`qnpPWh0;v&hD*A5=B_``eRDU z*Y;6I@~tL z4fSlHv%T5WcakP*@_Vsa1-F#sSyi8AdG5MFoe9EWykGnpqMOFpX4;<^Yn7_wku#&5 znn{Ef%aQII)|ReZi{;#dh2W>WWmT1iNhzxReSMx}Luq!lcI8o;Lm~?Ko&-0efwhu> zYFsFsS1ghqAC~m_mEQRk~&;t&ok=?_Bs#^_7I%(<;fiDEagtz;Lgu=tKa(hdNS@Kk+$v7ad} z4R&H{sBeYAB%JMJ#i2zoSnH>p8P>ZcbZtmFpNpL7|DX;4>?Y@H_`m3W~oTVsBE*!A&wH*ehS>h~dNn-f=L`NrY9SW$K-?4D{t;ZW^e2X4)XbEg5+q=?dvr~}t{B&Ny;7oy(AkT{AM z^^6N8eJqR1MU2|R3O^5Q{laoPfbI8h0%uT?es{5>NQmM9E#V}2)IycNKnH!WmfB8o zvuPCX`aO{gliL=2y@t@iE6&>8AwJ-)We83QWSc(^9H#p>T>puZ2?XTQ>#luV zepc$QK3Y8b84PuvQPj|vGdb$`rr7k1i?46J)VD?d&8o3`k8&?0ft)&b)J2i7^T-+= zdWQYCkT}EYh`?Avi02_862nDyWSCO)6%XCws7o)M{GJ<`C5B)uJK3rX7(2y{PEg6@Dq0%!H#Vej7eF$< z-}Udpc@?s2ltv7*s!>Z*Nl^h>bEofCF^c=t3$XyLynuZFH zL7;|?T0GWy&_+!f{elCL&Qh!Yp1cCHOi1K2tzljiyJm#8Ysl1vG*dc2Pr>wgYu0}lLBM!%| zm2)}y9YK0Gj-a)6o56;weSYO4D|?A)KVB&>yDK&TRO@j702w=6x`wY93^-=WT<dmodf8MEf>ASq?BB!y`?KYo1sRE}# z0h^WiVB&t7RC!KI*&Q}7*{qxK0#0*^&x13**CzY`*~}Cq?SmS@9kI>HyM&o{&^YQP zeI@yf=iR7sEd^QHd>5(#?r$OqbeYBo94^|I6GKbT~Z`LOhzqgZ-! zpysP$%OUYaL?=V1FJsVG4=UDw5yvr{$-tZ+3qp~fX{H76c>1 zAkq!VTdx4e&pVR(JDQ1p&fhGIEa}?TF5C!k-LxHQ_BJQ#vKigEs^>Too6TVj!BPu)nf(cW$91w#9@F`A`NyVybW#fc%c zNqf@RQnGlpT-@??gp;kEKI5Sj=WhDVs22N9OzT!5?q&n#rcPBQx6b=vYys)F$$ogt zxdMqhZ z-*_J-cpaj{R_l3}wF=#evd0*RgQAd)Yg>v0uhZ-?t6f6Fm_BS!Wgg_5d~(FPm)6vc zfUG>`S8z1es^$%fUv~BjeWU+>5FsNVoTph2a_mqCJ%sr9_$ZC_|K>sGrx$CwUqo$REv;_d8Pw(q9d&wU5%Nu;HygVZ zR+pBXz6O~5J?$+zrjYJ0z9}J9L3JF=AtgO$f2bc8Br#Kn_5a2P{j)7a`T71ppVjpb zXz;{)0d>=0H9SNfpBiy{0}2Nxri{aR{I55nE~_B>HZ(VJl#-6WA$W0AXqHJsNaUdJ z2pU XrFls*1V*7bJ>aR5nqrlL$-Dmo`3&ti literal 0 HcmV?d00001 diff --git a/examples/Example200_LowLevelPoisson.jl b/examples/Example200_LowLevelPoisson.jl index ff5e438..3c46e5e 100644 --- a/examples/Example200_LowLevelPoisson.jl +++ b/examples/Example200_LowLevelPoisson.jl @@ -17,7 +17,7 @@ for different refinement levels. The computed solution for the default parameters looks like this: -![](example200.svg) +![](example200.png) =# @@ -61,10 +61,10 @@ function main(; maxnref = 8, order = 2, Plotter = nothing) sol, time_assembly, time_solve = solve_poisson_lowlevel(FES, μ, f) ## plot statistics - println(stdout, barplot(["Grid", "FaceNodes", "CellDofs", "Assembly", "Solve"], [time_grid, time_facenodes, time_dofmap, time_assembly, time_solve], title = "Runtimes")) + println(stdout, barplot(["Grid", "FaceNodes", "celldofs", "Assembly", "Solve"], [time_grid, time_facenodes, time_dofmap, time_assembly, time_solve], title = "Runtimes")) ## plot - scalarplot!(plt[1,1], xgrid, view(sol.entries, 1:num_nodes(xgrid))) + scalarplot!(plt[1,1], xgrid, view(sol.entries, 1:num_nodes(xgrid)), limits = (-0.0125,0.0125)) end return sol, plt @@ -82,16 +82,10 @@ function solve_poisson_lowlevel(FES, μ, f) ## fix boundary dofs begin - BFaceDofs::Adjacency{Int32} = FES[ExtendableFEMBase.BFaceDofs] - nbfaces::Int = num_sources(BFaceDofs) - AM::ExtendableSparseMatrix{Float64, Int64} = A.entries - dof_j::Int = 0 - for bface ∈ 1:nbfaces - for j ∈ 1:num_targets(BFaceDofs, 1) - dof_j = BFaceDofs[j, bface] - AM[dof_j, dof_j] = 1e60 - b.entries[dof_j] = 0 - end + bdofs = boundarydofs(FES) + for dof in bdofs + A.entries[dof, dof] = 1e60 + b.entries[dof] = 0 end end ExtendableSparse.flush!(A.entries) @@ -122,13 +116,12 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES, f, μ = 1) ∇vals = FEBasis_∇.cvals FEBasis_id = FEEvaluator(FES, Identity, qf) idvals = FEBasis_id.cvals - CellDofs = FES[ExtendableFEMBase.CellDofs] + celldofs = FES[CellDofs] ## ASSEMBLY LOOP loop_allocations = 0 function barrier(EG, L2G::L2GTransformer) ## barrier function to avoid allocations by type dispatch - ndofs4cell::Int = get_ndofs(ON_CELLS, FEType, EG) Aloc = zeros(Float64, ndofs4cell, ndofs4cell) ncells::Int = num_cells(xgrid) @@ -152,9 +145,9 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES, f, μ = 1) ## add local matrix to global matrix for j ∈ 1:ndofs4cell - dof_j = CellDofs[j, cell] + dof_j = celldofs[j, cell] for k ∈ j:ndofs4cell - dof_k = CellDofs[k, cell] + dof_k = celldofs[k, cell] if abs(Aloc[j, k]) > 1e-15 ## write into sparse matrix, only lines with allocations rawupdateindex!(A, +, Aloc[j, k], dof_j, dof_k) @@ -178,7 +171,7 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES, f, μ = 1) temp += weights[qp] * idvals[1, j, qp] * f(x) end ## write into global vector - dof_j = CellDofs[j, cell] + dof_j = celldofs[j, cell] b[dof_j] += temp * cellvolumes[cell] end end @@ -191,7 +184,7 @@ end function generateplots(dir = pwd(); Plotter = nothing, kwargs...) ~, plt = main(; Plotter = Plotter, kwargs...) scene = GridVisualize.reveal(plt) - GridVisualize.save(joinpath(dir, "example200.svg"), scene; Plotter = Plotter) + GridVisualize.save(joinpath(dir, "example200.png"), scene; Plotter = Plotter) end function runtests(; order = 2, kwargs...) #hide diff --git a/examples/Example205_LowLevelSpaceTimePoisson.jl b/examples/Example205_LowLevelSpaceTimePoisson.jl index 30c5208..bd45156 100644 --- a/examples/Example205_LowLevelSpaceTimePoisson.jl +++ b/examples/Example205_LowLevelSpaceTimePoisson.jl @@ -17,7 +17,7 @@ on the unit square domain ``\Omega`` on a given grid with space-time finite element methods based on tensorized ansatz functions. -![](example205.svg) +![](example205.png) =# @@ -99,23 +99,19 @@ function solve_poisson_lowlevel(FES_time, FES_space, μ, f) loop_allocations = assemble!(A, b, FES_time, FES_space, f, μ) ## fix homogeneous boundary dofs - bfacedofs::Adjacency{Int32} = FES_space[ExtendableFEMBase.BFaceDofs] - nbfaces::Int = num_sources(bfacedofs) - dof::Int = 0 - for bface ∈ 1:nbfaces + bdofs = boundarydofs(FES_space) + for sdof in bdofs for dof_t = 1 : ndofs_time - for j ∈ 1:num_targets(bfacedofs, 1) - dof = (dof_t-1)*ndofs_space + bfacedofs[j, bface] - A[dof, dof] = 1e60 - b[dof] = 0 - end + dof = (dof_t-1)*ndofs_space + sdof + A[dof, dof] = 1e60 + b[dof] = 0 end end ## fix initial value by zero for j=1:ndofs_space A[j,j] = 1e60 - b[dof] = 0 + b[j] = 0 end ExtendableSparse.flush!(A) end @@ -311,7 +307,7 @@ end function generateplots(dir = pwd(); Plotter = nothing, kwargs...) ~, plt = main(; Plotter = Plotter, kwargs...) scene = GridVisualize.reveal(plt) - GridVisualize.save(joinpath(dir, "example205.svg"), scene; Plotter = Plotter) + GridVisualize.save(joinpath(dir, "example205.png"), scene; Plotter = Plotter) end end #module \ No newline at end of file diff --git a/examples/Example210_LowLevelNavierStokes.jl b/examples/Example210_LowLevelNavierStokes.jl index 232640b..318477c 100644 --- a/examples/Example210_LowLevelNavierStokes.jl +++ b/examples/Example210_LowLevelNavierStokes.jl @@ -26,7 +26,7 @@ Newton's method with automatic differentation is used to handle the nonlinear co The computed solution for the default parameters looks like this: -![](example210.svg) +![](example210.png) =# module Example210_LowLevelNavierStokes @@ -162,18 +162,7 @@ function solve_stokes_lowlevel(FES; teval = 0) @time begin u_init = FEVector(FES) interpolate!(u_init[1], u!; time = teval) - - fixed_dofs = [size(A.entries, 1)] # fix one pressure dof = last dof - BFaceDofs::Adjacency{Int32} = FES[1][ExtendableFEMBase.BFaceDofs] - nbfaces::Int = num_sources(BFaceDofs) - AM::ExtendableSparseMatrix{Float64, Int64} = A.entries - dof_j::Int = 0 - for bface ∈ 1:nbfaces - for j ∈ 1:num_targets(BFaceDofs, 1) - dof_j = BFaceDofs[j, bface] - push!(fixed_dofs, dof_j) - end - end + fixed_dofs = boundarydofs(FES[1]) push!(fixed_dofs, FES[1].ndofs + 1) ## fix one pressure dof end @@ -200,7 +189,7 @@ function solve_stokes_lowlevel(FES; teval = 0) ## fix boundary dofs for dof in fixed_dofs - AM[dof, dof] = 1e60 + A.entries[dof, dof] = 1e60 b.entries[dof] = 1e60 * u_init.entries[dof] end ExtendableSparse.flush!(A.entries) @@ -400,6 +389,6 @@ end function generateplots(dir = pwd(); Plotter = nothing, kwargs...) ~, plt = main(; Plotter = Plotter, kwargs...) scene = GridVisualize.reveal(plt) - GridVisualize.save(joinpath(dir, "example210.svg"), scene; Plotter = Plotter) + GridVisualize.save(joinpath(dir, "example210.png"), scene; Plotter = Plotter) end end #module diff --git a/examples/Example280_BasisPlotter.jl b/examples/Example280_BasisPlotter.jl index 57b3aad..7421fb8 100644 --- a/examples/Example280_BasisPlotter.jl +++ b/examples/Example280_BasisPlotter.jl @@ -4,7 +4,9 @@ ([source code](SOURCE_URL)) This example plots all the basis functions of a H1 finite element on Edge1D or Triangle2D -as unicode plots +as unicode plots. This is the result with the default parameters (dim = 1, order = 3): + +![](https://github.com/chmerdon/ExtendableFEMBase.jl/blob/master/docs/src/assets/example280.png?raw=true") =# @@ -14,7 +16,7 @@ using ExtendableFEMBase using ExtendableGrids ## everything is wrapped in a main function -function main(; dim = 1, order = 2) +function main(; dim = 1, order = 3) ## generate two grids @assert dim in [1, 2] "dim must be 1 or 2" diff --git a/examples/Example281_DiscontinuousPlot.jl b/examples/Example281_DiscontinuousPlot.jl index cff2fc5..5769628 100644 --- a/examples/Example281_DiscontinuousPlot.jl +++ b/examples/Example281_DiscontinuousPlot.jl @@ -8,7 +8,7 @@ on a grid with two regions by region-wise nodal values and plotting. The computed solution for the default parameters looks like this: -![](example281.svg) +![](example281.png) =# module Example281_DiscontinuousPlot @@ -35,7 +35,7 @@ function main(; broken = false, nrefs = 3, abs = false, Plotter = nothing) ## generate two grids xgrid = grid_unitsquare(Triangle2D) - ## mark a second region + ## mark first two triangles to be in second region xgrid[CellRegions][1:2] .= 2 ## refine @@ -53,8 +53,8 @@ function main(; broken = false, nrefs = 3, abs = false, Plotter = nothing) subgrid2 = subgrid(xgrid, [2]) ## get parent nodes for each subgrid - subnodes1 = subgrid1[ExtendableGrids.NodeInParent] - subnodes2 = subgrid2[ExtendableGrids.NodeInParent] + subnodes1 = subgrid1[NodeParents] + subnodes2 = subgrid2[NodeParents] ## compute nodevalues for nodes of each subgrid nodevals4nodes1 = nodevalues(FEFunction[1], Identity; abs = abs, regions = [1], nodes = subnodes1) @@ -71,6 +71,6 @@ end function generateplots(dir = pwd(); Plotter = nothing, kwargs...) plt = main(; Plotter = Plotter, kwargs...) scene = GridVisualize.reveal(plt) - GridVisualize.save(joinpath(dir, "example281.svg"), scene; Plotter = Plotter) + GridVisualize.save(joinpath(dir, "example281.png"), scene; Plotter = Plotter) end end diff --git a/examples/Example290_InterpolationBetweenMeshes.jl b/examples/Example290_InterpolationBetweenMeshes.jl index e2c0fd4..44ce750 100644 --- a/examples/Example290_InterpolationBetweenMeshes.jl +++ b/examples/Example290_InterpolationBetweenMeshes.jl @@ -8,7 +8,7 @@ this P2 function on two uniform refinements into some P1 function. Then, both fi The computed solution for the default parameters looks like this: -![](example290.svg) +![](example290.png) =# module Example290_InterpolationBetweenMeshes @@ -62,7 +62,7 @@ end function generateplots(dir = pwd(); Plotter = nothing, kwargs...) plt = main(; Plotter = Plotter, kwargs...) scene = GridVisualize.reveal(plt) - GridVisualize.save(joinpath(dir, "example290.svg"), scene; Plotter = Plotter) + GridVisualize.save(joinpath(dir, "example290.png"), scene; Plotter = Plotter) end end diff --git a/src/ExtendableFEMBase.jl b/src/ExtendableFEMBase.jl index 677310f..5c911a8 100644 --- a/src/ExtendableFEMBase.jl +++ b/src/ExtendableFEMBase.jl @@ -39,7 +39,7 @@ export NeededDerivative4Operator, Length4Operator, QuadratureOrderShift4Operator export OperatorPair, OperatorTriple -include("finiteelements.jl") +include("finiteelements.jl") # also includes dofmaps.jl and feevaluator*.jl export DofMap export CellDofs, FaceDofs, EdgeDofs, BFaceDofs, BEdgeDofs export CellDofsParent, FaceDofsParent, EdgeDofsParent, BFaceDofsParent, BEdgeDofsParent @@ -47,6 +47,7 @@ export DofMapTypes export Dofmap4AssemblyType, ItemGeometries4DofMap, EffAT4AssemblyType, ParentDofmap4Dofmap export AbstractFiniteElement export FESpace, FESpaces, get_AT, get_FEType +export boundarydofs, ndofs, broken export AbstractH1FiniteElement export H1BUBBLE, L2P0, H1P1, H1P2, H1P2B, H1MINI, H1CR, H1P3, H1Pk diff --git a/src/dofmaps.jl b/src/dofmaps.jl index ea6a98a..a200b98 100644 --- a/src/dofmaps.jl +++ b/src/dofmaps.jl @@ -6,15 +6,77 @@ information with respect to different AssemblyTypes. They are generated automati associated to each subtype can be accessed via FESpace[DofMap]. """ abstract type DofMap <: AbstractGridAdjacency end + + +""" + $(TYPEDEF) + +Key type describing the dofs for each cell of the dofgrid +""" abstract type CellDofs <: DofMap end + +""" + $(TYPEDEF) + +Key type describing the dofs for each face of the dofgrid +""" abstract type FaceDofs <: DofMap end + +""" + $(TYPEDEF) + +Key type describing the dofs for each edge of the dofgrid +""" abstract type EdgeDofs <: DofMap end + +""" + $(TYPEDEF) + +Key type describing the dofs for each boundary face of the dofgrid +""" abstract type BFaceDofs <: DofMap end + +""" + $(TYPEDEF) + +Key type describing the dofs for each boundary edge of the dofgrid +""" abstract type BEdgeDofs <: DofMap end + + +""" + $(TYPEDEF) + +Key type describing the dofs for each cell of the parentgrid +""" abstract type CellDofsParent <: DofMap end + +""" + $(TYPEDEF) + +Key type describing the dofs for each face of the parentgrid +""" abstract type FaceDofsParent <: DofMap end + +""" + $(TYPEDEF) + +Key type describing the dofs for each edge of the parentgrid +""" abstract type EdgeDofsParent <: DofMap end + +""" + $(TYPEDEF) + +Key type describing the dofs for each boundary face of the parentgrid +""" abstract type BFaceDofsParent <: DofMap end + +""" + $(TYPEDEF) + +Key type describing the dofs for each boundary edge of the parentgrid +""" abstract type BEdgeDofsParent <: DofMap end const DofMapTypes{Ti} = Union{VariableTargetAdjacency{Ti}, SerialVariableTargetAdjacency{Ti}, Array{Ti, 2}} @@ -384,7 +446,8 @@ function init_dofmap_from_pattern!(FES::FESpace{Tv, Ti, FEType, APT}, DM::Type{< if FES.dofgrid !== FES.xgrid ## assume parent relation between xgrid and dofgrid @assert FES.dofgrid[ParentGrid] == FES.xgrid "xgrid is not the parent grid of dofgrid !!!" - @assert FES.dofgrid[ParentGridRelation] in [SubGrid, BoundarySubGrid] "dofgrid is not a subgrid or boundary-subgrid of xgrid !!!" + @assert FES.dofgrid[ParentGridRelation] <: SubGrid "dofgrid needs to be a subgrid of xgrid" + SubGridAssemblyType = FES.dofgrid[ParentGridRelation].parameters[1] ## lift subgrid dofmap to parent xgrid to allow assembly on common parent grid ## by constructing a variable target adjacency with empty dofs for parent items in xgrid ## that are not in the dofgrid --> this allows to use the dofmap in assembly over full xgrid @@ -643,3 +706,36 @@ function init_dofmap!(FES::FESpace, DM::Type{<:DofMap}) end end + + +function boundarydofs(FES; dofmap = BFaceDofs, regions = :all) + + bitemdofs = FES[dofmap] + if regions === :all + if typeof(bfacedofs) <: VariableTargetAdjacency + return bfacedofs.colentries + elseif typeof(bfacedofs) <: SerialVariableTargetAdjacency + return 1:bfacedofs.colstart[end]-1 + else + @error "dofmap has unexpected type" + end + else + if dofmap == BFaceDofs + bitemregions = FES[BFaceRegions] + elseif dofmap == BFaceEdges + bitemregions = FES[BEdgeRegions] + @error "do not know where to look for regions" + end + bdofs = [] + nbitems::Int = num_sources(bitemdofs) + for bface ∈ 1:nbitems + for j ∈ 1:num_targets(bitemdofs, 1) + if bitemregions[bface] in regions + push!(bdofs, bitemdofs[j, bface]) + end + end + end + return unique(bdofs) + end +end + diff --git a/src/fematrix.jl b/src/fematrix.jl index aa2f46f..f784a4a 100644 --- a/src/fematrix.jl +++ b/src/fematrix.jl @@ -177,6 +177,9 @@ Optionally a name for the matrix can be given. function FEMatrix(FESX::FESpace, FESY::FESpace; name = "auto") return FEMatrix{Float64, Int64}(FESX, FESY; name = name) end +function FEMatrix(FESX::Vector{<:FESpace}, FESY::Vector{<:FESpace}; name = "auto") + return FEMatrix{Float64, Int64}(FESX, FESY; name = name) +end function FEMatrix{TvM}(FESX::FESpace, FESY::FESpace; name = "auto") where {TvM} return FEMatrix{TvM, Int64}(FESX, FESY; name = name) end diff --git a/src/finiteelements.jl b/src/finiteelements.jl index 08dc019..f9544b3 100644 --- a/src/finiteelements.jl +++ b/src/finiteelements.jl @@ -53,6 +53,8 @@ function Base.copy(FES::FESpace{Tv, Ti, FEType, AT}) where {Tv, Ti, FEType, AT} return FESpace{Tv, Ti, FEType, AT}(deepcopy(FES.name), FES.broken, FES.ndofs, FES.coffset, FES.xgrid, FES.dofgrid, FES.dofmaps) end +ndofs(FES::FESpace) = FES.ndofs +broken(FES::FESpace) = FES.broken get_AT(::FESpace{Tv, Ti, FEType, AT}) where {Tv, Ti, FEType, AT} = AT get_FEType(::FESpace{Tv, Ti, FEType, AT}) where {Tv, Ti, FEType, AT} = FEType @@ -89,33 +91,24 @@ function FESpace{FEType, AT}( end if AT == ON_FACES - xgrid = get_facegrid(xgrid) - xgrid[BFaceNodes] = zeros(Ti, 2, 0) - xgrid[BFaceRegions] = zeros(Ti, 0) - xgrid[BFaceGeometries] = VectorOfConstants{ElementGeometries, Int}(Edge1D, 0) - xgrid[BFaceVolumes] = zeros(Tv, 0) + if isnothing(regions) + regions = unique(xgrid[FaceRegions]) + end + dofgrid = subgrid(xgrid, regions; support = ON_FACES, project = false) elseif AT == ON_BFACES if isnothing(regions) regions = unique(xgrid[BFaceRegions]) end - dofgrid = subgrid(xgrid, regions; boundary = true, project = false) - dofgrid[BFaceNodes] = zeros(Ti, 2, 0) - dofgrid[BFaceRegions] = zeros(Ti, 0) - dofgrid[BFaceGeometries] = VectorOfConstants{ElementGeometries, Int}(Edge1D, 0) - dofgrid[BFaceVolumes] = zeros(Tv, 0) + dofgrid = subgrid(xgrid, regions; support = ON_BFACES, project = false) elseif AT == ON_EDGES - xgrid = get_edgegrid(xgrid) - xgrid[BFaceNodes] = zeros(Ti, 2, 0) - xgrid[BFaceRegions] = zeros(Ti, 0) - xgrid[BFaceGeometries] = VectorOfConstants{ElementGeometries, Int}(Edge1D, 0) - xgrid[BFaceVolumes] = zeros(Tv, 0) + @assert false "not possible currently" end if isnothing(regions) regions = unique(xgrid[CellRegions]) end - if AT !== ON_BFACES + if AT !== ON_BFACES && AT !== ON_FACES if regions != unique(xgrid[CellRegions]) dofgrid = subgrid(xgrid, regions) else diff --git a/test/Project.toml b/test/Project.toml index 1658796..7ef8fd8 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,6 +3,7 @@ ExtendableFEMBase = "12fb9182-3d4c-4424-8fd1-727a0899810c" ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8" ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" diff --git a/test/runtests.jl b/test/runtests.jl index f4058bb..5c7e8da 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,6 +2,7 @@ using Test using ExtendableGrids using ExtendableFEMBase using ExampleJuggler +using SparseArrays include("test_quadrature.jl") include("test_interpolators.jl") @@ -9,6 +10,7 @@ include("test_operators.jl") include("test_febasis.jl") include("test_segmentintegrator.jl") include("test_pointevaluator.jl") +include("test_fematrix_and_vector.jl") function run_examples() ExampleJuggler.verbose!(true) @@ -120,6 +122,7 @@ end function run_all_tests() begin + run_fematrix_tests() run_examples() run_febasis_tests() run_operator_tests() diff --git a/test/test_fematrix_and_vector.jl b/test/test_fematrix_and_vector.jl new file mode 100644 index 0000000..cbe6688 --- /dev/null +++ b/test/test_fematrix_and_vector.jl @@ -0,0 +1,39 @@ +function run_fematrix_tests() + + @testset "FEMatrixVector" begin + println("\n") + println("=======================") + println("Testing FEMatrix&VEctor") + println("=======================") + xgrid = simplexgrid(0:0.1:1,0:0.1:1) + FES1 = FESpace{H1Pk{1,1,1}}(xgrid) + FES2 = FESpace{H1Pk{1,1,2}}(xgrid) + A = FEMatrix(FES1, FES2) + @test size(A.entries) == (FES1.ndofs, FES2.ndofs) + @test size(A[1,1]) == (FES1.ndofs, FES2.ndofs) + + B = FEMatrix([FES1, FES2]) + @test length(B) == 4 + @test size(B.entries) == (FES1.ndofs+FES2.ndofs, FES1.ndofs+FES2.ndofs) + @test size(B[1,2]) == (FES1.ndofs, FES2.ndofs) + + + C = FEMatrix([FES2, FES2], [FES1, FES1]) + @test length(C) == 4 + @test size(C.entries) == (2*FES2.ndofs, 2*FES1.ndofs) + @test size(C[1,2]) == (FES2.ndofs, FES1.ndofs) + C.entries.cscmatrix = sprand(2*FES2.ndofs,2*FES1.ndofs,0.5) + @show C + + b = FEVector([FES1, FES2]) + b.entries .= rand(FES1.ndofs+FES2.ndofs) + + addblock!(A[1,1], C[1,1]; transpose = true) + @test norm(A[1,1]) ≈ norm(C[1,1]) + + fill!(b[2], 0) + addblock_matmul!(b[2], C[1,1], b[1]) + @test all(view(b[2]) .== view(C.entries, 1:FES2.ndofs, 1:FES1.ndofs) * view(b[1])) + end + println("") +end