From 15971114992e852037774dee9f02642cd057aef2 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 7 May 2010 12:58:08 +0000 Subject: [PATCH] Add Markus Buerg's step-45. git-svn-id: https://svn.dealii.org/trunk@21092 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-45/Makefile | 154 ++++++++++ deal.II/examples/step-45/doc/intro.dox | 63 ++++ deal.II/examples/step-45/doc/results.dox | 20 ++ .../step-45/doc/step-45.periodic_cells.png | Bin 0 -> 3855 bytes .../examples/step-45/doc/step-45.solution.png | Bin 0 -> 36656 bytes deal.II/examples/step-45/step-45.cc | 275 ++++++++++++++++++ 6 files changed, 512 insertions(+) create mode 100644 deal.II/examples/step-45/Makefile create mode 100644 deal.II/examples/step-45/doc/intro.dox create mode 100644 deal.II/examples/step-45/doc/results.dox create mode 100644 deal.II/examples/step-45/doc/step-45.periodic_cells.png create mode 100644 deal.II/examples/step-45/doc/step-45.solution.png create mode 100644 deal.II/examples/step-45/step-45.cc diff --git a/deal.II/examples/step-45/Makefile b/deal.II/examples/step-45/Makefile new file mode 100644 index 0000000000..a2f3f96ca5 --- /dev/null +++ b/deal.II/examples/step-45/Makefile @@ -0,0 +1,154 @@ +# $Id$ + + +# For the small projects Makefile, you basically need to fill in only +# four fields. +# +# The first is the name of the application. It is assumed that the +# application name is the same as the base file name of the single C++ +# file from which the application is generated. +target = $(basename $(shell echo step-*.cc)) + +# The second field determines whether you want to run your program in +# debug or optimized mode. The latter is significantly faster, but no +# run-time checking of parameters and internal states is performed, so +# you should set this value to `on' while you develop your program, +# and to `off' when running production computations. +debug-mode = on + + +# As third field, we need to give the path to the top-level deal.II +# directory. You need to adjust this to your needs. Since this path is +# probably the most often needed one in the Makefile internals, it is +# designated by a single-character variable, since that can be +# reference using $D only, i.e. without the parentheses that are +# required for most other parameters, as e.g. in $(target). +D = ../../ + + +# The last field specifies the names of data and other files that +# shall be deleted when calling `make clean'. Object and backup files, +# executables and the like are removed anyway. Here, we give a list of +# files in the various output formats that deal.II supports. +clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk + + + + +# +# +# Usually, you will not need to change anything beyond this point. +# +# +# The next statement tell the `make' program where to find the +# deal.II top level directory and to include the file with the global +# settings +include $D/common/Make.global_options + + +# Since the whole project consists of only one file, we need not +# consider difficult dependencies. We only have to declare the +# libraries which we want to link to the object file, and there need +# to be two sets of libraries: one for the debug mode version of the +# application and one for the optimized mode. Here we have selected +# the versions for 2d. Note that the order in which the libraries are +# given here is important and that your applications won't link +# properly if they are given in another order. +# +# You may need to augment the lists of libraries when compiling your +# program for other dimensions, or when using third party libraries +libs.g = $(lib-deal2-2d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs.o = $(lib-deal2-2d.o) \ + $(lib-lac.o) \ + $(lib-base.o) + + +# We now use the variable defined above which switch between debug and +# optimized mode to select the set of libraries to link with. Included +# in the list of libraries is the name of the object file which we +# will produce from the single C++ file. Note that by default we use +# the extension .g.o for object files compiled in debug mode and .o for +# object files in optimized mode (or whatever the local default on your +# system is instead of .o). +ifeq ($(debug-mode),on) + libraries = $(target).g.$(OBJEXT) $(libs.g) +else + libraries = $(target).$(OBJEXT) $(libs.o) +endif + + +# Now comes the first production rule: how to link the single object +# file produced from the single C++ file into the executable. Since +# this is the first rule in the Makefile, it is the one `make' selects +# if you call it without arguments. +$(target) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS) + + +# To make running the application somewhat independent of the actual +# program name, we usually declare a rule `run' which simply runs the +# program. You can then run it by typing `make run'. This is also +# useful if you want to call the executable with arguments which do +# not change frequently. You may then want to add them to the +# following rule: +run: $(target) + @echo ============================ Running $< + @./$(target)$(EXEEXT) + + +# As a last rule to the `make' program, we define what to do when +# cleaning up a directory. This usually involves deleting object files +# and other automatically created files such as the executable itself, +# backup files, and data files. Since the latter are not usually quite +# diverse, you needed to declare them at the top of this file. +clean: + -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files) + + +# Since we have not yet stated how to make an object file from a C++ +# file, we should do so now. Since the many flags passed to the +# compiler are usually not of much interest, we suppress the actual +# command line using the `at' sign in the first column of the rules +# and write the string indicating what we do instead. +./%.g.$(OBJEXT) : + @echo ==============debug========= $( $@ \ + || (rm -f $@ ; false) + @if test -s $@ ; then : else rm $@ ; fi + + +# To make the dependencies known to `make', we finally have to include +# them: +include Makefile.dep + + diff --git a/deal.II/examples/step-45/doc/intro.dox b/deal.II/examples/step-45/doc/intro.dox new file mode 100644 index 0000000000..a799344358 --- /dev/null +++ b/deal.II/examples/step-45/doc/intro.dox @@ -0,0 +1,63 @@ +
+ +This program was contributed by Markus Bürg. +
+ + +

Introduction

+ +In this example we consider how to use periodic boundary conditions in +deal.II. Periodic boundary conditions often occur in computations of photonic +crystals, because they have a lattice-like structure and, thus, it often +suffices to do the actual computation on only one cell. To be able to proceed +this way one has to assume that the computation can be periodically extended +to the other cells. This requires the solution to be periodic w.r.t. the +cells. Hence the solution has to obtain the same values on the parts of the +boundary, which are facing each other. In the figure below we show this +concept in two space-dimensions. There, all faces with the same color should +have the same boundary values: + +@image html step-45.periodic_cells.png + +To keep things simple, in this tutorial we will consider an academic problem, +which is a little bit easier: We solve the Poisson problem on a domain, where +the left and right parts of the boundary are identified. Let $\Omega=(0,1)^2$ +and consider the problem +@f{align*} + -\Delta u &= + \pi^2\sin(\pi x)\sin(\pi y) \qquad &\text{in }\Omega + \\ + u(0,y) &= 0 \qquad &\text{for }y\in(0,1) + \\ + u(1,y) &= 0 \qquad &\text{for }y\in(0,1) + \\ + u(x,0) &= u(x,1) \qquad &\text{for }x\in(0,1) +@f} + +The way one has to see these periodic boundary conditions $u(x,0) = u(x,1)$ is +as follows: Assume for a moment (as we do in this program) that we have a +uniformly refined mesh. Then, after discretization there are a number of nodes +(degrees of freedom) with indices $i \in {\cal I}_b$ on the bottom boundary of +the domain, and a second set of nodes at the top boundary $j \in {\cal +I}_t$. Since we have assumed that the mesh is uniformly refined, there is +exactly one node $j \in {\cal I}_t$ for each $i \in {\cal I}_b$ so that +${\mathrm x}_j = {\mathrm x}_i + (0,1)^T$, i.e. the two of them match with +respect to the periodicity. We will then write that $j=\text{periodic}(i)$ +(and, if you want, $i=\text{periodic}(j)$). +If now $U_k, k=0,\ldots,N-1$ are the unknowns of our discretized problem, then +the periodic boundary condition boils down to the following set of +constraints: +@f{align*} + U_{\text{periodic}(i)} = U_i, \qquad \forall i \in {\cal I}_b. +@f} +Now, this is exactly the sort of constraint that the ConstraintMatrix class +handles and can enforce in linear system. Consequently, the main point of this +program is how we fill the ConstraintMatrix object that stores these +constraints, and how this is applied to the resulting linear system. + +The code for solving this problem is simple and based on step-3 since we want +to focus on the implementation of the periodic boundary conditions. The code +could be much more sophisticated, of course. For example, we could want to +enforce periodic boundary conditions for adaptively refined meshes in which +there is no longer a one-to-one relationship between degrees of freedom. We +will discuss this at the end of the results section of this program. diff --git a/deal.II/examples/step-45/doc/results.dox b/deal.II/examples/step-45/doc/results.dox new file mode 100644 index 0000000000..63961a384d --- /dev/null +++ b/deal.II/examples/step-45/doc/results.dox @@ -0,0 +1,20 @@ +

Results

+ +The textual output the program generates is not very surprising. It just +prints out the usual information on degrees of freedom, active cells and +iteration steps of the solver: + +@code +============================ Running step-45 +Number of active cells: 1024 +Degrees of freedom: 1089 +DEAL:cg::Starting value 0.153970 +DEAL:cg::Convergence step 54 value 7.73356e-16 +@endcode + +The plot of the solution can be found in the figure below. We can see that the +solution is constantly zero on the upper and the lower part of the boundary as +required by the homogeneous Dirichlet boundary conditions. On the left and +right parts the values coincide with each other, just as we wanted: + +@image html step-45.solution.png diff --git a/deal.II/examples/step-45/doc/step-45.periodic_cells.png b/deal.II/examples/step-45/doc/step-45.periodic_cells.png new file mode 100644 index 0000000000000000000000000000000000000000..3e19509c8720cf56489f37c36ff9eb2c82f7b062 GIT binary patch literal 3855 zcmbW4c|6qX9>+(rCFDeImWUG-rU_$6l8z;5P?k}aLD`ZS%UH%XzG%KH6Qx*L?ijKuF*_FfY5u)#h1sQh_#-S zX$?y|DgN$?Y}|vV;{n}d`k4g9vZ`9a7>T?IDB;)NfZF+ z?)s-d|9WHFVUjI7yS$BU(I4alxDC>z$N@&wTlalDq^2P9xw7C}5iy;hXpt^4e+3r0 zN^2RcAAIP!Xi+2SQ}%W)E>KufwLGUuUrNlTbmUpd17>^f`2ZGA$AD?3v6t^`AWxP@{ z`P9_1K6R(M?>6!|1P#M8>Ofh1jfNwc9V0WTU)qv=CB=>r7A$@Xm#SoNGc~Q8t#h!0 zQE`7PWhu008`r5c2$3pI+&GGh;H;12s48ZzB>@#1Qf@^SJ!AD`Xin;C_IS!}P+1iC zz&N2;`zOehlTB3HVy`6MV^bvO2JZ)hdwegC+h3Iaio(}=XZ!36Q4pWUa{7jQt>%bC zSNzs?yo|LK!?)*(3(c7`<<=qT+`T{QI^X!4%$;k?t~}m?M;_mI`Ia)k=TqCyyN~uu zV&tp(n%4dm^z3YSwu8Ep&vtXYFM8Q= zji}L`6B)t@%rbc4xecyb*owjW=7%qC&aRvf4u8`48q$#B4!isNQT)7XbZ0m(F&^1l zS`LQIH*m(Rv;Y!41&6S*noxMDy4QyI8jASW&laS-ZsU$2nGjRnN9a~c@0k9udXh?Os4>#j9C4gWd~{#B-Y8?H`hzrB36 zPZ>)mCOmj6h6V4%rIy> zoO}}R3GL4GXI0i56Yr}!i}6Yh7OXon-g1;m?W+>6+{uyHPaRk}Ul|%%K*{rQ^ZsB< zJ5KmuJ{d8?tWKR?Uk<+>VK$P{FPDCZQ4y`RU-@6=aR0*$~$E_Wl=m7^L)n= zMs9NTIlZ1FMs_drM;3cqB4uW$hea(roD;pQF-*87rb@R;$>5yEiEwZ+lwE=K2&6w= z5-6AF9Y3eenVM7KZ2Ve98x(|W^{L8(p{=wt9j(%`O-hy6 z#5B>!^R4=$!lYh$j$GYR`LVilIwS22r)y<-GfyfgSB70?9GQ&h*ndph0>1NZK-9`- zki*vEa5(avtb_yu_8OixZ2@gtZ}XIgrt!14|G^@Yh?9{X%+XNq55`sc+lGU{GhGOF z0a5@#t)gWX2&?)}O`EOsoMdn0;w!+j3meF#8Wlc9MbuV*uA1{?lj(;>m#uL**n0CT z+f5%ecsV6+)wy?br^6Bxush~q87x=2GaJg=tnjx4dzhKcTH#=8CK{c1mgSMV#f2&- z@L*2_oHE1{Y&Uo~?Z>Y|@n5e2cXSODGl`uH-K|A99s^=H!w`lCCVHpsV~S}m^Ph`; zXomhN!J(LHk{BkeBoq&5nIeqcF}8I4Zr`^9TJQ*n$7$HM)RsGC#zyQz6#z5u%_rFD zp>?6J)bv71@JBrY`80Dq&0mm*BvfNGT2$(`=p%~u3eSD z>;z^iE9O!1bYa04ENhiB5wqNbEQ$)<%32FKGcg}kPDULB>7B-q0b?&dD!Ay*c1 ze6%g!mm@tr8Kd5#16^-Fnv=6HN_WV=z{}wALpAO49_=)w`U-Q(58WxS@_I-IQq5xa zD(3m*QD%(7-MI8^gGTPe$cS^!%2F6(hpk`3@<{`|ZvAv#2&GxNyRgQD(KW6qK*~5Z zTJn=<#Z_y_3Y?>MCW5}GF0#Fq*~Fl)KPZaFM#|j8i;dU(;qdXwjhg+>rw$Ljx+RPq zhA0csCNp>*^_ydbJ8_bXxB7pBNq=Ub#Cu0gd!NbQ1#T}P3Hs>S*j53B7V|v z+4oZ58(H%woDO)%>i^q)D}ccltJ981CBaFHX0^oX;{1F&LLEi)8||H&2Ta^6g7qej z5G*o_6C}EH!sd5f0xQeVr=y{yuS4Zzv0zMa zR1l&XRYqaV*UZ%icAOzE>;gl)a|(m(d|_?VuM+tK zl@UYIZYU6jx$S+p9pWL)#Q$;WudQNiE#w^Rn_4s+oux~$a#l+P9dWb6HDRwwKCH!s zUnD_@2HvktYyUySwWFmi&}@-MunF9oQ|bPjU946=!?Ngb)ckJ2L$Rx1eCS7{-M-}; z-E-N0%Y}nl@0d185h4rAGKH#K_TjBNj~v0_4yL8?7s~N1HaCSbnbk8nb>Hv*Vrc)o z>$-sFX}l)4l{0-Yq^0lQ@o6Tt>`lAWeK1D(cfR3~K(k}9lCV&|)-gp}f__Q8{M@lG zWKIcAdtZK$JUokyD$IJ2P~)xS=EFrmJdxTb;@oUgi5Z^RhD2?WNRaATI=+q;VkaJ zu{hiFvp)@sc4lXdbpsm9U+K9l@o@QBta6B6gqbOk3pP9_m9txni`Z$w0(b$~OgU5e6HARJEYijIeUkK|o`j&7K<>o3(?q$pKZvAtf?X+gImc zft?(FSy@?|CaqMJF4(Dqb_V#=ev>XW5p9d-^NuUP6#Q1BLw<=}{T?|NWFEWKH$IyF zJ}`Wi0A)5df_s}ytxrhEQxK<Px#24YJ`L;(K){{a7>y{D4^000SaNLh0L01FcU01FcV0GgZ_00007bV*G`2igS% z5Hu%(Ne@Q=03ZNKL_t(|+U&jAkt@s6rRN?f%0mTODN52xwE~ofKf?KtQw%;vL?%%C z6ztlyeX4*23dsDhmXAn8@OS&W{oVd<|L86Ky#;@_zuVvKk8yvuzuVvK@Ak*IzuVvK z@Ah~5vfEGd-~adj`rm(V%HQqp_ILZexBvJ5{68Aw5dW0j-|g@Acl&pB947pGQ~qv$ zx4+vz%Q%2P#{J#?ZhyCbS;zf()n9+ZU;91vSFHJW(@u|^wbNc;*Ux|I68geh`ukG( z+%2ou%-`+XzS71N*w0#i{THv;Kji=3FPPWsPho!WTjZHb=Vv{n{^@J}4`_e>{||%y zuY_+9JwJC`#J?Q`e?iH=+sJ97T=M@7Q+!eWn-H+rtNvb0;6JC(zVI*l&klgV!W93L zZ^8!g@4SX?%&UGaHGUV){-ugVTvVAv`4?m2JpFz(?3x$<&6?)tf$RTEnDT|Nz`yX` zpe(|H>+s8o|BJv)M8D|NeS~ozh0A9@j6aHf?VJB?1U$k0&1l|0{JYNTCwjR)^^Sbt zdH*`${sH6sX+5=_w_iv)BKnPylb>g}@HB3o`XthC_(J>mdHKu2XPgnvj}zxYaQ>@F z`3;PFj62Uj=fm#%xAhbcyj=v}bz#4X_xx#OzwjfF?J@sL11~>MlIU&#-G`nJVCBDE zqUUKz}@lf9^guYBU#d0ymF zru1vGXCIg!cZ()kN{L&cgMZVjo^g2T^#T8zp1|vV3O&Ve0{YcOnm-@70#^fMA(FnI zzYG4x2fvWNDmln}xc6Y9!RuXYO%G{R%)k3v`RDlya;-{QMx~{H`18@(^qoWM-ETvt zPej*VPvH1I+4Vq_Bj^0vL&HDMziTVC(xMl;1W9*O=HK45@v=V4+Zl*AWB0Y5qzHJo zWBmt$#lP6A&Vy;5N;h)nuR)@DKVFOzeg?ZMS%<6X`MNV8kpG47@dv!>H|VBWHmwHG zU(463PvMRW{=(~`gTEjxq+Gs7%MsDPQNZ-`k&Eq<0&4CLT=o8>&n2oJc;b&H8mhl9 zuPrEBH!@@)H`fyLZghTuj(a3}{V5qWOsFMDdgj~XD+HG78hV0+ulfz4rGMfYfFiCT z=28Jm`SV$6n;z!We9Z42l>L5htIkKdsdWjLLiPmc^lQa=TrMRXJ(YM?RdaY+i=qEo!0y%kz3T&Q1!O0YiqqVD?F|ec2IU7r4aYlwA5{bU z9F|{Z-%P9YxE-#a9{z6K=YC@a?`u_{?S(VA#5r}jD#-gx7u^0<-~Z(u&N(a$mVYg+ zMmRN1rRVME@-MBYyBc2W?d^Q$>1%TO+r*TSNkpC^?$(6ZDK;OoCO-oRTlUhv6Ikiq z*Yk`C{fM_UpY&D;HuiQ}kxijt`sP!tFz6-PHWJ^l+_N@ZAr+JiP7BYe|&;Gvm}( zMbN7ci>3}ThL1r?z9uB(sD30`JVPlHeY~D&!Th*!N%BjX6(w_;_p|UHyuDt1YW2ea z-Om#m`w+tlYr?S(zEcZFjBJ1GHxxN7>T?+SH16|;H!&+ptEtcT6DLRL6Zyc316&5Y z*R=%eIZz^Eh>tjGJfDYScc3q-c3RV~I!5b{OD#XM!hav)V{fDD<=jw=^76gQcg+hy zTmszvFk>e4g#(WZkui$JN#(WV49s8F=(=pyIo^%A)oLs7dfj&&Wc|7m;tQulNPpeh z$&n{;2a$Swa8XtcATmnE0bBjS(_nc-?`yp3{bBgUlT-i%Sn*vnK1*Vf4ifr zpAbh=Jp(vQ&+Rh`QSt>CdzJ_mqkeTD&6P7vA&$QH{QQ~c)a0AHPy3S(ty5Gbv*DJ^I zcl=%?eeoA+`v)!|9!>xhW7kiLlyL-l%Sb}FOrZ`fbjwD!PV@rNO+$<<>!N|>TZ`Hu zHlKcXQ02#uo3ZklRq^`=5Vqhc+{{Y1y_2hM(wp&22iL2+5(aRYN3J^Z9o{d~TtFo6 zY2!V?ZZXTlOl1CEdhC>r?N=SCbX^s~pF`LK{BpiGL+{~Ic){ef81EUp>;=CD@%kE| zWWpS~;O*Ds0@G3y`e~ShRl793>S6mVp6tZ<8Mp$YbzV#S-K_|lr z7;uoGNBu{DjLsLFle1}e4GOdG4<~ZBY(s1c!%UU$W}lxwEeHD?Qv7dQlkB zWtnKOM)(W9foJ<17FL1Xq{WmMID;^gU*0y{M><{?eqCLs@6DDr-d-A`BqEagC7CbA zMMLum9Vdgy!R*WL2g+;WE;T>Z0&%}rja<#Q#`VWQrjik7{^I4Sn{^FdH>fY>6kq9Q zL#zuKSU(m$Ks2k;Ob)u-v2|*t<=K+xIYZwyMPkj|5ocZ4Omu(DRNg_cNqrae_|HkR zg|}`K!<_EJ5p*6~W?Wcw0@0+Qi9-9i%oc?eVvw@8)F}r~LHfQKBtK)c8hT^zl4HtY z72=Fh%pjv(W_&-PzA^_kkM;=&x-VBh3PUzI49z*c^)myk zA%iv|$A-CsYmbIg{BlpqyXy<)Yk(I3g7!{2>)rgl_YNQ;8dG^4j@8nG#(8E!)|20J z)EiSQSR4hS^E+Z5nm@|ba#KbAV;~#UcQ`^Z;}7>1Wd05-y3*Co_}R!RdM6fNp00Zv zg{GdTVI^os1}zABo3SubSwd5}{u27jzp(8-hhV|nY`vI7pM@Y>OfmC(UJ9(03vL!+ zk`^Udkvrr4l=LdK6*M$ZI0?ZtXNPf)uptiaZxs!nO559bYb}(kc4Tdt{#dYoK;lWy zBY{kNT{y@uyT>d&hoqa7AN#_latwpTy)Hh!oIpeeN$3^`ns#_b26|3Pe^csA^XcA} zf4?2^_T(=tLntVW%i4O^aOiT5q@i%Jeio`EMoB687co*7Y>TG-P z(9lF-qJN5jUI+rvsIi!zTr%%fpD5p7CCVTkCX2FcZB>e&M0ihA-LK{yMy8W{R+ze4 z9nyG`R8xZANkDLqlZ$7>X;H{!-mA{?%Vk&NQU4_@TLm&4cCV}IOCCfU;s}pIuRqqw z;C?-xdzuPyHjkM{9+yYQQxIv+CR;lgh+GVcYO%wde&#kZKMtm=`4^@y=UEm784gvD zETr2Kw1umVtozUOZwf_BxOMGbtAAr@^pto_CrdX?8c`5>8VnJn#y><%Uyc);|``#Y78* z^Bjy{yK(<}1O7Y&1h<0oqtGAsn8D?bba})Mii1&rn*wJc63v4zjmAe{&hHO@vwh_D z8p~Z^fR82%$6U-0>Z!b>Poi*jI0IFyUeX4&YB)LQmq4aa6_W1SE3xBj&O2$ZA*X@TwT*m*<)s zDf(S4#QcvH>sO*?d@2j!Tdi_GjP#hu#PJD&@w%fH1Wh2?{RuP(L9H_5v)>-1{YJr` zT?pd%!<-ykz{CY!bxu~SV%Z(O(vLI$3r&M{Bs7n$`LDiJO%yZi!}vkWZq7qG;W!6n z2r2}80Dm;(1dVSJ^@t_`_u98evNfcf?)R#0gqbXR5@+lW@4%UCifIgC%Ud9plhHgT z(fC1@_AronX)7E|aal8fzX#Z-W0BOUy=KE!s2Dg(!`m1)-9|pp18@i9`m@T3vlMt< zTJ;XxO1Fnks$$dg$Yb42sSeE9P4zTDz!Cg?hLOkL0nY)J4*7rh!2s(*XWH#lAsr>v$lMw{;&rz1sX|_(B@1 zdYxz7LV@PQ8Q%5>`rq>av%ai+@0wRL&^2(r?=BKrsXB zBcn;fOyf%MVfNi_cnM}klY0O+vyyzduQl?jlY#IB#k6;7G!MMc{ZdnW;e&p0xyXf+ z*CY*5YW6ftk;PiuT{C(97#xFOrkn8>^SsDXKr*T!(WjqR*~~(qjjDka-tqc44B*n^TeJF7l*KA(Kmbl(mSVWd4ZZ=< zr@Z5&N#oJEbQLNUij z@V&W(!PkafIHI|OL3gZi4JfI<82e%~*QOL-6nw+%J87KU^HV=kL-I2c0~%^XpQ)`~ zBGbl*=2d4wmdVzBF~2@onH2J{{y*H_siIcIMfIAb{zB}V1mC>kQ=V^-ec+1AQiLiV zN)Ev4XZ)RCw3FeJ#?wKxVP$s*$Z*jZn1^mSQm=24(SG-IRT$K3GPa1n>d}7$*JUL| z`a_D31V1E+KAU}D&zl>@IE~`W=T#%9OC$XZgC%Kq$h0w{^^aP`GQ>65`kpyU$HL~& zlqm03fk6%m1XNrwi$6d3_#-GGGshJ|(V62x@PRsl!DUWER9amVr|LKi<1^HkzKUr5 zp$71@ifAF0Ar2z!w@iqef8{R^7A+JKg$XFSU{dN$lSa0}EdnKni=bcbNZFyy-6e?_YUQm92= zXQKYLC}J=N2jQ{jz$Ew*u;ik#$4YW<75Kn0g~@ZGo}b8PYe&Y%H_lENzu$YBGR-6! z1HoyQ8T&c!JR6?xbBO&p4kjubl-^hoe}exdS`0omHQ@~K@o5czbT1A*^H(0PN%8CA zHpWf3FbK7;$7PxgwV+H}LNx!!CfzgP(~z3%P%0lj%sb};?eI3xLj1A*H+vPE8+=nn zb8a`6Rp470LUCWq{67WK^2m&8k86B z#GKvR5SM`x{@BlVIrw_&k<+M(*qW7lsEIFGzO5L+q`!eT&K-e??d~%@o|E?Td8qMM z#fpJQ3p)0@YyMVKkm5_3@`_hvwBkq9AIrWO@O?Y=;^5OopPzv2`*TN}&dxd5r1PI% zj_TRFV``I78y09{!-EVmt7$)2eD^@)1Ht4^Lx3*#?|_}HJ6m#E`-{EOiM!(LBi8VR_yMITI>as3s?NQu&ARj$0wh# z?DLDh=c$*uvJ7&kVi(9gnVwGkm2r)a=2d3`I>XaA?u1&DYt3)d&=|`|(>@JQ%+IV3 zG<6}xTDTlcBt7*kzgr>v`9&XVb$jr6sF!`KY#E{Fr=B}D$5PM3J^&kklV*ZOzvF)f z13i%z^d;>C_ldPhsJZ{-7|WE`oFce^Pu!A^Tz)3%@{qLm3Mbj;RD8Y9b!xpnu9to5 z9G{_HA<#&Cb1d%~jUBV0iD5pfeGWqsW-PY=bLoN)LQVL4X_94Wnp~epq>Bj z(0XI`b~vIZU2cAdYpXMdGc96s51pbv^DI+d@j(0MGm*80Q1c2jBtMC87GB@doR2H- zO9KS-vyZrZByk^{hk8qd-ol6;M7nW9axdUaAMC?dy()}vOA&UY&h|^)A=VZ`%_-1? z^lOx2o&f=0HGPE0J;dA`FYV|G_py&J=^)SJS5eO;^oG1GU$8PqJso_?Nn8DyOXfZ; zU>th&!WhxMRIJTF%^}b*Cdz1rUkSeJmfg%Hd%o|$FUbi$p>Kgks7J`yQV3ITiMN&P zEx3bvirfQhS(nq05@6`Vo?@IiYiNLIu)4cojB9?e7J?dA?mNY@nW_L?d7g-#iM-`P4~9mKdhpjA{MPf<)#os$|A@;ez=ahRU^?IRadfNfytrdF1;PVh8KRZS4 z=MJ2~ists^VlCi6BT(}Y&6#G^C}y;FjZIh@MBdPokq4AI4K)agPTY1y64ODdI8@2sY}gK9D)-5_{v`9Qg^=j2K7v#m;2M^5pV0Xw-6kz zr-1KdU&xj&r=b{+OF>%ALv4u(9TRJQrKYDNJUH_gk<7T)g{p0>8+C@~7FLR|?-K() z&h7d}HeyKV%_GnBw!=a%NMPiKxq2?eO`Y|-o z(eVf>oTNKV4$=img2aX|H&J{LyS3CW{ z1k}bx;C{-4u7{dMG%wE9S=Xnvn@XUYJgKDw$BRY!&tZUiov7o6Sh~Y`nXCN3Y2^8Z zoB$4-V%=sy_$i6k^UY7FV}n~)agL|P#YoIeb)7B09PYy%9AP(Bn+v)!D5jC z>a9hdsq$>(@!>pM=s|FbvdO+E0Y9MQW_A>;FIyWip)RpD47DY`6(%-3Ea`u7HRZ&P zhBGLXASmLGb*5hMg{eniWKnNbc`oF^U>in=dL|;SR{Xej8;b_Y@wgT1qa6FRW-B}< zx#kpWfe>EYw?bqjx{0u;d>g)h9Cj`>H{ppLZ2%9e|y~v|2TBqLf27`$n8;G=X4x+Ys`a zAKZaFr_fV^dJ_1CwST%huR45DnbW=Jx#q0!h;nTnY9`UV-CaQU0d9GmE(ALhbNU@< zfdIsNsAoW@*|%-w6ULE8fDEcsM>J`By?+2;fuH3m~s~Y%B5V)0haGJuTfQu0FK~d5)xQhpt*{s5g-GRJdbU%ZRHrC7zwT@F5R?A7mna68(|+3aX2-x9FSz8=Wlp)M1n@O-ndLt9y`Ok0f{$;JVoe%ndQ+F zql6)5&FZ-G)o3dPJfk(fkZbY#Yh}YTj3XE0#0&K+Mcnz?%%GbA!pLhU#8rPR0fL$w zr-(f?v4u`nc6N9&iZ_V72=x?`J3_!4;&BGfoKS-mP46MJ-OV)|W}Jg){DnF!;y&8y zU@JKqQi3%83Y^MLe?seubz`f(Rt!tKG*LT_JYecez?idSbb3A3?%?Wk7+cF?&650p zWSikyfN?}&q6Ku^7mK*-kHKLGiW7Zt0(^1m^{NegJ369~IK|S=gzy{^k8!oEqZSrU z%Q*)bd=M|_oosRdbM8n8YQtQcLEG@17jjL9nuloqdSGK-=9lDBC(I~D*f$P7oDzB+ z^7MxfM&gEf2d=CIn66fO>GHy9aT=&OTqMtdEE11X9kQo^+{iiJDcP1-Q7_|&yil`= z7O4j|Ms?(`KpcTGI!v!136=YL!ACLV&AHk>^yZ4h`FUqNaKzFY5>F06iqzBO-Vg-W zF%Bt2v>C45D%m7#nUc+ELaTqPP#dV8lUQ~K!jxnBMjVdOHbrqm6X1(duSsPPg#ObD z&tcvfvA;P>DWP(GwvU+WY++@XJdHh<$fZ{W?NC7ak2@mUR4BY&ccxa zoJAkvs3Lw>S4v3QU>bd?$7)ul%dVqQ+yn zqn}Zt8Bo)>2IIAC^ODVGTwJV$pvI+H^E|@hkvNJ2F|wY9-kKp*oUHg~CJ%i))LV!= zUDGMNdqCfL8e03DgAaKi)fei6NpSAe92*?_Ge$I1tj%%_79;8zUNf zAC+SF4pu7Of)M7EUxp2`@V=!1>S09X5n%bw6ITDJydQqvJ+JSa63>G?m3kWZ7~i@k zIvVh*uC>RZ{RR)a(KV8d;-z!dVe>kitr@R0(*j{v@cOtrfIpUf+`DH`&#v_)HI_%- zDBc+F5Ok5aBz_l`cnb0W^lcAzy7}&Er=O?;h2SzOZ)f|b(Kf@iE^L=WP~(Ue>(38i zYo@-05dF`kZQBguj}ZI!%|kuHgzCg1y$XbydBdv28{nN2Dv`wMS zRBchoCZWxO+GyK@4zdWjhT|RNeins>3jPrN;A=Hhul1tH>q|y7i6`?82t4sSGQw~) z8N3_FvGyv<-`Axqe@O1Y#2uTuv_1#VqP#%=)_C9O4XTYI`HEStSx_VI7!n0!!OD&y zDTje!5Tw1-6k;Fxsz$i0BZj>8VnmqaoszFm#NBc;>){=AcQ88~v>G?kQVnI2;qK+C z;Hi(BZQ)wc5<`9wJ@n zg^<^}S|9J4_#HN?O)#@-#hn?vgE;TN-^Z#_FM9gfGi(MJ=kMb3Z{C|no2}qJfi^{~ z8J!Q!wdco~CXVh;JKzEuqEH$_;`@+XBF-)u(U7-_cTHCG`rLhHc8R!yZqgQ(zc>EcoG!4N1i6}hIkh?v(JmW@<~K1s!BYN8ymi-olEnouAYaXw8JX20h>lp(I&*~{9KcL zsKPZ9YRnyiX@{Hp971;&jjB>n!F~iR_?Y)9886mMCbUtr9OYdj!CYqcIdK;b-hppc zu%g|}i7@!F+7YLA7^)QUqm9Mfm!M57H<87fDbu{u&)2m&#L^C!Vlyg zUg)XDmgH00&$}vm?htoM<((w%lt}?(OE%ujf@8k1OX3xlQ*o@`$-qnD_cm$(03ZNK zL_t)8DDGaa$V5zaF79fy0bBf`8{0e(W#eJ$pKxQcI5UQwyC=eI0j|Cqi z;<~3J6?v_(RVCg4?<8^8sTezQH&%1kxmQpGrpdeRr2;q$N>u8#b~P67UFB8h(y3BM z=h)Z}J=}I5+HA?z$}z7}>)_~{k!chgN8M_60&5KrgX}{G^$3s=ajMEIjctT?qMvn& zJK4?nK|LS)?q!@FUxV2C=^ryn=NX9ut#2_@~HiclK2^*+I3LLt8sA`h20WyH^w_~6qJA` z!Rr6c)O#M>U%Za`I7g7qv%hM#jno%z(=()Av`KE1MVpgrCe)lnQ_3<%Dndx42aQem z$|3d<3V#s;J}!A`_M(;t*Ck$!AthMmF|)V8PIa>7N2ZHAP>V}gSFQ6-iMV+N#)(C6 zi}I~M{xF`kC)8ejV9qsq)|1Yuw(EzZTx(k+Ttu_bGdMJ9kuBisO9)W{nHPLa(!)R= z)0pMvo#V~_kI^Bt zdAVjPH80Wh%)=+v)}S^)Y!ZAU$WtX=fOk#4LNK#W!cGs{wQ5d=o&P_byi<_Zh&4}c zA-;DgqT|pt24fs;Fju>~Qn;yUgbB2fSFG71Nu*HS7j5U#MvNfNf{)XBtt0M{*N?qa z3U>wGdClyzuxnv72Ro;lvv_A958SEMD8HuSah`2>!)W_P;WmIarSZNhvI}v|OSFY! zNG9Z4QeKJ3J~#Mytyfdeq4GqDr}_AsR>4 zVr_2lp+o3F4oEzOcTO|=-LUi6*ieols^;9t<3Hn@G@Ur>9D}6u>{njfcQq*7e4ca= zZH^XK7t~llh-%JpW~?HeNw*zvGzmVo&0ZIJy2Mj?=Q6YZ2J9HU=5w+GYR=(RKT!jW z<}%Ru>iD(7Z6)27pv^fQGdf<$+uapG*d*y9%N(Y_=bJ2|_tkrO=M#4yft@n(>2BCH zrzmS?L7YKC(EuDx=NJ?XIG9%rFOZLqIl57}aaFq<-Mnb?aLs_4lW1$^jtrzDlJ_Zu zk(M&x$fll&yiwkH#hsF-xD9rStUIp-Dsny?h8(a1Eg>uM&Li9i0?n$oW*&Z=_7kkx z7?Z;7VY;nG+t7q1vt>0T<=}=!uJFOr3cyhac{9ALYkehgXQwHaD7(5QW8A~WlpTOt z;jHoQJV`uiZiZ`cK8$jC{uK9?j0Voe)Zt|Aof?mkua&G{6io}!kodi2*^0t~omboy zLd>0@pz?$Xde~N`1+x>Nb>L&4Rq%^D5AvK`bEY6hch;zGbO3EIH5B9~y16254BgBT zauV7!t~sEVzpramYfuFMWL|TUE(wRjWgykbjrwX7drn1h*OrMHPIku%ITI8Y>^!pW z{Ojgq-4N`A^0(&gl`;&R#FM_N32_EqqdJaNDC`seBD&on+|Iv3;g^WH9ETc>fb3HR zz})66mU?_;Trcz_ob|6~Ox$^3SBwjFB`69H3uV8-fUN5}7T7#&ADHW}GOF9}OdG3( z?IOJ6)tvAAG+8;wI1#869QON(l)^Gws(Vtlm~M<~9;i*#08h8Cp4W^)H{fG!9_l42 z){_IgGyXG=!!D4ZsD;?Ptdkq>**vV%!%DE@`V;v2I6v>4)Z$LJHV{ zosG44Cx$QXN-|iCbA;P|5#2^?DVt2t=;lG2zmz1q=(Dy&Zs3;G%!9j0%SIk@SCp4` z$ql$M*wu2$E9;DCoyWk|&npg7(ko6p9@a`wfQNT3>d6G0{8skttQkk!2jXsa&~3eN zv(Xmfn$xuAI~F_moPqxA3B6*?gFLGhU*ert+(ls5`9wasX~{aqz?Ka18d&Ltfw$ts z2|KW|u_o`r7M1R-9r2(oB;5W*y6I?hmdFKrR|iVuR-Y?YO`L$gG~SnL*Tgm9q4TLz zSO39^yAbRgvTi6^H_|ep4V=H#W{)ATrRTWz@t;;?C)0P1u4Wh^sAa^dz5Z= zR>7T85dYy5ZhUHDjB5_4I~Fz*6Rwu7BV)`8iXz{$F( zfyItKV|fXVbbb%3z|PORk-*#nedp>|&7Qw4Ot;I5@^k}7pQ%w1GOgXb!cFgjFwiy- ziCem=W*}yFowCAM5f8c8_;fA1IztLd-$Hw56LzV^$MNSRhi_$~bAs6o8F%fPtVA7q)KDX^pVh;R1q>9cXJ6h(~yo;$vL5>v+0fsl7z-n zv#nR2EQF&V&gxN(!_JU(snI?y;MyN9BAvUcP9j}pV5eo>JnVoO05#IlB|lMoFP9T- zAkX-|wUORwP9LJ1=}DcE&BL{sVBAcC^*sOddJ*27X2aXs5|7h&z9uaZQqA#sS;rYz zF+%qU=>Ysmi2>5-vaZv-UfA_W;p$W#)p7Z5x?S(t-R4vZy_s2s5ZxqCYUU+cYB->l z*L65prRT0AOZXGXmE*jV{j48$qO239QccvxMsd64IPl>cb9 zj?)j}ryJavIUNkW-9)z#+Oj+p(^>0#M=w|Cv4{pJUFI}vQM?qZhI4kLS{}BF*4YNu zcCz_bkj`gd4G$X?cOg!#)NxaEgIi`;E;gq((amgc@T1KGHGSjEZD$N(GU9P;hl9(U zW-ss5YBmzm^xK*TC8y)_|Qq1+g|H zXV=GR-CdJ(3(GgPtXE2yI@6{zbmJ9S4ofx<)Myl5O9mpfPD&#{ zRSKuod}y+cF|d57&PcgBbFcD5leMI?LUqt3QY{Y)`~VVr;`xoiIIw1&`?hc?p*NWA zTQKCCw=U#4e!$HB4*Lbshe8M#7tv zaLPI`H#rc0F8K&`N6x`s4!weILHVZV7Gma9Vxo_3%+CkWW^k>Cnn*O?DU#ieASK^Z z=*4k)mUoi4bHXl3b-$dfGYza1s#_AL>(db>8AGjL?u>*5%DIkVxm0b$-|h$n$#+9CIf>lXRPuZ!y5l%Qt7}&32|aD%`wi z^GsOEkJO9Ex{Y}Om-JJ6uyo$l`3f6$nq)AgA%E(KzIn~5Vp_9JH&(a>(B|WsJqCDO++m;5 zCm;0~mB&TaIukSkJBf6atTRbxh3aaQaKEjIDOC;|I?f8xAu|U_+j0%uAQO&sWt^_q ztu;bA5q1{uqTut7D)||so7bPhqTVdGKXpalobt^yr>ZR#g`0>r2iIsR%Sa1`D91#p zb_7k{S>i4VJJ!RtAnOb)AFAVUDbt7fZ8Y&}c-TDL+aF3tXoAj3IuGnD-VMNOlyPjy zn?$z<{AmL?k3aR4hUoIm6MD-O6K+n`=HQw}G#Vvc#H-3hzb5Z2ao5hjbimFd>l(`V z-K#!ENet2jt?OXcwX&{MZu_>a>r&E1Vdv(Z`y3v~ImlZG2Ibq0fQzRWJh_EC=;lY8 zhif#KH}Kq=1S`3yXY(#3?qp@>l67~G&JNNAt*e=fU>mL6f{wGU1&b}x>2bOc>}=kt z(UZWPwVZ>8Zk&93C;;d8rylucm{U2iP`?=7mU=vC784NxsKs$c{hDz0yHq-I(!8q1 zI|Fu{xRY(HOV&wcA|cW_gLGl*dW)bNw5}i~UiFNG0Xub6T&()nK{vM4PRThK=5$cL z%?IEn{izj$hioJWnIBcyz->uU`JA7%{Lm;2|lE$oP$8OoBb&dxIh5TkzN>(Z?klZ zqs_rJ14*k2$Lq6ARDUe;$OpTSlkF1}lB`=oI;n+AO%Z6;bs$2=PKa^{Ixt%!mXMBv zow!|9<=p_hhI0!CToxR!1+fQts; zTLZc4GjJ1)M(sXZzE*d<)Qxv!xV9o5?kLY$I3F zyef}8#O9ue$i&?!4!hpNO3}I~>G*oOB>XmRU1Neyc&Q;l7s*AeBAqGgT)Z0`qRj3MnF2k_SpR)2ztqWm7Z`)6Nbn~EXkZVqyg>v)e@sny)W!{CHECV~!!#1+6 z3)P9F;}Z^Mkj~1w4hp(@oO)2unbx&0qho?}*mUN(NN318i+54*`7^R+&OvEeeKLg_ z2b|(h_2ell-$=?TXhAnVFQKTmcHq*Td%2=V&{uHv+6~(7s|BzVWu0YUSxv`QnY1i( zSXU|Qx=PS#bnBjlZ~vMRbc3XWzAMX_bAZLSA1B~Y6dE6?5t#>^9ewl4w;{UCqRq&w zt|$O*&O#vu$JFi~e*y~QDB9mTU^h#;+P(5sCQ<9E2)g?)2XTx|VcxBPzC& z&D#OLy{uEicGF+dIgoSk`O|eVxH|x+buVy$WBsXLzWL}z+5})4jw<$NaX!mXbh-ldgW?te{Iq^Ogpn(zgzT=X3+hCL)|ESKx|EBGy$WoejGI z-uce_pZfqOS!`Du!`*=M1mK+VZGdhP+FV-)^IG+0KB_wL*fUaDzQw|Bz{3uc&SPDp zg3gn2^)y9%COl^gI;OZp7jzEN*|1yXKOA{gPXG=Ul!+*S>&rw&0LP{sf&n;PzD4P# zcv5}$AotlD=fl|P!hMmlGt6z?$*|KsYzw-`!^0WXPpEFhIM$Sbx9H&HZt@;9#j#IFY_ZEG4r8(&fUg#dS zMPW}m#&0U5Q>|-M(0wiD0EWVIcD~gTbS~1}2RqLoCePdvVOqvQz_riy2;ipuDI?$J z>E_{@N;F#03e)yH2wINeUE5V^%R1e|iU!uJ=~U}#3c5kx`W;yaJ3MFSTfyjyaFNar zJN}dg3_;Ciob69tfKzRk)v7fMIHd^A1YF3UI_NflHYeAjI3uC-1QU9WjEIVsYn2=e zJ5JV#u+t5!)4l4})f9AtzE!RhUS1>&@pg!~Z?#gcUecMc8y0Jkc-%6;IVwf`fNSr( z0&q$IE+pU9(XC6>1X}TFqP@Y%KBfYWbE!P+9I~#}ZlBeBEaZtWAT>BhYz;Vc1(9D2J=eZwnN|Tm10OyczPG`!GM4v*N zgJ@&OL;YD$QeHCbmU&o~bV^1-x2~3;yBBj8?~0m|fgXOCD-s4rSxGEH$0c4tTVbQ) ziGufcGmhkOwPl4c#Q~?MQ_GvwY`+|Ed;rcX--zGjFhn=Xqg-1!i-vUem9GI4^5U+> znI`Le9#$iroO`vbt0CxIzLgEn4QC;IrE}ov;n#zJW4=`vbf$H6z-}IAKE{!}MWqKE z^HPGPaL!uaFyM?D5jz0qh`wp&l$pz9qs_-P8gQL4Etm8Sm@V%_aaW)Eld>*AI^DVk zBXnUw*H;1$WFeFSV`s7z>|O|PRL^+@okqF^u#0?6liRE6@~!we12}&v+}(gvYDCNc zoG1Er58cSIIy!(e^PUOKBNyN!C+q^UPKTXEI#-bH{Gl14Tj*P7lC9Qui1}nIILb;) zLFXWy4|a8WYYb=pc$`DNEt}|jGu>8cHx6wcu3c}QY(+`kK*$g$uWnX-19x6_FQe6W zgbW;upQOg_R@vA#H`{LH;+G0mQt`Mfs~TdEZ^S9#i(fobUB&eD<9szq*R3M$j3Vi zc3AJ?YTaus>mDMVv+GMqwt}Ox6?`pUAkUuj=39X$T|eUdybchgG?X6wxGA7)&yPb= zQmDZ;hCceXoNi{}TPCezA&NE!*CeP_KhH+SSGmnM(~EjZ1w2aSA#DnN#$MN_jPi%%FIIc!frs; zEg@YQq_YJbmxy2^bnJ-Dnbi>BC@TTAji0fuz!#S|N$27nOcKql<8F~}adWzaZbHOO z6mCJZDO~eFt!{YDCd=&JhA94$>=fSAFS23SS8xqe=6{uTwAE)d&%4UZO0X-euSeFc zA)RSmV}g#0&`H&EFuEGTu&a)gYn23$gLF;<%e#h&DRKilYgD>{y^X;?soXZFzPOu{ zZWPHRbcU2go5r;$)bjWc;Hi?49X%*~hxVxkz(#47T&u0?g}m+O3F{;*@w{UrcP#9j zvTllWCAKW<8WMD@C|j9`7#DOQ>v~YkmUSZU?wP8@GwxCOW`^Dxy7dvUhHg@$LR7c~ z(Z+C%gPNCU*<1Zf>(Wp3uBRT=5Nz|OHTX*85vSuWB$ZCKLSe}}FYFl7ak6fTbi>xw zBj`Md2-osfDG>p7WeI?{hEBGwhOA>&N2|OWy&;9XYL#!8igz(`79VOi7S<3HLw$K-(qjFjcaP z7k4V`deX77&QH3&K|->wWhLE=pz|amT+3T05)t64p;N7^k92}3Wx?bDAG$AsA1~b+ zz}c+}ZuzFAO&xUG4i%@HA>5p3>$t{>H3@2zDabO{2dh936FMHwu&;e!4)AUhcL*Gn zxR{*;J0@V~G_WS=Le_PSpz|dn=7utXr?SM6a`ll8v`}3n0EaVQ;@#%n9wyOwlp z9(I6qwso!g!0cQEvSlBwa8pW!t$ZseZT!$Qs0AuZqzGM2Iw430-ca3Oc6D*ifk`hs zB;R_vDOm+yN`bD2M}?aYZK_!FL5&QMaX(xTjC*&pL&}%Yry8XC-ozMNOeflES$7)P zLfM%%wkMq~>jI>+f^^|rgpzVqN`<9-tG|sOU_23_dss6_r;@Hlip4uk-KDYA#E&N5 zSip&`3qko7F{d8w7NFZYwCP+k%YZ?z1QtNHJRV1!aO(?4$icfyWYwg*$dbEgon&M2 za&raOz;9qp(s_b(u3Ur=p+mOK-t?_acy7@!VNfbdJc3RM()mc|;@zlTi_VX;{OLmZ zCYjSlzOju8Pti?NZE>ziQ1is&h;A$WWzU}^mW=4UYe;vU=_QX!u?D%5b)ARZI`ait z7a&~$8qvC*&XfU`6`?EcHt1WGEW}t@iQ-$qnpq+Xx(MljZP=O_lj$G3^rJ_$wfVjG z0j`O@InAjOcWa#K1l^h_n~F9HL580Io`!UY2Z?IZRj|`#okqI- zC+9pntgrt|d^G;<08Q}_Vsa4-u=)sH?OPSho&N}7t_<6?U~*utzM^L_2YpX|)7{0+ z{&m4lJ4F%4#>q*wBLtq2ylVi*H!X<%R8K54@{NnTP0_6uZe4Z)^kQHsBA4~%n)a0Cjln~;C%j6mv4DTY4!Vv}o_aOZc-UnyNp}glB(x>en$9GhZe4jn=kcu#=A`i4Fy{K$Tn`Hh6HFe|agi{b6v!I^oLJdQM%CLX-(2R@MK{)#dixxNz?@hp=PUHR zMaM-Waa`31HYo*xbdV}a?um{*-6@c+vYf?5OIN{4WnC*l7xJwQ=F0F~RR?z)=JbM+ z=P_3)I59jY3OZ26UXIDSz=iOmxm2bYRL?VTS=1T8b^g>V-?-euskyw`O|zw>OE$6! z<-n+qdq~GgB2F3&fl1m>VS}Ek2+l#esGzGNbROT@U{1KDsse1kvGgXya>2WT;uQ z&b5zlCJx6xQ4)E_q#b;uD>2tb&B`1?0U8(SR6*zRtqSIfEW~=um9Sv?y1t=p~gkxiiR!dNr%S;Tw<<-T0^=%BZ2kigfboB5OfaT>dQiSG1q7=PUc+w zUY;~-3lw>0Yls?5d9Z*vf3npn=v=aH0op{yF)r0v>&s|8FW|ToYF8kl_)}+YAy;=( zvI@O$W9ddp88C(V9NJHd#o5AO93rw$(SlFTr81Bzx*9GTrc3ed9#SphZ8H)B*0oK! za+M`I=87yt7<0{f!@fKi9hfVL<5^}nR!y=Ba}B4u#nKw-24H7h1sK`OC2g6p5~+PB z!vx@5{?wUUIPPZVOqHM;LpL?gAabqvd5@dhgK^0QU&A}<(+fV*34#uknh%3?te`vp z#bA!hw-(_!J=uB_=CsVIb3(~vHHimvUZZOGR`8Vj&V7Xk`y3jV$|h2m1Fo+TsR7rH z4C(!;U%rLt#@W*9v8rSvc8zq+wWot|zA&7`yJirtNV=j}q{rNqYLHR{LXgfCp(}i= zly5D=b4;=|j=9RDN|-C|T7J!?~KZuxvuu*jgI+v_lhc+dznq7H@`vBLcP<{T? zFW;nkH@2*Vx20Ubj_=|+e%OdLzi(Zu<1mec;m#L-!LD|r=ht70v!W)h^sQ!!Knl{? zg04?EE+?wQtfGY!V z5&71XhA{F?tal@WZi=Q8B-=>`Kh*A<`ip{sYr!Y*?u1=Wy0$^J48uvj73?D3@_|+6 zp%ZjsB0>(&t-%~G%^J)BlV?%4nFhvE1%LZ3)2vk>KF0LMT$zn`hSn-g&76&?bd6n&Eu3+H2%pj-bdC)yOQ z`FrnaG_{H^ln&ch!!=>x9y&Vdcz;^|ywzfYsXTEe`qq;{f$&@$bH|5*{eRXVYo%FJ zO;Vb;SYQqePI*H0$S(0ys!1aGO&xO<>29K%s^i#;LF~!z6REWT?Cy3^4C;laR4fB* zEb1KK)VdH}zA4UBqMLy>S+DU)hWD(XJFiLhP!0IXC|D0WC+VD+W3-r~&HnPOLzrU= zJ6duzx`!23Esx~DoKUe*f0~?Q?H(yeIY(g5%R3m9a-pXXf^_!*4x~N@?k#Tc4~$D? zL|pFKSFAK;%0$X6Y88X?0In|XRRE{t79?}3(9JNU&JwvB>O~~s);}Y@(1x8ry7Vec z+yLePyfq&Km@5#_hyU1u68>sz94p61KFY>2&m-jwShx$~`luXxTPx*s%&ii1zLA+f zv1$|{ODH<8YwV3TV z)u~6``AQ_%VPTybI#hzr!rb|WtjFB`p$aL8`gz&^Y#pl~Jx%D1&#}6t9DmQum;Ctc zs&Iv?$%lxP3)LY{3OePIaMj0gVL0q%8nL`1M*U@~Tq@Kq^;&W&hv7UH6|@488gOFr z)Fs~vy6K*jZ|$nl#&zQdL|M$SW(|1~>?*Pd0l4iAu`yRPX-O3)ww;|qYbg17+A=+XwNUXf)_c} z)?Su1=2*wt0EU%vovGxK51y!;rt-Ke{3g5V%sglW7?pA+=0d(TBHw(BBYF27pKYi` zT*RgFVYo)1Ed-pypDLlZ9(Lm-@w-lzl3%}DS1iP2|5XL48^>d z^E3!5)2TYDt5gEMGl%*Wo$2xCU{Y50BB1?mXK8p>!p=jwK0+s9uH}e~xmJz| zm}{eH9Oih(+CG+Ztiw`H%cycvj=jI9?IC~va<;AzGq0KhBCe?y zVX!E;)UMpO%rj7_?JU6c-Sz!LJpR<6+wrTuEuWYygqLet@%C(G9D}@ucXi7>Pr9?* zq{1Aug(D^An)|mUBQ@sa6a=ri+TU{d2dk9CoqRB*99t`=6p^s?z;@v1h7vP7^Dv!J z9U$Wv-2y%NO$WXr;@X`onTft{AnRNz6Nc;Z3;4 zw3%F+DT|XH+oVl<-XR7^S7NRS&tc#B)CIhvm}6n!$v@~hRx|l9DdlX0DVc}+?5Zy1 zJgr^#C0kbm&W^`1B_kmb=X9xM7_Q4REC5^|Yg6T$PPYbdy(eYSrgM#M>%%{TEvmmmPAV0!Ru1xJVOX@N(!|HT$^49%D3Y_z@<|ewCP+^a5ff@ z-H}uu78CZoI~}M(x<=5Qm}A0oSd|C&fjJRmLzrU&YMNtpBp=44oSk0qXC6lEYWH)+ zle#kxu6p3HjOtkKBA#%54)3d9Ly$|HzX-)%KyqmV){@o^TKU0s+N)q-Eas5|We`a(WU5PnTxhU&8 zYpVrXT=PRexx@oFQ@%A z`hIStoTR)u%~k(aF!JVoRKc~ntRL)P*Jo}~xit2+^w09(Vwwo!%mb0!pU*y=`hr)6 zeNCzX&Nx{n@!2cwLu}%;A0pzhs7`-6KsOU@HqlPQ`!Y_T^H8ZQXvmC{Z9AoE>)K;Z zip`-^5T^#bJl#=x@*#jcCDLymWckdD52DRYpSpbHSN-`x9jAh&bDmJkHJd~Sild8w zWsuNmz_Ic!nRY|**Y5@BV@GhR()&L3YDk=5RC~YKi8$i%R|&T&Y}be3SmcRW)F9xl zVeehg?Tml10kVADDr}Yi<6_q3|7^uqAE$rDssK^$);RvaTz^A= zdymgx_GWR;=l*LKS31qVc${HrvwOVySu&C{^?HyIaR;;HR0HEu0C3RCKZ>kS9N)`i+9_T z3Z!cybUNlr$J$5bdg%0a)fad3hRkN1=}ZjnQAcdy>R$%cM9%tW?rDT_70IY_*@jZY z@fUU`$)(mEsnu_#Fr3q(n*MZ}_)kGX^XJmw@R$F@DVVrq_n zokQ+^M{4!!72RSz_ze1$1q_2;l^b;yyFXOfhE9%&BF^PfyArtudCLJeOt+k4P*d)} zYLI$|sX;nP(Dj%ziaI$dXNsT+yK4DU<&``o(0r{6Lvgpi+Zr90ZA|VhToUa12Gt=E z=W?lB7!I1slGT7KQwn_0t%~{^L1bsyviS*FJOmVO*%j3yi!iEs~Xn=g>rJ$ znJmzjRe<<3AmT3$0@8h2>_LP7~3bE&|U`<6KK z0L~1+Rn7b?-PWN^A4Y;NRE~zl3I$Qty~&9=`>N7byac{VwsMQ z#0f67{H>8^a08CA3vf=l&7+Opiip#whc7fq=fT|UsTj4ZL4?`L>-%#HcpI=DAllst z2gDgEm%Kk;^`iT-IOV()%H`VbD4|%i3B&c2Y3E$lb3}##R~G0}({PsC1zqhhk~i}| z9e>K*+Jg(Uy~O5Q8}U@v43=*ZzjbVWFjRd4%F72#ekB9|XMNmuG+Ppuce6oL20^Xh z9ecFqhwYz??aw_I>~$aP4=Z+YHAfBs4srkA0OmYT#2xu(IW?Ukvoxx@sl5g$IYk^> z4OfV`5_#?XjB}2NWl<}@S@Nw%n_sC#Oe-Vq%H~t>mEev^xh0}+#%wwtyT3^1YsD0} z@Ad+cd2hUEq8u8LV_9<11Uu55PlPTAdfQa}fLkJxdNi>=8HvZ4)SG7_+kr0#3Xh0` zx{D~k+ss2>>)UqPb2TC*3LC(&p|?r8MY9d%F&olJzSRuOeb>#BY_W8~rp25OQAC!& zu|t3Igp z>*;`#L_+w6sIK`!U zMBGh)D+Yeufw`Y$=WJ%z(@?cRHLUF_qgtcQDlRGh=KSWC*gL{-^&NP?newfVxJ{x> zM_y055_9Jrv^f1CIYaZ)H<5phx&PI9Rx;5_mzL^sa@X|_GS zCtZm-R>4V9&SO`bE`6qaq&%Bm>!BDjr?FsM|BwjOy+*aVvj#kud&Nl3!Qi92ATA0x z&-*_w?mUbmL|mJE5UPI*G*H{ykh!LPVC97>I~qVbYs!xftvD=hIehM z;mw=Y9&@ae`{2!qu#sd%6^wh(sD`o)P7!C4f?=QTQhmWePYO7MHtA74>D*E-4gKl^ zx29x75piO@i0o3?Fq{W)#O{$O2g%h{9bByWCwE4N# zKguMXggN7;+Eu$R^~RVQx{#p*D^xY28r7I!bY>eWYDgl^rDE&{DlXMjZHoXG2#YnC zbC4

~u()^<7;44hL1lF)pbWacZgo9K|Be7#t%E3v=2yQ`S!A{}}|Fd4%%0ga{yC z#Frj8R>dZCEbid@j|$@@5opgt+~$`FHfZ(O=9uuPTkkydf|L zy{zj3c1^uV6)Y2qMcB#aeJ4W$IFEc|(Z>8~|6ky}!~Y9SO;;v8FLU}9i{=!&%C<4w zQa%#L)&{0B7`MO^%^OwdjHM#RkluB~p-8}Cz$Xb$qf>M!36@ZT^$VH2s`MMb#|

qQ(JoRTAhRSeDqoD_5$2Am$DJEfdwla^(+n1@RVhw`#7(dOlj#J+(Pfm_5GfE$j) zIcwtjbgGeRh=k!}!14VI9>BG*OM|LS979h58|x5sKPEvj0p#N?`id$jkh=G&USF*n z)m~&NkvK)f4Fhfo-JXMxK=-?z4c4Vs1&n zN5dmSBx0n#!6o7}zzuou372Yjr>ZWM?Ud+U>WzTo%0!0bTR1Zt4T9B|vn}YRyy}vq zMC1t4ny$^anq~KrnrmK2^D9Zj$?-TXi|PlQ^edfWtJK^~DTb7AXpV_G$Z)ERBP!y; zfK#J1gJC#EuYEV*R?>}K4eo^xBj`BH{fv@z0(1T*|9aL4GcB$N#wjDkb)(8=8%#k| zia39gVSsMUZ$IcvTUudDGzkP#Iaarbn+KdzQ}|r!To{h+OT7bdEXaKFEl9U9>Jfi@ zNzieaQ#)@F#XBfp_%vcL>jPv-pczMK(+7r8m9h;{5qCG>ZZf$dRkqLkg>WPs1{!i2 zAg0$ft2!X!o?XH)=2F*$;WS>Te*ZDR6&<5Yu`WNv9E<-AekT0G>n9>LO8CQl&bc8a z^4vl|hzU9p5(RZ}oF^jY3dSjFPR%w9i?{*6t*VUcpFQ2yi+l$4`e0lY??Ysb>L@Pn z6>+o0B9mdbl0Fsn)Jl;hfHO&_%Qt=s7jYi|oKxufUFvuk?h637h;E*0 z5+CpCKv@x?lQ8$x+!1ADk3+Rd)xRgjNbYtpjVhmQP(_^XQaylMZaJSRA}K;|TC9yR zs=w$qo3xIz}55GRKPLuvRHS$ zHC>Kt*RP)xQ32~ue>hPkvQ;CVBq;?6r?QOr$-ptAYGfNCB2LLO^cea!Fjde=2y)Q; zB2*WQyRNvG5q(>^$Npx(x!*gtOO2KLa(9!yMOb4S1>A^y3u-t0Tfwo+p`#O4;kj5m z?tG!l?!vkZxR1wq3QMr7F1~o^tqmq^X1ezwV(!A%*w8-ee;x$Wt0(d`gox9o+d|R)gZ5OMSY2AdpGpTy1Eg|9HK}Ij)I_ja8QSHFD%%$FJ zQSSm=!yp}S&N{eB`R1bA0<@7iUy5vrXOm9XafFDF0CEMgO|RjB$^9<%9i+=%cq^Da zRq_5VUWk#s?@EN!n+8bgC(#nS+7L|wT-2^Am2pxFm!45AMcgLwei7id z0q6CngXVOQZsI(pA>o25@`~g8A)p(ps#$7jzu8KhUJI1%wm6McnTNTsz%&9&qf?q#kgqv3Jdt$hVH^SiwfOAp zYV~#yL%^n_E+u%sCncMXF(s+$7J0wm)&p*)P|hjhI#5oCxPJg}%?_Fwz_Gv@2i&wd z9ip4awW0nPoqF4IS)_C7xUS*h!WLR-hfzJINi9$4P0PmG_s@<+2HmFA{8fN+ia1uF zu_s&89`6xxZ4Tn}r!@`#1;91!sDpsx3Vln!ss1z|-*j{8p_>P7F0Lsj9+|;7lXM2- z2oVWn2Z~omBM`?;F^lL zDB##;d%p2Y^`{Z}CYaMvx-BcJReSR5QE@HlVzZSHX*-iP-NJUPHfe{+VI(ST0--lw z;fR0Mh@W)M!jXZrg9Nzp=gvF>)`KeV6LDJ;zQjA2sP+A7z%@wn1CFf`nenGL(5*$= zB-eV)#fTCPc^&C`L8mgV@u8PqJCs_h;99gHrhRc~h~n8>>=sgIV~K+Nr~w_6Z`!Zi zF4YuqJ8L)*$BtAoMcjV{a05kfYzo!mPsijNXHF;Rwqi!6H37J;5h*QdXI*t2*V@$s zYS#b~X%NX{w`|OBl989c!H{lvGS!?PaI8yprW&f_6sZO?o!Y8GpNOl%xc@fb28!VL z3`4}9-b%L?AG2JumTyVO6G?~U_o`qOf-D?*K9)lu^45HWqH=G0%qDqFin#v-;37qEN`^uAr-SlsDcx|kC%?I~hP4yOL$V^0dVsNuT zxb)Me9h~Wfl5`=!NntpKXw6Qn13j#WtB7NXxH7cI&snNPTotSOZvoC%BccW1CejPt zNZ}yejLC5<*C-N^Z|*E152bo-h({SyJT7cO_N zI^Fu;DlxdmrAoO}vvX_jg(yT^IbsC4_)a6@n94X^!HrH6`F$x=UjS}h^vz8-AKIp% zCe_fDz4zUrHtJQdV@c=9Ih68k?_H($ak#^-iu2D;{GEB&LemxYwhZWD!d|su!Y^u{ z5}<~!ROT^QDF+iHl_c<)F*pTqdRIf!e!pKzuO!xs6e6x%+EC^e$|eVvP4+yvKlwIp zPPMq3Alw3IBUm?dgniUmBLdjN?jqd|J8v&n0?Ji(wUBact?w0~D~n6uKYQyj*PzPp zUQN2)1h}Rk&Ty%nUc&}z+|1b9{psCwgLsMWhHnt~!Sh!V2Ri<89EOI0(7z1m0 z)lSOcHEO9^)71y&#u(>oO5$AaFHA>5z3zv`%H*9!}$;a(b-IpIW_2Or9ttj@@4J zq)&GuVv=+ov`qmn9)@$eR4dg`zc9xfk-O3h4!Vu|(>vswk8anmCg>JG8$Um8R?e_T zxJubB`PRw0!oV6{wUTngX}2b|s%&2^pRXGs$9)fcZw1Opw9~GdQqIJdKji13+dSa# zg(cw5Fr3tT-{ex8R73AlcMZ78pZ?qOO*f|_bTe$JL^l&{!qD?+?^X`cXu7^CsC5R` z$~h#Wyi2)KYF(PENUeUM5L=Uf7;MYm(6Df*9NT=V%eNW8HF<{8u7Xh*^mNo1OD^_bc0Cq>9l;a zOWyXoEudQfZBj4SMD(3~Y?b6Z!|heoNGt0Muez(Va@+_BQ_kt^dAM(Q+!xUBfSo`z zu{?PHI3lp~pv~z|J%H=;3}?A-zcqDs5wyur zo3X3J9v|z9-MVVT?NzrV#gU6GLUR3(rEdn)^w{BnP@NxJ{v)$0auwk8JcHs=)5~a6 zZ9WSTaBKbPSIrJtCEuKBQwQB_Tk1xez%>`tq`XA8ADRu(75USm=r^#>fe=f`U(dcJ zSP^&nm1&2qe5Cn*V%oVb?p%}Tw31+d!Ekz^QX*G2>bCo16H+E=PV%MTIw{v=9!eNZ zwxF3@GwBAx*h}4_HVW>{GnmC9;3$Z@Y^>j(-V}YiX||2q zl@ctz<&Juw!OqG#>`DzW5fBnmN|<*>A|yZpMA}Z0Tz?)2m6Y4%8}4Zb+{;H0krKe0 zimWW|7Wp`hFhil>x?o>I%m7K4zf^Yqm)!^l(G)2sVQx#^oc}0Ery!}xIkWc(3A#R| zz@klerYzl>41*-!nue`@6U_m*)W_gJKr|PTFgSS zqH@QDa@f`V45Q;;2tY?0Y!!6FgehOKvvRzWR%r}?m|l?IP(Y1=T_feBTtSIBD9Ys; zL3h1rnU)w;w!VxX`lh07hHkw%g}O$h-v%87T$;_9K8d!7E#+sune}eHaBJwc7Hzy_ z8x(7_sDDca9Q998Qjh9c1L&CpOqjOb{iA~m0RagQDThU62`lBcsN6}qu1b#Uc4&fb z2#HeO_+q9Vj*uLgcD0vsB>|h{!-YA3UeGmY)A#tzV~*vTf;Q;OdS$w`p*MgMXN`PQ z{izDL*Q*^mVg}uQX*KV(E!92gB-$j&b{#eXHHBp~-J@O(6lpF9N#(IGp6}%XN{D%X z$0;M?A?1j+ST9i_<)Dyq#bjVQ+-S;>2@nv{b~GZ|yZ5D#oJ>GHB-bF1Q*exw(=Z26 z3c6h%XP!r9_r8aQ%m(LIBjpUVL0{7=(@ip`r+hQh3&#r(0QX7_0e@4_?IX5S6>igL z!;lX(p>|JkoEBfJ51dS6`y>2Gd2_7W-WCne001BWNkloHVzeMJRblD;xu z{6ZNg+Mus@lj)|IQ>Z6}DE<`cZmsQCDWxI+cL8p2lI1gKi`mk#I$pF=rHpV}g|;Zy z?CF@WzS{xzT3@u;*Q_VR7n(CNw>r4fv0e$tUdrvL+KNT!`ZhlQeKvJ%$djA88bf+j zt2Jo{5I3u0!m@(I3dSK3m5#N?T=_C8L3iFb-i_RZ_+2~MEeY9Ghj4S!tv9Eo#AEgc zHt;$YGTkM&;=ENwK)=oU_!-CfU^3SwV76-#j4&v8@iQtgF+hC`q{{aA#t!5N?@aQ~4s} zX@|s-QsP=vxQTSDS*FRis6XWZhp)tFx9_*5i^|;8e1a|9?m*iZ*T}bAbweepX;WD)&w42Ria4;STR>HHjtCn)byN`5c8MX?NJ)KJX>nZ1+E1@@psuGx&euyX79Nz0YB@ ztL1dLV#5OwPpC5OkhvL6$0*I#$+C%5*{(`buEE?UkDLC7s_2&+J)ef0s*<_#aIHsM zNVqkn8@f4V>}n(5P6%>PQ+Xw`U&qU)3=?>_}JkuAU7KP1>OZSy2gm{_scz@F6)&DB4v)%JrCI zGp#*$dfqLp5ow?%aji!i^j-83x<%w0F#Z&vFM_Kw3^m}sGUzrv*6kkQ=2vaH8zr-L z1iWpk{ASg4MI>^`$eA#p=J=>i-(Ynd0$@(Bl;4V!E4%cs)~(&H78NATdIuvUw|7Te zmB3fp0qgHJgbb7bn;h$23{02YI})sj2rnqyXoWcd=rgT-C}347oVf8ofEqMw95~h1 zpluo50L1iy8h|SShp(j3?JwdcdD50mgKm_s9RXLHYF8B;*d9F)k|T3_Y4DKj&KjgG z)`enJuCE8a?CPFfMQ&u~6_OL|YEkfUR08if+Y0a!5VTx7GY==!B+PX|w=Jl4jCL=4 zD{xI(h)bI7(GzQER|_ICwIgkzZwzgv&@<`AM&IC^;wz<4d-;ZsN8BboX|QMB^rU{Y zjRo4gToa(C6OA{jtIJ1p%=N$3r@gg|%a~xBI6`>Tk*HGcu$;S;+kZ<(b`>FQ?^7ZK zNZT|M(XQ@D;;Q7TO5mdac$0Eal!-K$YXWn7ws&hIeV@{77jyssCt2z;mrdNQkZ%(B zlIB)G>AW+c)Jw+_}b9dY*d6H#SI6K$DHX-uNd5GHpQTrT!m;|c#fBBBWME{qnl(-;b=j{1dv|5KgIo5gPwG5j@wtD z%@J#p(B_eBJ=Z*9jlljlL^>Q`MYC%86fC@Y0!o+@87E3P%DW-mQE^B)qN|7mOS{TV z%{(&?MBK>C?S)jChf>O2=ge6Hd^pxK$0}fMkIiivQkHb(q=U=DPNh)2>}yxjATuov zI?`2!)QvVq%_(&2GN(}204w=ZOoe>g#_)bBNg8WQ@QgWHrrMrEtvv^u4dEnhacx%sdL^PAON;Ah|@`=2)TT z+L2?GFt?rR`?88+!S-ogk4RjZMZGJ3>Lr@WwE)`ch+B(1C*1^d3g_IntKTaz{xoeU z{JQ#&zZuzmgKdO}>2W8I9ch`D3RCY)<|j#`jhI@1<*%TuU5_$sKi;W>kA zEZW?{4NSUS2LEMoY(9=2(eD%mJhq0Kn;6;cQnZtKL}H zEmd%mMq__#axI275I0~|9Sa=YL~{zj18zfG;ZJwKl|A|2jJ6xn2=}4Q$u*}^Gl|Bt zj4+6+`c!vFKmu&0{ZjWpa`bZXQVzgZMAF>O@{zXZ;j*i3=Ar3%$ouq7Y|2DV$C{a$ zOIOcnt|D=zZ*55@z-}}Kho0gQrE@F2(BoBmL+U|W`P(qv0KD?81YCOU(*ESpV^zDK zMw^~SFwk}z*JN>L;fw+$a(*+*x>1&-4Ad~#=Rg3Rl!Gesa6A_?RIahBR7knAMSsgY zK-JfEPCws@NXjuJ)ufOk`EdACECX{;^y-(GJFkM)x7Mu-JnV)foQHai&|{&-a?OJ_ zCF0iJ(?K`MoF3^^0dTKZ#bj!idB5y_T$I74rRW2L^7 zWXIh06DWSz`PR+478!;Pc8la&?_1MLbgcY%xc-NnS8b`M3b*o5Ub=xO-}Wz@7w1oR zi<;UJkhS$6{~X%xrf zWgY-Ty9%)PJb>K3T07QjxoQ)bgDW-%`N}i2z7=*&r@(GW04}}I?^w}NaKmUTgSDf*_-vtAn7smx}<265;!layg|} z-f|2_f^~1C7iOTBcC}ja&=ty^eO>+D)qJe{+N|g%<%+-@LfO@D1m^050|1DKtu4R&z+=NQ+# zu{K|O{OIJ5`Jy)6jW#rl>NQrZNlGo?Sbe6okV;E#MR?kRprZZ+uX^m<+H_pTl!(qE z657gubt!jdSLw}U9=bxgM#{BOIZT)`WCL@%jx9lUkgoq>_vW#@bHPjZs_4!$fYYTq zxz-ud3T@yMZgtvJq}v`-+bVvMmv38({=<^(hTeEqvbng{9Z;8S9Yt zK><}TZg;6ggbrWVY&eN`lNQy?Fn~0>Tvc-}t_f(X_Fzey>bZq7=h$AIZk>F4Rpu0r za|TgB+Xp1uS4G(*z=0NnTLCpfKh9)|`LhbDN>CFH$Fq~lKC5;l1H3kB5oW|!kzn2G zNC*>QHi}42n4&7>&i$>xte3`=DW?sFNAh&X902Soh+~XJk%GWe&@nFcpc|fRswLhH zWf&OV4fxa2qNeslm&)R7$)=!<6K>G#TL{ptkZ-4~yXc04lFeK85u4Pu&WiHoY@=4x z<{HbXkXLPf<-l3EY9!kBa2^ef!UgJHsy`k_i7+3&7W3XyinNn!|HzG@avO5CeO;`S z+lolaX-l+AxvN&LYI~?KN4dzf?vw8JOzW8y<*n=aNfu};dAH>pV#3XlOVuNB9MMFm z>0Fc0CS@U9!VO@AZl!!H;&G#B3(VJzcEz(+bR@*)=8n7xUyumjxYv;++Az0gSJ7JaP!9)2sH7Ye)N7_Kl~S%vK3p{-Jmz5U^g1t} zEs(cAWRr_qxGi-baw+Aud+v82P!p_1OP1+RRiY_Sb8<~duvSf`s&E7F(G8#gT;j~> z0XTfEi?Yd*P01d7G%G5F*v7av47HkQ;d&7fVfV!2imwo@$w8_j{^4p*Fau#B=*oIv zB59~xQ5#%k9;m7ZJ|~f!Qf`~Jlruu21ekN1d^j8{)Umk|b9-W;{4i4`=t}p>q-t1~ zstoLvJSLkxNa8o7SvA7sRV}WGXk)k2>%t8nNVh`1ZQ8ApZ>MjKNj7ixU};wL$q?HZ z*G#YKfm*m&M4;F}BDMFzh?cb_Y%W~3yUS5KvS(Kz18%Xovb75(DYw0VdT{36CUNe0 zsN!*D`|8!DU)7OxjRI7GMT?(f0L=eFv*JRPM|?oCbNfs8CBefP~cy^sTsfD>=cVo4~$Q z$wF*`Zi~?2Sq@Uf>3U^XGpR7DL7VAOdDR2WghK;02iFX=se9>}4KnF0^%$XBE8kKA zOt(umw-x<{5Sz(0-GK^FixaKpRQ7?49GYuB7{`{p^_#oivJC8V4oI+PPHC&*rNgBb znTM*i3+kbA=b3>BN1a`GOz}(J2Pj&L4cdQrYt{M?V&{dkQNTC+6lhZMLBo6&~ zRS9QeQAk<3B0x>A9Z@73gSJhkS;Fo3twpyphmLzZt~IBKJKzwblFiPFE*kQ&wu)^b z*I2Q}Dz$-zEi;`;5q{<61LTV5Dla4I2q>*AWa4o<5MZwmp<@_{tF?>jJrAXnyH<8V zDdqOou9T0I!&U*;I@S_%8>g-sk(_|?Hic;ybcYYE+-sebXipCH!|}H$fV0F`Cz@P7 zqUkjiZQ4e9Nw}TgRbIN)@~tQwIq621R}Ewj99DEqi0xBcyMEH$A-nCT2WmRe$YPmv zmCPhkNgin4x>V}!RySVy<^b;-X^wCk>@E7#u6lsBh~(5)tzF&PSId=MXJ6`j>k+3t zcvd7II>!pfbh*9MP4KN-B4Q)mm0>u`;toUF6HN-q&EqUYG)bvd$44$OXG%65Z2+xs zE9`BUZj5|8|KB0mtODI-A+}(n=wRw#4c8VxjR>7MV_8P=`nZwgX%L7E_;kj-313Kg z9~yV4HY@Yc%s_{8z4W=uX3I8QQjXfqg!`%oP#o()v8`uT@U5@Ii{`YW5IETFJuVdf*Ic8&f3CKrcHS&b5T6lxw>5%lW`yI@Z&) z9`luY-+F2~9d=!Hh?G(Fs56sg5UGX}iB@JIb-h+fHf5h_kG2bE-8QULWkthu!@VAl zQu1@nuh$0RYd7UXL&?yNYvBqw7u3r8lMRbuc-_&jTF%GlvkJ9QMHk~PDR(4TODVV2 z10PA!;`}>L{f#&LhG|k-yN-6Q+ObmYSmBDz0rbALtq##)=aO%xK&uPWhH>U0TKz!P zhgOo!T4vgBGHuXy-mU32ZRoap+m}Q)KqN#0^l?p$uMJG~xT`MrX0g^oZH8zeob@1^ z;aBLfN|juHoPFSAXb$i@<10bV+fbuKm`f?Qy_ttyg{DyM$h4MrwOLC@{CwcNUo@3t zEqe9$dsl5s2mSVs1CSg>T&i-5#L>qwW{r)qJX0 zj!|kFBLq#%R;9JBR`Huq4`?!pCP8h6YfY>Tpq*Y>qOJICH6Wy~={@8dZHu@#x*d1( z2G{(>YilEGGZh|^STh>t32msYBbv!F_6}LQKiqM9w_Iy!XzmC7VgE2f$h0qU$a~Ax zcCrFkE4x~|N@uMcR}Z{NIn0$~Ez0E@%t7Z{H$iui4i8kdu)}j)uv{tf(o05^sx*d; zm{@hz(a<2847D+?0l#^^nO<3zlK)s|dwx zrLu^Gk#cYiLodtaI?V0AZv-8_fatWM0_>oif?2~$8$QL0wZ7U`5NN$R>+%>D(R8Rq zxdy-((lSrVpbcbOdc@uKpIqiNV`AxcPLnn>p#x&g8Cm-t2O5Q%hiHDB4YN!q*Sw&6 zGp(!!3Y8}nat^hWBg9-p<(+o**b`G{9?Bw;M#^oSTW7f(hdBT!1>H3=BwyxTgq=7& zN$$Px?59+c#s%PPf#xEb4z(!PfLns4qpfnLtu4I@dV6y!(9PyrT&#V{ ze}8$Bl1K3kLM@>MQ`GRM#^Yd@J$v2A@wX5clG}cWgA;Vh!Be=D9@*I@(&{ z*4ffMb|Hk`a(bO~yMl2Y*H{x85^G@xdRt+K;y_nHZGdPE&Ri^`Ft%0+I(BVTMKNkG zpyQ-6y(;Ze4vt0uET)f8akiBo5>pHeN3-jmpwOip?Q`We$0}oP`#B)c9&6j(>m}z0OsdbZ;~@S%4{g^7l`X~iC_^(z#_8hAFAatmS%wAI25DqFgrTWaVnr`I{U zrso=KLPND3k{0`hT!AxDphGQ6G#zK%3?&!E)FrA@ILo97H|8~9i^86Ldd}gDSz+Gu zqWj{(b0Axp8r$sZb$(RW=Vn*mRS}6IkHUbyn|-guj{Ku zFSxC$(t<=|a3;#MPN^wk4b*j|$`%D3ZIy5XzM5AKy&d&NY5T2dLT4PPnJf4$N{#sX zIGbo8oMopijjOo<^!36G)#LIV4(Q1_5C9!hY&>x-}3bn$%Hx_=xA8&*Fh~xv=Gibwp37a-sHf7JcD|fzKO#cd2$ZN z1SO)&TzUcJ{V6lS zHuJK_+}V1$O$^9%!S3urrAt*&ky5CtZF1(~{Vd^BsC89oVp@O)mLkvuoK^qOwv;k@ z&E5-hps>o5f)j1~9eTQ9I(AohC@&-PHw7s(Jx+TiT2I@11fX$!0H*A$eBbjT` z+v@4})L52Oy(SZ|?*QW_gERIM{L8Ki8tPoz^GHz2y;HTpRk+&y(<< z6_(9_OntbxyANPu7tDUQGKWxnU0ego`#Rpb9HgaNhE$?Wf|}EKKBmzg2&J(CEsV1W z%i2(#pKBtFaGqIVDRlI1&0&zW z(yUWwBG=R{FwS}6MFOD=&z-)qM_W}3q6)VwDNkNo`n*OPwwz`h%`eb$Si#MLjNI_D zO6`_K;uJxcka>}pl=JC0M4(N|0YM__qjCh0kC8*N)YYzw&Z>e8fZ%r z#IckN*V?YqvzfHAoOQOeq#Nxbu0k6L#&aRj*s&Z9(PjkN0MSfZMRBIHj3`VoiY8BK z5>}R!^NPs4asL`+0yf>n}oTd;0(g9MaQN_o$x?&WeoBw~(^D&f2b+9c9#_t|G{ zmUKHjDsg)K)!PM^O>Hq3XcaI)rq+uHEF+nf!|SZr!EMbHkfx=Sb)&~UuX;&IxUnLx zPIr0-v8@WCEzb$e3p5YTLM$Vz=$A(rjnFFJGjzlVCXvlV(W{x>3ul=$A+8yrErmCRDA&~Xz z3>n=UE0_pu&QZXcXUDPulXJ#2{zBus^4Z81lX5vdV$urDBd z?ZBd=RG5fhi^OFtjcOHjNv>F{lXTdPDVN%w${g#YB9lP|+`3{7XFj6oVhy@pD}sh@ zm3*UuR;P&D?vT$;=Q0vUy;uFDKpRkJA)I;4ri8N)%gCY_t>D@-B2JMa{gHJ+$JED_ zQVs!LxVJsLGj38Y)3z}%5oQ3~q#SmRHE(Ix5_4yKEu~%1!F#I#-US_`y*#fB)omBg zHpS%CM%j=(3!;6&MZHLCJXxrihTw}R9BXY6KRAjDPYJ_VK)Wk(`iIG?M)4yMsx2iVS5)bcA z%Hcj?X`3{c2;WGU5gW&Pj+Wancl1x}voH@qw};?>&SeX9I9aS<2W2hX^=DC26xfdP z&8`fCB<198Cn(A?#uZtK!u}_nB%ItUDotnWM}3qbHb%0MNIFKi5uTE5+XgZz(4siw z1={1WGZ)TeU6vz>)Jd*s*;JoAqkIeO&NH1v#fKgpM|+<`Hl!RNykKLL=vXrXBBs3c zyk4NTO(nCMZT3TNHtCVLqCTTdWj3u&PB}IN zvYuryvBO0gGud~H>Rf|;%TV+1%9GL=dZSATP71V-tdME% zJ{7QzWe?j-C(GhBy*9$cf|wyx%_z9+Xk&;1G_bRVjvlH_%0Y&Nyp&4^_3&6uxJCQ( z-les-D|VtH~Mf{qR$cI+J>n`I)-fa{%-dPqwj%EsWTG{O;UlhBrZ(oCn7pJiX-F=2$1t;h{B?NSZ_G9n_KTf6olOMBvAPd*URHj?B}Ka_F?Ovi4A zebnX|p@VaLI<@oNHrN#^wezQ##@@#ON-U|F7+O? znyyJ8)LgsBc%Y3g4pu}{q}eK#Nyu{5Qo`Trh^%LZr z5)vV7jS6WSvbhVlO%m;}JB7kgp?Bp`TD^wtI&Jq`FTa0+lf4al*jB zZu=>g+AWQLp~wQaPVi}#v+3D{0mKbE@~&Jqs=`dKs%^9b?>EXwRB%tyZ|@O4FP23rZ<>Ky4%4ER!F%_xVzM?h7JxXckP)b z1bh>em^+pTq-}~wg8WXeEtzoasv$ULaS{?>+J4IdXkSx30=VRh$BF4nwyUYhIp`{^ z=a&;<{XN4A3iST=f@n(pnk3emDBFp)PQF!j-2NV>khOG{k3H`z&6yu2#TbUV;U$)d zsKD~}Y^~OAA&7bc@<076KBe3?eQh7JRLrbAH`?r)ir@_g+5j@pi!dMafAd~WMgaK@ zF9gDKH^2eA;|(V+mQr5k0~rQn59EpcD=kg+Px=74gLDGIG9OW5G0Fom+S%R&)VN@l zmNzK9XQo!NokXLHbm&xMdquXXJVN+d_ zE;>SvAL%VrVVsCc=jKr5w{gbja+Ns?}2yzY~;N0W}VD1=NlU?(9O& z6tE|2Uwgf98 zy|zyp6HrdLx33nM+w%}LscX(i{gI)a+G5$vCj3W>e4*ITQauofW2piy63!oa$ z5S%q6^%Z1ZgLD-hGSn_nd7Ld0*GRUO(SO$gJFlb6>g5`Wiwl?+*b6BaKEv*iIHz>Qqibu#M5@ zo_ke@Weiw1zt}hc)drReQm1Mgdj2?+0}vFc-$1YpuO@_v zufS7lzhchMl_lDXbFIoX?l(q>V*co46T({EVIE(YM}=uGzp|_$+n80;ZuRZf5jN%Z z?mFaSiS+6?fQ(du20}#x5UD`K+#~X#@)Mv+NsW;Lu{T;%r%O$B_(&BFsSq_e_5@(a z+oT+hw)}n8l6c|F^7mMdC1%;5z5ZF8cB9&3!flpiFF5lM?fCD;dZlc706JCJ8Gg0u zuM6QN+0pduCYQmit5e>F=ZIuqOYqHTXu{~h9~T+IPRd11<@a-|r8J!IArb(JOpvHD z1t#BIq(}%2{n-QyaZGZO^0)vGx#V2& zlaZMAnv!-n(jeMF>v8Q`zb12+5axH}4~Li2`Q=ewUfJa~IXfD`t!2t50@)kX+rs zs$jy|$hQkM`oL(0t!m*U({E;a5I0Fdx(`i;45W(uJ~p1mbGCyN9W zLfu)$(VAXf^TENq8){_Sw*xBWzDlja+-6o^znywAZCQ)iD$l4AT@BBt0sFEq*DN{i z1f3L8J#2V&$gcO3tFbb);b85Z9OEIBmq(+i<#7VSVX>;a4O2$P38~5{$cvPd2jvGM zBFXm#P#7Q*WRi`0PLhtJq&n24(rbD0oF!iS1ub`jf&2P;{r7UHPud!->r*v6Torhk zFTu_=waHbrQpN0WYs3V$CB#`(NN&7u2PSTJ_;Gk-8@YVGDkruI$U7^zGFq~oS;d>P zN+=E^PT#>-bMS;eP8AUXhM_o=3-0)W`HQ>Twf&jTk1)I z)NRh-H*FsEDQeB+TU2CMDKD&~blTdJ6>dRiR}Ash02H2H-UfUXWjQ+U{c(Adb0_pR zkG3ABC2~3F{#6b~~NHmXldWY#B#t%T&G|X7LhKxd_!lqT!OQ zo1v9RYwXf-qOz58sqA5yzl7zDnp==hb0RN*#Q`S*q={TJ(*hRJH|swduSXp|dmCi* z-Hy)QWZTqMk$g2$L$Acf#{7Aojssg#4l-o(FiH`C1K=!U^9eQXZfjBu%o&@?X&4p@ zGualn&u1Mr(yC(W7JgnwfDQvwGj8;A<*6wnZLFL|BF_Fa0bb zkT46n9q(Elihe7sh78zewmV+DWVDK_L6A#TS$D7F#nvp@l2bo_Sjq8!1euyS@`eNvi3)&NGJz$E zOcZ?v8P=~C-W}~z)QTBs8mj7ZB(rt{IJ;Sp8N+ZJSd`Np5T&aLVpG=&yT&N%EX#r{ zBh505E5sSJSze72R;{Y4n^UJw>mU@b1kfcRh2iG}&N5C&DjBIzi6MF00OV+_M zJp>%slBh8pQdi3;525OKx^jMkjzAarJFL`#KQ5ABMCzj?hyYY1MBT+7^ubxXh1tJ0 zBTq$^vaMGWHMG5F{M8d?Z*im6Y; z`HZ_t;hyza#gw+ZiI;P$H|Zni({9So%Kql0pKJSHHN(Kii=c;e6*{Oigh!u&kzT@} zQ@9yQDhmG~<5V?a?k7y#q{{&c4XSag5|d;LfGnk{>-QCfxfgOJ;wagB&o?vztc?!Y z8LE2Bh--OneI+>E38I+`=raMQv1-Oe95!Q=(q{s z#BK25@y}}pNXi11fJ8dGGS<`E(-U+Xcty0*+2a*M2{fZT*-b8_k!GEq;>ZqIKB0Za zV){tQ(Wq@UcDr6Tu9W->9=ArusX`$UKQg7@!x3;Z=EBk(4#@Yb>$t72NPw)7sMc{% zy^r`VNiy6n^E{vJHanSR*6|Jw+=^p)@k`qcjk}n?AR3$V#{mEo5S92(v93Ra07134 ziMu-WLSJ&mu1ALy-jF-}O0amHL!!~P7Z z4@tn#Tjx#hq`EuL4Hgk!X~K-6+ZFr&>-P4y!sZVrxD94pj5=W=`S|y7UZ6P%I2>)? z>Q(v8AO7T@-)w(CEAb{>m2uRlQ8dt6TTz>2ps^IrG5*WX>?z6AHW*S&7W-RoZWy4TIPd)@0^_xiH+A0Z~I$2~5d QdjJ3c07*qoM6N<$f-nnA{r~^~ literal 0 HcmV?d00001 diff --git a/deal.II/examples/step-45/step-45.cc b/deal.II/examples/step-45/step-45.cc new file mode 100644 index 0000000000..962fb7d18c --- /dev/null +++ b/deal.II/examples/step-45/step-45.cc @@ -0,0 +1,275 @@ +/* $Id$ */ +/* Author: Markus Buerg, University of Karlsruhe, 2010 */ + +/* $Id$ */ +/* */ +/* Copyright (C) 2010 by the deal.II authors */ +/* */ +/* This file is subject to QPL and may not be distributed */ +/* without copyright and license information. Please refer */ +/* to the file deal.II/doc/license.html for the text and */ +/* further information on this license. */ + + // The include files are already known. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + // The RightHandSide class is a function + // object representing the right-hand side of + // the problem. +class RightHandSide: public Function<2> { + public: + RightHandSide (); + virtual double value (const Point<2>& p, const unsigned int component = 0) const; +}; + + // This function returns the value of the + // right-hand side at a given point + // p. +double RightHandSide::value (const Point<2>&p, const unsigned int) const { + return numbers::PI * numbers::PI * std::sin (numbers::PI * p (0)) * std::sin (numbers::PI * p (1)); +} + + // Here comes the constructor of the class + // RightHandSide. +RightHandSide::RightHandSide (): Function<2> () { +} + + // The class LaplaceProblem is + // the main class, which solves the problem. +class LaplaceProblem { + private: + ConstraintMatrix constraints; + const RightHandSide right_hand_side; + Triangulation<2> triangulation; + DoFHandler<2> dof_handler; + FE_Q<2> fe; + SparseMatrix A; + SparsityPattern sparsity_pattern; + Vector b; + Vector u; + void assemble_system (); + void output_results (); + void setup_system (); + void solve (); + + public: + LaplaceProblem (); + void run (); +}; + + // The constructor of the class, where the + // DoFHandler and the finite + // element object are initialized. +LaplaceProblem::LaplaceProblem (): dof_handler (triangulation), fe (1) { +} + + //Assembling the system matrix and the + //right-hand side vector is done as in other + //tutorials before. +void LaplaceProblem::assemble_system () { + const unsigned int dofs_per_cell = fe.dofs_per_cell; + QGauss<2> quadrature (2); + const unsigned int n_quadrature_points = quadrature.size (); + double JxW; + FEValues<2> fe_values (fe, quadrature, update_gradients | update_JxW_values | update_quadrature_points | update_values); + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + std::vector > quadrature_points; + std::vector cell_dof_indices (dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active (); cell != dof_handler.end (); ++cell) { + cell_rhs = 0; + fe_values.reinit (cell); + quadrature_points = fe_values.get_quadrature_points (); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { + JxW = fe_values.JxW (q_point); + cell_rhs (i) += JxW * fe_values.shape_value (i, q_point) * right_hand_side.value (quadrature_points[q_point]); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix (i, j) += JxW * fe_values.shape_grad (i, q_point) * fe_values.shape_grad (j, q_point); + } + + cell->get_dof_indices (cell_dof_indices); + constraints.distribute_local_to_global (cell_matrix, cell_rhs, cell_dof_indices, A, b); + } +} + +void LaplaceProblem::setup_system () { + GridGenerator::hyper_cube (triangulation); + // We change the boundary indicator on the + // parts of the boundary, where we have + // Dirichlet boundary conditions, to one + // such that we can distinguish between the + // parts of the boundary, where periodic + // and where Dirichlet boundary conditions + // hold. + Triangulation<2>::active_cell_iterator cell = triangulation.begin_active (); + + cell->face (2)->set_boundary_indicator (1); + cell->face (3)->set_boundary_indicator (1); + triangulation.refine_global (5); + // Here the degrees of freedom are + // distributed. + dof_handler.distribute_dofs (fe); + std::cout << "Number of active cells: " << triangulation.n_active_cells () << std::endl << "Degrees of freedom: " << dof_handler.n_dofs () << std::endl; + // Now it is the time for the constraint + // matrix. The first constraints we put in + // are the periodic boundary + // conditions. For this let us consider the + // constraints we have to take care of in + // more detail first: We want to identify + // all degrees of freedom located on the + // right part of the boundary with the ones + // located on the left part. Thus, first we + // select a degree of freedom on the right + // part of the boundary. Then we look for + // the corresponding one the left-hand side + // and identify these two. Since we are + // using finite elements of order 1 here, + // finding the corresponding degree of + // freedom on the other side is quite easy: + // All degrees of freedom are located on + // vertices and thus we simply can take the + // degree of freedom on the other side, + // which is located on the vertex with the + // same y-component. Here starts the + // implementation: First we declare a + // vector, which stores the global index of + // the boundary degree of freedom together + // with the y-component of the vertex on + // which it is located. + constraints.clear (); + + std::vector > dof_locations; + + dof_locations.reserve (dof_handler.n_boundary_dofs ()); + + unsigned int dofs_per_face = fe.dofs_per_face; + // Then we loop over all active cells and + // check, whether the cell is located at + // the boundary. If this is the case, we + // check, if it is located on the right + // part of the boundary, hence, if face 1 + // is at the boundary. + for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active (); cell != dof_handler.end (); ++cell) + if (cell->at_boundary () && cell->face (1)->at_boundary ()) { + // Since each degree of freedom of the + // face is located on one vertex, we + // can access them directly and store + // the global index of the degree of + // freedom togehter with the + // y-component of the vertex on which + // it is located. + dof_locations.push_back (std::pair (cell->vertex_dof_index (1, 0), cell->vertex (1) (1))); + dof_locations.push_back (std::pair (cell->vertex_dof_index (3, 0), cell->vertex (3) (1))); + } + // Now we have to find the corresponding + // degrees of freedom on the left part of + // the boundary. Therefore we loop over all + // cells again and choose the ones, where + // face 0 is at the boundary. + for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active (); cell != dof_handler.end (); ++cell) + if (cell->at_boundary () && cell->face (0)->at_boundary ()) { + // For every degree of freedom on this + // face we add a new line to the + // constraint matrix. Then we identify + // it with the corresponding degree of + // freedom on the right part of the + // boundary. + constraints.add_line (cell->vertex_dof_index (0, 0)); + + for (unsigned int i = 0; i < dof_locations.size (); ++i) + if (dof_locations[i].second == cell->vertex (0) (1)) { + constraints.add_entry (cell->vertex_dof_index (0, 0), dof_locations[i].first, 1.0); + break; + } + + constraints.add_line (cell->vertex_dof_index (2, 0)); + + for (unsigned int i = 0; i < dof_locations.size (); ++i) + if (dof_locations[i].second == cell->vertex (2) (1)) { + constraints.add_entry (cell->vertex_dof_index (2, 0), dof_locations[i].first, 1.0); + break; + } + } + // Finally we have to set the homogeneous + // Dirichlet boundary conditions on the + // upper and lower parts of the boundary + // and close the + // ConstraintMatrix object. + VectorTools::interpolate_boundary_values (dof_handler, 1, ZeroFunction<2> (), constraints); + constraints.close (); + // Then we create the sparsity pattern and + // the system matrix and initialize the + // solution and right-hand side vectors. + const unsigned int n_dofs = dof_handler.n_dofs (); + + sparsity_pattern.reinit (n_dofs, n_dofs, dof_handler.max_couplings_between_dofs ()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern, constraints, false); + sparsity_pattern.compress (); + A.reinit (sparsity_pattern); + b.reinit (n_dofs); + u.reinit (n_dofs); +} + // To solve the linear system of equations + // $Au=b$ we use the CG solver with an + // SSOR-preconditioner. +void LaplaceProblem::solve () { + SolverControl solver_control (dof_handler.n_dofs (), 1e-15); + PreconditionSSOR > precondition; + + precondition.initialize (A); + + SolverCG<> cg (solver_control); + + cg.solve (A, u, b, precondition); + constraints.distribute (u); +} + +void LaplaceProblem::output_results () { + // As graphical output we create vtk-file + // of the computed solution. + DataOut<2> data_out; + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (u, "u"); + data_out.build_patches (); + + std::ofstream output ("solution.vtk"); + + data_out.write_vtk (output); +} + // This function manages the solving process + // of the problem. +void LaplaceProblem::run () { + setup_system (); + assemble_system (); + solve (); + output_results (); +} + // And at the end we have the main function + // as usual. +int main () { + LaplaceProblem laplace_problem; + + laplace_problem.run (); +} -- 2.39.5