From 8d9bc9036e14a473ef8ec1869604669cafc62e04 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Thu, 14 Jul 2022 21:16:08 -0700 Subject: [PATCH 1/4] Implement missing functions needed to use multi-material simulations with internal vars When materials have different numbers of internal variables, this pads the array so that each quadrature point has enough space for the largest number of internal variables. This is not memory efficient, but it works for now. The while multi-material paradigm really needs a review and probably a refactoring (if not new implementation). Maybe "blocks" should be another array dimension. --- optimism/Mechanics.py | 71 ++++++++++++++++++++++++++++++++++--------- optimism/Mesh.py | 6 ++++ 2 files changed, 63 insertions(+), 14 deletions(-) diff --git a/optimism/Mechanics.py b/optimism/Mechanics.py index c6fda849..9dc98fd7 100644 --- a/optimism/Mechanics.py +++ b/optimism/Mechanics.py @@ -106,7 +106,7 @@ def _compute_strain_energy_multi_block(functionSpace, UField, stateField, blockM L = strain_energy_density_to_lagrangian_density(materialModel.compute_energy_density) blockEnergy = FunctionSpace.integrate_over_block(functionSpace, UField, stateField, L, - elemIds, modify_element_gradient) + elemIds, modify_element_gradient) energy += blockEnergy return energy @@ -124,6 +124,13 @@ def _compute_updated_internal_variables_multi_block(functionSpace, U, states, bl dispGrads = FunctionSpace.compute_field_gradient(functionSpace, U, modify_element_gradient) statesNew = np.array(states) + + # # Store the same number of state variables for every material to make + # vmapping easy. See comment below in _compute_initial_state_multi_block + numStateVariables = 1 + for blockKey in blockModels: + numStateVariablesForBlock = blockModels[blockKey].compute_initial_state().shape[0] + numStateVariables = max(numStateVariables, numStateVariablesForBlock) for blockKey in blockModels: elemIds = functionSpace.mesh.blocks[blockKey] @@ -133,9 +140,9 @@ def _compute_updated_internal_variables_multi_block(functionSpace, U, states, bl compute_state_new = blockModels[blockKey].compute_state_new dgQuadPointRavel = blockDispGrads.reshape(blockDispGrads.shape[0]*blockDispGrads.shape[1],*blockDispGrads.shape[2:]) - stQuadPointRavel = blockStates.reshape(blockStates.shape[0]*blockStates.shape[1],*blockStates.shape[2:]) + stQuadPointRavel = blockStates.reshape(blockStates.shape[0]*blockStates.shape[1],-1) blockStatesNew = vmap(compute_state_new)(dgQuadPointRavel, stQuadPointRavel).reshape(blockStates.shape) - statesNew = stateNew.at[elemIds].set(blockStatesNew) + statesNew = blockStatesNew.at[elemIds].set(blockStatesNew) return statesNew @@ -144,16 +151,46 @@ def _compute_updated_internal_variables_multi_block(functionSpace, U, states, bl def _compute_initial_state_multi_block(fs, blockModels): numQuadPoints = QuadratureRule.len(fs.quadratureRule) - initialState = np.zeros( (Mesh.num_elements(fs.mesh), numQuadPoints, 1) ) + # Store the same number of state variables for every material to make + # vmapping easy. + # + # TODO(talamini1): Fix this so that every material only stores what it + # needs and doesn't waste memory. + # + # To do this, walk through each material and query number of state + # variables. Use max to allocate the global state variable array. + numStateVariables = 1 + for blockKey in blockModels: + numStateVariablesForBlock = blockModels[blockKey].compute_initial_state().shape[0] + numStateVariables = max(numStateVariables, numStateVariablesForBlock) + + initialState = np.zeros((Mesh.num_elements(fs.mesh), numQuadPoints, numStateVariables)) for blockKey in blockModels: elemIds = fs.mesh.blocks[blockKey] blockInitialState = blockModels[blockKey].compute_initial_state( (elemIds.size, numQuadPoints, 1) ) initialState = initialState.at[elemIds].set(blockInitialState) - + return initialState +def _compute_element_stiffnesses_multi_block(U, stateVariables, functionSpace, blockModels, modify_element_gradient): + fs = functionSpace + nen = Interpolants.num_nodes(functionSpace.mesh.masterElement) + elementHessians = np.zeros((Mesh.num_elements(functionSpace.mesh), nen, 2, nen, 2)) + for blockKey in blockModels: + materialModel = blockModels[blockKey] + L = strain_energy_density_to_lagrangian_density(materialModel.compute_energy_density) + elemIds = functionSpace.mesh.blocks[blockKey] + f = vmap(compute_element_stiffness_from_global_fields, + (None, None, 0, 0, 0, 0, 0, None, None)) + blockHessians = f(U, fs.mesh.coords, stateVariables[elemIds], fs.mesh.conns[elemIds], + fs.shapes[elemIds], fs.shapeGrads[elemIds], fs.vols[elemIds], + L, modify_element_gradient) + elementHessians = elementHessians.at[elemIds].set(blockHessians) + return elementHessians + + def strain_energy_density_to_lagrangian_density(strain_energy_density): def L(U, gradU, Q, X): return strain_energy_density(gradU, Q) @@ -184,21 +221,27 @@ def compute_strain_energy(U, stateVariables): def compute_updated_internal_variables(U, stateVariables): - return _compute_updated_internal_variables_multi_block(fs, U, stateVariables, materialModel.compute_state_new, modify_element_gradient) + return _compute_updated_internal_variables_multi_block(fs, U, stateVariables, materialModels, modify_element_gradient) def compute_element_stiffnesses(U, stateVariables): - return None - #return _compute_element_stiffnesses(U, stateVariables, fs, materialModel.compute_energy_density, modify_element_gradient) + return _compute_element_stiffnesses_multi_block(U, stateVariables, functionSpace, materialModels, modify_element_gradient) - - compute_output_energy_density = lambda U, stateVariables : None #materialModel.compute_output_energy_density - def compute_output_energy_densities_and_stresses(U, stateVariables): - return None - #return FunctionSpace.evaluate_on_block(fs, U, stateVariables, output_constitutive, modify_element_gradient=modify_element_gradient) + def compute_output_energy_densities_and_stresses(U, stateVariables): + energy_densities = np.zeros((Mesh.num_elements(fs.mesh), QuadratureRule.len(fs.quadratureRule))) + stresses = np.zeros((Mesh.num_elements(fs.mesh), QuadratureRule.len(fs.quadratureRule), 3, 3)) + for blockKey in materialModels: + compute_output_energy_density = materialModels[blockKey].compute_output_energy_density + output_lagrangian = strain_energy_density_to_lagrangian_density(compute_output_energy_density) + output_constitutive = value_and_grad(output_lagrangian, 1) + elemIds = fs.mesh.blocks[blockKey] + blockEnergyDensities, blockStresses = FunctionSpace.evaluate_on_block(fs, U, stateVariables, output_constitutive, elemIds, modify_element_gradient) + energy_densities = energy_densities.at[elemIds].set(blockEnergyDensities) + stresses = stresses.at[elemIds].set(blockStresses) + return energy_densities, stresses + - def compute_initial_state(): return _compute_initial_state_multi_block(fs, materialModels) diff --git a/optimism/Mesh.py b/optimism/Mesh.py index 7aec1b2a..4720060c 100644 --- a/optimism/Mesh.py +++ b/optimism/Mesh.py @@ -136,6 +136,12 @@ def mesh_with_nodesets(mesh, nodeSets): mesh.blocks, nodeSets, mesh.sideSets) +def mesh_with_blocks(mesh, blocks): + return Mesh(mesh.coords, mesh.conns, + mesh.simplexNodesOrdinals, + mesh.masterElement, mesh.masterLineElement, + blocks, mesh.nodeSets, mesh.sideSets) + def create_edges(conns): """Generate topological information about edges in a triangulation. From 623430954f2317c07f7716fd49b2f4c54d4e2d79 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Thu, 14 Jul 2022 21:38:28 -0700 Subject: [PATCH 2/4] Add example simulation using multiple materials and internal variables --- .../uniaxial_two_material.py | 167 ++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100644 examples/multi_material_tutorial/uniaxial_two_material.py diff --git a/examples/multi_material_tutorial/uniaxial_two_material.py b/examples/multi_material_tutorial/uniaxial_two_material.py new file mode 100644 index 00000000..0304be64 --- /dev/null +++ b/examples/multi_material_tutorial/uniaxial_two_material.py @@ -0,0 +1,167 @@ +import numpy as onp # as is "old numpy" + +from optimism.JaxConfig import * +from optimism import EquationSolver +from optimism import FunctionSpace +from optimism import Interpolants +from optimism.material import J2Plastic +from optimism.material import Neohookean +from optimism import Mechanics +from optimism import Mesh +from optimism import Objective +from optimism import QuadratureRule +from optimism import SparseMatrixAssembler +from optimism import VTKWriter + + +# set up the mesh +w = 0.1 +L = 1.0 +nodesInX = 5 # must be odd +nodesInY = 15 +xRange = [0.0, w] +yRange = [0.0, L] +mesh = Mesh.construct_structured_mesh(nodesInX, nodesInY, xRange, yRange) + +def triangle_centroid(vertices): + return np.average(vertices, axis=0) + +blockIds = -1*onp.ones(Mesh.num_elements(mesh), dtype=onp.int64) +for e, t in enumerate(mesh.conns): + vertices = mesh.coords[t,:] + if triangle_centroid(vertices)[0] < w/2: + blockIds[e] = 0 + else: + blockIds[e] = 1 + +blocks = {'soft': np.flatnonzero(np.array(blockIds) == 0), + 'hard': np.flatnonzero(np.array(blockIds) == 1)} + +mesh = Mesh.mesh_with_blocks(mesh, blocks) + +nodeTol = 1e-8 +nodeSets = {'left': np.flatnonzero(mesh.coords[:,0] < xRange[0] + nodeTol), + 'right': np.flatnonzero(mesh.coords[:,0] > xRange[1] - nodeTol), + 'bottom': np.flatnonzero(mesh.coords[:,1] < yRange[0] + nodeTol), + 'top': np.flatnonzero(mesh.coords[:,1] > yRange[1] - nodeTol)} + +mesh = Mesh.mesh_with_nodesets(mesh, nodeSets) + +# set the essential boundary conditions and create the a DofManager to +# handle the indexing between unknowns and degrees of freedom. +EBCs = [Mesh.EssentialBC(nodeSet='left', field=0), + Mesh.EssentialBC(nodeSet='bottom', field=1), + Mesh.EssentialBC(nodeSet='top', field = 1)] + +fieldShape = mesh.coords.shape +dofManager = Mesh.DofManager(mesh, fieldShape, EBCs) + +# write blocks and bcs to paraview output to check things are correct +writer = VTKWriter.VTKWriter(mesh, baseFileName='check_problem_setup') + +writer.add_cell_field(name='block_id', cellData=blockIds, + fieldType=VTKWriter.VTKFieldType.SCALARS, + dataType=VTKWriter.VTKDataType.INT) + +bcs = np.array(dofManager.isBc, dtype=np.int64) +writer.add_nodal_field(name='bcs', nodalData=bcs, fieldType=VTKWriter.VTKFieldType.VECTORS, + dataType=VTKWriter.VTKDataType.INT) + +writer.write() + +# create the function space +order = 2*mesh.masterElement.degree +quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2*(order - 1)) +fs = FunctionSpace.construct_function_space(mesh, quadRule) + +# create the material models +ESoft = 5.0 +nuSoft = 0.25 +props = {'elastic modulus': ESoft, + 'poisson ratio': nuSoft} +softMaterialModel = Neohookean.create_material_model_functions(props) + +E = 10.0 +nu = 0.25 +Y0 = 0.01*E +H = E/100 +props = {'elastic modulus': E, + 'poisson ratio': nu, + 'yield strength': Y0, + 'hardening model': 'linear', + 'hardening modulus': H} +hardMaterialModel = J2Plastic.create_material_model_functions(props) + +materialModels = {'soft': softMaterialModel, 'hard': hardMaterialModel} + +# mechanics functions +mechanicsFunctions = Mechanics.create_multi_block_mechanics_functions(fs, mode2D="plane strain", materialModels=materialModels) + +# helper function to fill in nodal values of essential boundary conditions +def get_ubcs(p): + appliedDisp = p[0] + EbcIndex = (mesh.nodeSets['top'], 1) + V = np.zeros(fieldShape).at[EbcIndex].set(appliedDisp) + return dofManager.get_bc_values(V) + +# helper function to go from unknowns to full DoF array +def create_field(Uu, p): + return dofManager.create_field(Uu, get_ubcs(p)) + +# write the energy to minimize +def energy_function(Uu, p): + U = create_field(Uu, p) + internalVariables = p.state_data + return mechanicsFunctions.compute_strain_energy(U, internalVariables) + +# Tell objective how to assemble preconditioner matrix +def assemble_sparse_preconditioner_matrix(Uu, p): + U = create_field(Uu, p) + internalVariables = p.state_data + elementStiffnesses = mechanicsFunctions.compute_element_stiffnesses(U, internalVariables) + return SparseMatrixAssembler.assemble_sparse_stiffness_matrix( + elementStiffnesses, mesh.conns, dofManager) + +# solver settings +solverSettings = EquationSolver.get_settings(max_cumulative_cg_iters=100, + max_trust_iters=1000, + use_preconditioned_inner_product_for_cg=True) + +precondStrategy = Objective.PrecondStrategy(assemble_sparse_preconditioner_matrix) + +# initialize unknown displacements to zero +Uu = dofManager.get_unknown_values(np.zeros(fieldShape)) + +# set initial values of parameters +appliedDisp = 0.0 +state = mechanicsFunctions.compute_initial_state() +p = Objective.Params(appliedDisp, state) + +# Construct an objective object for the equation solver to work on +objective = Objective.ScaledObjective(energy_function, Uu, p, precondStrategy=precondStrategy) + +# increment the applied displacement +appliedDisp = L*0.003 +p = Objective.param_index_update(p, 0, appliedDisp) + +# Find unknown displacements by minimizing the objective +Uu = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings) + +# update the state variables in the new equilibrium configuration +U = create_field(Uu, p) +state = mechanicsFunctions.compute_updated_internal_variables(U, p.state_data) +p = Objective.param_index_update(p, 1, state) + +# write solution data to VTK file for post-processing +writer = VTKWriter.VTKWriter(mesh, baseFileName='uniaxial_two_material') + +U = create_field(Uu, p) +writer.add_nodal_field(name='displacement', nodalData=U, fieldType=VTKWriter.VTKFieldType.VECTORS) + +_, stresses = mechanicsFunctions.compute_output_energy_densities_and_stresses( + U, state) +cellStresses = FunctionSpace.project_quadrature_field_to_element_field(fs, stresses) +writer.add_cell_field(name='stress', cellData=cellStresses, + fieldType=VTKWriter.VTKFieldType.TENSORS) + +writer.write() From 9f88ed3c1831e0cbe1367764dc8a9b3cc1b35eda Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Thu, 21 Jul 2022 00:37:54 -0400 Subject: [PATCH 3/4] Write a tutorial in a Jupyter notebook --- examples/multi_material_tutorial/.gitignore | 1 + .../blocks_elements.png | Bin 0 -> 42278 bytes .../boundary_conditions.jpeg | Bin 0 -> 52078 bytes .../optimism_tutorial.ipynb | 914 ++++++++++++++++++ 4 files changed, 915 insertions(+) create mode 100644 examples/multi_material_tutorial/.gitignore create mode 100644 examples/multi_material_tutorial/blocks_elements.png create mode 100644 examples/multi_material_tutorial/boundary_conditions.jpeg create mode 100644 examples/multi_material_tutorial/optimism_tutorial.ipynb diff --git a/examples/multi_material_tutorial/.gitignore b/examples/multi_material_tutorial/.gitignore new file mode 100644 index 00000000..87620ac7 --- /dev/null +++ b/examples/multi_material_tutorial/.gitignore @@ -0,0 +1 @@ +.ipynb_checkpoints/ diff --git a/examples/multi_material_tutorial/blocks_elements.png b/examples/multi_material_tutorial/blocks_elements.png new file mode 100644 index 0000000000000000000000000000000000000000..f4190d0d62626c8cf96d428a046c2c7f497ccda1 GIT binary patch literal 42278 zcmZsDbwE_z+C4GCAfOB(-7VeSEe#6N3PUrfq;yIT41;tGNJt~8NQ%_Z-HO71AR-bX z&F{?n-h03Me!ma?;Be;5KC}1P`+3&0*4h(qWS~WQlkO%C4i2e~wz>%p4sHt$4jw(2 z05}ut(Dn)VgX?XgrHWHK#;^_iW90Zy$4Orw=Pqyz#=*lS$HB*50{r5{y5fNTI>y1_ z0{-IQ;1}Zl&#Q%a|9%y>r4aw$$9VME3%eJkNP$KyT<<^hd8n@^>+l37_!#=c-cj&5 z>?!sF9Qo(6z#+`h=P}E3n1`pg>~jUSzpju4jRUapS7ow4mJipc+Upt>A*qi^? zhxxmwf1L%os(4dg=s&}zc$3!lD4Z{Mg@C4p|z}hvVUr z;beLu)g%@_O%um6DFl~kDSS6}I%5~j{(K~Ngnz^)wylB+5p|ux333Xw1A8rS#Jy0a zGCD^-XpI@$*B(F8gGL8Y;ojm~lNI1&W4{QDX4YEzYIE_KN0S1U5YAQz=D|x)HXtj< zK$&DoHNIG}|5CafzkVwp+9X66DTSPQu_

)cv$Tmq*SU&4ky_=h@3AlR>Db&0J!;WU2IF~BhW3kr zGZ}pQW$mj$EZ$!HjZ-gW6&okxWsu#;Yp!ZgEor6hS?&h@u((x3%6iwJ+m5o6W&DWV zONLYaA;tmI&B4REH7c2_y=%hHQ)OVvo4A}Usse8fB#9%`T-=%~kRt*HR~9ebSl-cT z@eq-oSnUYJ2eBDbqIW`B-15QWg+=#7|&Oz z5H@fqO4a<>I$D!f+vE0nq9|#$rdqrJj{PoD4D@ortEEsjzmffE5(^}UDg^E{xcj?B zZ#<(E&(4yER`S++Q3=zGcS1(-7j#y-cMl$#VWPG=_xolaAxK;ezt~wJLu6c=RNuQZ zYp|GudFnUZWPC+e2_n+c49Ikhci5x^Y^qhQq6Xl88irzut*-?rprXqwQXOjIEpq!o z!VfbrWSZLKk0rtvTw*TtYSGJ+rmhjy&AfOTYRd67E_}6z`Wl5V5LY~3sGpIHCT?ow z+*wqGZlx8EyhyKP&=;U-MGjEXeyq|=TZ^4Ek1LQD5uL>mNKusOTzs(D;|JoD#Lv6p zrNqNkT37vAlqMjIH|jLto*?ct=q!+(9YssNjJKdVl1X}k(JLepmtTJLvgtKN%D}Hm z-=dHX2u#^dS1%zwJUngpxxwd5LaA_!xM7j{m#fF^!kT7f4_#sPCGb#!*6{iq&Kagw)5^)v znlSgnKYJU#sNUY|HzNVkV9q6_+ebH@DYnk?nX*G(E|6xP2`XfR})h{=%*{uqSt-oBnUCGAKqyAUNRP`Ty;Yh{fO7I|MRM+ zu5aa&$RC7oXk6C3i}G;dqe$o-3iGz%sf8>U^d_aRSRjEYe*JGE0lqKMimjLL1g~jA zBCg^|(3ZJuM&iJXNoZMl1s%#^%<$p6z)yAWWD=KbVstU`LdO zH0LT**zw$Ms^w93kSj906;)fSpY8v2z)(!XUMyk-zM5{Qff5bKOl@4Uxl{isf@Bm? zwb0d{F+1-wp#+3IY@fG#zBkzl;sHaILbjO8jl91!O?QM=u`1?83aC3#1ky>_ zcij>zfhi(i8UEZk5e@eGLlC{0R`0N53>l;c9&dn^asD`;$|L&RstS3FG&bfn_KFBY z&p?;TBlOc1Enzj}%kP(DVYY`bQt6kpjREJHt0X$KRU7-(z%X%OhpA~vc3zCVdvt-? ztO{R@gRr4&ut~-Lt!9W<*(ze4D?wy<@W~o8ol^Er)M4sO)CsxEo4dpd9ZtD!9(RZo z?vMg8>Q&9Vu|G55BQjI3;JhqKo5!f*S;g&aa%oSF5_)VL!Au*ITz5$saTpU$z8qjB60*9-5K;45f5t|w)Ps_O6gI~#6{U{$Jmm2> zkhV?VX>bkpsI{Iul{{*qYr$W68tnCB4{4pz{@E0C5Tkq|B=JOS{*(ug72I9dUzl2| zD^X`S%IHb2p$#{@wWz6SNzf<~u#Z~L9b7|UUNEE%pB#`uwHrB!U%HusC~>b0^CdNu zA=(;+%Q`)~&gbV7CUVK2vM_1%JCVN-#u>$^YRbL?Y1BYE@9plsh$(7U1~7(dStrf+up7wTq+7dQBjSxW8w`axc{>dLW1+-{i%$${$o zKT~7#>!qEBf_Ws<{Zl&$%C=*h_ze%ad^EX2663rI|#fJd!jJ}pvF2Kn{K<9Ixg zz8>aquy*cck?~3`$aO81E@zqP14#gtqD=W5w%SFr-XmUSx^U}2mq=lk*EaSh26aWWMhF{nT{cWmCZ zQ}vy#AM_TQq|Kw${Y+imJ;f~w@#dI*-ViS{vy=qo0_;jT{QiB#03JTunR7|!#&?P7 zOCw+_diJrIQp7`a zW$O!B@*K)Arc1fJv3Ft8V-ZuadkiIx1WX=eQh3?};gfMuH7*)b{n z${2JrI-J}p;*Qnp?^@wk^$v4JYBkO2e6xDENLqTLhOqermH&J>qj$X*f77@d=7>b!B#Y+RZbX;8ft;(Y`f8Ng3;mNOc&SsfM+0{-HWZmPc71xunLTeIq_)I30`e%GUCUa*5G z_YBQnjHuVBP`W2(Y*!nql+#iI^$W84s|2AY&;eEK%7kgV!hR*nf15a z1m)MuX_udp8Ge-d(YAglgphyIDRu{we9wGq?7+@C_YBISO?9;Co)oh}2xt6q=oqX) z<3s*x+jyA0)%70@;qcjx3vnqCrE!x7Fz=?S+MG@Qs@HO}4O%#Cd1oSa{@KGWD38OT zg}|d_gU-aoMk`{%^<<$uMPJQIOe?QZginV8Md1RIS48Oc_*g0>{B#Jyp?G%csO{Dz z80OMH4P@+8G9!DwDYqSKF-RfNvr1b!gq4}In=O4x4mmrf%Onx{TjIJ?3TfGhD~2&p zqme%FJOazGDMYYL;OTgE?RGfoWtzW&o$4Gp&Q-BSaN=tXaZO;1|2*@v>!<=_$Z$!S z&|hDDe>s+n(_Dk-ufbiW7bxy18uT+NUWNDR3O%@$^hB_i1Re2q?Z%?oltbBI&QEur zPC_raZl!^}*oRsnYTF?iUshPjPSWK452vWjKx5I$JUTE1jc!5;rgdscy_kuE;}+Ql zsYZJ+H&lp@d!NTPeb!|>i`hM*w*^0HLvYXtfCTDTNbv5v?6>=%htlFD5=OY;HF`?c zU+P@G(RvxFx;NI4b2=>Sr6rFr9W0-~T5pCiuyecPe&Aw=eSq(&v6B`vQjFJyZ}24>|w8*)jo)-3&TuLe_t_oxW=nD4_3N#wDvjwXT(4%g6r{i>xM0@D`w z-o3n;@~Q-VAhj>&C9{-n!eYGaKU2(xAOt?v!J7ox*vuH^!Ob58mF{BsW`vK=q3_6O z39{cXjP(q?>CdgZPs;t9!!}h4`YC!>thbAy6QGe`ZF(a&$cP6wL=+Cc_y1FfA}=;ot!mAjo=*ULNz5=R8@>-j+xMb#aOGX_kg z#T}M5rx!w*yBy2+F1L0fh^EY=M$a(s=3#f}Z!`flggQ_|#8uUL8)$TEd5}i$Op~Bz zEasTMrES}fOGqea8ZRAL2S*^a9F{fjub(gE=nKFD*cOMF;Bd&x4W5h?-9IB`PuNSs zC?g_)2x0c-#{9MuJHpWnC?9+&DIK$vJc5mRMUGD6cDICCKQ8&Ws$*gFyANme)mbG2 zB%<|@m!08hz1)9h$R5Wniqew;Q!>vG81b=h@3*bG?t(<~Oaw|MUS+f-+{E!vgA3RU zBEACGI}c7g{!ks9;=lIijYADv9c!*qE@uV5eKG9;6S=Q|R%{3@=&g$>=w=+S_tr z3)WDIFo#p&v!Ru-?{Ye}b-Z{1D$4N!LVOj6&r5+Yx)d7IAnN_Gaz;?6RAJ==gOJWX z<=A#9`%^p^JU8Eb zo?MKB`3pO*LyWp2nu_Q4UFrjv;ZwEtUE|0t=nk91BV2x$k=OGshvPFi1k$I--?4KB&dyU`U?@SgN)|ym!J`$rS zx>_K0hS_p(cmZ3s4hYf*P!0qeU3%{BFmzGj@f$a$ymz17Ha0tOwOeWHS4@B0010Mh zbhzldWvm8wZS<4|CL?;P^PSjFmV#6i9ctw%zJJEZiTCd|HUu|zp)}>lW|~KK|1_fy zybqOp7^@iWG{~66YJ*~MY;=XA6NQS7TprxP@EjR0#Mq=cYZv=6_UV4gxRbYbJn5~$ zl7g*}pFQw35Tc5Z)JyO#uVzF(?o~VimFec(-JQb;Qb=5MiD^&7^t~-M?8tpfxc(I* z80;kwo%H5;1r|aVzEH~r2%ZwQi746wFjrr;#XxJuQqc4DlX7k`2gTN}pD;=4n=P;{ zR|)8yR)~VPkaK~zyYgDvxkoQ|CU1!pRYV>Z4;qYnh|h~5Q5%>GP3R8LU^&n9w{9VX{;ZxfI~1di_T!W_QRr1_?pnk<{<5jY;{O#Mcp`*6ogKXqpxiQaoo z|M<#clv8p@Y~;F0v9)^tcL6qd$t!ip>luG51V^)fLC(xgm!UCA3b~hJ&(H{`!Q#u> zkz-RAM}@`Ng{<`JN{=DR7AI$k%mqK_XG%)ORXPM`-5<&bHWtKXq5vMAtd>T%N;jwE zW{kC}it+fIm?^UVR{deM0*|pnT`40(%GL9k!*A5Ol_<55QcYrKi$a`Gj%iJs6J!zy zq(Bv0c1=!+ONd()GWk{$NxvGq=_B!{_3IMnN3%xmL_S)@)}m4hR%Z109xJ#qn$C^6 zD8|yVF+?teb|=kYLOaKBOv1MeUNuh@AyD@-!J&2pa?x??Ny&;>P*SU*V41T&wC-J7{^+!6lVuhGh?MQLcV;1G?(cfIj^?J@f0vhR9wYGa5 z(^mJzT4U(k7LDN74pc$)liN&!av;{4DArKJpumQ2KZ?57D5?55ZNaxl_Y;dqf z^om1FvUuZbp5^B$jTFOLAvGB%;?t0eJ8#iw!b$V?mrDSQUMoIAU-E|H7qJgiTn2m` zg8EGXkmVL+@M>f8X^$)?P^Co(y1i~TknJb2s7>}PWXfbGYS>@cISTPQ{P-n#|EjD` z&Et1D?{bGsfE^QwTe&{+t3*`@fM?lxw$sj~KXQK?DuGu31KG}2NM~hhZuZTW3pA9j z-7>e#H615J*D&rr_2LPPm?$L#V2UP`hd;{Tt1tocv?9AHvn->3E_^>#v2{N6_}$I- z!)_TQjb~d5|FwFZRl`fBcUU6shEfS}2o(swnq7E5G50VCwn1rye^TAmC#9z{IiE-8 zM#1DDILSOfY27Bu_Xi^z>=o;3LpTouZG4F=;K940qO6~K!@?!Mj2YL>x-C}#K=us3 zh95bRlT+$vFe0oJ9+|e#%lCM&U%R-K|NeMa%r=#(LC2GDkR$%}P%IOQA-2KNw*O6_ zE@Y4sud;r_%fi>hStY!8&NHGAk*-uyuhiTem7h$+QZdsI;>9ox0O?bTELJ;aKAxeS zG`c9wWW=u%)`aPH?v!>#nPe5V8Gsta#&5RT_hT|t5u4A<&F0HY_S&S%>qQo#Ef}3#;m9AH`XxhxiZ9c)NLU|Eh!yc&jMN+)@|)Uk&t8k=-wRHO*e}A>nz`{EG}~ zO*S(JhrXS#ge!vg2f3`L4b_#Sozup0Sqty-mHz@sM|PnQbDMAS5L|-fnbdUD`6LY^ zngwY>s};PU4c){STuXJ*cKU_tFv{?2OdNj`MvKvwnZBKqpI%CG>6N>|scxJRz1@^` z^899WpB|eekr#!!roCIgyqR%PjoU~HfI){d_vuu%bMR%@Q$!X9f>}PUi)i*119;&D zmUDcnTT5Y}7GBnO344>f`oJm1g`K>?V^9|wgl}=GhUpfEWWg%hxb4{d2SZx(+LL1= zTHY`Q(P4a*I{9Oz>E-XLCAYG<-kk-8^4|;&5Or@_yRWoq!URfulUS8RFZl{aSAISp z&{Ama@p(8grMmoB9&Y2Fik1|m&IiTA`VJdx@7=`#5xH<{uyA9O29(v8Zj~v*trvYm=2z(RZ5tlqaosfE~xl5f`-ZQ*6!%4j$ zaX+;NJvQil(np@xQ7A#rCBY>FZTS%``9W{ces?vI^ThB=~)ZxK5XhHagA$AKWQNs!!GKPY8l{LWNP z;egTr!2P%Zw6i~Ud^~l0RaXgci_>qAJoO20XZqE*xP``n$V5C5NQO)p*XI0WExIDX zX{G3jtwt_LCqLB65l~yoKe1^ghsFBLDgwwFKnla3kl&*_N%5-VcHN!!mO>l+;VSIn zrZM7H(c(WFojD_?$EN)r=|QEh}t)teDyNH>X~hIb~ZZc0+nII|gP^*4$S zogN3!)IB(WGCp3+UhUn>8q~>QU97~lWWv{{d0>&a9P_-3isf!KI8O}=Jior$=IfFL z1#?h+49v=+9OkN0TIGl)Nm23e9(2QNcAt0{$5jPm*}cQg5k?y#`K#uEJTCVU`qPWB z-SsUd(UjNA#}tiRPSgPg5-6zVJ%*H3b5)xeUQaoD<) z1b0sZZ6MO)*s1$%cG^E|8w|~6Ym008X4`ZYOpc?9qP5B(*?XG`gFbF`)61w`6xDSI ziq>QQ8|+I*FMjNIFzeP?OA9hDuMwH(b#B2tUjjIk|Kwvv+yO#K#lgwP110^HA4&%A zyaMh&FWobTlO?!lP0n#z(#*Y|&>r#9Y7XQ|G|I${uBW_J%=P4BSO2D%l>wa^7G^60 zF#C2@x%YEda7i{*S~aVHMXNkNkDEkGbu2xK;Xu#ePnB2O`s-{q*O+$hc}E#p1`=3v z8-V{a{AqdLmz{MsS9eB5{8onyvI0etjP?EesbQ_c=>+}n)n^yN)~lQqP25b($Ixx0 zS_~~+5###p{ak7NIRLpxVUcaiD_N9zq%uBlQ`RCX zw*httv7)Z0e_8RTSz_oEx)8hJ^f^t79e~(#;RC`pQ%5Jyj(B5*N8JLB+nA^#UST=> zm+5l>1+)Okv=UBx&|I|H_aV4!5Fnj_(xE1-|1_I4;Mv-FQ&E6I#)DQm7<}bV>Wmai z%MK6iqYy7jGdK*SdTGJ}%CobMaH=tV0H8u;;`D~>+YO3R;%X-Hy60@POpo2Whv;6h zGKzXfX>oPnQ7}zcOB;Bkoh`Ymp;jU~pvj7&4BORqlo1SAMyI%MrlOq=^v)lSxvzsf zAXF$>jxkJzZk=Fwww|sKuM2A)U~#^tT@sKKx1{~P-rlF*97-v!)BIU-<8_kHln@q{ z0AT@3AlNeD&9wVHit)LGvKBX1kG?Y;8*~(==YxttuLGNC58_l%uoKK@#zQoz=XZb> zR0c-%=DRJ~4Ihv;)Kl&Qlla+`SzYuu8*#-}-xEx2(ivvuS-=UKob7P0&XYUDKq|Jcc{?^2+f7r?SS-Co@`|n)8rE#Mw8DIXIitNd75z|Z zRI(t9(Jmte@k^R-D{4_=#qaH(0$$s|-ItaHaxn#bcmZmUj?9s77q^2cBTU)4K?5xS zSN|`ZPodQtmqRY`h5rugs2f;^c12Nn$ex^T`mrN`l26JAMcJ7%I%U%h1>J*zC6B?4 z_$tJAyKz#62XB=eps{JhRZHErP5fpC!L*SuU~owG(XC~Ju>hG&ZDb-|)!I83ZJ1mN zGuTNRrnB4o=hT1K+mC~=m@)FE0LV)OpEuq#XNVUcuA)p;ed-sQ%`iZC2}w<|+461q zw0({dMAvqdXaOY<4lv`-#+4yY$=)nck0P|rZL%4O=l$R3nEurDjvu9D)#zEDLWF5* z#AFyy7V7Aw%kda|jC>{a>RWLLcks?pSi=#z-b6AH;dF}OA}PwGj99@|MCY7yL|9A2^#h|AfORSA zqo41({@wS#A);nD@0$h7(E?tQ9LJavzXABN0)vhl30)C+OS_paLH4hCx9E0_Fj}Q+ zy6dM#YOu(-R%mmZUmh8KvxFz%MpWpA<5a!tR`MyE6PK;G~?J=o@gGt`- zRlvjDF$NWv0hShv=#ynx7=>N8@#506K2R8k#7qo2ksG~Om?}A#ltc5UJtXP9yM>|* zpO~wBM;RDH8}Y-=v%tC5+Z4pbGO^oh?^Yi1rR=;}_=_waNyjoA}kyjFl`5y%{_Oq1oqW5hg% zW?fJ9{itSTKtcO84>tsl+af+EolaTnorBO8*Mqf*-jD z_N$Ygcy!c(F{BZbJ{tkK`D0oT=^J<@6&1Yq6%&`WZS(HTNi)0|cUKjtiyH3{6+I2Q zGJmG*@qEEnoUimRwk9U3EKkHGIu&sa7CF@oRtqHyiFy)}E~23sFK~mVBQFT37eB6y z+i{nJNyr{etwnv36o2#3jVX=oMR|RyRKvl2g|7h*4YB=lpy2Q*B0S9D?L`T&6 zJ*8%h=euiGx4ICTkaWq5L0?U^hH6rY-WOXmPTFCOmIYm}kBvYFF96=8bI#QHPymt@ zNhInX)hrkMs=m)@+)CS1d&m+F=lGrSnuq5e9In{Ny??XkqVX-kfmOjREN&kZ0 z$gvl*l}+_BZPWLfL~>|qtdg`sr9QB|8~65{C{s#?(r`pMzX>sI1&UZJbY^&$!1J}T zQh}~GMJW{yT(&7ty<==)$V*laeB{yT zy*lvDw+Y2VFhy4#SSMFvn^w#U7Zs8IU?ndfkyHnw%K1qu3_r8h873+5{YSb&w22^ zVJoM@7W=c#I!au#-eU)KYpcTy`0o)Xni@OjoYegEni_IHOLECD!UR2qjuxN>) z=S`?I;7BnAfcnqwyumLB0v`BfD`)roR$5RT83jR3G+FZ#c_1$3NrAY6+Q*mqk3y&- zGNYBX%e#106QT!_lAy$?U;S{@V)z~`{^X}=C6BvXB?3jIH&oPEe6J>O2>=G`Ri!`3 z6gqu!&`?kWP)-1#q_bi&fAG!62^36s>rc^c(Q!_^Up1G{dfM0FH~sXvx{>_qX|n)y zp=wx7nuD~XqK{&K1g(E=zKw`F0Gj}ylnnr-X1x#QWnQy=In)8d-uLMPI=!gN4mUvy zeL6t3TyPoYcWZy9lnYaY!>6AYuznnFU~eAWOBtrq1(N(XD(usH?{)UM4iv%P$To^ti8b4$W-T308rAZcIdCEq z?@|8SapG4`l?y>gy>jMq>CItj#4{jAt!Sk2WM5W~@ZsJgd;sbk06tkr3bmHqt`ZI| zeW}RTNxAVMO=2>Y~Us&|% z@+R_Shh(UHPnpL~+Ik0Bw(|mF^wRC^SG&ehEt8p_TXt=XAvZ1cS|BIT9}8Hj8>`xV zHoO1r@E+NM^^=_Ffobory_#ZBhdRF5im&7R(F;s{Rpm1+V&MI6(ZN>uS*Iv~nkw2@ zbRrrHoiAoHq6Ecsy?ZtC@LL5q4>w}{V*$TSCEUe&jxt+v1#kjXu0&=0Oct}rYCrk2 zpW8sjs76RH$*l9t9Y;UwFa*WG`M`mG?X`2CD!_d_+^K0)z{j&uA8oYv`!qQxon(0K z)@YWyo~8@6ejm8+$TYzB3;COD^Ljf`%Cr_Sp%}lKnmgIw%{@@7D5^(mAjF6FjR_bR z@2+Us2y=KhK*;hMnyR&CPt2+1xAZ^Oa$@+4D%Z!*lP&uOC@vrv58A(J`|F8e1}2@^ zKTSuPYRYa?Yp}$EqpwZ8r33PjY+uyEE!2_`x7-~`ClQKfUIfd52?wHxeA?-QhAhBh z6ja!}vDHnB4!eh1;at>s6E?kS5d^SnK$TU4lUQvkFNarM)c@o1_v!Qov2!F~A2CL) zoLC7XCt+dAq!KO6T6V3OV&gh7U5Tq1ry?lpGQ{WmAg^TLxr|XieC;m_Dlck@ z6WM;|QfycokmSLx-H|?<#oG&8-fwT`OEE|ni7wN1qIsOv^_S5Fm|@dlL(Q~b@hC@+#?g{P7L5B7~OLt&S zoi(~1;%cFHr3`k{hJ(Ygod^dyaqEup>t6?CPoI+-!*yi&bH9n~@(Kczc|d)j(Hm%E zc407)u)OH_@it!_?G$H=(s78FQhcXgkVHoo;+lGIYRep8T=Lo@zPso2t+KrWxE?@Q(!b{9 z!MlYAU_P09BM2;;t?k0^?GO(i5I`9wf4AtefbZTSXA%n#Q9GW_*Sy*B!~h=}^R+>rju`bI>-UX8R? z_zbpUyl;?nZkd?O1?&8DykR&0qpJE(yB&NACzysNdEFRI5H05r0S?m(~BwnyRcs`eb33PQ|c-)hnp8Zr;2*BBz0M4Etu1b>= zz)NzI*Ay?YQWJoY0M?b}_B#m(Q999K?txBx(cSClK$ zSVK-a`g`%;8R$O}L(DPrPrnk`-D}sQ0R8uadS35hqX)zsc}b)0aCb&=N-hbf7(WBl zEeqb|`RXU|WxR*H=hTzmvH{Ncb?>IygFPVPU)M?1iOmP_;yHjPNOWMNF&Y-u&)mx4 zzZ>mp`up1vgmR4gyZ*&_fP!EN|27pHa*Jk6W$Y=Sw{{wFX zf!oE^mLKlm$+o?ZHPQjZxXb&Jn1Sd4Rt@!<3~JHNV`D&)RPHtuF`S-l<-uD7Bl5}y z+{;Vjj2w!?-~-?9FBqa{7l3L5ztXKVNii^+j$X;G*ZgtJ(K&`k|8G`2;*FgHBz01T znml55zARGgv@ws=zuI?MVaE3G90M+&UFfW}O@9D5g;w_^Z-x=GUd@IUqJzT?lc$pM z{^r#$GwrbX-cY?TBuUmySmURcx1^*+CWO^Bh4PD8_$|GZ;|0$T9|)Z6K^u|>_q^2n zr@EkcIy8Dc-Q=O#p)Yrv$OIO_;49N9qZ4Tz70E?A5J>+|ib05_>91}(!ti0CcjkRDaz$#pe^ zS(W2~X4$RJT`FA+rcct>b1v?r*SQ-9@gRFtR$Mxm$HQ03;JgPF{85JI`KC~57qQwQ z192q=yKz8V$FfI=7jv#>kgN_P)!pU}G+9*^iC*3Yn>B!I;ND(ps-!I@yTsbWxhnId24Q0r-I=bW^a#C74iz+ek zYqdl^cv2q5)<|rUoBJK@UIRNVk{{Ev_h5W_lpYdE!nyLazUj$WRbeIK zSiHYh;$&IdV>-^4#FfPT^xMfLqIO9OGRO$P|5huL{E4?JEH6u1=5j@YiaZX~GE~uI z^<3x6jEnh2DI{D5ps47kC(0g~uM}vgm%*$=>@E5h2>g?ce_9~IU#Nw^>09Z~wL4w} zx+F!CaD99&P{A5GrC3CMDk14~`!sh2F%KO+`I==VTA~oC^b?-MR>JDw!mRNH(9wGR zY|RZ;&(m4RR82_!NM&CLnb&jT;BxT@(oNAE+vO&g%#Xihcvb&QPWeX8CviMoHwaRU zkc?hGLcm1I&rF~Cb74Tn^7zzpz0e`p9Y#h~tc`VcyGy?J z`k+(GfV|@_SCD2^-q|1GPO@JIz-YrYDz)$@!)aeeH#=o^08@@UE7>nGMZ@Vye*h>l z;whSPj3NN%{)|nPPP5(PW$uqsOTjtnM(?R`$hfml8}y2P8wqBsqgU>9XFL#DzI9B~ zd;H$xrj$$UHqy)l;Gh+;$Z$aRJHIVCo+359l)09$!r1Yjw*7irpkQu|61?m65uq(m zw;D{h{ONPsR21`GvrEvU$Px?Ds!8rY<|mj8L~SCi{y)bZu%2-6NP70STuGxyOkL{@ zLj)HuV9Ys1Vse3qtA-dcYeS5uEV!Bm<~#`$HI)=Lx>*9;@Ag<#n*T*=5O>Zc8i$yIc5t4pmxc;07>t* z9^#${k>J_e=h*RjGV0s!tQ@{Xr{lV)54L*Hy+yP1_KUm zqp!I<{nZyob)S+I87HTLgFgRcJ68(va5r8{TmPacns}1?_57YH4ZthX zZ);vz-nuD?I`$0TD~G>OakzN<&XSZ#a3B4l5KC$QFkW~R(+(Y)9zofAOAx0b zSFFW2Q0aH_dV;VL#=mW#9h4xN#dmM485g5A>9dyT^1bvrB&I~Ofc&n6Mf2eK@ui!> zqD-9#l5Gbo6>W_o#BvP9L>PZp%t|wRa32){s zVgULy7o8YJ-6e1c8J;eY%HjJ2lv4Co_z&kavBnt=aMIO1N-B-QFW}zfp*fBE@0)-K z`jpvH(_3AyZVo_>x>R+54?`m~%0x^rw@>)GdXt=foOaU?l}C^OEh_u-KQ7Doa9tCt z9UaoMvWF1)3yMj*0&IL8eE?2X!;4LQ#>4a zi(Lz>K>t%;B8ft-St{<4kndIId7Hv@-uHHD|f&1NU1p3FVZr!ZC1d#dK@+Yplf( zh=P~oJ0Qa-9Co@zI_;v_oqiJ)(&nZN1dDd5(DcW-W+!~7!6fj#%*bksXF!yeIs>@M zf(&K7dY%uJ0_}z1O;>f z{%sG%-E++MZ6&@c$iar}vqg-5a2(6E)nwr2xh1)WTFkCU+-{BsJ+oC#VY*IR-CwQ_ z#1-+#PLi)}f-(V~h$QrFuJBrs?D3D__q0g0hZH4~2g)9B{1eQayrPVjqXWTp8=y$6 zZ@=lcn{=}y*?n#Jv+fsKAOJ9pI|cJ!0%jyeEG7Qgx$e*u)K9kMXCHmS-O}gP6|?j( z^Mg5(g+7IG&^V~d8Q3Lc7rXG?*!0#q?J3|=k3P%&U~=_vcLyZ9g!+?v4tq{d9K0>s`5QX(PL`=yK?sB&VCG|Gjlch;5o+u=_F~q;Z>%{zC895Lhye zRfPNCQoWzDkvpKe1attx~B8X+aGNo*y;D#+{#04*;F%|WuI}l`Gdt( zz6|2(Wp+rbcC>t>?MSmS8$R@#xs{+ihE1Z2k>7m)$E{B6`9wc}li`CS*Bs(p@=jIH zft@7|Zz+CV5ZxS&CZjx`@(75oS-4H6a$3CAInS*{t7-G!Ty$On4!%|f7HS`r^=sM} zFPS+jqK^cgL_w_KFanE<%#Ei{g?6Ip=-q3djpqtQP9>}TZx2do;sf3rS9fXnRZTIF zo(uE&3q={l88;6@W!G!pBXbNOjk5C2<5f&T(Lb);Q&+h&2YKotgs?s}KeV7{6(Q@CNb3NDLnj8}#F1nRD zs7g>~P_XxsdZWZV?R?*nhGtom>gz{#MUR=W5$@U{qOsG*$Je(X#azgDay7o|0bDS{ zm;?-sPFGG|uLS!r_};Ht0CYda<}**UZ;vbZ2;qlu>5BSsVI23sc*|j3)AOgp^0J#+ z|C7WLUZmw`6CVy#-~4$5ySQ^%2UVR`p#UhkuzX93NArDz0AADug4Ti|<<`dDa5Hj{i`INz)5y9~2$K+BJ$z{7vuLy&Z=LKH@@X_g z()D!pzmaDieBhH@Qsnr(GLFf&NH8 z?QsnG9p;dQWwG(8|2B+wfG30`E8g(v=qj7dj4jZB8rnVys}BUFK^KHxb&c|P$FGbL zT(~5Pu`R@y;fK@=DYn6zqA~AY6uQ>#^%QGTaOfwBf!T>tzdHTACJ%cQ{ER(nhZ>8y zetoOdhZup*z=$~2rJr1Dlto|bjT?zwk7{7A?n)8CZ|2uowlJnS4Kiz{Om;aycYjx} zizMZxsKoq}mwhDOIMeo#z9x->ZAx94;3{D&&4ER`?j_AVRncVO)%8y&HWd7-uo?m{ zV%OV7Hs_X2v{MfN{<~F;45{}`MzU-^So3H@?%sS+1D2*aKEq1apS7IAv= z*=7$Q_edkcQJ8>SWM>`QQ$LWZb^`0UZq=+xd=6QSi@(_7H_082d&iE9MYIQYYs^E| zK6sV)@O0kkEYEiW%DaYMb|{}-<`KO74^_k=pj9oOS_*|MnhCENitJbGf+R|Ew6m!q&=uh}%nj^28swIdyKQMdg$X@$^#GhMb;D(>QuHeE z)@)|MYtKwKRU&PE_y&2J`0D>Cq5g+MHbM`H$1S03H0=Ih-Kz1bN{?1PlwMKh$F$}b zutZDIGX#D4G~J)-%Lkq;zvpnr0`2kjgF|J`FxsKj3`X#a2rk$0`5nM93lKNE@^uBB z=ek-f!<6Wf0gn1JYx3GBI(KzrgtSpBY{39hfA2KR58pN1r|(xfp5EGXN}hY0fkobz z1d1|q&m9DH_ztI3rf{hu_GExmuzGl7YW_SKLUZ@Nko>y{5L&O|frf!j3t;*?lHacZ(}Q$xwo?MMV35$uB8z&eg@1SK_;Fe4&M)+${)z>M%l|fk zPk6FTgH4=g^PgUdlQYOOKnK}MhBL`g40S<{6W`K4xBWtE^q05{9Hdf)c!ecDVUxnt zQ;Br8exL^SP5)=RR-hkZjC;Vn5ZcTuuqE9>1--bOM$B$+|`x7Cky!n=igh`CDm3f6&e|EdM-!9yaPXL_kxgul?N3`mL3 z1f!|U*4_pWCC-WN*4Z!*A8xT%FnaWN^6JRY%jFpkH(p}B?WDNhs59SIV4W11QOb=G zKN3$$t84Y)BmiQ!K&@1+rxm}UmlqX%cufpgQwCb@=>PcueDT8H-vzL%_g-P2?mSE& z^NqxprR^ivgY&D<>UIkqO_>Fb=AT+}CwwVy6;0t4lh?^rQ54M`-eO%3j{ev|($hoJGYr?QYG4I3unJ!|Zn&7pQfYjt2EXHICSIjm zUUO@q^@=9>kO=G_WWS*?{kcy<0}$6o-9IlrQqURWY;icjya$W5Ca)ltC~-9&0ebi{ z8L@&XXoSmySZ#+R<3;I^_L_arxWap=Hp-fa_QJgO{wG0-@tS_k2K%gQ{gR0T4N>Hv zM-UFI=EYvs>tARK3x&IQ*u^g}URC8SCjqeEq7pO8CSr|hY}w*8e@fwG^ub%jKr0V zsO3=PRs`k1wX&()HOSpK=jh1OAsKfk`yTpCK+j2TJ{tN0|6*A?r~uv$Z^a2Hvc zxD8lyAYS4*ET^Y$Phb{tbO_*!GmZkx+S~_->O)9T3~HC(N=23?Po2;cttCJu!3sb@ zzFyHC4mAkN{+k%f#?&mzMohJzb27_2zeeug#MQ`r){Up;e6@pV&u_@h%YqNbCDo8P zi!_f^d(JpQwL%?R2&D0o73PMLY~}#z?k8mMjA%3TKvqTL5d*FhW>uY3!p*ZKX6ZTe zxq^tsR4|cf@9PNR>&1FX%a%F?OWIK9vQ)q#-|5r&N?j0Zf*O5h=YKXlF|dHNHaqSNQaW5bST|O2-4|L z4lsa7HwaQvgGee!cMl;dLy2_vx8~VpTwjpn+Spa7*eOeNa0q8w1U$m9kX0yfPJ*2yV$|x(|f!c{_*E z7Lx;lKurX4D3%LCCsXt|c`t8P5dTn)N8`UEaRsX!MNh!5xc*ES9&1>CGdIhIHbM`% z2AQr1VDe+KKwzs1D~A5mQJBWYJj~pRt3{FEunF?G8XcGWq@oCaZjdJD)JeUcD-Z~k zHMqx=%oOD_qWo}lnD+3w?CH;_`Ik<_w?jzkV*0o*d!f1BKHAopZSF$cOVJEC#vB&y z#|X2>^}fQse|=wD)Q8U-dNPPfQtaVsEye2^;Z)b_>}+b8rjLM@0Cs@ zq(iL28CCkSRai9DK+xnAWXrj0DDQ$TzwV)!5N60PwEum?w zaA7v%dJqD?@7?JB-9D=1XhXI1{klK|^1|K0a_-_C@OqHkN`7tY(m)rNt3y~C zjZB$%7a%5Rr4omkA)wnYmd>4*5}Dqw;LAeubGY5-s)UB!ru=N{Od6f|BdajE-J;{d zXM?7F2~Pam`@Xl}@@htJ5aCgUYm9fXYt_!=TFw$~{TVFml)YD_E_frTu-x|OSR@7n zSI33$IPzOeDY{}8UKA8O!fy(3T9`}fJNOjx*&f;=h_6_f{aW~V8N>>2{(p+`(6~_n zUl-NpZ_mj6Weu)Q<4-|uMI6_-C_iaxZ*(F^5-zAWbZ0K?=5Uv{+^qH`+BH?Ilmr?qp!HaXVG`A}J0 z{L9ONmupL-@WU6FT0NS^AI4niQzANe6`8Ur6#N1*uqQvl0v6luWfL= zk)fq_(pivln{vbz_~ryDr^3gSZHoaSbFYANP1G4e5@5Ns*0%@DEV_Zf%yJ>{*hXLk z!p;sfR)ftB6ain!ldH~ZIo)<7I~TB1`v>7n?B3gL?xi?G%e?3Fk0pWRsqOFH?UCS| zp?|0TdLf6xx2%{kb55y$F5r8K8VJ6Tz{HN0gwC#ph3Qt0KtFq7$)*UP!l(6FiM&YV zqdFhts%4&7upXk}#0JRHH0D9DI19lHUxoUnFy!!_XMJkt!{CcV-f_azNAOiEMf^!M9&K0`q1 z%FGNsdA3jHE8B;CaQQb-%{P>#H63uvu=JF4-u5e)NlDj$zYt2NPK}uan3299lII9_ z;2J?XQh6JYLOm|0(1fU#5DET;R4JeD+17G*mS22(DL(}G`ptoi*VGY}TfWa6D6L5T z6AsrE5__N$*{7&!h58Au)luce;UzmW+cku1%Y zV!jMhU?Ia@>GPW}uk(=x-CK*&KcYd=4BHulpuzh1AXmpb_kUwlPC3uB@5tJ94cK|= zx9w+vn*i2!p5pQqT(T<-XLcWXJ!aczM4Dbl^Tf6@-2{?B*1}4^o@&=Nvs{88#h670nkYPOlsU~y@5f{ZQIvD z9}WE^p0_m4qT9af>k41f(J=+@+6}y*3R7c>f)ckU1YaqZh0Xuya(^K?GJCpKqLlbS+Bzr_zbxyehxQfq!gSmf)v!>jdLA5QP znuMTco8tjWb`{}})(w(+t95L#N(71l9A0VJ{VF)iZ$<0lK2U?+g!1*sBGA-rC$>)l#m6opueAc$tn{h>3S641o;c5IODn!$`PWdN~_E9$wf4 zmNrF~2sYZybdbV5V*7>bYBJ_Z40(nl44wkDeSGS*yYXFEkGJBC*B`XCPFarazN+rV zv1uCNboKgu@f{+FUfeI2B`bYNdtYE)DZEh9URu(#T@ohrNf@*4BqFPtVmz~R@YK>w0N$~jduAh<2L!kG z9^;V%V#O6bi*z$ld{rZf!fg0kifa+9n3xp07%Y6JU)K0mHc%>`Wn#4@)X3*5%> z@Q;N3^No}xu7yGeuB_^EV0Kvqe}*$G)lj3uw`=dvfCwUYJiNdJc%LEzmvoaRY*!m@ z1a(2&E9ySa{!ux@Kefm(z1sB#X<$-?m&?|hSo|H7W^Q-%V^X$NPNPq+y`1p%VgzOr z>3Cok9qvHiJpst8fKt_N$Ho!}mKtwBL!@?RJmL;M6-6OSD?~(ZfG}_HHz>>N1-hr) z^P)pLhWKqPp94-|VvK8^A=R*G6*lmI_FO4_RL{LhNj!k|?GJWY7SMj4cAOy#fX8|$g1j(#- z9{ufndGT-OOH{e#gr=L!O3>)17;q+Fl@bG{bjREV;L2M`R-9SdY`8K&ud<1>k?j0S zheysA&P4?Up7-#zi?80hOm34Y|Md;Kuu38~>dEHMp+&zu4C=i4x|9Uk5VCc(YOLU4 zCpXa&3Y+8-8mh3?H&L=u8QRMl z@MbQys8I5T%1;)=nsdjSu(Y^=dHYJ-DktBQ~zictli(&Jt$^U~1QhR`3Me|V5 zGuuv>FSA{d4S2~rzj!DmPS^p`-nHTF<(O#`=PqRb=Z34y)3<7PeXX{VqvF_MrE;;9 zF{6{$H;;_Uy5?UgoEPotCSU#b+Qw@jN1;ypxaEP%S~dhr0!pVeoh7}hn4?x`n(@7`W2knv_k`d!Q&HeXx$0Kx$?tS=%n zE!G#ys2&rTNCSsUl860uz?jH1B#<>AVaiRc1x%G2+@=rj-3Q-(Ii#fIvH{Hs^g)(XyWc%R5}*^jCaZfh9CgoF|DnJD;buT=Bv{N8s^+ zWZ{BZ(4=>MTNq6aZnGe!h&4`D+_ktV-zkA=QWQLOJVhJD&HnJ*K^EbZZ9uP>_A=n{ zo6wC81mqOctRU5M)%RnH7Er@*)vIL!cHT6xVZ-Q>NDO*FT_WaBtcDsC>h>^h#&=Dt z;V#G`*5-1_7?S(W??wdetafr zJag4ic;6xR<6J5JNPh9GrL6-O_xN1;+QSlb!FnkLTrFh2#bB2xnP^)f7TDg_BvJ-#9ke=ZI1Cl z-;J&kR#S8BeWlyC?IIjbMQyzS!9}$XT|AUn=V-~fH}+3+s0lrONiMn>*=R8$x7G_H zn^#>n=ue&=Yn`(SUM>98{tU&yFE^EFiMC*F*JWxG+SWeA%2waMUj}M{zXPd_?gSpS zQ-is&e-3}~tG$f@WVsnYO*Fjhu+#N%F61yyK#rHCx!F+)^r$%;GrE2$KelA1o2uCR zOU(wtUuj|g(p#26gbGUkV)h3aD64&b5UIZVV~(Q~;yf;dt@$&rRkq;0v1LWWwM%PG zr_rpfQI#UJ6H17XXz!eVEIAA8|HhKgfs%J|qzT|n6~QxbdzCQelR0LG2TS`Llfp10 z2>+e9tV%Tgn+GS)=eRt-Y;_?IhBp2|pVB}|lli2raHl4dPK+PMvO^zc7`KcE zzA>e3J)T7~@gsRC)*HQ|biY5o3&lM&( zBHSY!LBz%XZl=-Dgx*qS0PYcQZWnLWX893!g zYf_JGW@C{u)EB~H?ZW$SwTyiWE z$zJ#-$Vl!Gyf#g!ox-|7Y+nfYfq6l*3Z~P5{Dh5k$m|A;KHOtT2X-)t^0NN5QHb7J z&$wo)u`<~uRveY^5eSbYd2qiaVPwt#c3YZeOIJz%Ku}Wk?XT_nIY3( zc)X6I0*u$9J?UgrGG{t~Tn_ql(FM9Uf_!m4s*SxdOgbfA!2p|>V`(i#S&W~}N&_cU z=KF?Bpdhj}dV9b>6qw)dznud$5d3*2IIy~BQ6e3)svdnz7V())h^ zhk&J06P>;K<6VU-$DP|La0sY(39LIeRRT5owJ<4J`SrJNPbp?6Pe`;`)$r}-(`|Fl z$KO{qVOc1Dt5R3`?q<*oTkN`Y=UT!`vq))i{8{yTTs*JdccT&{o3G=HEaP0d2aeqy z;iXLJo&M##59CN;rkwu8e(5aNw{h!n9kUg1oP6m{>c{2J$N^U5rcK4m5B@lOibn<@ z)+m4S&&>M_*?KS)Mpw!3C_nLeecSw1119|Yg?^6RS4K4%sBftG$l2LCD1`^j;FjbV zjcC8Oli|VWzY|Lu@`2|tI|sf|&D6>JL!+;Xr&qc_;8e$?Cmiu7Mt}vLcFNHNkfg`M@0kV$qm1YwhZDzYxL zu<|aMBnajguME{letvQI@uiH#>26m3R_C9KnGJZe9;z`N{REZY)~P3Ut7wQ$?AD!$ zpI2S$7p{aJ4|LLO*$ZAab`bx1ReW5yA45AHfG>DI!d^M$ZQ2OU0dVV_^7u1a+NH$W z(~EWQwE7p6qmqMpMJ=mUgtGn#+Drip>+jK@PH)@|{4Z`QRZ=84CKO>5232$tDNiiy zUXAW1%E$ynIg?ygvKi>LftFk@>S;Xu%-S48!!H-Dgb~ND&AL&(sl&L2w1)5U`PSs43j_Fy4 zo^R{u2nUOOiWJ&p`Y0ua8(+f}{g(|>uM(}Z%x+mVT#}WD&0ikmmeXI(+`dTg-;Nv{ zehymgd_d%wo2wZT0j>5z7KcpC>3eJC$3#VRIrz^xzj^r%qx_FZ+dTP~dYMjuZwVMeNKMnidHe`^jg|b*7VUu2>ge z>3q?Cfumf`fEbT`_gnIE7f{q|uRI|a!*Fvjyq6P<9?UyT zJ!E{+{arvS7Frv-&(Z`fzlzNw`FFrGs0HUENyPa`;5bVPt$*cwZkZ~(mS9ArR86c0 z34mSD956qQIAx7u;Fl0gxCk3E58Y?iieGicO4VsRLs+;jT=^3KHr9UtnSgo? zyRM~qgQ0qbyQ&E7E4bSQ`w?$0mIin$NB97;;@`mOoJANoJFl252zRlr0n}9GO%Hc6 zqMEO4x2xJmaN;#+P_9Xor8Mnu%K`CPS<-J+V>hu3%YE+e!`T)Q`gb$adH{4*MaeNH>hQjk14(7C`8|b(^LA$>yZh61&D=c zzTPmxXZ-h$ipk|vJMJiwqeLVwllZW6a4j{z;!MhaM0K-uncJ4;`);IHZ4I1`2ku%OxS`W1FYPGVHAMhBObgNs@u!0r&J>bXB!4so&3# zppA|HOAJ4toie!dyQ?E#qt7w(fzvQ84jgtQUDG7vZnKpiag7An_?w=-fJt zE~@v?+E>hYyZ2hEPRZ4!&AsbfO0bt+8}z!DaFyClgZJ_J91(s8WWw6g(UUA?b#=AKB6GU!CkP)-hTt*~$8;ywee)6l9gcBfH-AgSRM%@qZ% zo0NUt*B*^RrRv-{&83uT8Ciw|LgbUKb5{>>)skh-XyouI^~PpzPif%PGVt{pJxw1; z2ISZwN4O2(-B7`Pc_h7X@^uBbeD9!4ZP|=hOU`XtZdB7gD1)~X1*|I#5d;SDzOKo# z-r5wMI46~u>9h3jcPMDt(d-W08O8u!bE7#So_*P@Q?t_5a+FKtoO4o||5KYA;%F{f z>KkVA5s5dEpu%;_S4SYQ7*E32RW|ZS<7SW9iC6X9Ta5oT z9v{C`(XxdK$O;K1(*dz%j3=|Nry6P+XJL_eT5fEX+urdI(=3s%;k1d5J5l`|Eq1_} zSN9P*NAiDB zZv$JkhVAw(T(0QLNim#fX&)*{e3yr69PqwtUIW_ zQO(QdYOT6}Y$$sfkPUD9e40kxjf^dTg-s9!v4G{;RRJwT9AZFep{|Uc4btHC1JDr_ z9~zKM#oi3sHGO7cU9$6-K#+WSPvtLYlu{MqO_Wysyh#mHG_)ptiRV>=nR(N#4f|$8 zWbE917Yy|NIgUQkRk0?xWz(1F*!fWs^uFjZMhK1my$-coPe7kv4A-T8Ph31A&VVm) z*VuCkG9PFM%V2Qzl`o!6W`~sdQTBYT?BH4|bbts@=rr&ufQaub1G7K2p5spU<|$z8 z)uHb9aJ{FWhJ+sHOM~U{-$nH2+t;Q<<+shqjd)t~I#HC`(CFnJvmDaUqX{9^cV#34 z);R?;F%#IA4i#m=WtAx$MV#KQavr3gzXKZw_-kv4l=pu>`ODQ{G)i!A&bz92_J*N;=#GgMl;;=b~%2LvoMseGRBX+sH*y#;Pf{M*SSAwb6rjsi4 zjqp6^yIr-NRV>p6@sWQboUN=Lb?zj2`(P16lhRHz35>sh^>1u)KQFq3HIoo%-%3d2 znzR+R1MrOYEkA=0&41y+lN9Ye>ey7HzFnd$`1VWSF2~YLLtm#-3bR*#7R4gUf-@jP zI|RfP6W44Pa!%iW&%By_ax-*82s?*3i|;&K2~0B(MZd7Vru#jHsdMKr_zC-Tz^gn& zVkd_ZkX16^E}wDPqrayG*ivNaS70vHU#a{;bet&TYe^7uYh6Y(3pJ@6yw@LbD)xDy zPt+{TmryVjzfWN1Vney*=63XScsjNTC;$HghEkZtYZI12Bed&7uWijir<6n>1~3?Hd#_UuUpGTIc zx$vdZI+hSPzA>PKWilxb-1_2^?r{%j-3Olqjrn?J%*0kSs5omxeJA72WAh6T8pO7> z=%%txHQVk7@?>gfuoT+&xzuJ%N4D(Bd3N{;>ShE{Khi??mkz1TSot~`hP3OhrI{ZGsviWBnXN{f;I=|BY*3NJ0|4Tb3--GdFwh*DORHsCU3D)dS zO|n2+moy9`6qTEHuwm7|JOQPXnv%&p%0u?mzG6EHdf0K;_l60F{luM4j*Jb{F?8gB zwnELC{uSY1=PE3T5Oz9;qY}=hM4@+0K|~pF)V>7}7IYiT{aKE$pZ3w!$(1w}`Tsam z1Z*sp^n3CK@3!{RK17sk!7m3G&CnFepR?SwVtDRDcK*p*KG#BTm*b$mqTR5?-FIZ{ zB(kLr1Zmf1c#@koi{*iSP9L6v{SGHg`4|x4FRTG}6#L6fRmYoLKs^T=OhG&U1!Xnb z?ox$7poLr_wj;k>+Tq-Dr{VK=z*c9(Z#Bi#`MmUduJB`ly@hn5shxDjjd(Hv*{O^s zPqw;zbBurgYb5CdJXUX2R7tEK&3VmRn>hmHral<_v!Q;+poz`d)M5a^h&gmuXj`Bd zxOoRyTVn(MyZN=LzwVKdn|1r#Fm%!K;6t4S%=X!8d~qP8G_k`U)ZcjBQ=0avu=){7 z*9}Ow9+H;y)WZf?;Bc zmz8^m<~zj2muWL^bq$WKR|=!8l|L&8Vb_*3FwItKvHEyie^E$~AMl`R+~y{a|A{~# zm1YUagHNP#)?c_cKi=0ge?*f*R~b}o_%v_D%q=25MbE$|)Hl|?^75QPk9!3!#);j3 zRt=?Vy3P~RXN^u1M&Pa!N1b5O(OTLZChrvFhF|xobhp(rzAtaiaER53w^7C=LM_7d z<-DC)7O3~L_pGgOfE6$lDTRLtWgEv+>>Qvc2~rUM!HRifAW07~%4S9z_V%}MaKoh< zln3Z`%UUKiK8Q3s`n0?9OpjRwG$uT;oZHFG{l!l`K;^sO43LrWIuI4gC zKWqccDsXD1g$xmDwwgY!R{y!So^5g0qxw3i-(PeUtIn98@m2prHsn8jOAUs{Xej(h zfbiRQQVc&YEn%e}o+1=LfJ`W9?OwHG<;8%g!S@z)5}UjbRcZa$B+^Tu@V<2~NB!p*+l%Tsf!DI!FeT$12{6wl zh0S^~%t|@(Q1RK*HP3~iBH6v>2wl|x#8EvDkLb6xTwsg)8=?U^5Gvch5?{F0Xvn!S ze|}wxSEAlJ3}@#F!RX1Yp^OrWhsMLhH;5xZ{x13u>nn&uo(bWbH-DDUc)C;oNn;Aa zK}j^qKBU9xyZ!|w$(5?xI-k*8L)T9l%?9CZP2C9xj24Zb@rHCeU3K%n>;7R&R3%V; z=p?lAvE@PU4lPbuf7bP9rclGlD1#`0O-p?Uk<}=`zN_#B2Hv+q)qXU^9-|KCa0P8 zy{=W-0fWkWU1*D@xmOf1PQu4pKI$IU?Da!$^|O>WvVI#~(R}@(x1>A1MD%k^L#?Ui zsI#EIe_tG6={?i@*VoLP|K{Hvsb1=yS{Fjv2+>Ilme^XE7eo-;_ zxgHG8BZZQ9<;+(z$Zuc6XG70eU-umEhyFKe@lA!X`Hg?wwz! zV0anCPUwaYJII455L zoh#~}tZHMxWsiEnUs059)A{jjY!svgVX$jeh}bM6 z+m!RMp*g6`y5lP86SPiM1eWUeW;-%hOv+t_~L z$ey(L2Hb%2lz*)?4xdMh6-J)Y2ZW?KGQ(!EBr2vCFvqZF2=L(tnN}ne$y59i9F6= z`eJLUw(JR&JI3nt+FS|TGASE4e0j!05E&g9Ub@U(a+FRKbT@|uxBUCSMk@(vH&a%E zNf;&*jAH-t_~M@2Ys90-b1+0 z=l2@cJOklpc;IJ0+GYjLrgZG0$kLD-2EC;E=@9|y(UasZX8xJqKYn2rG9j&bqm+Ez z*M9UjHOWhKw3iSFy^Hb&3)-t1)8H#8u$*&Hzg$dtM@}o31&)DQRn|JZ->MKF?aA)r zAGZ2Kzb)*21YLpCA1KDWRs#Yp%Qcc`K3u(mu1WnVb6M;3u!?Wt8qsG)1o;edO01gH z1o94+`nG*SH27hI4Zp>zVb>J}b2Xyr9ww68#9udW2bKcegioGXBn&&Szuc@eF2>HqA)ECEvTW*4r)%t%5iWfKp+H@g9FdYjTfmHy zzivS{*mIfd84pLrB^ijinv=QXWZ0*b^I^z`m$X^LC>aO5s{Nz`CH5zomxPU_R_`S( zbzRvh-t7o?b6}OCe4x^kv9mxwIX|wt6~aD{Oyhd({`XH|YeE)xyGFv!#9w8|LlQ{w zC=4+=Mmv3g88ZXpLw{r@&AV|kR78+IP72N`Wu+5P+@v7>xW|0jO59mW2{mXiEVY9A zB11xrk@Gp8Y*Hyty=ryV){!_AW_USzBod>ZsNN2O+_#Cq$`p66x`0mU+4EklXm`JDQ=)*43y< zRBUu^qg9|bqgst|;M3>YD<8b!6^Zl_PI&q%Y5X?WIiClW&(CS3^&hA-xz#F0!v=pm zbbZX470(u2dGag|c}|XRDb14EQBR_56D@lQ^3lwLDSO-lf(wvlU%NDKE)cw3iK$qr z{gCxRV*7s5z1gN3eeCW4JG3?6(cHku6UR{LpXA^bb5en2Z-hvj#Fz(WJK&qx#`lw{L^sCIMxd^B!PR0cqMGi=J1sh$^u6^_+_T6`L%w15Ufg+LU0a$H!AKu zS?--85)S5UGRwZ=oj>--QkwtD!Hu>_DkF4Vm0e!Sq&LX5GK`nGaX=v!UH; zExdnDL`e7EPyq@C+=yse3cbuwUv3`6Zyg3Swk`?U~Lj! ziiL<`+B0emC)e6b@y|p^Boz+8E`Yn)&T6pSjUegXZ{v0Vdfz+smM&F4G$zv4A(7-~ z{Jg)&`-KWqRr3?*8HKFE-O1=i5n&k}HiL|vu11stX;spuK6vvvac{nXV;?M1`pqQ$ zRbwGtC1EgFZZ9lq&$jj>$0bz!UVe8h2}bD2?Sy2x^@4Z3U>EsnJ`Uz*_GOMxAlHNr z0mCcU+i{@pP39O>+zm#K|4jyLd1*DHswQ~&tW|6wbeYCn(Hu4(`%aLk=w12c0F?tCHJ# z6ZHqiU$VaV4t9X_gPWr#%^)B=HNgtEjL$-ZCd!$ERn}k%Os=>QaC%Q0>P>q*A*J3# zpTfKF9GGaJOXrz^WH$hC%ZrWPT5JO znDpvbWiQexRSWeK!Enyue(CUF_aHgQz|wVNq|3uTOpxqkDqRZa5?G6RTe(WiB_-yu z7oyaVjd(xMcf}rYdi%)x6_Sf}t>*|U^{|QV7;EOINI1bPoAnG_&oy*bN`JA^vsIVr z*qKG?$6l|kl&M)u`=USC{mm)3GG*M7!Reo#yJY4szY)wOt<3;0>KzC6d8bd~m5 z9V~Dtmbk#{cLpycND#dG-18so`|Bu{a~{4K8~;SMM&jxV%=iy2K@{|RI;Xs;o*m+3 zjCL_ocRBv##`E*tE)}_+sFbOjj9?=f($(9!Plex0{;P)?3?zA`V^Npf@O2-)H=d;7 z)bJpe){nJ9fW7%q=E7W$_?+;gajpX#{-wCFj5tV;CHBIe=Q4dtF8O9S5;`%tnlc8u z{kImfSkw(;o+IR!k~FTbB!tZmk{rFQa~hBHFF^f7R#4g5G|hlnM>$?DoSA8Jrc1J5 z*54a{4u!PrMJPdkp!Dw|2jfvgbLETEU|oBwuX-jlKJVzjHoUVVtG@N3IxZjxL0z?cOrkiCCjheOFTGN1bUj{^>I>L1ZZS zHOj?Bb-vC7w<@+ii4|s33(0+!@;xlr=^C5Qq((bWA6Sys_&~O;No$hahsFT{srDtf zu+`#R^4o&=Tjk}QRTDFzEwuRRCw8@Znm4REzc~m&iI5CsVK|oTUS^uLfII#4 z`w)7&iAHxl^YSZImKWOjM=06MnRkl_eYNZDsn;*|)zVv|#6bsH8Z0-#3gb6+*jtce zHz;pd(t|=VPb2u*2boOut5RaZWDZ3V7a+FZi-e8a)J-A~pwVzz- zT*nH_HtBJtgEZy-ji(c38jgf>V@Kmt%W@d?r9Pe=7Fon`LrB(dtLGGh&8J^^OF8!i zAuoi`#o1_hVRj-cT1Y1Dhq9g_0JUUvD?T=#g;|l-y%9EGb(x}>9_I3-ku(E2RxaVP zbs|8+ja`$z{mTKy#Y!%dMeO#A;?&;FP1(ek^1VQX?lsE1&B-t0HAO3Jqwl_ECa}xx z7jyBr)p~ZF_7>w%j%~oMG8uh_Y2i9p9Nke-3!5;?{COoE5){Y6f-Hp+m)b7MW#MN| zt?Pj^K*nVgA>&W;*G=6^#Ny86euX3Tl?p4goS@Hx%@~#p@3-hUocqSkTnBSFl-d+; zKrQ9LSJ3`+B8`hp`hs3Wx2kvQgQ=XO*=kckyJ5P=z^%Z^-xQg*rmshka+wy0Xv&r1 zlNpLmhx-xw+Cy8t27j_YF(y{ftd#twpt*lB?98gk1oou`efvJ);7e!rqvOh&+G0&_ zn1o@JW9~1z`Pp10K{-Bu_5;PWQaNjs&gCU=W>;Em=l`zB($AHlOWgv_WnfaBO%DNH z=h*ojp2;^bM%LIM9Is{KgV2!_G0ySQSc*6{aut<1r=f=2na`6gFQa;@&+EHiG%+_bHaeL;hE>O~!gcaxgb1 zL&BEButp|(XuMDHeN|SjN-1qxno2(9UQX6B@!{f!8X4vZ8l{PGM$CJ^8KL^F;G)0< zQ1!BIbXFuk+6-wS&(KbRtU7<)omBNIktr{hwO&m|TG05TYqfYzI??61rnCF@n0;Yg zxI-r0pebW&2bcO&n%C0VxFeOf7p$&|WiVvl{o{0n{+?pXrAzErRpe!!WS?ZeVRuqZ z;aRb^XL2uf=G9;Z4xu`QwtN1l@1g3iiASgK3EFIw3udRTYR#(fW#3p~6I8q|PH;SQ zL&Yos5>$z+e5WJro01`;R_`6;p^cU^4r=jGD%PDe*%9`LH`ulvczEfERZ_sW(O5-p z4d%rkQ(hhpWxMdcdVYn+tsjE0hpm*t5|3>Zp#oc-nM*5iv(5WxFEie(bAt@2nclMv z?EveJ5x##ELE$aJkT)+!7s-DrKQ?Dlrlwy&JwNN^AAS(Lxw#+_@B;S%+7|Y@KRPx$ zXOhs}lV!7b5XY*9I*Ezys}z&ZPBGr`*oGdeH>hZ1FCJcmaeQV7&pEZ-Hm)9kSATx+ z-Z`mwwqpow$AXK%hJ)$vAC7$a2nt_5{lU#Z^T_pbAM)T+vd^k*=H#a(9{Z?EJyRri z`?4rhBkKppt;kgIzrCS1V9GB4-ev!HQ(+)n1rRUf)%kAt2k@I zv*OzNLSu%X)?fKCi_K0}`IyfVb`-Z@C-uk~BkAlFAwF^w)T|1l}~vTqNjsK_I0uQ=R@{YS}?3~2+v z!-3dx$$bC12Hf&X{@s9k|y2 zoMyE^39bvMn3tu~oOxn!B(NJg3k8SsVGqY|wZtCExSrg|20xyxSP;kb*vZd6^X5(m zM)8;fHxD$tG%)JJoknj^e9!m~uTzam7rW6LF58Bku*hMP-0#%-n)L`Opnq`qgnqD} zXbjF_UF(<3R8XFiHD0F*FcSSeyf>nS8(jv z`fBto&X}R36$9a@Cmlyv&2g*Sc+5=GANQvu^IDI+wv`x!e;t1V$cN0dOhxXP#HIEL z@KK1sM?vXmyMd1~2||XkeY^TiZK|j+05_o*UEA9UQ97-d5nHW7&3xYlb{QJ>BDvN0 zdLc07B}N#Y!*b_82WMYSP&I4aH>?*lK#J#8kovhkD!lY#D;qq}jO((7_X#;lQyL=< zxQs1Whr^a#*4|i%+TJq;QXsB?3Vk|y%Y+-hul0PE1AiU|D}l?ct!I~ z{MGz#fw`T2nAtz?^9@XxfzjQueL#%x1bKp0uVG+E1B39xzSVnTObqatiyLS{31()K z0OFuwT6rkd@LyzWE zGYI=)9hQi}P(^b?kuwMObIEk#ph`CWZ^1+ydX}O(FTOAH1bGd- z`%ZIQx)IbYduk{kAu>Yw&b_}72(NTG;rxP8Yf(|6B-c!O%YZawfujIdn&UX46V<(R9%JO{_u|BzCiW ze;+f0LAx%Xi^(ES@ny&~drJBEL7*XB8 zk0%EvmrAfMgPF9XHTPEHS(GgY!`Ae3tiEex!^Sc1S*&$@A!hYfXfUy-u@Dz)r! zbFGwgy89`Cw(ZZAj>v#JHVK15a*d@ps&TzH(Q!9>MoC6XXCxy;K=9Em7pvKcm#msk zz!RlltlzhqPyQ?-`MO;p>f>+XpxH6n*<_aV#EFBhAQCy#x9=+jR4;aR9anALz`?;} z3U%Y4ZtGNJr5+8j?d0M@@(`;T^NHhUYkXXGn>Y<^xkK=Ic2sy@ydcc{UL^{u{IlY! zpacA*s8*fyxA6;$&PRKBmgo^Lw;qX)?hCN{FMlg{kJ~-cbvwe|qRXIHr2(gITEM>K$mUAECh|lUO@>2znx8GY~z)2t)oCMxY8H=O~ zo~)E4KlJ8YjIZ(|q)1rDh+hk;WC3qx&qMnREdphS5^rNJnQvaC(|C`;M6>T@J6+!% zaM9R&efK=+Vh=F+Aq>qMtXPOHb!3rKxh}YiW5$5`zkNCP9UW;`3#!zeC?YNs`Mvir$`6!C|B+qX zIPq!5KWtW5erKxt(y3zi>}q{Dg2bVF|GpT~ks;!ve}pRgh?8gSoXuA;7<`*EPhh_z5)=*^F76!Idlw|XP&HNuqiBCgto^BJ&x+ghK= zi^N>AfL`eC88+H0g*VmK#^*HE zZ-#6-^Ylp{#qn;YRoEz~QR#7uqrRP*FyfMmy)&_eK|yUBzw^GGw?F;PcL&%MnT^ZuSa zTHeS0a~tYroybL{v!DF~_4Vu=gC09GA(i@oP@Q%bij(;Lm^@5?+5_P*mfDN!%`+gZ z0h6mF1FS^=|CRe?pFKyhe(P7iXBW&z1GK`DaTjA-1Dm*wFU+HU_vyT^qZ9L`gL26x zKZ`TR29Y8%=5|%iQ&nq|JjweO~lY?q; z*ck207x@QdR${YmG2>d01i~QS{?Or;?1)6hAOh%HrIf42_K%@o1UjKHAJ;hVUK+K7 z2!1iA<&RTZ&(2yVi|ZKt_IVII2tdLNB=FXU*b0F z2=y7-M^BH}wu6>V*GD`aXD6n;va5`DRgBE1!(1xxAN%lmJuAi=d{0fanDGC)I`eoa z*Y}UxI!34nF?J)$aukLvA-k+6CuPaLk0FF)Un5y3V(gSSwk8?bg)kUo-$jZW{4@Uz&&>1O_viXt@B4Gj>M!@-8SbfdQ7MkVEe-YVm~K-1uQj@e^!+2`#^x|d6dj6B!dH0n<21f zI}hI)(yG)S(hS=N+VC)bj*hmt$4m`VH9oz;wCVz)<^k)2vX*PrZIPQTkYu58(ri2G z@4XR%H-|hIx-jw(%QB{0Bc%7tlJ<&`7 z;#4d~tZ>0;OmeCJUz0~g7WJ)Px$;b^`{1j6h<|EHf=YdQLM=-ePrI)4LeyAwSng?} zr21`1{26r0in`AwkfA+n9ob@L_$5N55c9|wl5-C{yvu&y_U5b4M5LIf)(dCjTl7-T z8~&gu%IoI{Dk@No;=GARPE&mmk>V zxqvA_it7n`CFX{tss3jpTk4(us?Dz5@7-#L#7F(FRh&~VNTFE2^%dyRa_u&=-mU_J z#~YoZR|Am8J!h#dfn>+|zLDkKJ;@E#^{jl6g>z36QvU|I=v8CmIU=0vGCWqI|HnOE zzqRHa$n9R|BJR#=}gT9o}a%R@aw>s~5BSF?e1Lxm8B~p2YZaeIe(bLeaKZAXOoy{ITBdaKa?lQ$q>%2RbW6eb)7x?3% zQMh&i~5)ox7_8Zhx?rEyUc*uVcze5tnBUBo>j z5C;wI-R{@|F*nEYPAS$IPc=8I?7M!AYpQ9SEXdh@npc z2DP=K%jxxl-7L_9sR13_NKh1s|Dr)WanWeYVh@lvb}w~G91@j%J6%8$CC>xLipW8e zpRN;W;Q-S&^{r+G+PcRzgUlw+9T0*|UH)Ug$EhtcN_=@lv4FM1(a_g( zYgvb9e&J|rAbKgG@8{K=lp?WRvX<)+mDGIIw~mUIynBN)=A~N=PQzu$p5&XeUT@}r zW!c_<=oYB+t64_pC@M-T`)T}^q>y0jnf&?^QE+mq?a-|%)5h5?u0@ zMte9Ct{TAzpZqN_U>Rl$fVs#zDmnep(>eRGWW%OS^inVV>d{wiYxOW-yZ-NG7fcx9 z<-jq5(lrMVP_6^L;x4zYDns5z)}^`$Z=iP5--)OnXEi;BcFNH^IO^t*v?P)ImF452 z^5`ppO|9!Pp$|ZLS48>1!?TD`Fq>Y1@1TZ4dqBe-wF_f^+BeeP6uT~d31N=k+!y5f zQvGCDx>auh?xk)?5%h{mn+cAO1Kz@IuPy?Hqsn+ZlhVp~1{Ts%@$$fCHkvUDd$eID z8*Ca%&;K25$1+U%YfwVsUL+u!xmdYx{1C{kMNf3&?lrNq9gcel3xgfQ#vwn&I4LTA zGBRb&If5eRAA7@VYXpa4xOQyo@^N7$d}LW*nY$n~I=<%B)$=&6#3g!>j`l>{OOYD% z>zX$2X&0ci4QG6#Nt>m{Nwf`+-jzjQ44#cuJ6@NZjHdY2F;pYAyznQ^`*+}Ic(lV? zb-6Ib&Agk$&pSYPFGD754BJ7bNaH)It8iuXpJioZFz_D(d&rXskWA(jb>#FKK4=Bu z5nA=;U46R)iBECv;pq9UPVu!L%`iM}$sflXqr3v4>?MQl_T2I?A+Lf4|9$|Jm;DDQ z0Uz-n&)U65_D(MDr}-WnJIdoQIK$I^yc;(%S!&pv99eYtar+JgE_c{}pE|s(ZU=%h zp%{T`Em9wz_qYyx)1W$6vo9X!zCo{~y*yY+JzId^Aj|lsUL${W1$Vv8)HV6Sz;Eoa z2ZL>O{u%G&B&@_JTe%u}r;pMh)ufx&iTgYFGhD~;MGoeGFLl@Zp(8z-m&^7`(#>PM z0c`ZU9DgnT@0iWR%@I{&k05N#$qK4Nf!5IAjs2;u=VaozZ;zmK!cZ9Bfrj0E)v%;t zt{_Kq?!$4?;)A1w=AX#M`~0Ns#gKoCTO?gX)1-FlP>w$nE6>j+-Yq|7apW?KyT|5N zJ_uqy*v@nMVc3jRo9e7Y#}JgY&VKDj=Pfn<@fo|n2}bWSmupA0A#7)Xy;GQB ziaJn%^uy!{LZqV4?|b+;&Q_W*cBBq}BO={xZb)~A@n^sLiMeTuEOvIr+!PppD*X8Q ziD!0lQa^snblt&{?!oowFbYgUKBf}(5AQz_mdP{KKccau&3jc&%r6P`%J6==dC*}v zuY6~+g_?PAKPb*1E#Dx-c6$}H(l3GqO}Y|GgldKhYe#<3Ue}X$qSZeP_`r|+X+EQ^ zU{t)Y^9lD!_e5o*{vd9OEKXcYkmMA~7gQn>%XCs#IMP(ie~8c;x`sLm5f}?9^Mx03 z*E$$~gCy6QCbsVOTh7sUv{G7)k!nW8`sRy6%6&;{H2D^Mt2}RFYbQ8&@+ThQywCqs zB>FTsETa2(Sq{L);p#?T0H@8Z$X(u-;5xND{Az_9)BRs}?!X#CxLnp1+WmJ*Xo8m3 z0ZK=4QTXr+iD-u0HE|9NYs!7}xnG!9Q;SbSyido&ph`KLfLK@62VGPMkR{LFK0|ki zYlO37C_KqH3Y%4KY{_rsDXmPUP-v{UVxspYX!66Q%@%W+m;HYRyGeE_bXgD1hh#ql0cjIIg z>>KNhR;l;u_uK6t{8bf@koe{ZkSP?H#7n&7XGkR^TFkkTexy!&7go%3HSiFG*QtPt zl=#pYM{+pXNSX9FjY~SdR*!^Z&uEMT9z(C_QOrOQVO2uF2{N(FI96(C^HId9E-TR4L9oNEc8#vTB>!vJCbHlT1JDWod z96+P8WVooh|J#cS&G*X07s25th>5gKD|V{c#+57oY7^Ucj;kq3WW3R<@Tw$#xXbgt)MN_MA{OnjDPcf-1|pb01V4Pwl18xrmc zT#y#W+TC$%kbL8Jy+-x&=rW{4A6pEt0c6J2G`-|inoUf}i z-r`2YI=0^Dm&YclgW+Z1cE=W>2NAbyj~LoF?aUk(<9;h>on0qP)F_nj9%d6BJ$d7g2Gax?T9ciQ@$P& z(`G{-QQ0kNndj`ULSD8H@hPpi*?PlbMzZMbN2YquDxrlQg8_+COpCLYV1(i3uU1y6x0SeS}Nf@nqn?|@Wkr70UTy9j%CYu%-SH+(od|5^= zYwU&CeKeZcO(nlU@-SJ^j!XIq+4vm z3IMVxq}s4#y$Fwbo9IQV#bq9Wn%n6Cq`ls-7CMIl- z_xqMm8XS5vmCTAkwq}26Ap4k6fg9@jzB{0wC!a!^9CP=b@B<}@4p8x=Flsaz0-zB^ zPxbuv35@^i{HDjRJ};h6!_nff>ad8oAbgQxA!$Wl;+g$ zMXOz>2twZn@bZ^x1!fwoVhBy`s_o`PYW2z9-L{w+vEc&9Q~GH`S66 zIRE#?wHzi2QO$m)kQ^Jjz3S8Bs@(Zwc9;n*rdH{J@AJqa5Hg}KWTw1zxA^;{=j-pd z%RMHnr>}Up-Yo9Pnu))-lkD0aZnNeY=iV zU1qr{O_N^R1G!$C7yKqjBMpX)=}(yunW%rmz81OaViH}fBF^9JC2MjUmq2-r3XX^@ z3M5rA$1_*_I@&9#EW|oOJ$w3LC*Dl5S@gUr8i37bK_=4vRi2#{e+@~z=%?}f3uyIS z4%plyR&!Fjd(aOxab$?YK;6LCk(0tC;1)I6rr;a`3qpi*cbw-7d^U~f|B?BoiRm|R za#|fYleVO_+s>R7_$otQM%2EKPFMbHLbxoJLi_le(Gk^XbPTL>WKgQKhKjkP?>$)V zq}ZJc$N2m^i8@>4oK+P~KIxL4xHFtSJ6@X8PNQ z7A%?7_bT%%@?K6BVD=YnzH&rN8?+eXCMG!<%m>XWR5xxjd9~K`j8=`B-H`%|C5Li@ zwJqIFNaC44)J!kZoQpV!LE~l92*+Z&?t57n(sb|ZU=h_fn7OR!Yum!fb|3e?kmnK% zZusa^nYvt8!kt<13?y)Nn%Tr<>+C*%4bWO;K>>bB@LTgbLd&+~j@D%8pthC~Owp%7 z+DRt)d(Lpo6lP`#<1)x4tIpc=)&j0HcndbUh%dd3trqZ>PTuJ7C%Fy=mBeh)b+0PB zid_%W^%5AsIKk@_qRC3kRn#>OGP$?z>ye*Axv|NxVnWLSYY4{$&33GL2zsT?%NcK( z;9Jcu!~7@RsW|BO9aVo2d-8@X1+fb9#T^fZwK$Q^C*hZtFmAHbB~u%j#m+PDEgl?= zcsk_ITqjj7^&s5N&(Qi0Z3#S$LS8F3YIGsmwjMt^3_%?B=Ty{v4JQsnz4APDtjY?i z28vmvsMA)^V)6=^NrbUSaj1hsRnz8k~L5wcQdX8)qdedynj$l|7TW=DMzH&V4j-Gza3kt#70c zVqgG)41ove2nYJ02X%h}0-2bA+kpd-72Ra4ivqIp&KmY%+Wq0wz~3rj0&8(TZ)$1YD?-P}F=Uib$DLW5q0 zy$OFC5g8Sol>9y=HSNR4^qkx;dHDr}Ma7j>)it$s_1_xW+B-VCFh9C`hK5H*$Hpfn zvAB8spM}Mx<&{<9_Rj7eX`g&>_}4DrSNR{afY<*|ySM?n7@3)wm|6eY#lRT&*KlrT zma{6ycy8Qdedx=3?$Ybyd^Z!bD_T!TUcOJ@fAoCtq=1wvPMY}Fw0|u7zcVcK|CeR| z8uq_-O@VX)SNq!-fhQ9qBk)%-0gd?>^WS!i_1Hg+^*`;{QVrK%}nwyCm1P0ORayg0d$J`hqG(j)+@&&r|P2i`< z33Rp_^KNu7ba>9}6`a4b2py(Y(kEcN)`wX4{Y8qu_vxI?uij6O;|({`$t9$!BT!f( zT%O)->F}VxRF}@7g>wGEfZ5@djeR0uuDQP+VTfWl! zuf<9Y3yap)t64wk<;Gv<<2E2&r3ixm5vW7Jq>plSvS?T55{#a%B2o7I9inww=W<5J zFw_z3st}MR)DeroUAm_4i=l>I9u_N^a^TIpHpRZX16lM9d7iK%fvuHTxjv&QI!??4 z24sB*C1HR6nyqS1mOGt5wR(DyRkl2yb|h=^+%6RR&a1R30pZmaWA}V6!2E zb;2~__DZxxPLGb?U%m3U@mY5CX+vX>CRX!lz%S3kB#!NLq91z+|LgRE=&%@%Ra%Sw zK#6ckM+o0Nr=W*=e0LfC?fg4t-)FJybEpl$d4vF278gEK=#OA>xZzc1>v1TPmt*F= ze*=5$L?ISiUGbx$N~*#q?)_sGW)aY#UGxE(^qN-E!c6ZaVQ(FQg!jH(P-@cEO2wLH zIJ*QEuLf@4%zB)I?P^nqDWX=l8$OB(Wz^uVx3(r2<7e9-SdH(!MA^CWZ>Ji+oQ_&w zfR4_m&a6(>ny_?SDmG>wyx5U=JL+7ZQxJ$`Cztc;qyi4igJ?BEKk8?t-fkva&Vxmz zMBb(p$;1`Tl!dRSE5_$Zo>R=5DSfA&>E+uzs`at_{`HRYr|Zje1%Nr@bP_`6iW^k|Bj)-H?dAR*)fOX2kN<6A; z#Yu8-1E!0k)T2t3vyC%6GlB$cp1i&HO-#(lO8>sFp%m8*=IGy_8KdR??Z=I-{BrAF zcTaeKe>y^qn&}B#gWs<_&!TA?sLz~Y-)_LLF;Aazq!oX1@1t*bYcGk4y5^m~(xznY zuC&#teW6aS%Xty^CVxq)d}D%MsjD&vf6*z0a?(Qxi^~E2RH+q$OW<{P63ig^L~Z&< ziC+U;@Xhz~+w91LHxl<86R3J(#E9-41TK98F%2B*2If%yQ+Q3DTk*X&#YS@dmc`=( zyXOh{o|>i?vt{t;v_x7Me&o;nP(7zzF(R1IiM(fwM+t1++B2n3c`a`)a$`46|%+le_CJgtcEh>=!+CP$2J6lPcEj z*PJpA*^$vwOuwB`TBu_gj5g_Ceky=U}0s;K|mrwM|h^ zViN=f4{TiNuNhkoXkp;FW)(PAqQ)Qs;WF)vGr@#yHR<^v{5aOuP6Ug@9B2Qfg|51Z zxxzLjyfn}+PpbU<tQIQ`~_QP@19Sf5Es3m@pk2C zaYpjsoS4EqvMuh(c!#*?xYyX4*~rt!yf=J#5j>Sk?`|;rTd^22TjvZg=v}wLYhGDw z5g}+$`Ep4L>=Xn_AdT5k943>#;E;9TP;~rQd#Zh4I z;Nx#5c_PX{J_%y|X z`W)^2iK_K6uhzY8x{*hQkaEd;^Bf)!nO~TBi8yc9WajmFMZ3MpuMDyv*Vm9_WNkoV zA?GVD=y(TrkR39*CdLH5aZ>2#a}UA9s|I07?9V>nTQv`Y{-jeN$Is!@=Zn#cm6Rip z?j(@|n{`6@Pr@IVu4>{`+c?ixi|fX`-$CrjPLgwvQ`fr6N*&(XsYz4|@wtjChy+L*|rh7S#OBoBP$EJdCG$c#Zr;ZQy+F)bbO z9IU$RM69k-k5!G)S5VJ8hbHT~ctX19zx6M~Vh zzqrEnKpc&(C#|H!`5u87_wxw4dFZ!d3uojVosx{cgUn=%sPY#^mGq{;E%VDEW|R{i zCbCq%j5+&T2)}&jPP_acU&>nFO=e6B?=WQ^JvuAIB(_G)Ld<2UZyQq{!#Ky47VAs8nDLE<{v1b$6}Dp*H9@w{Mx9HV&zs+TD_0{!qmJOULXe&_%ferNW;`oIg(#v-)j zx^e^>N6D)9>E!1YbRybj7vU8+1mlDHROko1e2;j4>hId|bjs`ieVDEU1dA>Sf+^zd zjBo^m;J*~va>k+i`8%D=&+c9aI)>QtW+Z+ zxD$QK{Rrer}^oG1--(jrB}wCxY?d zaH>X$1N64nrE#kmpU}_fhUgx=WqnOWC~3t@ac;*Y#gV%^Mnn8AxaM%5%DY1{Bnr02 zvEsaz!oP%>5NhbTMa$y?A}O@^f&xd6VPkDPc2#|!dscU;TkEY}JkBCi@U*36;0yGG z56N^0&J2P7Or2pJULcqjJlh{n)zESInWj-`D5H_=@N>&gdN#e}+_Y-GzrU3=3*Rfw zUe4D**f|1upkG01j%ChCoPcSNKhyb%f)Q%X`fF{o+)bXiGC%8d%5B3xbs2i;W;_rG z_l2a8^i|DIwp~MQlsSg|E&30MTAvL+yw4OGj|hc16AUl|%l+kKyNRh;&ld6e8p3-) z5zB9hGg4z;4y1p3e|*Xl@W9lfC%}4@^E+~4oa9f-sM#>+l3-IcPDFWpHE<6s^*GEV zeJj>3HEFo!Q@a`{XszHK8+d_DF3RkzpA!bx|2B(}P!L1aBeJ%?d@Nu$9^+D*Vbgg1 zR@!shwpv2`b0VN&W!VX#4pzQhkP-gqL6g8r$^y z2-Ir^hEu7nw6$e4S?bw)G86tO{OS~~R!7n?&S%5s($-b=l!b!>v5lru?;qPlY)d9F zg39Yhk|zQH*K;YRq4Oc|m*Q!7^ z$#%%o;JK_LP^HEO2R`Nq6p7$#QcyiCukz$bMadwGS4AvdHrLd<#EacNeL4e*!s|WC zRvNNA8Gv`|1O)|Bxe3Ncpi|$;0)jHX$p&WXA6H1)A8V)R{mE6Xlq;yAhMqxHMc#Hh zMdaHiP1j?{0(XLem#3$fNGCJbt!$nvgV!x?-Q~Bu4q|7B9{m?Y{RLDfhef;EYZa)g zEdnt92uHGi?i%Wp-QkdC_uh{F)*58FoSNVaX|Ihc0_M$0urr6`yku9U$wOR>uvd;; z|B|CIv!@T&$rCb0wpMWT|DbJ0imEn&-;V%td!FY?IU5CG$>|CnT zPJ@_#+=Uf3U-~(5Z`CB2UPZVGe%k>4jyg7-lDa+I5nFek1!Ux}lTz4s8FUfk9?0#K zOT@G_1jepTAWBh2S%mA6{1nKK{x4>bu9tt2F2 z7|EYI0rscG0rSvNabg&D7H15OMDplJ!t}QBN&SWibZ$rO_ASZvVwdJ?zTwkaq3W%@ zLdRa79OiyqzGH^E4Q78dUX}#;s9N8a-m0aczClM65fJMO1Oe>YoO0OKkhqd+u;z8l zKHiF6vsXaniD*T-0*ik=nYT;5*xKxC^O)APr#&Z2j8DY`X%8-K!4}LVV>^7aYG$bC zTLxN$;q7z*a$-rBIEQR&ddSY1t`*^@_y@0&jyb%`m0~o#q_CPmmW8xZkl~POgY}&} z&s)7E`9&org+G5jc2`EJYFh=TUAR8|G3LRX?<*FqMh3f74epr2n&c#ipBFV=Yc2X= zWN#{m@vXFlY+3+-WEvBG{&mRXF^`E3f*lEwEqr8?IS#-1i1^LEtO#5$Tb%LT`s3RN z>SK{ZJrqf+Ky0)rzp>~8Mrip?JZ%kdnyEw1T6hOoV8cMi^^IB$bv>W|QrZQFM#%ye zJuB47`;*nkKBOR$tQl0bLpqJX%XWX9X_%z!OWCs-Kp)` zxxbQcZx4Om#RisHb(f4kz8K=MWE>l`rNo_>X`vu|k|PftZAQL6OhFbLft*^>Xt7!X zspwjXw&1dQZgTGp@2s>z3A=-rF#V-_<_mV!2$Lu9;|oGhTWT*|iK-LQl+Is)RM7eh zHc}|JzW$)S*O5=0qRTiPKLRb5qpt^`J{Kw;fksdRmg$GN2mYIIS46Ak>E2PWn;OM|1R4sKq^q&6Z;*i1p5+QHk7aRogDIQL+}bGqaG08DXb?|_V08xLr}hI;Raz

B}?THMN5YR#pvf6Pz9S}aK-<^RM(|$$-bWD>y&da__ zJuNNJ%S+1lw*=SkZW-i-+ph5@`u;9SCx$eoekT^Jq(veg9)TteIwe@Nv?N=U<~7t&?OaV&L_YL2dY&2@uLEcNPG*R)e5+KUhT?V+8eJX0lH_nN-i1fqou z-nanMR0Jj@&MmOWa6-UzV-Qt=)c+mK-gJV{GizB(IF@$o zefn5YU>W2BC_#r#qjoPnP2K^3<7MX}R0<4m6Xq#eH1z~g(v+NuL$D||Z^jiGmAD5e z7tWgGtkmsXu5y@Mx#Xap>8Z0K8ug$yZlpq94+QswT?1B4Yx>zuawW|2Pmly526mr5 z&g;drk+W}2vCcpqxG!zZB0qqe;qz1hM|VJZx&w{*NgtqheFYBFtA1^p!?!{ve|at{ z3n5xHz>iQ8RQ+ah#(V-xOh%veeG}LFUrlZviV5Bwc_i2CAM^$DtigYTk|@-xR@Hcs zB}#2SfCU8b$$vWHXz(x7d1+C@{<~beC_9(B4TJAMOGw$|Zv^iV!y)MK{8sdf?vCd{ zPcI1CHr6+LeJ{dH<|bb0mbrdyJQdMuD;YoheSwSR9-?3_0OWYq+k+pJWhLf`>^PI2~u2{s{?s#}bge?jfnbPr79 zXbQV#9G{C8liD%p%LjnFT)VPbT8v?nrB)aZN+ta8gso<&XDuORjo;j4B^(Tk{q&nX zQSDkFQ|RPmKS90eS&}~~YG>=}bG=dk?37YA7o!m4bwXKboF;RC-Mo7W?jC=ADSs53 zHG$G1CiNAPNn$)=9kcdjjogWmb%M#CweQ)47@}sgPmjERDNU9D)oF9KTJ3Ofd02lk zD{OQ$?z;H6Yw!rf;FwD&7iSCepBu!j1u!}q-)*P?`Mq)f{o|SD_XKcU@EL~dCoi(z zCLqI2lxm6`XDVc7P0|;oL?4)Fl|p9EUf*52XXJR>C)aG}t8^~-Ci;=wd>9E!7Uq-Z zww?#SVJwaX9YeXlQlME$E}jV~p%#H3jZIRXQBO+DkW+BnN1zl>a~g0Ih|>3UQpjde z$gL4{`yuj13WDV@X;6>JPZ z$kREFKsX-ACLKw8N%mGP>JL)*|Jfa0H>9=O20&%=H&HZ8%alA~C ze_wWI?zP1Wz6Y*>zQ%eWA-*?EBlicdAHS~c_I)?+f9UsL3jW`IHh?ziJ!ehe#LY&k zKmiesrG-bob0ev1*m0Dx{;uL~)&M{0zC5=%uI7pv zzCH@>LZ_i7mhcS8tpPKaS~8DiWF zLH~C?k5$Mr$DPe|%0Pnkms$*1*T#-yydx+t;u zcj#o+c?RvT2Dho3z#3T#K-m_q;p$#aLGs#N;`z}%xol(4LdA1!XW56~i*l{U1Hv$; zBb*`C;1ds-NS)*b4U<_ODu4zGf-4|N5(LZB->7~>Cfjk|pXxfoEyJpwe%l8=eo}+| zFM1O`;9=*xnVyTf3Zn*^4)x z=><%`bp}*6HJPTYsrH7hT;f$v|M3F%r?|lB34Oi~# zags-%hFsZQn&*!u!GdX9Xp1_!A6qco*YW$oL6uQ`Li?4&lk~7$JF8Df=bbOV!tN6` z(XRR3fxh#`h4^r`Mxh+0lb6Tc=>CgHR;_C&C$crgj>Bml)Bi1VM%plAks~zd{Rd1)Laqat$c%i;zeHg} zp#3)eS9AUF=9c2-?pQ^I?PO@tzFHtyw+}&MraIw+D zz>v>MA9H6t1vRe5NcVlvuMX~wn(-`&?dA#E?v_fv@_YT~pX42xvV8aixMsEs97d#C zk&R{BW;v;sh!Rw3>;jEm^RP%MJJ_pWdU6Wu9_4(YvPs{wUvsI^SpUSSL2lO`+8{X_ z$a1h;q7`VFP_%x-lm~t2XFCa=-)B_DHg~tUVKQLJpGgPYubsfxn73dwvba& zdYB}A+}!(sJ8c9#myVQxpM?*BtC+z&F#qDt`AGe5W$uR>QQ0*D=jRngy0&dD)Tl7x zBq?_cPi(hg&eHyPR5TP4P%%#2w4x?6A`)k|QC3Vn1NcYuR%?Kh!hlZ`lR!m5AL21K zZlS(n>;)yJNO4^pxA8Q(TQvD5vHj(ba$6~D(iT~Ezf~d%a(ve;%1Mo!@_7oXWEQ2P z+DJ4>cvD&~*${EFHht*VNZYpV2ae94V?(nJk!-hQzdb8mQks*`-q{cT?AB3FY7s@W42kOjc8enwwy# zhaCtt{==nkKi^d4Gzya9`{&7&J{^%gkKn=+6)|y5nxAK(p{tv&jg*_0OP$vj_-pIi z;!RKKWR_s9dNJKze>nMPmS(P)t<{0adqMAswMFyO*E;&@FC3QGwu`VzxXCJ=m%MNzVLuXu^|T0hV@0BRL(8#l(7Mods(ld|V>_~; zW&9$(v_g3F?nh=U3z$b)Rx$z@VD)vxA$plS+VdlYTlT`{4rBGhjR`Q(4#-qYxx z_BhQn_2yu6JOXE$hF1VHt~__2IonYE){uXUFaBJieU`@gey#C!Rw`E;gvTDy`lRKU zjpxqSKVDBnFu4`$?V z{)!tVi@Z$$V^fu@>A$kB)UBa!j#%%6-rnz6BU{J72cHE2xdAPRz~x4RW{FRnHjfZt zV1olODC<)_jYWVD0L1G*F>G)4mCoslR4C!rL|0$1iBHe`RX~r&iGTqS5}av}1Mi>V zCf5+nTWhl})Q>BktsULSWz`*5(>`^_SoL9@QU0|3HSX%CPaX|^N>NDdpFe<^cj;z$ zqJGjFei6$SEZ_UX5_+>~52p~tkk6O|YKwOcy6kQEwb$=O0 zJ+?Y!nr!pUCa=kjNjA%Cj;C8xZ`A4f>-eG3zE&XkjnXo0nGkITM2p9q2PNVfN1%m; z7e}DFgjR4ALYU6$7C3wCQFSBH^%d1Rn~X93IwGS!hK)^E8|szLFy*Hv`?6+;K2LFp z*SC6L!40X)C2wp(m)f@k;6uo&f{nPj*~vP*V239CJlSs^60tXJin%-bVoL0q#ju(f zqwzVStI1pJjJ3<*Gnjxb!Ub@*1=(mQJA;aYX+I^Hd_P6Ah3As-IBefejfeDEGXIx83#932yWZwsEd}Y%87dGvE_153?p7$p4;z z54UjA@-Q1|#W3FFdA^xzDuMX=7!v`d) z!SQ=JEVT*vt{szekVs$Cu#Lm{$i*yDQLJ%%P>+cW8-G8`N*WT z_Zo#Vu~gqLm07-UDVSje@@Tbo7`jASA}14~0kY;gLRTU;9D(jl)-kplk12VJ*=t^V z{GnACd-<(4tJ;ydRo(N@&%`I58tr)@gXSU~194RqoDZr!P)&>qq z=t>Y07!kp=D}*X9nahE9+LSZtpXZcy>bLH`5}Y>^+nuGAn7j4S%I3>-?O0y~yt2I+ z(j{Vkupwbt>gZ`%JJ{qeQ}E(?yqMz0X47$5rnT4P@m*1K^tlak)lSL;cKF69%sHfl z3s$D1Ru!4tMe0a;ryX>IYzbk5*snO9svjx8yfD>|kaVULX%*~Als$#&nq=G_cxe2U zA~Gz{w+meIsG%kNIO7; z?gPe2NGu(JMj=Cb0}@Z(Au2ryo)m6$BZRxl;vm;}=}@8=>JmjV?>N||7sST_Wl9AH zB-F~ufk#6stD$JtBT&U@N)4sggLN`jVmj9-8NCcb$%j!pxjD*rwfuTt(3E* zMs<_;>77-dZeeK25PTqR@#;4M>BC1A8$0RWxpg1(e{==!`H@f3eA)^rHUH}Xe&);zbr-m1e zM|B*!KnssR&wU61>TO3L5%H28@Q^paK^JmY|AO-4AKo!+Io%@;fv?StN8vCtLz|c- zw((N1$Bs{|ZxHe%Raauanv!_XHTqT>z6y_BZbOLw9R!Ic8(;uYY|lEFkPGDd&D$?= z@B7WiX&2rgIK+lhz zXO68Kds8f0%UoHZP}ZypF|GRf4w)$$|1w2w}I~C)w^KKB{we#$BZi zj`0rX^VVtzndiK>!_#zVy=C!9bV2Iy@|T*T-o9(B<-FotM?WZUv?nNszQ{a>zaOj{1D*U@ZuE#1#2EbP6YY4TGA0{nvTFjDyul(QNkc8 zkcVdf5XBbP$U=&l>jGkcV#?VWENr<`w+g^I5|D2V1bjCR-fe%L%V52RoInhqkX{WU z00xB9CEK^RNbP8C)ggR5$3Mjcj~B={UCewa*fSv@TP=$FqL1?83Yq0 zTgD;=Zv%OOgSZmh`n3Y#77|Ublicrf)wt#ebXo#JH76I%fkPo2fNB{F=J&s{3*mr1 zYQvk4Y!!6J&kT*SCuyE!6P1&a5B>R!kH;4xNtFkFpVv_1dDady4yqb)za1Eedifg> zlclBXH&46DW*qX!9%{ER7diF&cXKmla>z)pX4T%b`*ezKMMWO~DjZbHtgV1iykU-x z;#hl3)gXLU)oS^1V%gLD_<6jC+y(Mm4QYMj2d|wxDb7@DvLi8g9=)E6n+=Db8;0qI zl)((cI$pk9NxL^KT%92%Xz9LX-una)UGJSc`<^Q!Pp?~v;qWq^b3g*ihQ|OGIGHv# zh9-{BXLfvD&JwSs-0Jqu^6Kk|Z9Q987Zaw~*-`qgl}j_p+Cnsa7sBpoO_wB$Wuj)? z>+Dqb$G-kHdKkQrQDTt7^x?_2lTFi)?g_d>e8y|+`bD$(x{W9wf<0-NR^7q|VjX>WHUZH@B&8x)?s?G8#d`V`_-vWtv+*d;Q&@7 zmL(&_w7UsJ%D?!EF?yaFoa0HqK+vheC9th&h>MQYD4dS5H?bSKl3)Kd%WYD!K3;Qh zJ}^8o(aPen^s}D!b7+Ppz~@-rk?SZI$m2G&WVrgWh(vAD+8h1Y=&=wQ_%$9Ot6h!d#`!3Pndr_YRis8yjk8woX=D``|A= zdbl5c-*_!*?;kxzLWBHphjtPyq^6le@2nYqEa9HMixD>8cT7w%l9U3OIarsuKbvq* zcMuRrV&G`F;Q0eAAM4+54$&lI+UR0tk#?LB)yu5!hm!qvl3(>qU0*`=m;_5u6%ed* z+0~kBnd@MZ2Q91gK~bCHE_$Y482e-X0>8gWb&P&T+#^Pz5sP{NKiSclSp=PKgm<3WzE!=ESg@zI@aqe*SI{T%`b2 z1u)mx|1j4Az5ivd9NO2*ndW2RKJCrVFaLV#%}Uqx*y^9*?}FVlyZpVO?QqiM2=rOa z)z;}#TC!SuqeaY)!JqzrgPuQ}X;=g-x{RvJqg8AjfwTo`;B82t%oUsx#2T%ZfX%%{ z{+NRy2)63*5|uJrre{JB9s4gud!Nf_#_IXxta?0R`+AUnx1)+l2+Kx2CCEZquA7!9 z#rs{uxr*%~5-C1H9nQ-N9bh*wmK8U?-b<+5pWo^fF1I$TNyiScd7chZ+CrgF`NqM< z#`e(GB+rfD8aVyaYT5_ifDDyzQ#Rfzu{ytymXhXd*s_gx_ zMv1vn*GKlvMw0uzOWnL?CmTWex0t_fL&)+@dCoYqw-W4dX=waCSKQV~m^s0_sBESu zU~1APW-qGT)wxZVFR5$D@Xv(@K2Jd=@HgCA8p-3kiE4Pj|x2ENvk&T@%AsdnuL$Wf-h2VsoM+w@^ zdW@$pRycoN#_GIJ>@8ePPt87KShkU5doK2?^J9WZkuVqWky%p3>~ zH^}F`a2rp0R;Xv^d_>pm@SV@%)8Uyfv9~T7KVO}4Hu__rvULXhe&IP+;*L*80t=Os zRtxyu7i$bYAz@Ej_7)j2=OqgvKgkj$N_V{$qoE`9evH&K22Y)sOI`BNWvSa-RM9&j zJyM1_BhkWVy;MNum7GV$w_Kq2a(m9GjxQ3!YmxKH+<GWe!E(T>ei?+4hMOXCuybn7vI;d($EUZQ1&szJYm1vRg>0Y5r4N z*(z@w#xi8YSMh_Z^AeYEZoo-m%lKVR+Z;6C1@K%Um<#OMA8$a?rj2pfNzEYrm-*( zj*}k}nDhx+NkPLaRqpP-Dgl1d-3?Jo$-!|LSZQAp!(z*9+f`HHRQamuzPy#-s&V2e z<1|)&Y9jSMShY#F1o3b=?P}3MtV&Wzs4PJ=-wjR};x^%p=3t3*XKBx}1>2?!_v0BQb5QHXEeeQF^$7rdXkK&U`a>97gLJ|8YmMDHiyY1 zF7pyZ>rMD6&u_=B0V*wL<(s9%I;+}C!R%EG@n}Uv>2r3yKE!*(I?EUJ7|KH!%h!}T zBZ&^5GkQ@-$vok_kkN7(3i9?=0wXe z3NpXtKd2z7tQ~R^Q9~8oXU3PTQ(KTgYRl=H3e?BTM{@A{d~T-E2#0{9U|lccW(C&)+O0QC%m42R?xhTH0>k6PdG(jssJfb#!E!Z zh>SLgk)nV?+^;Cr65HMWY7@ymDYOd|E4C#Szdm)zaWYCO*(qUemYJ-9>knUV;wJ1H zRt~&yyrjaBg)FeoD=PzP{SEkT7e;IIF*s`TfdD}Kw=;B-DXmyU2ZYfR!dYWn(yL$+ z(6(fjwPgA&>DRX0X0F@Ex%$Y`cRkrr@ir5eF0y~&c2b*zgw0l>)^P;*(ws)v3o}l# zbL&h0H{*g>ix?jJ)0%ZBHkFyb9E3WemdrGrCB2n?`XTc#YAX;f&!??N;e6s}3HN8k zV1PMqhVP6tY0v+ls^EPhb;Ls$J58aBEokeuXQdG(Qo%S)E0@zc>hf2dxN_nm%K!J+xLkF>X%od zTzsRllX%^v>SrQd>XVEtZ#JaU%>Gr+jw=~)t7~&064KrNqRN}&WXM2Lih$FO`GQQX)#r`g3`__@{qEdPhje>ffP^8IR}`%F{L%iA@f>Iu$^SCoQCC3d{>(0)iakWRDeH)VEJ&ElI4`7%R}elS&|DQ>~pe*FE@NW|VrY!0F^MbD5F_)S&B=%d9^*;5h)#U>K$t(@!oDf;$b|Mc(GlOUQ+qZAuQ* z@R@l(vnQguZF~x<$&Pugj{|-?+1X*+p24ZbB}nq|ZlKP|kAB?Jg(ZT@4=i3O^&J8G zSO0kKIx?B_&-+S`S;y?Wp)Uv?W(q;@x}$K{!}~xb3AH-p?q6pg3gYj9;xwRO3y8hZ z5onjm(fkvaXD0CjYnxp@V(Axx1{z9gsJIv4jsy~i&f;j`o*a>7EeX#AffZMJ*4zCM4Yva1OyMW7Aau(N2H9oKHlw=ywxq6x>1o>TfbL)=TtSrHoD!15&i>NVbHc=&^~(#NcRBL z;P66;dWmc%#}{|)uH)&4K@h+G^AGL}HyWKilPa~?@@~lAQNv+qW3O3YQDN$IeG*jZ zjIY%F&U0_A4eskIoM+}TJhO<-^x>rUh)WP9HrR*_Yi7h#2}H>@Fq@WnIqU|q9YKh3 zFj30JJ9_^L@+#dk#$BlJKu>aYSI&&h%=W!@skNg0Q3*Q?1EN3+LF0{%7`;t`Qxyz? zCSf|$wfD-#Dtro80#MEm^s54OK0mCfdG}Fqqbwa)Huw4$**VVnLFpF|=YYsR*kBAO zS4aZgT}0kVey=((NPpSjm=61^QM`lQCY)kUO*iNS3w!Q=&&SatN`A*Yz5usW)O!8B zRqEEWYvXZ$*|w|=jEiEb65@R5(1pl*32yg^MIUS70S?_b0=YXH%v&{k7I0#2 z;eKByy-pPc+%tV%veOjnJHq|j8o(O_%Oen00%z6%W;Jpq;GJqlHch50ETgL0dQSO$ z$B5S)T1N;4Av}AyWA0gvH}5@N!MiuSC-kizSa(Y$fpwZL@>-6mDf3cG%jbuVN$mV$ z{}VRw_FbLncWaG#FE=;13i3WkAm)zzXM%7qHCJozspT{)Mx1%5+CJo}qP{wpeo?CI z1>f$&0zcZba1%l)ue`@*+fE=M=(0# zMODsnK!fwUe)b2~nBG06%k(>>N-Z@9K_VPH*<5f7l#;@k8B>NMQ0lWgoi`4D^n?&R zYY6~al*&F^_mLv3=8a>Sp}z{nLjDbPgJJz8RK0$I474#K*cE02j|jZl_fj0 zn`$Dz3Vky;<>YwPUQaWs#gvCZ`6oe$sz@wYS0pCS1x9Ky5tX}UE>sei zVtHF<&Q}oObvkuEg?b_RGrtXj_ySpydKOP-7=BPpYGsJ~qAMtG{^M%@q=*6@HDydn z9z(eLgx3yGuY8ShB*)?H$yCj`G65_3I_xyIQQt)}rKToMs-fu}qNMh%zVH>+#Ul_| z>N3o>f|j9W4us4&rz@lIbDne|atvPVlvl9rhhD!I4e%l@`DaGl`s(vn%unz<W$!|I*(!tYSIT$btI8&&|x90{2p( z2va0YS_#4j8Pg(2YuJEY!$)apy5!cqR+ZpW6-;*!7ftUL*{1&*!>eyRKTN}BH`H6Y z=fn%~H7LlPk?^CpOB|=4hiMa#@j6m4?_sYO68r!&^Pi2`+P$XP`TOeaJ}=LQJ|3=S zNw%rj9@%fu?{Vegj^SyiW?xt@EKzwH?+>|>Z{Sg-Z4bmk!0l@y73KS_bk`LSGj?PM!$5v~)0vWWMFCN=e}Qfo>X= zao?#9TcAy}9-a%^Q^SM#sqQ$MA4Sr1$miYLCUMUa=f_V=^Il`57(K^u(Da54X`%Q+*H^34uV`t6oxY;F?K^?3CEyeZGke&QVRrF=rCAXYt zx{b~kxtwT_M{=FoPv{SPD;9_s1)BeW5mRWN?~GPR)``)HfXaIC0Y~xY@_{ew89Ho~ zJBS&Y8QI!kuJ_f9@=!$CABjUDlrnu#XLgBLUeK!1CmUP+4!lHSz(~>_5$PeGt17}T z;{;`bE_P~~NWtF^TE2c3kVOT68xZthEC9RYCr1z^+91cZ&d))@5BK`R_rQ@eR|s2S zd1Iw-o{-vlEgIMY6^jy`RO{8d&*@4^9-nRxm%ffl(e~RFiKeC#e8QV7DQD=@h*R|A zP~EPX^F(lHiyYLaqfW3PFh(sn6K=0hzU{_DP3(;)8{kNSBM4vhCvs-N|qI$2+9CIl8T;BfJ!QtJU39=-3+=l=ZeKfd?x`^P+v%pr5mxvuN= zdT!UPkqQrhyWj+4BYf7~(TSnBG z%mMt29WkpNd7T_m&+rvKYK$Ej*BrVquy^fY8Dz%c)r_@vwvBc98x;ffl*-ZvJ&NVg zx}>(x)LWf@%bC&tq{oO_M4c=Pont7;Yv{U&wkI)lo>BI)IJH@=Ac~&U1pHS$db!r{ z2n~GF2S&qQ!Jm9l3vnAPzw5gg?~NI`ia4@7l=*mMzDC#V%7+}C5df*8!;sSA_z{={ zVRV19Ar%J-+S`gU)@huvC2*$sxH43r%3ivRj86r>nyFWRR=>AHcT3Z+H_*eWOR=@b z2$4U5idRFiit-FLO&oGBAgn7Zzo%PjWf(lBm1srC3-IN&PITXGPSn0+`b9a|+?xiWv8TZ`~#|8e z@HIn~q^l@%gRi+x5r8WXwkf*=HAPR#uXp;80$$n)QeI5TwwtTvBSEbf@QeN9t2=+F z?u+HSy>r+;a11{J=^(!|mjX549bnRkk{Nz@9;i@sKUg{RG!={&L`{V~)mJ7}P)a^s z?>mD^NM1jU9*0zgG|h0Z-$Ua<0-pqkD+M#z{*qCnR^K6x9<^UT9RLXi+)HFll-r&f z2}3j(N&&drp9T+3CZ4Em*7*UfAmB$~Ifz2GbfVl_eaz=z`vZTdd`T3K0|D-k{nDb` zNo@YCeTDaVbyP+e!E07=eed`d%Ki^k^C4a;q~4jFu;dqVjMA+fYV6*RN{1!xhc0ZB zdao&Nv~Luq&6-c2113{ou6b)Ia~kym69lLHrIYR?Jc5V`Og9wk2Xvo!#DkK}C9ZqT-2cdBk5+GWA0f(6)uDJSnyM0|_uyeC> z8MrvrP`k3wDD!N{ZVBg`@s8US;!b1R9&~kFxFTvfcHfJ)S9e~vt&Tw-`oV#u!dbc!*bv?o< zrmvxJ6tER2eL4exQ72d_HJV@_ST8X@L%)hVLm2tsSmcm0oIS!9j7HI01GU#=MzqC? z@LNB`j6~17zcPIndz-dV%zS_{6Navj=u?2BMPwqtfnA^=eMztpk`JXzi zXQWv%zu2rOBf$KM71WM(D!T4+Z{}_hakZCoSNjycd0JuvW;V3>`M%Lwa zk=xwP9y1(_b}a_a&p)2DMPE>$Klkuk*^S++88vype5(TCS}ShOFjW6}{_;2KC+$~8 z*T@yP{J{X`M5nGlRL+}pdAH1@XO5?Le=DrHL7t%xS&BLkGJ_lIB)aS~np680VJ~05 z{$ie)E@8n%-M&j;h#TTakq0TIYUJ zrqe)h#fAEQPHc#j7SbB!1lr=;j`H0%C%;2Epu%hELXm1DX`tklKc6&epo?L3sPuWa z<|KYygvIrIr0A4r<2{q$^n@;XYZUO**kD`KAU}z5bGuUq6%w}(J1dU4RHSsz z8JQhRZL5P=!d=QukFwT5fV;6t2Kbn%ji(7(ALh}gDZOYWcN>^yZSb+dkff@3&V-Gn zja7C0m5d~n=kSH7A!TLO&Tp)*Wo3&S}cKJ#APyA2Pk1br<4UIJ7lC{x$nm(aVhUjlbOpT}E~V@Wu{Z>qE24FtOd~*KsV! zuQP8xhZzzsQ_lb;{&@KxDq7G6OC|9`Wyz5ErAssk3wIU7A=a{}Z)Mi7$UTM#5cfSv zE9aj_@XCN#mJf2NCO~Iqv9v64u-v}mYQu^BN!n8nfzGjeT!E499bVXXj*WCCw~A%r zy%#oN?ZuMvkK8d-larR9-Qe2#aDjuwO)-YORTj?W6c;os$Fa}REQ`&m{qv{o`=2}1 z7wch7@gJQiGkW}yW{T#?@oA4`%fFVy45#8bued=<9m%787yLGk&p4VLiEl*VxwWgB zbWx${%n$q{;|j;!_QZ=&g?Ida$K;?A71xc0Ark?!k0}SOF=tswh2sFSZTgNk)K3DB z;kbzK5GT(J{X?PvOc_@G9}=a$UV_Dm;;vtBpb@jPi&<_&Uhv!pWMdR|mi>@LR5R8X z-x*SoF8)IsSPaw+!hSzDSoFO$8oKI6K_-&9H9J28XXpLWbM-^m z{l2*{8KyC`%WL zD(#l!+Ia^T+B{8ImwEdG-No_+d~X>OiglD4&jVOLM#@#V`MxLd)NCk47;y)PG@z6w zdDVXEckc$@8x%U5S6g3ytTq|WIC`V)ysOoNZCUB&Oj+iZeIs)2b8@z-Jnk-#tdWP0 zi9j%64Cc@r4sJckPrHnQ&t6adQ5|`;mu6)2h5=)pVaiR~1#QOgQ+v+;!bbkbk@bNq zgZ*(Xwu62SxIy2J#6xgm3VDL{Ocf=ei5mCaMmMj3LNWVlH65`_eSCl*C?*6;xS2j- z_=@k^^&Sa_-Y@Y~J-1hygFu#KuV!=v_w`Qju}6%k!Z^$5t}^C)keoC~Wc01` zi-TUqc!ABQG-0a$*t z1WT(QbX+RDC9kvV+r2UWUD}VPP%lrCT8eKvR-KnyOkvMPMYeQ5pn z^P;6O6BwdM{d-_qsX?&C9FO`CNafdbNzfKtVpkIG3Hb-M>E8^5{s{w7`Imt>#lg;3 z-#h={tyEyT_&e9C`}N>=Qnv4bu&HfGOsO}R1_8M-{~`i8^G?cY4jXj_unb6Q<9{#z zp%TSs9i!5*y(Ec#b+)yq1`$HM5vmGB&%ih18;t=HEueAk5RJ#^6E9*UXLu18Szm6B ze4$%N}27G+fbe{j)$** z(g`YBPneqYOY*X4m0-URAXeD$`VE;Tv(b~BjCDY`%WW0juWZMA3 zvzAc!rH;mHM6Y8%N`Hw%`soqUeR-Df#%uPqCBBa{IE-t~Ml5x+oN&6S@jz6M-OoceL-nVLG`FTl{N5-%4UTg~H|&;z*48Kj4x z_NM5v@}#KggUYb6U8Q-Zt0wrc&XqQ`B4{D0_otwE8(&dUB2y;ZlOk88dDprv?ee$P zHBw>M*LuBzRn&Rm8fV&-@c||ryRqlyEu|i%$$LwmX&XPjpl3+v^H_N!(L&;v8_$c$ z9a+3a`RKBfSw(Z#TFef#`5A1y2I)KRq_rpBHHdZFX&YL z+ZsJkA2YkyfP&^hAXdwxUTBxP#Gn1%aE9tRZ`Y4`UdX*oqzpEF^bTAy)wyl{S~(5l zT+Y9E%^FOgvb!?mZB}Wh^WExdD!+?Q$XA~zX+@&1Z+I|=453Rv5%Rh*6BDx+1IMpT zmBR?K%~zklYc~OH=t4GYj_|(bN1bFbEa%RC2)O^a)GG7Xzk8? zlv66af2f-VasmMU|Htm_&Ekn0Dj@r_rPGGsA1^KE751G7HeVV7E!XCE|`ZdyAK3!m&@_OB|xwPI%{Pz+1465lN zQh>OIt6FCuXf&ai-1GOkJ@rI4Jbj!?x2o&o+buThL(N6qOV5UD)f}ezk2KyT79Bq# zQITua6)C-77SwYHYsb8r!7ak1t<)&4IQ-5ySMrnW)yuv4!M9%K_(p^Z>w51$=oLSE zrv_BTLs!AWNu1jYnFnnEl-iki@`pwSICrxTsyT)S?f|dDIk)H8&DCMJrvFeqaqyJL zw*F!7hx**>`PieXJY7*USbAydJI$-8oSc{yF!2&OW_=Wfmo**++e5;YA??s8h)vxP z0x~anoia{xN*Ak%=BkM|cGzM(^JJS{;+gR6m&Knt-*oJXQv;Ubp3A3zribdJ{ zsKIjkorlmpEQ44k0nv2`-kZd9?UpyagN_^GI z(R&lu;D@j)Vt678_axm|J{1~6W&!8CnikXeLV7qLe<~@&pq{HFR{N-lrF1P%RlY+8=i@17B|}%9E5Lmr1siS)TfyXp%AyK#ERR-NuA$^KjPaNzVG=N?y7;f&<@uT=^OCySSwl`qkg_b^& zKqD1;5=~rYh7Mk3hXFqw=h<*~CaJqH_*iTZvLEQ&&LIr*S3!7--nV>y({Ub88&RS6 zyel8i-x-1*TI_E6wIMA?%KpN5uDH{YryN-pO!&wTfQ6I0a~wnB%9ULFpLm}|}x0IoQLvR_Sl{W+rhna#C%M%BKAlN+McL4oO3p$Gw^Sg?| zE_OghHRlu6unuZ&`)}~pK&xXq*UQp+^fi?|3U<3@_)~T_or@@W2$oyN)%9w1TZ$giG(hTZNwsk9yUbrPfn=k7Mi?db{wip&VRw zfGaOXRdm*XctsE$?%vfa04S%OKQF(Ux^onjliadl)V}I#`s#T!f0(NI>)7uri8_|& z5_y89tN9H+0#{i7i+n05(3X~nvCo@YRfe>nLy8G5z{UwJt$rriB#c zwwqI~5}4!p*7S-KjFn${jug4_y&qp0n8s=OGu#hj;5MC1twv3-P#De~t2^?|1cB|r z{#dSh(ODQfNdr$thfTDrOVtNF7$kw`y7g#cZ98B{yxL-hs&lHYBvsx)vWVy@G4lAg zEGqYGy6=F3i4t-w2Lwr4;n5+0d>F#^h}iN?q2D?uy2kFacA<8@wfxWA{YBmHAL5~5 z7{=aLvpXta>R9!f?C}Oy-Tw#?{a3#U1PA2$=hH&($O>&a;hw>5fbKJ}+P% z2HY}}@(DqXxu2FV(0pw$s8_BAyVTU(=vCnY^*8j#82Zb{lFjwg+X=_X%fD9PL8eNx z^Fpyvl#Q5<8GaiBuHb#HXADudlLjv6-*_R+5!c&Uqh4>O_c5Cb#A-J3>(1R2%Ysv5 z68@I#)P=;Ur!0&rfRL2GbWQN&Wl9J)eQtJ7NPDvChDMvC|cM)nFSo(i3R(x9M!Bx=_7Rz%{P(sT6?aIu_v(r?{ILpt0-{a zkraj9eh(an*=48%RU3Swm{Od<_yklrSen8~iQhlmHYWFLN4Zv5b7un%kiNOl*9B;) zp%P8A>DxDz+P(KNY($%5JxULRe{aK}1$3rSmJ1&^Gj{OJO=)ue6ZWPi8_(c@Hx}k}yw~I57iF#DEmXLEyfpjz*3n)4Ymr7{WYg zQ(m9XmoZ@M`9+t0u4B$dDw4srA|yDNWj}Cjk#_!uW(e}*?QceE1bu>?^N>5pPaDWs z(bLJcHzmBt_#*Iry<;qMo>!DfWNGqq?}Ko1(#E(!Cd%&@?l6(0A~eIZHa|T{ac~R` zupkKKf3|lwyWtlujZXGukul!;O*eNq6kuH0Ptw8J++X$=V6*x@xi`T(>j$TvuB?~M z?aCpVYA%BuKiJ|t>CF#Hf2ZEu?NNUlwsPHNMQZ1AM2jNc5yf0izQ1%q!cX}MzhHCw zIkLoF1m!GSp7pPjS)-9w&ktQ;Z9{yELu_mGW;i{#-67z9bP5u!u+5~=gyg{wY=Fnv ztin!-y*T?&Wp{my-fh$aIRX4E69KT5ITOWLpxwSe&a3iAI6`mpH@y6)8HRsBgv@K7 zsNwVs|MsHyPGp}P_(}5~Yy~NEX5W~Y=$iHt-&j18_@!dkcn|Xd>>L&dD%on&!Pd+P zI%6etpm}r@g26D-p|k3rhLjiqw4EQw1`E)clcG@$V3fF99ArOIF87b-6)l(^0xp4b zipLr|pgD*ZvLpQh2W1q*@Pttv#qAI=j<8b)lQ0v9xApfH;jhR+|QNe*}<+N;D(DC*#jnl9=aAkY+3N5 zOC2Q64~Zc z7jkoLWXXWF=CIGj(3PoRLmeWyLG7Z^6CMm$KDCuN9=`cx=_6{yNJ!OdT<~ zT-WKKve<6}D|~G;W5OpBRzG(K38DEfepsM&oWtqD7A)vS>1#LDtya2vx+gsQC7 z7Jj-Ja6f%hKaJNu@!KD&MiWKQhPmS_b+HzgXO3UEpT_^>45nK3R+F}Ugnv8QHURgu z(QVmqKCkHI)`@#|F(O1bA)rJ5icy}Eo0pC9+h1s4(I_$AAXg*LFmDjVn}J60>VDpl z`B9PI@F^$>fRCO=V0j-=X24~;K9iN{Kv~Q|67zpzX7tlT*R5S=`IOSi;pkCMMk^G3 zpTys~{CeJD4Ip*x8xm1WM=>YOsqA0%*;Ty+hG?P^W(S}Ef^U(lP|LCKd0_ah{h^}B zLY7Hgv9fKVcEy~VMtE+#v0hepAgnGU{{Qc!^+aHi33^v0;lIF&b&nlCog%jj#wKfS5xzkJ# z>&U5Y*YfN8kxiL2pQ4A9qcBP)!rwc74IsnKb^QFD|Wb>#$EkblhaJV=CJB#TWDkyAi%~ zXl8`Q3FQD+G$3pQ@>O7i)VZ9FDkGQ{o;qxJFm#*sO#8l1Y%)##p&9X$+9#{@Ep;(C z^$!&r*3sE_UzE5-$eT4-KX6nJ0=2S4_y~Q-`;kTTJZN>yP{^25^6)d;9MvAYkOU1| z;EK7>8n=FYcf<9D^p_g43jw@NJqKcSY~&e|qzJ=jbO}`f5DPy0sG_~gmqwF%t!fk% zxnvuL%lE4ko&@Ti=jH-w~&h8_IOml)iT>B^4=ZUFD50J}uK|`xHme zHsXmA$6Q4FIrG*7;`}@SKA>^-QJKXFIwd2BN3PD2DIS5F7)h?AeT7Q1*f2|_pQWLy zhj%W=Ytd8f!>AFAcgfx#I3EH(^p1(sf7u)&)&Q+SlHn zSeh|s5F$ilKzt-3-nJFmDAOUGzP90h+SEBoapJ_rG3m({4konES@MX$H)uRHz_pW7> zzS&c{U8@l=DCXaB;CBX9D&HYp44RUs$lFqt5Ps(ceaRSnR?2m(DD3`{)FFnS(pg%A zV2H-{d)>}>)J8{4D(S`{C3v(=@QIv!|DHH6o2}d%c5aoYNU#8|4onY#IMblY(%Qqk ziNy!43m=}Tq%Dyb0h0?GN1IpUMLZ?!%a4~kx}ML2?2Kp)4J2h=D>Dk%r|dOxx0Z0J zPB0#=1DKGi`TDS<(Km2M)#twRerEffR%SJUq1i92KQ%4!+!Tl^>U}6tN_a5CDk;a% zRrup2>3zm)ilqYX#%o#w<38UM4@X5K*4lnysEHv_Q*@2*CEoR1N+pE@wjF%li0AJKqH<0oj;dH5+;BtMF7{bZq$1&O z>89`xK=W7O7#WGS2AJH+7g1l~=zR@h&Kxk1Ee$@edNs0-e5wFc!@23+3+0JMO7joO zR7&UzPKm|R8gTBNQk@_9P(q)$W@+?!rAMA?CGUXVQ2)AAJ_QC4H9)}p7Mbl2l>ll^ z?K6syvJ3ll0K|XY#6y4xLsuM-qmDg`6l=;}ZWQa^Q&~$jI0LPHtLf8UsOEK4ufId6iT0{#k=4Rl{7kigq9JqNZ2jRLgUk)R;59Klmb&+Tq9^2cJn!>xdH5J- ztY9|lJ#bz%Ck_OnoIJut<3^s1F3gLmEo}?xYw!D9#@cdfyqE2fNOfhp zT+)^E3(AoT&p5A5>>$NPx-+mL57hU}aJQzE>xmcFB~&aO zs{jFE?|JpdD(v3BBI~3Vd4?1`8GW@|AY;ZY>2!qlJOZ>G2__`R0MBy9B=L2N4dKk* z=i2_=u6CVE(M^)h*RExFi-mK?+cD^GF#0Y${$wT|2A)OR@(jr{F|-TO5fkRIbi^)Q z{O^z;572o><)AOM>DauLZ-;cI#~wLT6>ZIBKXD7U5(;~-cRU=U0!5-&{!l?7yvW|9 zBvkQMQ$&MXv!&on4CP9)&xLe#7TX(j;wZJLlygmZ`oTHvfJS~I^8o~(n>Q;_xJKq^^RdlrX=R;W z{W`EJ+vTDBGg%FoJKl`8CG%w3T|haNhxV$})*D^Lfou(R|m) zF)^sQO9LLgLLGBm$6GW`>8>*i3_4maVNog2}Nfz$y6(b?W+ za4;IL1$x`4Ky=HnaIEbO>TMBp(;h1Nq*)hZTt6(Alm7d~Lor#~bH8;%BrQEsZv76Z zM7>0^k*J6sc)R-qg|OLUZ@$sCalquMc3&l5bua$F{9TvZ8!gq!r>{OzDPf}P0c9FZ z9Z_XY1;apeCE3sGSW}nEvKE_G<{so?wo1+ z5~sbOQ|GR3M01t%7I}|&^UIph1)|0rgca`pYiq0Amp>xu$=g9)df9jt!PYeMW_4-c zDq11<3H}6@fP{8+nE@m}9BWiNEI$)6!uRe5(d%If1){#k!~dc&czW9`ZqvC4Bl6Z* zGBvBE7@XU$!N_*4o}OLql&Q&MDrG8_Q@7Yp{AMRz_CKjJ$c=yR#nXJyEa8W()b4on zX0GuPycGyjvt35rKj{YER$J@Lqrzn6*prL#kw+dM9;!nsXYR*$_1v5|^~OF{G5xr3 z&#q}&bf8}JSRdg>P^h{5wpX7RG!HEdW4RkHXt_6reoX7k@4HZWL8Un{@MguAIegr4;*#K1Hg@}r6xfH7nLv8k z!GtwE{IkAsxUKGo?-p~_4KKt?!+wiz=O+%U1KaR4Fl6k8tXp>j{$Z`gGZzHfvH`F< zcvhubz33(KgQ?-@@i3+IYCY31#Q=*p(0vrM5rfNgNKC zy|+KE=3-e0w*iYpfA|VNQXHEIbXP|S0Ge_bwh{zXE$80&^!uM=+YVu!9N)@UmHt z%T|ww{<-j|rtSHmlpZQSk+7~CyIq?CN9M)E08>7pFt>tOi{_xjSNe;~S=8xPr>paK zW7aKHrSha1@wU#dZ%&EczT!!;5?9|V!-)BtZE{Tj>UUL#ioH7$F;l_VedeY#ZXVlS?}QWaHatP5kh+j4Ncd5OlQ^&w+4 z9ko=kTL*0s(a&{W4ROS8cbJj8@@UxOnw{;-b?aX(>HZ?e3x3wcO9N>+#~q5l7>c?m zR56`+_1G-`eWjSiPot~iBv$`aQ((LR@!ccmrZ_7`8tmA5;k(x~Cy!qLG>=AFk_!4g zO!_Fi7EWIo4wc;{3vNnITTr>KzLMCRpqv%YYK#HcN+Tkt{mpAvt=HPk)D{vnhJZPk zfP*&Uq`V~Cb?Yx-X@50x6Z`fcA1KLkhQQ8%Ad{0kDW3E1El+Nva({&z8u*GT`+kRf=lv?ps)$Cql(tiYgL~*zcAD#v`|MUZqJaoB+Mf2 zr-tNmpt61MV z?xJgLZU5Oj8po{kHw6#E-E=*lK+nV^W428bY*ddEZkvgI|C|PS0qd8?KLpK$_~RMe z`c>IrGH^e6hKdPSWVhF4f!VI_H0emwaqR8zIp=KX76AQzFT)xReJz@BY+^gN;Ef`D z)5eId+5Xh&C?*v)HL{njSpq5+5poY&ZHRTjBi@nWCLq+!5?ahmd2N)sJP2&5XY=St z^q4KsV8UzmwG&QcF1UHepHZml(CE1&L)SCNo};@FQa^MkU(w7QG_LcjhL4I zp)vsi>$U$Yus)q)N>W_a(L8Uy*pd66fn)#M*SFOXA-xAs63jhMxdIcsD<_uyP}^Hk z@UHOLcAwXzc}&*^WIgB=>J>mY1EF{L;Hjm)9;||w_0`*~PxUWU(L8qT=&e6~E3vml z?+)+y&2lmUoP^#nquY9`4Eq38)F%L3?PrexrpZ1VbIYN6-GlzLN8eo0Q%z~H1~hU( zwVtPIU=g)nY9IK>CH>kN-wid@(!FQ*>W%&6XpV!WwU0uKb55(*wR)$z7`JGHp>#cD zli0kD`Vmq*qHpKdUQV9}hGP}Z*styK5aqr9#GLSwSq)m=znhT>w-uEB{vudFh52@D zKttK+#ju+P80qkH^00741yS_}2An^?+yqrqKPFiR5c&46nQ1V8aFYgmii>W`=Q#Xo zmnbm6)Aem{@b(mC5Y2BbV1D6h-`mGr|7CI+jF+k28;9-#sW>oyU>v&;#bXYSM^U>i z`J+IS(R4B4A$}Yl?h?ajR?KL})_hJ}eY5+Jo?Hve7y`v0ygB5pv^&^ow3lx1&Jeb} zUj3?BSII*g5it4UiTU1Y>#|;8w{g0m_UccXig{1>!%V5J8z+@7y zY^)*f5-_0%2*U>eF=274Mf$(E6PtT)DO|U5DM=lg?LGVAx{Q{=vO32L`L!3sAVqld z+VkZKHFebtdUt$U1U#2)30BU+%a3O;70qmdH!hRu7wfOR3Cp&UQ>#O}=|}qp!{<&CP#rS>gwYH)nEX`g$Z}enV!AKP=P; z>?@Z#fyq$9EMGN!?)alI1@!S!_eZ;CH^1FQmYr}kI)!E@% z#}Dxx`-$yY9`<6%9Kqq24C4g?(lEQ2M}z=|#S(mH+{Zhwo-OePi!>a0Dlzoh470R^ z(g<;HcTjJbBmBf&$ycrpmNw1pY6PI zz)2}J{5o_ROHa_Oiq7iBim&<}0!ee9j{qgCw3;yB-uRcE?0-9+EJlXWGynhGW*q2`wk{J38WJE_qJXPHnGNJlI zb-oml-(kEXkmkvGY~xhwlwwn{SA(l6nWM1af45iQY-}w&tg`QApIgj7l_qjk=9Tut zdumNz%i1Y?bpr@lkIzg`kM6{6Sd10IJVi{Arpp%o>ArhKQXxkVEU(uUI=-4MDzLiC z!)Sk23nlQPT=j9pm(@y&;A(8!#r1Co;|f|VE|02w5RF4imbSmFo+V8DIJ0+?CAg8H zF1ab?IbZ}9-nrKQ(^ZBlHKcD8z`jYWsrt#9f9nLM!3o2fI!FH2*D^nEn;K@UQn~s6 zxJf#LdB03Y*L_NF*B1(drrOqo(0x=M+T=?yKtP{(Y@d3Sh-SvdA$Wi~+ zif7xCP?)|B<;q7MGx9$%(QPMTTT6dqqEB$F|Lq0(b7HvrCw4a@brY148yquR$^rtU zVvOqQ$bK-Z{Ov`JOyWWNLG+t9y}z9f0Ak|?m zWwSCj*NV|kk6aLeKfVjcjRyQak-{!nsIqoYGP#%zZSF$15Z7|nUH6_+uKd6z?vf3( zTj|3yn=h6In)h_fL&ASmFie44bkc0}7Q3U<^%{AW2Hl*kJr;(j&Jf=jdrRYM3bsH4|5EK>F{%pvpM6*`9&a+Ob2Dg}zE^TZ;N@?^U6 zTi@%D_De{u%>;lHfpN~P(Or_b$Yto-zA2=E{@M?ep9Sf>*Z^%jpfGcz6ioXr!Hl@P zSwOaEV=D^hz@}W(>FJ)GOY^V~DjaE@OKza!l_>A@WbysnWc#_mAKaFS>j`!gc$10> zHHC0hX(>g=&Bpz21aBSCsvZn}mD?IWo|OS~YaKta?5Y=CrJFzvsgEi=XLp*END^d- z3F*_Z`JoB2b)#y3n#jk33R7nT9a@xQ_Dsxv+6g5R&BIC;I|fiwWCJe=3jaT00UCel zniI18|IbUf2G|>CW)mnn1nGvvONh`Hl06ft5xTzcvnmFuGJ%lqA7^eu)}OorHl(eE zl3WnV=2i)J2DcfNZBlF=oc+GTi>Pke6WbPl{+jyhr9JsO$i>+?`yWd7-B9`Ry;ekSTGkSVs0A}8#RGc+}B#DaxDXW)eqWx zd2s3AOsV>Zx4`2Z=qvu<)eGeD*YG9A8NB2HX6)nP`59PSS=+Mgt@-v8fQrQ!11fDk zSV+u%yi@wcmvKHD0-~VpcE*ASj;t-ENSs!)`ibWlD-ZZ4ZIeOi44evfIxGTie>g)ptCzW?&_X?-%? z?`u`h9#ce8fiUBLVsUQ*Zn`QGhsMw}G;h1FLK$C;eo)2<+SH>kw~zf-yo_G|zs24E z%2WEk9Z#KDV@^UtmIlyzsA?n-DMmp-%665lL3l&(xN#fb@Fd%mVk?@{aU@PP{Ht(` zX>D95hGr|Sc}a{EAYca;LNz@@80P&V3Rzsc{xrp-p@KYcFV`T+#Kq*bbachu(pw}I zf{!o?xT|-p38O*oQ!l)MwXPvzdO~de_pnQd77eJ?d=Y(u2Ol({LPp2ALjcsu&R_s0 zCB=@LD|&1(WF;2HYkit(hwt{=!hQ?a5$rvDXC+^ww`p7DaDx7bnGJJ~LXNU~aY0Un zPN9!AF{YXDXvG#(RSta@mu0lA@Ox29B@OS|#bNshr?c(>3#;@zy(n@nO31ocIQs#F z_uy6GZnuvrin*RK?ocOt!fyIT$IgX%m?j<;0Gyrtx%Sze5v6H4qk7Qgji>6EYCkF< z@l(l7-)Ft|t=TYIvoYfdpEpL%DLv;KC>`2@WNN6GYR@aORThhIl`j`Ao+|t4$P;bA zXyM5_VUriC3YrxZy<0sbHsJjJR%9oaNKI;Hr2pisIRhRwcM$v>@9xFd9n{d)fmTj- zr>w;dPMps_me+y-rZd&HoS43tu`<9xJyAHTS(O$v1M~v4laXj*s+I5PQP~|gAt(tz z1MCc#i01c!wI&eR%es8z&Ss5m#nRtj4=ev-xSCX>`k3WoGR784b`s`^q$hz|i5H_X zM=g+F01M6RzjDa_p9+U_z&?OL{r7B?IE4^7R)*WE-SIZwP`s-!?4e>L%tj~ltqSJ{ zCCAx|e}|oP8TH)2AZ#*%zi7w#T_4vGIxpGFqo$p~d9)SeEcI1<6qsjE?3r^Ry}LIk z(okyUYv~`9y>}oqj9<>_eedKsW*o(gDJ17P7F$~Bf#Ocn7+2dlMKVK|V))>?%>xKc zAZkoS-(TdBoBQcEpr7Vfl?iVGcQ+-*V_JGoANpG9RCHJ$+9?SUG2;?i^l2hrM@a@i zI_FN5y04%vu3=FAr>$Piqu!r42As#zj%n>|8;VvH%xd)pCmmUVv*f)xrWdpH&QOBzbfH2zd6G(`#0U4j-9- zzF!cRyQ-zU@{D;@gebBPBi^+e{>saaTNj={e*iAmZGe$>ubPKcg76y7cTOtur6vZ4 z9Bnsn=sEgt>2lY!i!p7G5HJfuc$$9TE6jbu+dk)kp6otVQ zuRfzu984QQVduRRFW)5eBJ6W$hrnyjA4l|fcwYP?q&y(`YB@>iysH1VaQM%Z&i}1w zIC&p*Vo~^K)4LfjJPl%n&bEG?VT&5)|9uo!u`y zx|sO&W3Mz-47!+#Jx4}jUhh1Pc^5=?39C)IJwK&Ax7CgwmuO94bXLy?wSswX=e4Hl z{_yiCve9t)bd&nSM;}6r^<$Df&~LHh1u-oke`RfdXX{S68uvoa#2T6J$C=1BP@MJR^hTFBn7PQh@_o+$xa6=-1 zV)^6Squr%=HF8TN`yK7vP~N0(7Gxc-4{B7QHC+ga0$(hS{Ny&`yo&+P!~WgAtBn`yl&oLCGu3G0;7jYTBQ9z!rZ#8 zLGW8C68Ud1p#hW2&|Z;o>)q#kx9wTzHm#OidtX^sN-s-@-WwjggXf(7a_bWBRwW~O zN%lttitAgZ`ND6T(P=?RV0&-c4l+@C2%Ps?k z{aM@%P1L)Xb?WR(0YEY#&461{08T29GjwywqFfTi(Tyr}n!QW+3DaC}fZEuNpGY>S z|Ho111oRPde;sxITmSH{_fPZ3?6ClubLv7Izvu$x7GpOw%frIdsK(A=_L+e);^B$x zF3QbS0--vvn`+c`J&&m&Gj52^u%Tx6R>vBo2>ytF2*9Ezu>a1ZBUhBvH}==--4J?_ zSMUSSfU+LSNInLxJ1Zz3bB9ngS-k9Q=X88sYRdV7+F8eYXR0>~lSo&wp}6Sjoff>7 zf|%EwiO(sTvAj-H1?gq%TSU<%(2Yos6huMHYgPqU@+vvSXwfxc8NHN&JN#hFQW9c7 z(64~H)QN^#wg>dzVWk@sFiBZ=E`qtt&H$mV>ZR!~%Pw*7-d4ti z%csZUeOV+IW9TI&9Y)iI(OXXOISHhh}|-Bz+uMKQ<~^6;G7 zW3aMPig`77H6Lus3l2gKpqdzpgiC-)znH{aR{1cl&Eu=q8&3uxKW2__{e!Jk565x5 zz1Dw$i0@2JKt%Lk5b=M*hECo`{oN1k-?f4xkXgU$&VNtU4|wdZ5|d|dn&?6^5Qc=j zP#9%DkV2RtN<^;Ey*)?6vU>~y_$M-89pgtlA$sFODD*7H z_pRaMRzQ8e3OTQ8;-67cU30f8z_@noQyagKkha?xOenYUgu*PM^dYH?PYFx!3i^okbVi# zSr!I+W9$&(-mfRFigC>P+KI5lc_BBP#ejm zeNV-oDIO#z#41Y~-;5MEs-Y8vf9E-od49+7J)ZA-rBM&ya+;KMD8q36HV;+Vq0pAgi#CdM^ggFBgHp+SaOY z{`b~;Q~X2Q(bs3Kb=O=?Ql8bNEZ*0@&1=``zq6bj_C*IS(J>Ik$EvAFC$?BD){cNA z&8ja(NXAUvlH_ZMz4EqfoR;x|6>Gs49Z<1aA&+x2;LThI*7JrT86!Z;LN{WT|r3o)^q0`OK{vfcaV%X zc&hugh0|NW_lmuM?--LR#+=L|+DU{3KDPCPo+Uj%O5DdK_1rd{>s%Z3Q9;hExzEXX zCrX_mBo6*EhkONUv=o`Kan_!zSPqv{sHsJvi@=+h{o7B%yF`EWiUuE5^l82veH!<9 zFjIG4~(Uew2tfO1T+ z7VTOOkv95y+eZX)-8bmCr?TGb6A5z#72~!yA3Hilud}u+;*Vk8I9BCc+;b$;{LEm? zX-^H`)eH9UJ{F=z{aA8x^yQ68MXtlyOKVoh=1Ea*sb+(*(yxZIGj4ABXA_2T#cO8V zmhDx_m%lBT^V@X$ADei4W7{rxNrX%hqls;3Off|t7zpYl8M{vA?*rE!Ekx%QSW0HH z0cY+z$~9CbPz#qk@aelxti>n`kFt!w7OD=m%9ZttnrD!#2gT_Rkoym7$qa7{pKD+rtt(d_0;b;0$b)C+7_18Y~?Cw{?%h(U6hhAo0J4jQ^xMJUZsIl?N5& zk~1klGvuY9dP>H{V*CS?xQ-oOnIwz(4J@*5UtPxTZ%ix40tZ($lh=^kw{6 zR2S=k)tw2d)jy-RZEq|l#Snc%PWFS*^n!KKT8!zvIiJM!j_KA9XteQbTbHx65rKV~ z&=Esvr7l*Y?!$R3PoyP&FjH*F19GZuICZ_aV49cdF@Kul`LU1oFCgAlOf9wONW)-k zq1>kvf&_f@V;i^NeBs9HH13}hIeorWHo)mP{kEERvW}Rw_lEFdbM3=L>JfL z7-kl?=yaS1Et4-|Q+@dXdWKsE!cfd#i93fn8=TjN~2m`Y;*5lWU*H3hkH-T*NCA~g6< z&nn~(m4KamG&%eaV7Rr2pN`xq=AZNJ{KnRjcz)+wz^rf7wyylM6fsX^v2_O*WqsL90L+RsQ#^t(H)KsMDJoevzJaIj`pK86}`|%(qT=%zn4^UxT zCwy_1`eru+-U{N<0Q-%^V!wcQlg%gQxwsRGn-}BF^{xzufrirBtxt;j@*&dPfu=;# z4%&yghU-href|B}lx`PkDta9jIoau*VTQFF)es3x#LOa;Rv@6BPkxyS9a#QRIZ;FX zY~w|YLl&=_drP(hJzkvggxADbpW_`&Co;8aKb2WNy2^fKA;}F%;fr8Bq-6C|^Z}~=OIu4Br6I?hAny3<;C(!@ zV&-H-db6$f&f~xLO~-mx`6suBw-S@=Q59eOeBGH859fzvMQ7sRLeOhAY5Akq@GURmv^4 zGfjsOx2h9Ht_}UixTy$#kLs_{B8uFbnmABZJm;L&OfI?|>N;S%B|hRP?f*Dib6^P8 zrm7(!s>`t(_$!PUU{`XET`UXd8ruFcgrypT-}17U4cs|<>&a1~8PzhZOy0@KdYj~} zViP9$+K&MhMN;>Da=r7c<>$^`@dv4O*em9IxHQ8xgA-H&O0v zX#1rcdeS5RQb0h|J{l0cM7{f`XA*iu77!_+yn(*uEnj1(T5v;if1beE4%~B25`7Ur z&NZIEnRG^g^`;%h)+10hb%R22Am8MUy)BC{duO3qB@V<5fxL*e>9cSHxeZRZ29hBt z?_kjGrOvB|hS*kWBX-`h4o%g_<*38w1a&9A!*aDNI;aNxl@4k$M`bTLc^~lDkB1JV zbihupk*++_8He9>J}K$ZM0CrRjcBGR)S^`JAWk(DZR^vO;fm9T^iC=s7JFAkY{M6s zel6c)T~9aKZp{A~8&*A=rgMHxf~2y_b*td%ttbi`wRPw3Zk4W2h_g2@KUrq%AG=m# zD=vHT_zefenR9ch&80c)iGq?}j|C^{pXb(dZmnnOgAsOJ{jC=7R0qV<^3$3L1R+Q0 zVZCyB4f{{Fp(@sTeMn|_qi{%7i118Oq`@birX(AET)5%TT1E7?8=Cz*`%jm5+Mv=ZjDwe)J1_3|%|H$etVZT)o~2Fp;dN zGLsa^FRptS3+($oGVMJG+2zASA|4(#G*XT^DI1M}tA<&H?{Y;%tM6`cH5W4))`A(b zKEXowo8sd1$Yfk0_#n>{jz6ex&kwwbi`fW+tlAH*C0v8qGX_1Gdfju-RVpfmL$SE- z0N85}YP2Ma!D0*^{Z|t5|2BO7pU=_$?ZwthQo`UBlCe5EViK&2C}%a0Lv8AD-z-yW zE*+E^zZxzdgdLa?_TVB4a7{hffzOTo0(p==`@{fVS44sw!2m~-C}cJkd3MLz-D03Z zLdc+R>eN)2+3V!`nKRWwHZwDsqkF8Yc5RK)r^8%NfWBuK&=FsHT1iFNDKEA)2pKm& za&cGdVk3TGu7A#$pO~YnX0iLt$e_Bb{&}gpS`PS=$1^Mi9$AEEtEvXAZgyfPbeNuQN8GWc*U};ZsXMTVdd6A@ zLT4>dA=Pc5piQ%`K*eA4?e%6$iko1uQBLS_!of3d`T59i9!8RHFSj5$)QYyC{T4_3 z1EZDyd20W+pa1X=L`Ag6u5K4P9eHnYJbFYvkKlw^ECXs$=atprGwMujDE2P?G78bW zGPt$?J}zIV5^k_Ea|vQ>`0v~XBC zjI6eZt07(ubRihLQHt9;Qm&FrvuKQE<$U%mzRg*8RUby1ye?0tiov)0v5{Eq^vZuU za_42_cf?y-PC1sTpA|MS>-#B4m-U#!%9inX&|)Y{!%MY_>2?uJ4I(zPG{VFvj#2C_ zKhjdFmjuUGwsmS+i;g6y<%*6K+R!p0fFqvtQ*`meF$;A^JcH|6$AGSc$<;+lE4nI` z23V){HC-M?=xSx>rf8J%F7kVQUprG9@VjJj4EIo>pnjgTPe*=^+LZi;PHQU{CJ*5& z3`zhd#n122oa+9=9gx%tNe{HpcQ7@Th7HE5ZJ}wn+nJ&gU(@e)<&y0o_mIshInJ1I zO(24PNE}OPXF92e5i#WBXn3QU4s`HYhq}o37RYxHFgOr}jXas+)E~E;u`d2JmH@g` z|E$%`*!z8$K}YByQ|^%lx~6O!(|VA({$GqEhkwO5y5T|2CWm10w2#{}EQKF`<`4e&qGJlBAsdwyc6`-eyF8w@+lx=mEJaUI-BI5g9v$#wJgwPvmKP z^mm1vYWhu4-)QrmchD$V z97d-(s`YHt%*7M@WacE0t1d`6xtG|=Af(CQNj2TK7rWIICH5B78hlP|&z9#5-v1GJ zW7gHspvE?A>b3*>IC~TOoB()p+F8q>*AstJ zNQ=VV$P7Z*4s_uJp0zRKRP1VK3}iSL@d1D&$i|}R{0}-hAT1t4Js}8XHr~P<72QwHLW<> z%{cg3(p6W5qaA->;!8bo(B}wMeIlIPdx>-zuaJ~urjOU@I9qGu=n4lsfl5i)C)Lu- zn@=w#^kGz?GsOo&o!`W6e(Z0&NI+>17RyWA;`GwtmS2M*)|$URd=YxSdB0|cQ5TVh zZBXn{PBY=Rx3}-tMRJyx#>@%nFnmvS6PYaD?#HR|;EX{vcR#}rcSa$X-&)a)mMJy_ z7b9GH1>dUgIcM_2HO(ajNI}n^dYY&t%e<}IxwC8?T7tr#qjybQsZtN8I+u?kMS);4 zWi!&%MnYb1NxqYtD{tp3Uco1=2{Ai#nu~Gi!y9s0>5G!~3lB3J6iBewy&L26U7OA zH9Mr}=5&a%Z>(FO>_KKqkhAOY!LOUj4K+SepPuP)BN?3>4J+Bdymm|ZDPIl+_xk;6 zdXN=LZ$yHfJqqqEbwcJM-*rF6fTsj8BNgx-y+aJhxmvT<8#9ws6C&J_Ioo9S5IW}b z>#^_Ut4#9bis<;K!!uJtMhhK3hW5f&dXPtX1b=#kOMmQ&{%M6NHUXLVunTZP7+r{y zKO9f}HrdRsJVYfh;punH!A&Ld8uFvczKkMI7^Wi!yBj_??oMwmQIhY-T6dGKbkz6E zr}?CC_O7287%mst=uATYcim8oH!|bkSQ;iw*fof=c;Dxbr%-!4W9p5K6}9HSQT5@I z;G)U2i0bV)U!}RjH_nTgh+SYdjyDaLv(~&*mGQIg(zU9IcO5^do=wYq-ezA+-4*O3 z^V7S=7MZ_RUDJY~U!(dk4~x`A@1SIdHjtM(0{A07F}w@ZsO=qI?~!BH>9`g)dV}i_ z_Gr6vz_eoJ3GH^*g6cUTXuSLu`nii5uM^MZA<4DMD<+8wo3BZ+-OYOEpK_vJ*(H6y zV(Kv`M+nJ-yKRzIs$yC0lzE@~n)^|vB5dnDcMk)5Gt z%K$T%d$vsKloOjY)E{Vu-767bEp28*1HvKnnZd!N^WtW&Hw5VJ^h5onVrOs~uUx!2 z%_MA&Cn}DOlncA3C2~J;-W4O$f|B8+YYOR{^Nf95U4iAK_(yo`kRtZc$noAC3hge{ zlT;wK*@lhPsL?vcu}IXj|HqN*>mP?4MPQ8s9CH4+z8<;c{NKL*a_Vf{z#17}3ToXx z8YW5y3WT_?q43|`Hv)C=0O9Ix7y*5N2q<44)VXs7YBtmchnT}}&FS~w`6a;lrEB)( z97F#+>a$U`qF^uxD}=V%NcPn;5Vgs@YSKf|^txi5WW{>4+vW2!gQZr%i6``=Wi?ei zQb~1(RVp&?AHm8e!wP3XgFaq6lw<$i332ZWeckHk1uw;Y%@^l26mZAva)U~r%T0LR zHgO<*lfrCkmd|n}hA~2qR0<)I>)5i&#NFkpPP+0M)n_roF>t;<6a3T0TI+I~(aQOd zDMwivRcTsSOC5#-QZ){QG$$DfZICqCHl<*>m@rbVEd~gxnGVIm! zb=C}Dr|o@@CfsgeM?17GaD?x|mKJ=YOpBy$NU(EzId#^|^r@s9M*9`=R<(16WycBpJJ9Ip2Vf>)ERiwO5ILkxO%m|TXjWMTQ80~u@n|3sWQftK^Hh_nC93jc54 zkDdI$D>R1b{L)-h-eTc?7;qCQ&q{8thOAf~_QU4Nc=9EqaYIAo;}XK98P8ibss3Km5&{?n(_F6| z{7l`-K5}U<9Tzez_(3ebv5)VEQ7ea2;58$UM-8qej4jw?ny**)2%96Bq-q$%bTC!U*D+7H+@b&Dr|lU)bF< zk1XPLA*KQWYo*<)QfAgPSl76cs~5SIkAbR+`Bjt2&ist&>NM%JzPdF%nJx#%uWKdq zY|h%$EiI#h`*sbbw2LFo%(gMj_hlIIE8QrG`MbACom;+pFJkO=ZN=7v;`tuM=S*Iz z0fAs7dOT(@=@yP(=ZY=;>UDpg*Gz#-LztVmv6PFh;wnCvOFJ(yVd17+#>rau9)dI; z!A)&Qg1_6EX@_rZB&4=K!DV*8oY>L&`RbD5Nd5%Jw}Bm6uO6O;!=8fzuRPz*H#LeR z%~flrknPo2Xe+K5AU9gC(jp0Cfn9N19aor9E%p8FyFEoDjjCQ8-K08)&Jv%7EGXoH zUB?Ntn4sr=@3g(N2!b`sY7AwcDtdEX5}{7D7k@tX4}UQqLSDL}nd`UnKHLyjVZ#CX z+O|%E+Q{ym)8gx!T$yT1Y~%$SX{ZZnn`kdd%aT6RHES~L{Fm` zZD*t>(P0$CVf!)`ojqM9)4A4wQ$!AG;=-Uf4@cNaDVt28S~A%fR%Z7LSVT~xHUJhF zvASIJM`bXA87%DsQPxt22Pp*P-JlaDB(}d*25Kvb-kzeunSPl*M zxZlAUSSj2jrvnPmz`ym>;o7sq^o5Vu|JE7EqE@xB_r1Aq7~nZ$-uDNkp5j7XJXI1@ zt+o~1$yXy4J(k1ZrSEkPFhZqOh^Al17ZUB|IVUzJ@*<)!;y(Ee4tHi8Xs~zVD z<+*$TPuPd=__t`}XS_d;hi8#?&T8Yu`s(99Pr8%|Z;0_4^ZcfeffW0uK75rb6xM&s z(6Z}kp!Q;XR>gT6O;SvYR|)o;WDSqeP!LRF^FH1 z{&OlEwS^_dwQm6tfGK-nj1wr$mU+(gYTW8{Wtbo$>yABrx;S*(W0+!it(tTiY*HQL zh?8L82ld||zw!n_9fS)kdzq}LHa7MXkNT##=4XW|?lG#t1TG~Lt> zfWcpIhZ)9DJk9P3VXI_?qoX>ndz!i4)+68Pl6fg9v?Q4KrZKs4zZ2Y>&V-|JWti+z zszO@PYHIFOc>%9jfs^c{=r!hxyN7c!6YO*JXsaCn%LHD0F0~ed?{L`!sR-bA^)Gwx zp9C*)xgK(w1e=Ajx3xGkFq!sO&?}fq)~zj?*50R~3tYab+a=V!7DWW5jOQtJm$}pV zVjXNcU*4$dsb`N4*^q;Cy;ulth|!u2XwGbLyYQ(IeQ7TL=Ej+n!( zVJ>^ku0@@N?|@F4_cYv!7qYqId!WxtPF5Q;c(XHDxkq!UYy(|}4F9#QPY{;+wRx|; zF4YnhImNthpONljtvB*gNk}EW9trn1Js7;`^9}{-|BJsVMjRnvdhslR#5UOME~ZaF z%@7swdpg_}gSl$?iO~epgU?wr0j5}laP`g#sLTFVAN(;zd&Lv~4+q9agJVkSR_V(n z49Hm2vwBL~DXx#TMAu&r%?zsAuN1IzOVgI_>ZAXt4^illVTVw;UPX}P9+5P@P6aC_ zm?{ZHFOU}jFMz$g56BL2Ok8k_b`|d|xQ|$Ljc!}}g5Sl5!U|Orp2MP}2LZz$n8H*2 zo988V2sGqO7zDVNXx$$m-N=R4(E%vQBwnKxe5foD>D$dCWje(Pp*OgB_K-yDFEPdCv4zPM=vRSNHZR%L$FUCPrg# zRUoh^z~#o%QOCp$zxCx7u3;*xI;Bs! z@2Xhs^>^VrMF_EoRgQXPsf%J}0`+Mh#IMEidTiSf=sb&rtj&pqQb2bAK8dhtruUXW zcsD)%#H=P7J^T7csQdUIPcLUMApX}xfxmPIZ&GsQ

cf!8YuF#OzZ6n^mQ zL%Uo8Smw$*=?W#pG2?9ag`2pN1VZlwq>}az>q!g+B(=lL3ingaiIp9>4Bs>4w@YLD zYI z#dDXiWgB&{mbfP>1~QIrAniJVoDmrFGtn0ixmw;s+LaUzrLvbj4&5EO_~L#%M-U;g z2r=lGa#$L+eOaY3aNiK==^hR=2Wwm)QYV}u5J=Mi9k8(C^jJ!Sw?xq9IRzq@z;B8t z;wg$;o)?iDsf+bD-fO;r54q}HtLy4NNvB*+^BEG^bTGVqEUvyN03l9VRtqCX*bN4n z6s|_ykf<7`jaVyROhWtEN-K5~O(0B829Ks_YX{TaEd5H~j;}^8@G?DtlbmyLRmm4) zvn-eQs4Kd27laZYcs3j$D>~}GXvRAY-L>Jk0dP!Tm8MhVpC?E6cO8VGT;432Q?*N} z#6)3YIQb+o_ggM`1WH+4qTvR`+QBu<)0DoJzPGrv14fr*6LP88-nHG*&Ah|N2-rAs z#tD5OqIUX7{{^6fkLevik*NRwE;s55b_K7 zeq%4NUsqIX)AMoV4_NTBIDITpIJ0bgPU3abGobRzltxLY;(T6`Jguto3$4so$rqlr^8MWS zRfU^*u%F0#C#Y&35QTTnzt8Hckndl+wppCD*{|yJ4yX9JzU~>Go<^W1^tvm8*Nj)-@*Q?hP?iX*@Bqkw4 z6#h$?_~d9Xkn=Bui35L1g(C-o?s(;yqf%ijkR2}+90;nN<0lStRe;<0{T{u#I>PJO zqAT8?Uh|n$RS&w2`zkKH^_@)m!Logx0;zMecF|()U?CeetLJ=&w2uP z703MMy-P8HRJa`MVl>jVC++|hQz6Enr%y1?l;$jPolD|%s=gYupw!ubutPT6210Oq zEf!Oupqv7pEEZmSmt`iWIGsJtjU$!pV&N3`H6J!{E$Ckx-bT`|VpkXCEq*?^tN*bo zm3zYVy7BTo<247C!3L?q{ofSWV*OTSg;3W#iNi+Qt5zHIDJtdO>l4En&SP;&MS+^- zp5Esxj7^ifu4qNPq*8?rxv0%>?6iXb=mZ~bwP}2;HVJP7#U6^?9BdS@5gl#cm5WUj zpx;CbtVzV7k-h__8x`qNN@BlydJ-oyJEnwjeHV4Nsgc9iY(&6hu@9SKG`JH`FzZ{J zDr|zL`m#mO+(qU?3YO~PC1bI=%-22*LiDQb=+1`NTKD`g4L;J^W1icb{*ltWDw*%f z@U>;uTUME>>B*Jg<>i~)^4Edr0QPFAa}rV}6(Q%09xtC!&Bz!q{0t(b>iIWwPW+~L zRa~7dPcBJ#l1#phENe{5836(5z?-dzf#*AI9x{@L26@(c&H5kBGTvWXMc;+|F;3tf1abgC~$R@yj+TZPcVA- zXqSgQE=ccy(%a|8$CoB!1=)9a&@GpCfNU`ZvXYJj1Ov}NHY|*>3v0=l7A-80sCB1Y zmMR_P0V^Sx3YV6Mc+Ugi=WnA=Y@;MRs)e*vrED^xw`t?@@GiAx4ZAPS|I}i!|6nEt zdqqsRSVBr-NZ7={Sn<5cT)qV!#=gvhblT*ZQ&d$jAH4SC{a!8hIJay9{fTH^j9+QEzY`r&-M1^7`sZ{ZHDV(9eL}c zXt=*8zoLwZs#mEsY1UcduFvlfph@R4O(rSz$+>!33`)7ZpOvmLzZX+tV)w*{$KPIB zYo{cCRx%pfaHaLqVQRB~^p9gN}L{zy~P6Ua7Dk;X*R z>iKR`lXdX3;pcTq;~3YgvhcZ|okSC&J!nQGhvi~5$#g)EcYUm5R_T=$JUBvw0n8{Q zG<#D?j|uKrDJG_4%Z8Jq&Emx_F<+O?Cm!h#(1Z)FC9j>qeO*8hx|d_2PG`yOx6UC= zIVl-l=wcMS^p-gPuwGJl-VY(aigbW>Wt`V_{!>M`9nlTn{|ma|Up5BIf5=GrSdf#@ z3T=6Cl4u+FI6hOqt@J!OwS5wy3J}n1V{@H1+F`8gl$-7<+75Zk`*ysvLkW5VUDdwKCte?DnyV zJQ5-33eUX3BwpTsi?D{HjgmBb-K5rry$Z7L8w_-lc=z@rNMd2QT@+9B+@p%Y`%`+a zay(Yd6esit5}mV#>71s7`{iVBNHnf4YZ`3^mZ9{_%QG<&n(>&kWvh80&eDP$8kgX) zo(94?wv2UVxmLz@)3PqEMg82J<4)vqp*dWpg+Z0f}Y~uqNpmo@ME-wZ2i}{QP=Ly{w;eCLvans==Qo z)*+2}-qGi^c_sVZL9Lkr&#h{;h?57iMTAU6ECxqSWZ;|%u?B6L2X_gIF&mdh0uf}4h zG+sMxUf#PsOz{pC$m@olD)chvkM6MHb)=c+YOG(9+g?e}eN;=Ib12+HT4g{%LbIoss zPCR_C9@s{Ng%*m#g4JnOoHtb+(n&Huw#Q{}9)uPbzdP;k9&^HO!u6{D5XKQX6S|8V zM>d-yXvW4%&WaU#7&|DM^m_A(wax0vezjb|LTfHM1GZrhU>oLwQwaN_IT-p6q{&BT z8}gT*!Cy*6_gbJ-Y^|j+Q)dYPH3fyoTi(G6S6>6n00>hm_o{VE z@boV_Ws%(d<{|B5+9~r+E4R{6dpd_J{_TT9=dYFeJdZfgZFQE*#ibTk%#wF$AqeAmS&UKLiVV;p zJQ)1z3D2JG% zfJvlvRPfV)atG9~T%keXuY85demr%B!!qXAo$;gGb35YnSruNWgv}6U&?@M=80j_U zzEp&T?&8`8j#*pB`SX0ztiy9-cGoE|p^mStenPWZy7L;f4u}T~vq6_nd9zB*ADVC2 zBU=$UOQ2}@y4~_`iUv2@(+WWVXy^$x-${#RAGw^IYu3&IyqxU6oW8Pl6}s^bppQDJ z9jT+rBAsA9n{nkfE~5c8BG3NtNndxDgV?i9iYM}gKgGR?dDS3=+0UfpEq9p-tw^;6hxL!ROEf2kSaIMme){r`c3f_ zpV`U&zDIy8I8>o^a=NggI7{?Ws68k~1x|UcZP;`^iFD8V9RIYC>ox@!1x0DPmp+Fh z8|18eWlBY3(iiK3LYH6N-?QvrFiPh@G)iyd+1hlKd$c~8zxO_0eUsGiqJWY&PDhHW zF~GzmTkVm$Wb{U_V699ePpxyf)~nOKKwNqw_&0^b;`#Y)PJ;l?u+s4!!i6t$WV&hG ztxHGHUr@=<8b4_{v;ix7l=PA59;p>KQjH#9A*?gXTxlvRQIk0{a1Vre$ zaG9>Hi!C9jT~BmN%{M9ikX!h9|De({H*)zd`mbppLv;+pCb$cuki2An&PBda5!e?>gHn+_2M9%iqZnx~&8*pNhV zJ9`O9^3SEx{R`yD>1VgtJ0Cj}YzUUv2DWA(r(g`3b9M9?bDZL=#-zu%ZO&GoRk#$K zNZHBjKPSIOVo23IxwfDAaRFSu7gb!QV%`Z0Z)27atby8u>tl;+LRc{|LHDwygEz|r zS|8Pksy?gF^j;Sg+pJE_2`S3&BFb?(w{T>XckxR12lLUJTi^I?#XUOcFP(i%0rb?E z8#p_IKDp@-|Ha6)ugaYfDyE)=(?v2q<{R6qt0<{-WRZzk6BR^=oog5aEaaKu%5s;G zwA;8@#lftbqm*LK9L1uE(Tl*I<}5jrHVwK;)cdOL#tFlcH_ z@H=N$m_nvcry;dLyr%vuhZ37ke)mvYy}g+cE_mVG936a%Q*NQ%3-MuoDpF{xuNAIbKQRPVzyO`^ z)6+dT&6X7B6M<`zL$~LW8mN#KL}&6zlNm=a(ItrQwa|Rjh#GF~TUGeT+^3V)KyQ#7 z&V^W)VY7my2d~;p;;33Igu=nP4pK7V)z;qCA1iPr=m5`4q#`coM~Ro1dnxLZGLzRK zI2~x{@D#^kKADPib$G3f0~5J(7+~3WC#Fttph8#ABqnS>I7(g%=dq{En7r_iq^(!Q z%*NJcA;$f1#n_b$e99*6UkU;=*vvb-fYVN z3FXCGmCpoc1PrX2qO!bZ7zelWpH^H!F@sLXs8RNsFw_MQ(1)%pJa27KS5r8LyTL$@_wb0 zUNTBXXTFl8-)ULDG9_RA5u4r5cbU->Yz%qbS)I2-Cy-|h&W+=RJ%<6GWbAoT=13vw zEXS5Ux=r>3Flp>(Xr;cXW(F+p+G$oi`#wRrCMKB9%f}|UX6|;m@(VD3-Ece-l z>+54&-dln-@gk&4N!z9R6%`@+AEUZ^siRoUdi{^j*6TD~Lb|O)KA1T#R~P*yY!TX+ zYHB=}cYQ}XRT|At_mwQE!cGgPA^$Rq4f^GHus z_V2oyvkDnIBu*za-o8LbA~cYXz$y~i8IA3Gxp{e!EI4K*;_JB?mrZ*6e35FeeleSj z&&8gH+<^V8AzltNzU)6nc((g+Q{~dzYuNGAv?DOSfIO;_Cl_KeU&qAhiV_O%J+72F zNm3y0=-}0o!ajPx8uoJa&7TZH^6zl&KEn&fAY=)8o-Ks>^iy8SdT*hxkELO)AAadX zS9C^@T%VlzMPhe27kuA0+t#NI`<4<%kXNqS(e|daOaTX&Fb(|i%d>gl$JtK7% zrc2BagE?2z`Hj8w72A*RXwaEKw$tFytB7GF zXfbtD%t0~A++O8&=0`oT4cGEDN8D|p=nPg!xapqjR?QCEvqt_U9wq`FpWi;#dc3^9 z*sHE9qT{B|P~5ke_)EbnG=VY0emGUhKpC9i4w5AsWr;ssJ1nm*`9RFGfLTQEG0lrm zemZH2+s8bT2&1!7Q8Mt-Ip3(jV*;90c(v?2`oRa{939_O_G`YM*o=lixB41!c>T_3 zd))V|a~_~9>9<4gGI%|#`Q{`UDtEY&p0cft$TQQUi&ejCRSV#;py_2hM1R|jc68Zd zi_|lVd-1w+x>D->B{2eGnt9{nni52}d%>^lWRUQRrxYNx9PxE*GP!519B0^Z>cMJi zHU-diJ7MNv-%NkS;(eyulo7$RuGi+^D*V`S)pKndQ|J5=>}wvYtyJ|~kd=ks>%4WC z4nR|Xk`}6C8fO3t@wYOhA*zzgA752!;*dB%{@N)KlPl)>X3IJL`paPHW?6G2R=EF* zs9oZywLgeY+Vuh8+$vQ0pHlgtf1)ejXQ!~NFi;?6Q(!nv%vgU5w34lDBM!zWuTNv+ z2>h8<$&FQed%jiGGR6J&tfS)jRs~eN502@bT!oz`PoQzvqzEExYwm9v!{4qT8Rg)q z@Z7JdLwRe&YK8&F^GmAATYZ>{oAI2`pv5DMIdJX7UtXMn00C4!dAyqE7#QC5xgS*~^|b(jKpVQ{-J4+yoeSHpW7*8=XuMQ<6BZ1Je-F zwR`U;ES9%yzNi#Bh;~(;dRbPoE_>P5;n~%nn2)%cJ0c#LoGdfaGBEknMaA^)w@$OE zu3}OTcij*erUa1a|A#-9`-aRz{zjuN*sI16X%m^z;LeRlol zirlA~0bW~4K}YFZqhHe9I)aLb0p6b%D=#c~#25vB%v&~&^6+tmij8lnS`{g^+YfSE z3G6=@M9mn@I*>R4pnqJC5bd3NSlW&}6{BoFHu*fr$oLRTXX@DT0`4+1m+OBE*<6po z-8G<$CF$bYJ|kCT-q~%EKbNk@JcYmIKT%Pv>@D1%tt>Mix^Ggq;y>_dSRvogJKJ(9 zxnynJj=LgXkBAI=YD? z#-I8X@=so8^-6MGu5-Yd<4|OVOmN$wK~jBj_4DC&B*SxQ_xt#j)a-}80*e{q4j)Pf z(1eSNq%#h?n8PHHbs}{N@dLvJ90gqUZ}ncl+^~Q9#aG~`PMx{LU~uY`;2cs_@w@C} zv4cpuOBz2iCLQk<(p4@*CM#}r@A(WHaZH%LHf}l8^G5nQTTIeZC8SyizelD0^P3hB zIgVJGkEo8N91HF>> xRange[1] - nodeTol),\n", + " 'bottom': np.flatnonzero(mesh.coords[:,1] < yRange[0] + nodeTol),\n", + " 'top': np.flatnonzero(mesh.coords[:,1] > yRange[1] - nodeTol)}\n", + "\n", + "# create a copy of the mesh that has the nodeSets in it\n", + "mesh = Mesh.mesh_with_nodesets(mesh, nodeSets)" + ] + }, + { + "cell_type": "markdown", + "id": "413f82f9", + "metadata": {}, + "source": [ + "Now we're going to register the essential boundary conditions so that the optimism solvers will know how to deal with them. This is done with an `EssentialBC` object. Each `EssentialBC` represents a boundary condition on one field component of a node set. As en example, let's create one to represent the $x$-component displacement boundary condition on the nodes of the left edge:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "ebedd68d", + "metadata": {}, + "outputs": [], + "source": [ + "ebcLeft = Mesh.EssentialBC(nodeSet='left', field=0)" + ] + }, + { + "cell_type": "markdown", + "id": "3ffe98a4", + "metadata": {}, + "source": [ + "This is one boundary condition; we have three essential boundary conditions in total to apply. The thing to do is to collect them all in a list. So the code looks like it would in a real application, we'll ignore the `ebcLeft` we just created and create all of the essential boundary conditions in one statement:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "70344154", + "metadata": {}, + "outputs": [], + "source": [ + "EBCs = [Mesh.EssentialBC(nodeSet='left', field=0),\n", + " Mesh.EssentialBC(nodeSet='bottom', field=1),\n", + " Mesh.EssentialBC(nodeSet='top', field = 1)]" + ] + }, + { + "cell_type": "markdown", + "id": "9f671f9e", + "metadata": {}, + "source": [ + "Next, we create a `DofManager` object. What is this for? It's a class to help us index into our nodal arrays, keeping track of which degrees of freedom have essential boundary conditions. It's purpose will become clearer when we use it later." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "a34991d4", + "metadata": {}, + "outputs": [], + "source": [ + "fieldShape = mesh.coords.shape\n", + "dofManager = Mesh.DofManager(mesh, fieldShape, EBCs)" + ] + }, + { + "cell_type": "markdown", + "id": "0c4b89c1", + "metadata": {}, + "source": [ + "The variable `fieldShape` tells `dofManager` what the array of the unknowns will look like. In solid mechanics, the unknown field is the displacement, which happen to have the same shape as the nodal coordinates, so the `mesh.coords` array is a convenient place to grab this information. You could also manually enter `(nNodes, 2)`, where `nNodes` is the total number of nodes in the mesh, and 2 is the number of components of displacement in plane strain.\n", + "\n", + "We use the vis tools again to check that we've done the essential boundary condition specification correctly. First, we'll use a little trick to turn the boundary conditions into a viewable nodal field. We'll use an attribute of the `DofManager` class called `isBc`. This is just an array of size `fieldShape` that indicates whether each dof has an essential boundary condition. Numpy can cast this to an array of integers (more precisely, the `int64` type in numpy) with values 0 or 1 which can be plotted in Paraview. The `dataType` argument is different now; for block ID, it was a scalar field, but for boundary conditions, we want it to be a vector field (one component for each displacement component)." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "3d92324c", + "metadata": {}, + "outputs": [], + "source": [ + "bcs = np.array(dofManager.isBc, dtype=np.int64)\n", + "writer.add_nodal_field(name='bcs', nodalData=bcs, fieldType=VTKWriter.VTKFieldType.VECTORS,\n", + " dataType=VTKWriter.VTKDataType.INT)\n", + "writer.write()" + ] + }, + { + "cell_type": "markdown", + "id": "1103ae8f", + "metadata": {}, + "source": [ + "Note that `writer` still refers to the same `VTKWriter` object as before. A VTKWriter object is always associated with the same filename, so when we add a new field and then call the `write()` method, it will overwrite the previous VTK file. Indeed, if you open `check_problem_steup.vtk`, you'll see that it now contains two output fields, \"block_id\" and \"bcs\".\n", + "\n", + "Contour plots of the $x$- and $y$-components of the \"bcs\" field are shown below:" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "ceeb6e96", + "metadata": {}, + "source": [ + "![contour plots of boundary condition fields](boundary_conditions.jpeg)" + ] + }, + { + "cell_type": "markdown", + "id": "f4ee2543", + "metadata": {}, + "source": [ + "The first plot shows all nodes with boundary conditions on the $x$ component of the displacement. We see that the left edge has a value of 1 (which means that the boundary condition flag is \"True\" there), and every other node has a value of 0, which means they are unrestrained. This is exactly what we want. The $y$ component plot also confirms that the top and bottom edges correctly have their boundary conditions. Of course, the top edge has an *inhomogeneous* boundary condition. We'll enforce the actual value of this boundary condition later." + ] + }, + { + "cell_type": "markdown", + "id": "f8cef770", + "metadata": {}, + "source": [ + "## Build the function space\n", + "The next step is to create a representation of the discrete function space to help us do things like interpolate our fields and do calculus on them. The first ingredient we need is a quadrature rule. In optimism, quadrature rules are specified by the highest degree polynomial that they can exactly integrate. A smart way to do this is to choose it in relation to $p$, the polynomial degree of our interpolation." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "68a56fe0", + "metadata": {}, + "outputs": [], + "source": [ + "p = mesh.masterElement.degree" + ] + }, + { + "cell_type": "markdown", + "id": "d0b6b351", + "metadata": {}, + "source": [ + "In the linearized theory of solid mechanics, the natural trial/test function space is $H^1$, because the operator contains products of gradients of the displacement. Since our interpolation is of degree $p$, the displacement gradient is of degree $p-1$, and the inner product of gradients is of degree $2(p-1)$. Thus, we choose this as our quadrature rule precision to avoid under-integrating our operator:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "cf915ada", + "metadata": {}, + "outputs": [], + "source": [ + "# create the function space\n", + "quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2*(p - 1))" + ] + }, + { + "cell_type": "markdown", + "id": "fbb9d8e4", + "metadata": {}, + "source": [ + "The benefit of specifying our quadrature rule in this way is that if we later decide to modify the mesh to use higher-order elements, the quadrature rule will be updated automatically. This helps keep us out of trouble with hourglassing problems. Note that our operator is *nonlinear*, so the quadrature rule won't integrate it exactly, but the accuracy requirement of the linearized theory is a good rule of thumb.\n", + "\n", + "With the quadrature rule (and the mesh), we can now construct the function space:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "6c36cd8b", + "metadata": {}, + "outputs": [], + "source": [ + "fs = FunctionSpace.construct_function_space(mesh, quadRule)" + ] + }, + { + "cell_type": "markdown", + "id": "f79a4bc8", + "metadata": {}, + "source": [ + "We'll see this object in operation later when we set up our energy function." + ] + }, + { + "cell_type": "markdown", + "id": "d44fd1a0", + "metadata": {}, + "source": [ + "## Material models\n", + "Next, we instantiate the material models for the slab. For the left side, we'll choose a Neo-Hookean material. Material models in optimism take their material parameters in the form of a dictionary. For Neo-Hookean, the required parameters are the elastic modulus and the Poisson's ratio (both taken about the undeformed state)." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "2c1f89f0", + "metadata": {}, + "outputs": [], + "source": [ + "# create the material model for the left side\n", + "props = {'elastic modulus': 5.0,\n", + " 'poisson ratio': 0.25}\n", + "leftMaterialModel = Neohookean.create_material_model_functions(props)" + ] + }, + { + "cell_type": "markdown", + "id": "66eab9e6", + "metadata": {}, + "source": [ + "TODO: The property names are not documented yet. For now, you can find them by inspecting the code in optimism/materials." + ] + }, + { + "cell_type": "markdown", + "id": "176065d4", + "metadata": {}, + "source": [ + "Now we'll instantiate the other material model for the right-hand side. We'll pick a $J_2$ plasticity model, which is a bit more interesting (and thus has more parameters)." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "24778eaf", + "metadata": {}, + "outputs": [], + "source": [ + "E = 10.0\n", + "nu = 0.25\n", + "Y0 = 0.01*E\n", + "H = E/100\n", + "props = {'elastic modulus': E,\n", + " 'poisson ratio': nu,\n", + " 'yield strength': Y0,\n", + " 'hardening model': 'linear',\n", + " 'hardening modulus': H}\n", + "rightMaterialModel = J2Plastic.create_material_model_functions(props)\n" + ] + }, + { + "cell_type": "markdown", + "id": "19011b92", + "metadata": {}, + "source": [ + "The meaning of the parameters is clear from the keys. There are several hardening models currently available, such as linear hardening, a version of power law hardening, and the Voce exponential saturation law. We'll keep it simple and just use linear hardening.\n", + "\n", + "For multi-material simulations, we must create a dictionary of each material that maps the name of each element block to the material model object for that block, like this:" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "aa2b0a4e", + "metadata": {}, + "outputs": [], + "source": [ + "materialModels = {'left layer': leftMaterialModel, 'right layer': rightMaterialModel}\n" + ] + }, + { + "cell_type": "markdown", + "id": "3f1723e8", + "metadata": {}, + "source": [ + "## Write the energy function to minimize" + ] + }, + { + "cell_type": "markdown", + "id": "94039b0f", + "metadata": {}, + "source": [ + "Numerical solutions to PDEs are obtained in optimism by minimizing an objective function, which may be thought of intuitively as en energy. In fact, for hyperelastic materials, the objective function *is* the total energy. For history-dependent materials, one can often write an incremental functional that is minimized by the displacement field that carries the body from one discrete time step to the next. We're going to write the function for our problem now.\n", + "\n", + "There is a tool in the `Mechanics` module that will do most of the work for us. Let's call it first and explain its output afterwards." + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "e20f5769", + "metadata": {}, + "outputs": [], + "source": [ + " \n", + "# mechanics functions\n", + "mechanicsFunctions = Mechanics.create_multi_block_mechanics_functions(\n", + " fs, mode2D=\"plane strain\", materialModels=materialModels)\n" + ] + }, + { + "cell_type": "markdown", + "id": "35578c4a", + "metadata": {}, + "source": [ + "The `create_multi_block_mechanics_functions` writes some functions for us to help write the problem in the right form. The most important part of the output is `mechanicsFunctions.compute_strain_energy(...)`, which writes all the loops over the blocks, elements, and quadrature points that you would need to compute the energy from the nodal displacements. (Now we finally see the `FunctionSpace` object `fs` in action).\n", + "To use this new function, we'll invoke it like this:\n", + "```python\n", + "mechanicsFunctions.compute_strain_energy(U, internalVariables)\n", + "```\n", + "where `U` is the array of the degrees of freedom (the nodal displacements), and `internalVariables` is an array of the internal variables at every quadrature point. It has shape `(ne, neqp, nint)`, with `ne` being the number of elements, `neqp` the number of quadrature points per element, and `nint` the number of internal variables for the material model. (NOTE: For multi-material simulations, this is currently sized such that `nint` is the *maximum* number of internal variables over all materials, so that the array is padded for materials using fewer internal variables. We may remove this restriction in the future to improve memory efficiency).\n", + "\n", + "All of the minimization solvers in optimism require the objective function to have a signature like this:\n", + "```python\n", + "energy_function(Uu, p)\n", + "```\n", + "where `Uu` is the array of unknowns, and `p` is called the parameter set. The parameter set essentially holds all of the information that is needed to specify the problem but *isn't* the set of unknowns. These are things like values of the boundary conditions and the internal variables, as well as some other categories. Parameter sets are constructed by calling the `Params` function in the `Objective` module. This is to help organize them in certain categories that the solver needs to be aware of, such as which ones have derivatives taken with respect to inside the solver. We need to cast our energy function in this form. Let's write it like this and work out what the intervening steps must be:\n", + "\n", + "```python\n", + "def energy_function(Uu, p):\n", + " U = create_field(Uu, p)\n", + " internalVariables = p.state_data\n", + " return mechanicsFunctions.compute_strain_energy(U, internalVariables)\n", + "```\n", + "On the first line, we use the parameter set to pass from the array of unknowns to the full array of nodal displacements. This means we need to fill in the values of the inhomogeneous boundary conditions. Next, we pull out the internal variables from the parameter set. Finally, we use the canned `compute_strain_energy(...)` function with these variables in order to compute the total energy.\n", + "\n", + "The inhomogeneous boundary condition part is handled like so:" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "id": "217248f1", + "metadata": {}, + "outputs": [], + "source": [ + "# helper function to fill in nodal values of essential boundary conditions\n", + "def get_ubcs(p):\n", + " appliedDisp = p[0]\n", + " EbcIndex = (mesh.nodeSets['top'], 1)\n", + " V = np.zeros(fieldShape).at[EbcIndex].set(appliedDisp)\n", + " return dofManager.get_bc_values(V)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "c578d78c", + "metadata": {}, + "source": [ + "We will store the applied displacement in the first slot of the parameter set. In line 4 above we extract it. Then we make an array of the same size as the nodal displacements, set it to zero, and replace the values in the DOFs on the top edge with the value of the applied displacement.\n", + "\n", + "Now we can write the `create_field` function shown above in the proposed `energy_function`:" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "id": "2755c6db", + "metadata": {}, + "outputs": [], + "source": [ + "# helper function to go from unknowns to full DoF array\n", + "def create_field(Uu, p):\n", + " return dofManager.create_field(Uu, get_ubcs(p))\n", + " \n" + ] + }, + { + "cell_type": "markdown", + "id": "34511c5d", + "metadata": {}, + "source": [ + "The pieces are finally in place to define the energy function for our problem:" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "id": "be482478", + "metadata": {}, + "outputs": [], + "source": [ + "# write the energy to minimize\n", + "def energy_function(Uu, p):\n", + " U = create_field(Uu, p)\n", + " internalVariables = p.state_data\n", + " return mechanicsFunctions.compute_strain_energy(U, internalVariables)\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "5fab41f4", + "metadata": {}, + "source": [ + "## Set up the optimization solver\n", + "We have an objective function - `energy_function`, which we will hand to a routine that will find the unknowns displacements that minimize it. In this section, we specficy which optimization solver we want to use, and tell it how to work. \n", + "\n", + "We will use the Steighaug trust region method. This method uses linear conjugate gradient iterations as part of its algorithm, which in turn need to be preconditioned in order to be effective. Currently, the only available preconditioner in optimism is a Cholesky factorization on the stiffness matrix. We need to intruct the solver how to assemble the stiffness matrix like this:" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "id": "cd2bb5dd", + "metadata": {}, + "outputs": [], + "source": [ + "# Tell objective how to assemble preconditioner matrix\n", + "def assemble_sparse_preconditioner_matrix(Uu, p):\n", + " U = create_field(Uu, p)\n", + " internalVariables = p.state_data\n", + " elementStiffnesses = mechanicsFunctions.compute_element_stiffnesses(U, internalVariables)\n", + " return SparseMatrixAssembler.assemble_sparse_stiffness_matrix(\n", + " elementStiffnesses, mesh.conns, dofManager)\n" + ] + }, + { + "cell_type": "markdown", + "id": "2b5da4b6", + "metadata": {}, + "source": [ + "We see once again that `mechanicsFunctions` provides a helper function. In this case, the function `compute_element_stiffnesses` takes the same inputs as the energy function, but instead of returning the total energy, it returns an array containing the stiffness matrix for each element. The elemental stiffness matrices are the Hessians of the total energy in each element, and automatic differentiation is again used to perform this calculation. The `assemble_sparse_stiffness_matrix(...)` function takes these elemental stiffness matrices and contructs the system's stiffness matrix using a sparse matrix data structure (from SciPy). We tell the solver how to use this capability by building something called a `PrecondStrategy`: " + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "id": "b1bc05ed", + "metadata": {}, + "outputs": [], + "source": [ + "precondStrategy = Objective.PrecondStrategy(assemble_sparse_preconditioner_matrix)\n", + "\n", + "# solver settings\n", + "solverSettings = EquationSolver.get_settings(use_preconditioned_inner_product_for_cg=True)\n" + ] + }, + { + "cell_type": "markdown", + "id": "d4b55505", + "metadata": {}, + "source": [ + "Finally, we can specify custom settings for the solver is we wish." + ] + }, + { + "cell_type": "markdown", + "id": "3e10766a", + "metadata": {}, + "source": [ + "## Solve!" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "id": "28b893dd", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Assembling preconditioner 0\n", + "min, max diagonal stiffness = 0.9999999999999998 1.0000000000000002\n", + "Factorizing preconditioner\n", + "num warm start cg iters = 1\n", + "Assembling preconditioner 0\n", + "min, max diagonal stiffness = 0.9875365254464693 1.0110599625088337\n", + "Factorizing preconditioner\n", + "\n", + "Initial objective, residual = 0.04170630570636782 1.5580526760827706e-05\n", + "obj= -1.47133479e-10 , model obj= -1.47130159e-10 , res= 5.10352716e-10 , model res= 7.43900588e-21 , cg iters= 1 , tr size= 2.0 , interior , accepted= True\n", + "\n" + ] + } + ], + "source": [ + "# initialize unknown displacements to zero\n", + "Uu = dofManager.get_unknown_values(np.zeros(fieldShape))\n", + "\n", + "# set initial values of parameters\n", + "appliedDisp = 0.0\n", + "state = mechanicsFunctions.compute_initial_state()\n", + "p = Objective.Params(appliedDisp, state)\n", + " \n", + "# Construct an objective object for the equation solver to work on\n", + "objective = Objective.ScaledObjective(energy_function, Uu, p, precondStrategy=precondStrategy)\n", + " \n", + "# increment the applied displacement\n", + "appliedDisp = L*0.01\n", + "p = Objective.param_index_update(p, 0, appliedDisp)\n", + "\n", + "# Find unknown displacements by minimizing the objective\n", + "Uu = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "2b1968f3", + "metadata": {}, + "source": [ + "## Post-process\n", + "The last step is to write out the simulation results so that we can use them. We've already seen a few examples of generating VTK output. Let's create another `VTKWriter` for the simulation results. Adding the displacement field is straightforward:" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "e14a722a", + "metadata": {}, + "outputs": [], + "source": [ + "# write solution data to VTK file for post-processing\n", + "writer = VTKWriter.VTKWriter(mesh, baseFileName='uniaxial_two_material')\n", + "\n", + "U = create_field(Uu, p)\n", + "writer.add_nodal_field(name='displacement', nodalData=U, fieldType=VTKWriter.VTKFieldType.VECTORS)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "da5415dd", + "metadata": {}, + "source": [ + "Let's also show an example of plotting quadrature field data. A commonly needed output is the stress field. We make use of another function in `mechanicsFunctions` to help us. However, before we write out any quadrature field data, we should update the internal variables. Currently, `state` still holds the initial conditions of the internal variables. That is, the solver finds the equilibrium displacement field, but doesn't change `state` in place. This is again due to Jax's functional programming paradigm. The following function call returns the updated internal variables using the new displacement field and the old internal variables as inputs:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e46298c0", + "metadata": {}, + "outputs": [], + "source": [ + "# update the state variables in the new equilibrium configuration\n", + "state = mechanicsFunctions.compute_updated_internal_variables(U, p.state_data)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "d23fc4a2", + "metadata": {}, + "source": [ + "We make use of another of the `mechanicsFunctions` member functions to get the stress field, using the updated internal variables. Note that we don't pay the cost of iteratively solving the history-dependent material model response again. This function expects that the internal variables are already updated when it is called, and simply uses theses values in the evaluation of the material model energy density functions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d3b3e743", + "metadata": {}, + "outputs": [], + "source": [ + "energyDensities, stresses = mechanicsFunctions.compute_output_energy_densities_and_stresses(\n", + " U, state)\n" + ] + }, + { + "cell_type": "markdown", + "id": "b4221fcc", + "metadata": {}, + "source": [ + "As the left-hand side of this assignment implies, the scalar energy density field (at every quadrature point) is returned in addition to the stresses. In fact, only the scalar-valued energy density function for each material is implemented in each material model. The stress field is the derivative of the energy density function with respect to the displacement gradient, and optimism uses Jax's automatic differentiation to compute this derivative. When computing this derivative, evaluating the energy density at the same time is essentially free, so the above function coputes them both.\n", + "\n", + "The `stresses` array contains the stress tensor for every quadrature point of every element. There is no interpolation associated with a quadrature field, so in order to visualize the stresses, this must be remedied, either on the physics solver side or in the visualization software. One simple way to accomplish this is to compute element averages of the field and plot it as a cell-based field (instead of a nodal field). Optimism provides a method to do this:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0d3a7bb", + "metadata": {}, + "outputs": [], + "source": [ + "cellStresses = FunctionSpace.project_quadrature_field_to_element_field(fs, stresses)\n", + "writer.add_cell_field(name='stress', cellData=cellStresses,\n", + " fieldType=VTKWriter.VTKFieldType.TENSORS)\n", + "\n", + "writer.write()\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "2558d9ca", + "metadata": {}, + "source": [ + "The above code can be wrapped in a loop that will compute the response at multiple time steps. Note that before moving to the next time step, you should update the internal variables in the parameter set:" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "id": "7e0acf41", + "metadata": {}, + "outputs": [], + "source": [ + "# update the state variables in the new equilibrium configuration\n", + "p = Objective.param_index_update(p, 1, state)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 84a6fb53999d205bced5757f7c26a6bcdefc80f1 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 22 Jul 2022 10:17:30 -0700 Subject: [PATCH 4/4] Refine the tutorial --- examples/multi_material_tutorial/diagram2.png | Bin 0 -> 14028 bytes .../optimism_tutorial.ipynb | 360 ++++++++++-------- 2 files changed, 192 insertions(+), 168 deletions(-) create mode 100644 examples/multi_material_tutorial/diagram2.png diff --git a/examples/multi_material_tutorial/diagram2.png b/examples/multi_material_tutorial/diagram2.png new file mode 100644 index 0000000000000000000000000000000000000000..9cbafb49a5e9bd01412fee45e5714ae85f8353f6 GIT binary patch literal 14028 zcmeHuWmH^EwKp@QF z{l0tW&RTc=&diT#R*wg*f`qf6YNyfDrbmzt6 z@n~rT66BW`v7UN6bCD?W>7zC16@+1WgfMgHdJ`LQT+DJ`d;jM6 zV&Q5-E8axo%~OIkE%V&>(A=jbD0u zxB6x#BbT715+&L8Mb|mTi&(C)-74HKB`!lSauXezI-@>rl&%=%j?*0J%-RsMQY32Y zM1at?pWiLLka@x5V`QxjO-12;tDr(i+&*?OUv5Jr;0=~NO6&GL4#RKyGwtDVLYn7+ zc+dOklN_pOae3y-=TCtk3MwSeBXuEt<)iXn`#wW|F!}OOVo%sRZ1zVO0)n68d<)*b zv%@*1MR}Ti(Zclp2USdI1lK(VIfrYX9=hxdFb;tZ#wv~)8c3Xg4MsvkCPG35EM(vh z37HBB?H?NnNfnv;KW$xP_J5Z_K|+diK?41|j45z@{J?=P;Qj9-N*T(3miSnP`X4s3 z>SIy9NN+sgfQf*c_#z=;6F+{DkqU|^kdRQ&UG$Cpj5XB7?YumoHuhe&4$wdk#G@CI zWS}@;dN}ymKmtA7J$=Okr5OJyAr9D&r+FD6{}l0alVUX1(1yr+`8Ys?p?pw2MrkYv z1S08U?EsAOITh1nTAE#LF)xCdSJrz$+lY1C-$Ledg(B z6UgJ~%k-~G{-YiR2VXlM7lfaSmnY=0UK?94e?KWk#>a;K^YgEM`nfp%&z3xW|9xA) z4e~x-;pK<&@&0#i4uLNJ8@II4nJRCfI|M3wM;FIM2Z!i6~t|tG*l~444as4ls{_QHs`}m~(%ai-p zu>5lt7#C?QN#6e)V`;1*{}WFnBmpNC1zG(-4;%c39(KsP^|kTN3z1lQWRk$tcBs>;V39?@!m?Kc3eJ$ z2J8z9Y4~f2=vi1;_-+<82VEUCZ3bI?`SUZ_Ad`si#WDvc?J@wszJ9J>t%yU2djgFgN}`jO^>Tp zP>^`&!~O5&HvwL~Uvtd5-P3fZz6*tSN@r6ePwJ0FMxxMNFW~s4VMl)&d@kev{bM>^ z*tJd2=`hCO<>8VJAzF0gWxd6lr(b_cdGsF;*0zE&$!v#2oabTCN^xGxzP#N~psXLz z#AK0ly_{eMV;T&kTg*+JQQb#hMZUg08t1Y%mHyqT!x6Xt>Vw($$@}y3v_C~K{KW7$ zLPVCeHvY^XOtOppca^9Zc-pJY7GmE9{`{=;8obXr(1(~jyHTi2R5b}Q1QEr+V!-Gs zL`kUL)&w^{-<$+)#JnmL@v6b26-NkblWr5RedrS_rz1g5Wk6VhP<2ivtV9BSn%^EZ zwA5-83cG(44Z37^@;@jJSp50rjY$y4@OyF#jRO|Y`y}f$zU%Rd8k6ePXiQ=X)V$B> za|5A&zKvk>?zR5jd|o=&w)cYt)(K4vCylsZ{MhS%{4IvuvP>nDb)n5}QgqOnr@4MNk{S*&rm2I+Ug5F2QgerkYy>i`G#I6B zgj{HN^8Sp;6ld&j{h3`@Xog-|Z_sd#N5CzR6ns?oGSGrsw}i-OE-^fpFg}P@y$Vym z*R<1hp-%6m^H@CL`s<1cMpZR6?9ZAE8?QY+f7>tgs&RZ*aMqk+&|pZoZ-?mG43lVy zvUx)a#*v2?i~9#Q#|f-(piIZZ^oh9jGu*%QX0P=8_}rz}VAxc_8sX1H0>0Tv;%!$2 zZw*}UU>_~l5!)!o6L1o_WD>eJ`F@p$RDg)H35c9WV{tFH|F!J^8Ar9!+{Tqu?bW#8kw>KM@fG@5icuXlS3I=!b654#in z5_H^jc|OCd=jg~v%~`<$0&`WN^^B$3mOH+s%;5SYDVF{?q>KbP*0y=ZgJDv*wBlK1 zc(yI&v*W-xOlY8}RIs4##3XSoa_iEh+Sabo!#Bg$gR+c?Wh&>P*K~JNv?4C|En2}UIAH_&O7U6@*dcj}R#95xf z#a$C5PY7nzh7L!=%MYyaj&wjI5v|~lRAM6&Pp@Os$kr0A-s6#1uW_JY%9fd9l6f?h z$-TWbsS>R4v7IY-Pn8Qf8x^Pg!Mce!OMI7$qwbI6qA}4HvI5(A3DM!kT!?JP!^gVI zKByRO?32742nHht&xaD5)s?~Mi)C8~RB3aAQ)pgyQ4!5A+vj$Q7Mv5E$sS_Lg&dPz zjszYf7fEbnpD)y!>iS;5VlY|v#Gjor|M=9&bzQIW3>Sm4%YBU~IrhtKV8OxHuZpv% z3bC(7jUvomIS-G;HnyFAQwR$`$BJJ}M!``b>zB9BWEa}}{!@lB8`M7Y+1;h4=#na{ zgz9E#!WYLUIiDFA(# zAmI9RTRo>(BH-M*W8jHnE}qui-Cf8?fjO`asqK4K-&!9gY&`U({WLRE$X$h(6Ynb? zI=5+PaMUb0)oW)?xYhbM#LfQhH>0R0p{$4ya8Tj)&qRn(i)lWKkoPmY7}UwE_rM|h z_(S(;J+^H*Z6LeU{pFJ1-W1^kY}tO=Y>reODbTRNFhxXht6_n{cOhBhjV^!naF^ZB z{@<2IOE2-Uqnux*ViDtLq3ieBjK*HhSp8jW_9PrVMSpURwDx{AM!Faarw;EcG3^|q zxb0n|fKs~H&e!}0E1j;|yz`{N>YWO4)@;HU^%iykJ{{fJMO(`)2r0E~jPD}I9pR&y|wcK00Et-p}7OdN4J48Lq zg!{{{!+*hrcbj#K#nke%v+m=eHLOuSDQZIcF?}zL1!^(biOc=BqSW64% zuFoldEt?h#$wZsA$lqt+ZMXGC6LWvSVpF^O_=eHMbpJfe?U#BLUJCJ$hcc}Y8F2Cqqm%C2)(dCQS?$@=)hC1qfp%n1M0<{o3!!$T!Kdcc*f$a3M6=`lf)9{rb3b+Y?d-@($QfW`-@5$n8qzYk@RzMaUZADa|mX|8zL)2FIu! zwehND9n73%4avMyu$2wEziN42jY}=eH+@QK=${B{XG|lAjL?TyVztV@(3`5JU(Tewn+0Ie0uFc{l#}q z>XrAd(hI%8uuv+=Kt44nTtJ81Uo!_sUhNnRGEcWYx8mM^S2+0qZcu9w@wh{sG2pNyozn+UK$8)8J8+LhG)EhQY zZ7NY!d2IBFSv*~Hph3SP%E0mE&$cgTKk8g1h?f|)55;~kA>d>zx_S`Cgg~&TqDAAe z*~UQmuy)|d>B-6>`=4xBtSBOD&Xp)Zeu7X{w%R*7@|jO}pYM-MOI^z{rV(C18LG$9 z#fs8!hzNY!fmgbHSe}J4P4_?%w4=!`^Q|K3x+pt)JQ|yY=NtSg%i`3&P0)iA_H;K`H2)aQ3h_qVHq z`~9(eGhUkm#qTZ(g*_?&0KCYQDk+CcCAei8DQHNpDlj|U+EibdpBZ6)*eC7y+Qx+V zu`>{q=DfoN=`=Lv28%bs_cwFTYk=%v>GtaQa3}&btFUn9*!o)Q6!HkEN91gk+qI-` zqQ(83y~rxtb}O4WUWPAmEsTuRikQm{2}PDm$qCek0qAp=9!}0#cS^s*B*rz($0qSc z^TUgQwVrHOLpQs>#eHj=U{q{lLtCeWtm+Lt&`u=Ggq9nsnDvgtK60g{wKdHP++!U^ z?pC&Pts{K3#>F*lO%u@z4p^G|ys++2y7bVHugHdt`z+f$If;%a3TCfE7$+`kYm!0^ zA)d-ePbcKYUF8?c(PzQhK086tcLvXn>C*DEQV6wQfV)-sbTG~e-c~3kfeS+}E-s!~ ze2n#k^Wfhifz6;+@Eno|crNklP2zgAyUwRLq6QW(YC!YK=x=onphU zSBPnsj+sUGxy74YQsFB$BxEBnYxEFaa_~ui?#8nN+M~4&r_+AkSJa$4ZY}*v#?d*? z@@J4Gs>&^*$c!AodpHqBNAj~LPU(Uu3Z}P-bk<|0ATT` zgb?&qLjcg@9bjLGl7KM}sX`k}+4n})vE6Kxo~jT@aiL%;4v}(IQ1ImNTD?z^0>D_N z+6W;yAv(tCd(C8|UoElEcAwj7a~a?!cl$~|l<*0ONWUgGiZ>et7#YJ3-BLLsDjikx z=UI$cIO3xQJ=F_0xM)#18PvQxEz|lp#8=?n5b7caA10H%txcap!Nl~A6 z0>T?YK$gtkMZ~rjxwcVHhSg)tvv(htAqOO4T?B}1YgBZcvItZRB57JAF2-dkQ3A}( z1LHSOdA=Y74eD8d*Mzx&3`I zXQ%+gG5nM82y@^Ji`2#JlRGTz#fZQ+0VjA=g3;#U2Poo48OX>+%oO1S8QV?vH>^90 z)oRapa5$MPC&tLDe0+Sa?W*8Wiig=;Mk5WrP$gOX>DX{L5GtRX0KDC;xhM8{o>k1x z67m(^ETB5o*W1-s{)&pPMsUtySEwdi(dOMsqDB6+>dSKLZjp~ zgfa{yBTakhbG_4y?FBn%Vp43>D3}a+;r137nS4cUn&_2r+>FQKxumpLY{V#-EP1rr z+GzBFTgZQ3Au$(a%;2%$s*!=rOQW(eaP$>Jx={D9*tXx!$e#fDpzM$@(8)AgRcH$i zVUw~V!71}_BYJA0pym=MN`ev)2h{}GO$QjaNPyt{GHp481eqELv%+Q0Y5)J+|Ji47 z36qY14WV|yABMX(0Fi@hSU04RN%t4B$NJ{TJDJ`T5h&WK~0N4_^xU?yDC_1VH0F zKD|oNES1&(fXeRivi=X($?#MLCPgXq);=;J9-b6v->IA-F+*K%z=@_+=oP8QT6+ut zBa~OIo=X{-tc)i9rg;?$a4&|F-Y>}XSd7Qp0;Di>egpB`1t3QBTP&qmsrKy=+N5!K zw5px1D?7*xXBI8BsaVnI;TQ!)5iWwvYFex~RkC4}e1V*X3jV~BbijM_Q%Q%?`%8@m zHLH>0xmHaZVQMdmRkBd+CxFMGB*tNu02r)Cj*d#ibCx+tkSsR)MWUe$pV? ze6if1UcVf`VvlSN1s=SSl*@3wOg8G;bcf}Mp=GTEo$le`;c}Okl`xGBCg_vAAhZsg zlbbu~F`@@TB^k8$fsoPRQCxo$l%$+HuwfXhv+NhY`nJFXmR+egs0%pR2;tG%+|sxj zOGZMVzo4X~ROe3n-g8h=38bhxkL=Q&idFQCBnLQFY2;mjP*O?AxymE2h_8*KyN`R0+-YEm)mlKo0}U!&vnt)IX06c09O~ei1a%tQO#g7n)e2l zt(d}ooh9#7_7|IsyUfh1KfAf_T=m0oHR>F>a#dGNwm_3j`kk({ks-bf;p4Ks>Gq{p zsI-!SoN>Z_FJ4%$S(XM*J7NnlaR2M`+5^B&b&NLxu{X=JDk>dK-(0^!JdtEO-jhxh zLn={SOYP-dU&yWB%}YTbMJtdF4QgctkT}lKxwg#mGi+qZ5WfAh@Sxvc0x&0-YGV-# zu3-MFBVf`RS0@PjzTA?N=dY6tx#n|9aN{nsedh$wr8MXj7GdM0%VWk*+Z*gDONwij zK`QN8mq=gt^s9HCAe*)Z()9~|B`pwcvd7i~a?GA3VXqDKA!1z$e4JVu&t0r7*AI=; z#DVBK_znc&HR#!4Cma2Ob`v;3C;h*#rAI1Bh4~5H7pqh*r4J`yIgBtsRDP#h zvYuN2F1)GB$jlri6Kq6_a9%kUdP!b-j~bd8O4$kbQg=$!oA!8traQ~H*H&=3jX}TA z{?J@GICw6(O6h1}%^HkKC>s3By>~>CWIn;O-W4M7Rsx`X-j&Lf76`db)Vsfmilx@a z#5*4vZKe?AJ{mRogB6xADI$*tbm{lcFGVz;eTcqVCrbP5fn{X*U@)lhA_SR_u>Y(k zec69(clqg7sGoS6JO}a<-b_rLJ}wWxMUu*Dx{q1P#Lu58w+gB!2YmlHyr_y|V;g7S z2S7xkth6VHs1Ce>X#Fv}7}GQ3c6-dZKMr5i=oEbSXUjdm^qvjOBXoIKBXKDOq6$IO<_-0IcFyU*^gRx$SvcPCTzfE6byDV+UeL7Dz#p%<0x!{^P|u(q1K z<|bW*Nt6lZpY{%nUmeCA%GhHHe$7+~y~$Z@CXGA%dDe9%Ez5ko-t)xktzs~&8EVlx zuOae9!NkO*c_*Q_4DZ@zRjSA=bUKD*82Mv3{2QBb@Fu_yXT*f)(<|W{cg~pkCMQ@M zS%vz`l}AJ%fkvS*@^N?zADYa1q>-!vMr*lFk8MCx!un|@wNyR^E;lXoD8W%4w^-aa zi>lLA{Y!IyJhb<_Z-I;#atggniLmGT)YT~TH^3D0Z!4ux8Px2G%Aw$$>ubUmO1O#a zFlJ-{Vt4nJ$7-O6{(CFPZX;wf^!M3>>lQi_csWFZoOt-V0D=iei&)UOq}^8Few6<2 z5Va})pxE&&iXTUvx`5Gij5tLP%%v{hJ>w-##qAZ z{DiXw^BRMKScN>}2^(B{nTVE>>L&`eU-$PH4U8;E8$F-Vk|x*EtrDb3WApIk4EVeo zA~T&fZmHIRkd;8TS+m>rsavkMzjFpcfF`Fi`7HhCLx@vhx^O-&bxXFK&gS&HX~t(v zow6s{sjx+O>v>@XKMC>)LmBA8veY>xBpRS9l#A;$&gcN_baIvOVM-a`972}!Kau?Z zLjP|)g1`*x@$qpIFen7nAzahkt`$6TH{gtGZj zx=DnMfYEoN{~|F~tpb2bt7YJp(=O4C#iSjNx2fwGG`fxzjqEM=T&@bBI2NYTnGSz+ zZ5n1AQfRn>7QZw~O`&z-4^}FtE&_tlwRwXfkL%~xJO(Y6#SOq>NGWXva0TiU$1Q^; z)9i;|8ar~T!TKQNRpz{G8v3wn23$JHyr5s-RUR`SI;o&}0HJBpKfmMeVv*{P{!J_K zw2IoZI|4{TC$`1e8sJ85pA;S$m5z+#t5vI)R(Ct(>mMdksP&#c75;i~u+;7_ljrxu zaTda0+#hsR8~W#$x&SkQ3M;20Bhr!pX@v1(wnZ)MZtc_Qyb!7KD&WMDlbhTA`AiOv z3lJ;(u2zgX?B$t&H(d=NAj&e?v`XX!xLoX9S|k31pk+PjpX06u9rThn`0FQ%c z!uD+jdZ`V_P1_qsD-B|n05P-DVpR%2z4yNSwMkRhzE6s>?Qt^cB$Hw52i9^ZY0dyI zg)@G(J9(2FyaS^c8)sW!Te`R<&U`!|WsBBIWAQU6Z-$Afv z&v%?s-yvV{Ca`%kOv3~y6c-m4^X$7ZRARfu**cJQqa-696V`fjB_;G>POwpvsgU%H z6a%LI7^gw~d^(fz^#F}RhK(k$%6F9=(*2IHXIX4jAkFrBR76u?>jnyf+WNylD(VxA zPa_I@q;>>kX=mTQtpMO^M%%erzA|k(O{B_=@3;yFB&IKY?LBovw zLfJ~^3Gg`yXqD*j*KVuOiJku68sy@~hQ7)_dRbpJL0`z^n6SrTAY(?DZr%X`P> zrVFTgn6dy9qjDUJPp=6)qFuE8lfj3QUq3!8i`s2J+a8Iz+@*iW<2#7#P~}tcVxc6f zWw>dTW{t(qOvFl>MYoY}?_ecSX>w8=vHKV@CJ*bbHi&2uTMZ z-rg?Qh$Zhu+Y7L9adX>d46{E*V;#p>7sIoBqR2=GC!u%nwIeCTA#ad56NM|D?a=U1 zZf2%D07!S-#NvhKT@|x}%bmEUdw+k^JanmZtW~2 zBVY51r-*^kJdjy3+P(*udhrmCmYt;o{XXy`qBro)o$AvwaIt8y(Zn7-2XFsN@E9tx z5hg2hWOO7|{0%u-u8m@}fp#+00+fzPrm+wZ-U)?&oUFf2|9kgRxmGoG<8qP~ZMjPo z>ua|97o!IjkMrmw}SRK$A=b z=$B=7lP?zbrZaE40?%wqO1-1`;@bBOxvNaSZIJ zGW0hw;ao&Gob*+TK*B@IOiq^6=8L*F(gQfHJLVnA-0Ttm`fNhnv^RoXcpWx*eYShZ zFW=0aiLaC1##nIm1g`#G5kC>O^Y|Jhyp6l!75&bUMSIaL4N%Q<<)^N1ftc*_zDNu# zJqcmlL*bBZF53&b)M0h80V-gv$@ddKvqic}Lx6Um!!RrYoJIJ)F$7#hpmq0A4lN2J ztXy--K;VE=2Yd1b?MwzX5rM95GIP1hl_?>d1YnKw1RJN|eJ zD1{2&!K3(K;=|x&DIp;u;U;`!1Tt3!s90hSw7TVH3klD%?=^8_!=Xhn?jMWOb92c- za#Jj>4EyE~@_q9(b4adgGsOJpjN*}*`owi?nrrdYr1i%CXxYXm90tK`O!x?mNg2Av z!zfC_esoJSuO~xi(AOe@Y+(&6ByKnDJ#HBQ7|$C6(nrc@i1H-KM)^LAL+_fF*E%j> zG}@fcIo!O0>{DHaN_dxHTL!~c0V&}kSmDia^PB58)aMQyTpM~XH645&5m40Gk1=94 z&DPF&Rj2`WQ*GTQFo)&%kdS;5-Ucjbh!9#1toQnJ6&{pBZIzWPl~O3U@}ste$nq|- z%`NHVOuOLVZY4f_IXJq^Ktg=!N5 zga|4?uUU$B4NDA111vP0D^3(mCP38z#3wDnRzTnB2rsu(9C87;C{3pN8x9mQz#`*Y zZUBp>0IZ=(BSf^WFsn&!`o6 z>LIJ+?{Pz8W1}@wM@Pq6X-U}ADYXLYD$&6AIb?COLTZXbxS(iuV2=&!9gLg61rPy; zcVzxv+kNWsWh)l*jm5&oJ_I(4)QI~Zjop_@aSlgXl$F>k!oac5#ke0f%|7XNx*8ag zo(%l)#H=X&d-+l@@N{d-Us6gOBZ;lIpo2(hxFV3< zFL!&}9Rbq1Se$9e-mLGDCN5D0EOc|<-A-=_~T<+RM=I6l}?Lgzvye~Jvmv`*2VV=dgn|K zPOBs`9%mqfUWXWC0Yv!wzd*1)@56DkArjeWGK?Ijwt2TH1cFi3h9Q9Zuyg)Aa#`yQ zfC@#kPu8QIN5uUR_kGa%8rzn8@qHRKcbqmU!}#A{U4#E*bS7OJL0ohF>sglU9|w=E z9iuNGlMFHZ(BFt0lGEaTP7xcvku|}68>;LM&w~Wc2qjK(WO!S0GLGd5u;`OSm_P2> z@*^Fl9vffhA;bKp7D;MXRES(+c&nD;)5>Jods9vIw5rO!fE%)PCZR_2vl<+yFQBEZ zhkRGz?Wys)vM$ML2VBEA9`yU?yEhW^(#M8%7aqLvWIWR$ihM4CPn8OWp|DNTJS&@U z0bF0zTop!JtcD-cC>6Jg<9@w?No&aS^8_oPEHysV?i((=ycjk4*4ajN)fx7c0j4x+ ztxv$Lco9IQ@Y#ltx=bVrY|S06E+30y<3lSQl81$b#c5v_-^vhd4j`reX&Z|UPFtE#mq%m*XX;zP{BK33?n75 zM{=tIY!#72XI|ibmnE4dbCZrw?fC+Tc$j9*k76bMxbwLhKdP?M+D+DRf}nu4PQKIM zbMl~5k-_aC9@BFUi@VVQzna!DZfq*WueztbsEZ0OXscW}q7UmR1^?`(eoGuST1;X3 zo60KnU{r6`MQ8SUUbY->UOw@r3|p>nKi|z7`v;N?8$wA&Nxl(+R5M@o)G-|W_+eR3 z%VrisJmB!X)3s5PxlYiaC3eH-w zh`r62kiJ)uE_WKsDaYXXPi0PC-GTB$xtX>%uhodDSXbM!cw90J6<=i>=R1{@ zW6+1Bg`GyxYaf4TdKQd>elZXl?DFfo*t%gF&Ijfy;%r%%XlX}yO5Y>>3Q)fl-y1nE zE9tLZt}FHQMNxnD&EYYdXwg>B7Z|Ti@*k~FI?%(=*bN0`=ti!kA3yHkdsp1)w(96c z2u>^f9Nz^0hOVxDhEhTy$GIrBDM$jCPupNa<$%o>u9xw{#_n_l;`06 zuUcjk#iMGM{y5koAvaCWr%pHJM>cmDw+;n^(YlhY3fr80I+{J?$=j%^HdQFIA|;T8 zfK}me?eAN$3Zxldo__XFZD;rXrS__dF^{&VRnAKp>PIjYZk5TZ@t%Pa%=o+V$1Gy6 zcPEjn79dg8q0`tm*wc_{Xzw@1+mq)H4H@c4^~?uV`B5g*SLAky4=McP0f&n*fchQ$ zVz2!dEj(iQ#pXtTFaW163NzQzv58!_h9YErR7}dhjt+U>Ub@`g9=9AWX~|SSszHCx zQlEb-n#Ib}x>sQs{yr76Iyl2beml!qR|b9CMIpB`4!?I(hR%L1BQs?oe~aZ0jsg+~ zoFm)KfzZR^fQ{{d%!)Q>?EcHu>6q{XOh*>-#+Wn!v~2ysF86Il%+CyA=$v|>!)Zmi zrnO%XOc$g^GvK$8?MT1gnkt&1+I-K;95ps-_Bsy`hh48*2`_J}$Oe4+42 z=V^x@&@{9t+Cf00pI`6sy?H81JdkaEaWIEZMC{wVzP>Kt3w0CQzMIloR^;Dq=%v8p zc^(~{WBqh5(v;$T8KGK0xHTdVco4I$Z!~rR>k22Jy0QhFJZ=^p3P)aUG(Hnly=+%?Rjrsq}9a0ZA^@ zLX}Ok2>cVpyMrig8A!yWsqH}cA~UP)MT&r$CJ}fFd2UC;+qKOi2B(6TlUoKxScYC1 z2=f(5KFfGPFZ14&iP-{3t2uaidD*}4i836X76KbJq;B>yC7n}gB_ki21VG=(Y4A;{2}hHZ20aj4y(}CNay4-_4`fh66le zrRz#P+F%ME=GC9ghmzG1m}6B`HsN9;0A~G0aYYam=dk`|77mmeRH?7CLmf?XG7mRe z?l$;&osHo-#~;0f=gx0Rf`4LCw6E0^N(Bo~@}o-Sy$_j48V&o^QHnN8Gb?Q2yp_SC zhBDDTe3mF zXk^Xv$ls@*JavTE2eR5&2NLh z3jXa0h!^6|i%G*~({Xk}@G(N7)^Xqi6Bw-9vs(uz=GF&jGy~6&+0N4F!e-0225yhI zq?w9^yg(z#4|f-H7om-X!wh&~wb=W9v3oSo8*Ay98-^p?-<@kwV76yW#SAtG&#Q5A ztJ6@iGPven9D zwKE;c(73Dvj{gIW5+cGg)M>~RhKIJPBRLrj2o>LE;zvv|@kabeACB#0`!1E0(ne9# z&Gx&=z}>+Tn-D|I>Ra4jmf|YXB>E3MTQG_qwws}I42(lewTcy$pvLN9n?G5i#?ohr zSVebSdedfL-kWU=Cs77Et8c^P;J_TvnOa`%qw3+)rK2BbZ!GO9*#l;e5+qQO8Ht}{xWTZMIdY2v# zhhDmYmZXl86}Efx{oBGBh^Vn|<+k2n-2dYk#Yat=W|)i(ciLBpr&1pU8r(y3>SQc<4k6%Eu*;M> zmz7GJqxR|>_SSJyae5keDzz>m@f)?3xq$XM0@%WLE;M3OZ**VoM3GlM&Sz$5MB>LGFLYCw*HpLG=y~fA29} zIg5mvuc+$3p>(tmdObP$+?6^ru$Q5B-OpE>_#wzgDd-dzMMajcjMjx{ttLcf<6EM literal 0 HcmV?d00001 diff --git a/examples/multi_material_tutorial/optimism_tutorial.ipynb b/examples/multi_material_tutorial/optimism_tutorial.ipynb index 54af4951..6f375ab2 100644 --- a/examples/multi_material_tutorial/optimism_tutorial.ipynb +++ b/examples/multi_material_tutorial/optimism_tutorial.ipynb @@ -2,16 +2,42 @@ "cells": [ { "cell_type": "markdown", - "id": "b5d338a3", + "id": "6356bf65", "metadata": {}, "source": [ - "# Tutorial: a multi-material simulation with optimism" + "# Tutorial: a multi-material problem" + ] + }, + { + "cell_type": "markdown", + "id": "3a2ebec7", + "metadata": {}, + "source": [ + "## Overview\n", + "\n", + "In this problem, we'll subject a rectangular slab to plane strain uniaxial tension. As shown in the figure, the slab is made of two equal-size layers made of different materials. One material will use a hyperelastic material model, while the other will use a $J_2$ plasticity model. The slab will be loaded by applying an extensional displacement $d$ to the top edge." + ] + }, + { + "cell_type": "markdown", + "id": "b60dcbd2", + "metadata": {}, + "source": [ + "![diagram](./diagram2.png)" + ] + }, + { + "cell_type": "markdown", + "id": "045899d7", + "metadata": {}, + "source": [ + "## Imports\n" ] }, { "cell_type": "code", "execution_count": 1, - "id": "24fd5985", + "id": "3f49e265", "metadata": {}, "outputs": [ { @@ -41,25 +67,7 @@ }, { "cell_type": "markdown", - "id": "9c7de108", - "metadata": {}, - "source": [ - "## Overview of the problem\n", - "\n", - "In this problem, we'll subject a rectangular slab to plane strain uniaxial tension. As shown in the figure, the slab is made of two equal-size layers made of different materials. One material will use a hyperelastic material model, while the other will use a $J_2$ plasticity model." - ] - }, - { - "cell_type": "markdown", - "id": "37c03013", - "metadata": {}, - "source": [ - "*Figure here*" - ] - }, - { - "cell_type": "markdown", - "id": "b1c62590", + "id": "15441e66", "metadata": {}, "source": [ "## Set up the mesh" @@ -67,23 +75,23 @@ }, { "cell_type": "markdown", - "id": "a8b6f7ba", + "id": "da642525", "metadata": {}, "source": [ - "As in any finite element code, first we need to set up the geometry of the problem with a mesh. Optimism can read in meshes from the Exodus II format, but it also provides utilities for generating structured meshes for simple problems like this one. Let's generate a rectangular mesh with width `w` and length `L`. We'll also specify that it has 5 nodes in the $x$-direction (hence 4 elements) and 15 nodes in the $y$-direction." + "As in any finite element code, first we need to discretize the geometry of the body with a mesh. Optimism can read in meshes from the Exodus II format, but it also provides utilities for generating structured meshes. The structured mesh route is often best for for simple problems like this one. Let's generate a rectangular mesh with width $w = 0.1$ and length $L = 1.0$. We'll also specify that it has 5 nodes in the $x$-direction (hence 4 elements) and 15 nodes in the $y$-direction." ] }, { "cell_type": "code", "execution_count": 2, - "id": "b7750f13", + "id": "4938c5c6", "metadata": {}, "outputs": [], "source": [ "# set up the mesh\n", "w = 0.1\n", "L = 1.0\n", - "nodesInX = 5 # must be odd\n", + "nodesInX = 5 # must be odd for element boundaries to contain material boundary\n", "nodesInY = 15\n", "xRange = [0.0, w]\n", "yRange = [0.0, L]\n", @@ -92,16 +100,16 @@ }, { "cell_type": "markdown", - "id": "f18bd7ea", + "id": "619b23e4", "metadata": {}, "source": [ - "The mesh object has an attribute called `blocks` which lets us group elements together to specify different materials on them. Let's see what our mesh has:" + "The mesh object has an attribute called `blocks` which groups elements together to specify different materials on them. Let's see what blocks our mesh has:" ] }, { "cell_type": "code", "execution_count": 3, - "id": "bb1bd192", + "id": "219257f7", "metadata": {}, "outputs": [ { @@ -127,7 +135,7 @@ }, { "cell_type": "markdown", - "id": "2bea405c", + "id": "5afe9f12", "metadata": {}, "source": [ "`blocks` is a standard Python dictionary object (i.e. a `dict`), where the keys are strings that lets us give meaningful names to the element blocks. The values are Jax-numpy arrays that contain the indices of the elements in the block. The `construct_structured_mesh` function returns a mesh with just one block. We want two equal-sized blocks for our problem like in the figure, so let's modify the mesh.\n", @@ -138,7 +146,7 @@ { "cell_type": "code", "execution_count": 4, - "id": "ec126372", + "id": "77813b01", "metadata": {}, "outputs": [], "source": [ @@ -148,7 +156,7 @@ }, { "cell_type": "markdown", - "id": "d88821e7", + "id": "450d807f", "metadata": {}, "source": [ "We'll loop over the element node connectivity table, which is the `mesh.conns` object, extract the coordinates of its vertices, and use the `triangle_centroid` object on them to determine if the element is in the left or right layer. We will store the results in a list called `blockIds`:" @@ -156,13 +164,14 @@ }, { "cell_type": "code", - "execution_count": 26, - "id": "a35d7cdd", + "execution_count": 33, + "id": "70c7dced", "metadata": {}, "outputs": [], "source": [ "# initialize the block IDs of all elements to a dummy value (-1)\n", "blockIds = -1*onp.ones(Mesh.num_elements(mesh), dtype=onp.int64)\n", + "\n", "for e, t in enumerate(mesh.conns):\n", " vertices = mesh.coords[t,:]\n", " if triangle_centroid(vertices)[0] < w/2:\n", @@ -175,7 +184,7 @@ }, { "cell_type": "markdown", - "id": "b87eed4b", + "id": "960c5019", "metadata": {}, "source": [ "This will mark an element as block 0 if it's centroid is on the left hand side of the slab, and as block 1 if it's on the other side. Now, let's make the `dict` that we want to attach to the mesh object." @@ -184,7 +193,7 @@ { "cell_type": "code", "execution_count": 6, - "id": "e245933d", + "id": "2f8ce025", "metadata": {}, "outputs": [], "source": [ @@ -194,16 +203,16 @@ }, { "cell_type": "markdown", - "id": "61954ea8", + "id": "24f443e0", "metadata": {}, "source": [ - "Now we can make use of a function that takes in the original mesh (with one block) and the block IDS we just created, and returns a new mesh that is the same as the old one except that the blocks have been updated." + "Now we can make use of a function that takes in the original mesh (with one block) and the block IDs we just created, and returns a new mesh that is the same as the old one except that the blocks have been updated." ] }, { "cell_type": "code", "execution_count": 7, - "id": "e2a1c548", + "id": "7b5c94f6", "metadata": {}, "outputs": [], "source": [ @@ -212,7 +221,7 @@ }, { "cell_type": "markdown", - "id": "bcb47738", + "id": "35f85f8a", "metadata": {}, "source": [ "Let's check to make sure this process worked. To do this, we'll make use of optimism's output facilities. Optimism provides a class called `VTKWriter`, that, as the name suggests, writes data to VTK files which can be visualized in ParaView (and several other visualization tools). First, we instantiate a VTKWriter object, giving it the base of the filename (the name that will be suffixed with the \".vtk\" extension)." @@ -221,7 +230,7 @@ { "cell_type": "code", "execution_count": 8, - "id": "ea13862b", + "id": "da768ca5", "metadata": {}, "outputs": [], "source": [ @@ -230,7 +239,7 @@ }, { "cell_type": "markdown", - "id": "619fe562", + "id": "067f340c", "metadata": {}, "source": [ "Next, we can start adding fields. In this case, we only have one field - the block ID numbers - which is a scalar field of integer type. Finally, once our data is loaded into the writer, we call the `write()` method to tell it to write the VTK file to disk." @@ -239,7 +248,7 @@ { "cell_type": "code", "execution_count": 9, - "id": "42104d3f", + "id": "19b50965", "metadata": {}, "outputs": [], "source": [ @@ -251,15 +260,15 @@ }, { "cell_type": "markdown", - "id": "607566b1", + "id": "cdee4b9d", "metadata": {}, "source": [ - "This is what we get when we open Paraview and visualize the `block_id` field. If you try it, you should see something similar." + "This is what we get when we open the file `check_problem_setup.vtk` in ParaView and visualize the `block_id` field. If you try it, you should see something similar." ] }, { "cell_type": "markdown", - "id": "50aa988b", + "id": "bdc4fb70", "metadata": {}, "source": [ "![paraview shwoing blocks in mesh](blocks_elements.png)" @@ -267,25 +276,25 @@ }, { "cell_type": "markdown", - "id": "f0637278", + "id": "da15ed2a", "metadata": {}, "source": [ - "Success! This is what were shooting for." + "The left layer is blue (element block 0) and the right layer is red (element block 1). Success! This is what were shooting for. (The particular colors may be different in your ParaView installation depending on the color palate chosen for contour plots, but the spatial layout of the blocks should be the same)." ] }, { "cell_type": "markdown", - "id": "bddeee2a", + "id": "6ddc1328", "metadata": {}, "source": [ "## Essential boundary conditions\n", - "We're going to make one more modification to the mesh. Looking again at the problem setup figure, we can see that we need to apply boundary conditions on the bottom, left, and top boundary edges of the slab. Similar to the `blocks` attribute, `mesh` has a `nodeSets` dictionary that maps names to sets of nodes. We'll make the nodesets we need by performing range queries over the nodal coordinates:" + "We're going to make one more modification to the mesh. Looking again at the problem setup figure, we can see that we need to apply essential boundary conditions on the bottom, left, and top boundary edges of the slab. We must create node sets that group the nodes on these edges together so that we can set the boundary condtions on them later. Similar to the `blocks` attribute, `mesh` has a `nodeSets` attribute that is a dictionary mapping names of node sets to the indices of the nodes in the set. We'll make the node sets we need by performing range queries over the nodal coordinates, storing the indices of the nodes on the left edge under the key \"left\", etc." ] }, { "cell_type": "code", "execution_count": 10, - "id": "4d67689c", + "id": "50df5e33", "metadata": {}, "outputs": [], "source": [ @@ -301,7 +310,7 @@ }, { "cell_type": "markdown", - "id": "413f82f9", + "id": "244059e7", "metadata": {}, "source": [ "Now we're going to register the essential boundary conditions so that the optimism solvers will know how to deal with them. This is done with an `EssentialBC` object. Each `EssentialBC` represents a boundary condition on one field component of a node set. As en example, let's create one to represent the $x$-component displacement boundary condition on the nodes of the left edge:" @@ -310,7 +319,7 @@ { "cell_type": "code", "execution_count": 11, - "id": "ebedd68d", + "id": "8df8bbf2", "metadata": {}, "outputs": [], "source": [ @@ -319,16 +328,16 @@ }, { "cell_type": "markdown", - "id": "3ffe98a4", + "id": "1b5f25c8", "metadata": {}, "source": [ - "This is one boundary condition; we have three essential boundary conditions in total to apply. The thing to do is to collect them all in a list. So the code looks like it would in a real application, we'll ignore the `ebcLeft` we just created and create all of the essential boundary conditions in one statement:" + "This is one boundary condition; we have three essential boundary conditions in total to apply. The thing to do is to collect them all in a list. So the code looks like it would in a real application, we'll ignore the `ebcLeft` object above and create all of the essential boundary conditions in one statement:" ] }, { "cell_type": "code", "execution_count": 12, - "id": "70344154", + "id": "a8f61e3f", "metadata": {}, "outputs": [], "source": [ @@ -339,16 +348,16 @@ }, { "cell_type": "markdown", - "id": "9f671f9e", + "id": "91ac4490", "metadata": {}, "source": [ - "Next, we create a `DofManager` object. What is this for? It's a class to help us index into our nodal arrays, keeping track of which degrees of freedom have essential boundary conditions. It's purpose will become clearer when we use it later." + "Next, we create a `DofManager` object. What is this for? It's a class to help us index into our nodal arrays, keeping track of which degrees of freedom have essential boundary conditions, and which are unrestrained (hence are unknowns to be solved for). Its usage will become clearer later." ] }, { "cell_type": "code", "execution_count": 13, - "id": "a34991d4", + "id": "01f306eb", "metadata": {}, "outputs": [], "source": [ @@ -358,18 +367,18 @@ }, { "cell_type": "markdown", - "id": "0c4b89c1", + "id": "d0291001", "metadata": {}, "source": [ - "The variable `fieldShape` tells `dofManager` what the array of the unknowns will look like. In solid mechanics, the unknown field is the displacement, which happen to have the same shape as the nodal coordinates, so the `mesh.coords` array is a convenient place to grab this information. You could also manually enter `(nNodes, 2)`, where `nNodes` is the total number of nodes in the mesh, and 2 is the number of components of displacement in plane strain.\n", + "The variable `fieldShape` tells `dofManager` what the shape of the array of DOFs will be. In solid mechanics, the degrees of freedom are nodal displacements, which have the same shape as the nodal coordinates `mesh.coords`. Thus, line 1 above is a convenient way to obtain this information. You could also manually enter `(nNodes, 2)` for the second argument of the `DofManager` constructor, where `nNodes` is the total number of nodes in the mesh, and 2 is the number of components of displacement in plane strain.\n", "\n", - "We use the vis tools again to check that we've done the essential boundary condition specification correctly. First, we'll use a little trick to turn the boundary conditions into a viewable nodal field. We'll use an attribute of the `DofManager` class called `isBc`. This is just an array of size `fieldShape` that indicates whether each dof has an essential boundary condition. Numpy can cast this to an array of integers (more precisely, the `int64` type in numpy) with values 0 or 1 which can be plotted in Paraview. The `dataType` argument is different now; for block ID, it was a scalar field, but for boundary conditions, we want it to be a vector field (one component for each displacement component)." + "We use the vis tools again to check that we've specified the essential boundary conditions correctly. First, we'll use a little trick to turn the boundary conditions into a viewable nodal field. We'll use an attribute of the `DofManager` class called `isBc`, which is a boolean array of shape `fieldShape` that indicates whether each DOF has an essential boundary condition. We cast this to an array of integers (more precisely, the `int64` type in numpy) which gives them values of 0 (for no essential boundary condition) or 1 (for an essential boundary condition) which can be plotted in ParaView. The `dataType` argument is different from the `add_nodal_field` call for the element block ID, since that was a *scalar* field. For boundary conditions, we want it to be a *vector* field (one component for each displacement component)." ] }, { "cell_type": "code", "execution_count": 14, - "id": "3d92324c", + "id": "6ea61275", "metadata": {}, "outputs": [], "source": [ @@ -381,7 +390,7 @@ }, { "cell_type": "markdown", - "id": "1103ae8f", + "id": "8a3cc28d", "metadata": {}, "source": [ "Note that `writer` still refers to the same `VTKWriter` object as before. A VTKWriter object is always associated with the same filename, so when we add a new field and then call the `write()` method, it will overwrite the previous VTK file. Indeed, if you open `check_problem_steup.vtk`, you'll see that it now contains two output fields, \"block_id\" and \"bcs\".\n", @@ -390,9 +399,8 @@ ] }, { - "attachments": {}, "cell_type": "markdown", - "id": "ceeb6e96", + "id": "93df5456", "metadata": {}, "source": [ "![contour plots of boundary condition fields](boundary_conditions.jpeg)" @@ -400,7 +408,7 @@ }, { "cell_type": "markdown", - "id": "f4ee2543", + "id": "fb692616", "metadata": {}, "source": [ "The first plot shows all nodes with boundary conditions on the $x$ component of the displacement. We see that the left edge has a value of 1 (which means that the boundary condition flag is \"True\" there), and every other node has a value of 0, which means they are unrestrained. This is exactly what we want. The $y$ component plot also confirms that the top and bottom edges correctly have their boundary conditions. Of course, the top edge has an *inhomogeneous* boundary condition. We'll enforce the actual value of this boundary condition later." @@ -408,17 +416,17 @@ }, { "cell_type": "markdown", - "id": "f8cef770", + "id": "c02948e6", "metadata": {}, "source": [ "## Build the function space\n", - "The next step is to create a representation of the discrete function space to help us do things like interpolate our fields and do calculus on them. The first ingredient we need is a quadrature rule. In optimism, quadrature rules are specified by the highest degree polynomial that they can exactly integrate. A smart way to do this is to choose it in relation to $p$, the polynomial degree of our interpolation." + "The next step is to create a representation of the discrete function space to help us do things like interpolate our fields and do calculus on them. The first ingredient we need is a quadrature rule. In optimism, quadrature rules are specified by the highest degree polynomial that they can exactly integrate. A smart way to do this is to choose it in relation to $p$, the polynomial degree of our interpolation, which we can obtain like this:" ] }, { "cell_type": "code", "execution_count": 15, - "id": "68a56fe0", + "id": "e6cc7a6c", "metadata": {}, "outputs": [], "source": [ @@ -427,16 +435,16 @@ }, { "cell_type": "markdown", - "id": "d0b6b351", + "id": "61bbbb95", "metadata": {}, "source": [ - "In the linearized theory of solid mechanics, the natural trial/test function space is $H^1$, because the operator contains products of gradients of the displacement. Since our interpolation is of degree $p$, the displacement gradient is of degree $p-1$, and the inner product of gradients is of degree $2(p-1)$. Thus, we choose this as our quadrature rule precision to avoid under-integrating our operator:" + "In the linearized theory of quasi-static solid mechanics, the natural trial/test function space is $H^1$, because the operator contains inner products of gradients of the displacement. Since our interpolation uses polynomials of degrees up to $p$, the displacement gradient is of degree $p-1$, and the inner product of gradients is of degree $2(p-1)$. Thus, we choose this as our quadrature rule precision to avoid under-integrating our operator:" ] }, { "cell_type": "code", "execution_count": 16, - "id": "cf915ada", + "id": "eb4e827a", "metadata": {}, "outputs": [], "source": [ @@ -446,18 +454,18 @@ }, { "cell_type": "markdown", - "id": "fbb9d8e4", + "id": "6b967ae4", "metadata": {}, "source": [ - "The benefit of specifying our quadrature rule in this way is that if we later decide to modify the mesh to use higher-order elements, the quadrature rule will be updated automatically. This helps keep us out of trouble with hourglassing problems. Note that our operator is *nonlinear*, so the quadrature rule won't integrate it exactly, but the accuracy requirement of the linearized theory is a good rule of thumb.\n", + "The benefit of specifying our quadrature rule in this way is that if we later decide to modify the mesh to use higher-order elements, the quadrature rule will be updated automatically. This helps keep us out of trouble with hourglassing problems. Note that our operator is *nonlinear*, so the quadrature rule won't integrate it exactly, but the accuracy requirement of the linearized theory is still a good rule of thumb.\n", "\n", - "With the quadrature rule (and the mesh), we can now construct the function space:" + "With the quadrature rule and the mesh, we can now construct the function space:" ] }, { "cell_type": "code", - "execution_count": 17, - "id": "6c36cd8b", + "execution_count": null, + "id": "fcc72db5", "metadata": {}, "outputs": [], "source": [ @@ -466,15 +474,15 @@ }, { "cell_type": "markdown", - "id": "f79a4bc8", + "id": "c570688e", "metadata": {}, "source": [ - "We'll see this object in operation later when we set up our energy function." + "The function space holds data like the values of shape functions and their gradients at all of the quadrature points in the mesh. We'll see this object again later when we set up our energy function." ] }, { "cell_type": "markdown", - "id": "d44fd1a0", + "id": "dffeee71", "metadata": {}, "source": [ "## Material models\n", @@ -483,8 +491,8 @@ }, { "cell_type": "code", - "execution_count": 18, - "id": "2c1f89f0", + "execution_count": null, + "id": "4f7d7828", "metadata": {}, "outputs": [], "source": [ @@ -496,15 +504,15 @@ }, { "cell_type": "markdown", - "id": "66eab9e6", + "id": "81d32097", "metadata": {}, "source": [ - "TODO: The property names are not documented yet. For now, you can find them by inspecting the code in optimism/materials." + "TODO: The property names are not documented yet. For now, you can find them by inspecting the code in `optimism/materials`." ] }, { "cell_type": "markdown", - "id": "176065d4", + "id": "c9a4acb7", "metadata": {}, "source": [ "Now we'll instantiate the other material model for the right-hand side. We'll pick a $J_2$ plasticity model, which is a bit more interesting (and thus has more parameters)." @@ -512,8 +520,8 @@ }, { "cell_type": "code", - "execution_count": 19, - "id": "24778eaf", + "execution_count": null, + "id": "d5d64680", "metadata": {}, "outputs": [], "source": [ @@ -531,7 +539,7 @@ }, { "cell_type": "markdown", - "id": "19011b92", + "id": "34971c24", "metadata": {}, "source": [ "The meaning of the parameters is clear from the keys. There are several hardening models currently available, such as linear hardening, a version of power law hardening, and the Voce exponential saturation law. We'll keep it simple and just use linear hardening.\n", @@ -541,8 +549,8 @@ }, { "cell_type": "code", - "execution_count": 20, - "id": "aa2b0a4e", + "execution_count": null, + "id": "e6040db7", "metadata": {}, "outputs": [], "source": [ @@ -551,7 +559,7 @@ }, { "cell_type": "markdown", - "id": "3f1723e8", + "id": "0eb8be55", "metadata": {}, "source": [ "## Write the energy function to minimize" @@ -559,7 +567,7 @@ }, { "cell_type": "markdown", - "id": "94039b0f", + "id": "47989532", "metadata": {}, "source": [ "Numerical solutions to PDEs are obtained in optimism by minimizing an objective function, which may be thought of intuitively as en energy. In fact, for hyperelastic materials, the objective function *is* the total energy. For history-dependent materials, one can often write an incremental functional that is minimized by the displacement field that carries the body from one discrete time step to the next. We're going to write the function for our problem now.\n", @@ -569,8 +577,8 @@ }, { "cell_type": "code", - "execution_count": 29, - "id": "e20f5769", + "execution_count": 34, + "id": "637b03bf", "metadata": {}, "outputs": [], "source": [ @@ -582,21 +590,23 @@ }, { "cell_type": "markdown", - "id": "35578c4a", + "id": "39c60c2d", "metadata": {}, "source": [ - "The `create_multi_block_mechanics_functions` writes some functions for us to help write the problem in the right form. The most important part of the output is `mechanicsFunctions.compute_strain_energy(...)`, which writes all the loops over the blocks, elements, and quadrature points that you would need to compute the energy from the nodal displacements. (Now we finally see the `FunctionSpace` object `fs` in action).\n", - "To use this new function, we'll invoke it like this:\n", + "The `create_multi_block_mechanics_functions` writes some functions for us to help write the problem in the right form. The most important part of the output is `mechanicsFunctions.compute_strain_energy(...)`, which writes all the loops over the blocks, elements, and quadrature points that you would need to compute the energy from the nodal displacements. (Now we finally see the `FunctionSpace` object `fs` in action). To be able to write these loops, `create_multi_block_mechanics_functions` needs to know we are working in plane strain, and it needs to use the material models as well. \n", + "\n", + "The `mechanicsFunctions.compute_strain_energy(...)` function is invoked like this:\n", "```python\n", "mechanicsFunctions.compute_strain_energy(U, internalVariables)\n", "```\n", - "where `U` is the array of the degrees of freedom (the nodal displacements), and `internalVariables` is an array of the internal variables at every quadrature point. It has shape `(ne, neqp, nint)`, with `ne` being the number of elements, `neqp` the number of quadrature points per element, and `nint` the number of internal variables for the material model. (NOTE: For multi-material simulations, this is currently sized such that `nint` is the *maximum* number of internal variables over all materials, so that the array is padded for materials using fewer internal variables. We may remove this restriction in the future to improve memory efficiency).\n", + "where `U` is the array of the degrees of freedom (the nodal displacements), and `internalVariables` is an array of the internal variables at every quadrature point; that is, it has shape `(ne, neqp, nint)`, with `ne` being the number of elements, `neqp` the number of quadrature points per element, and `nint` the number of internal variables for the material model. (*NOTE: For multi-material simulations, this is currently sized such that `nint` is the maximum number of internal variables over all materials, so that the array is padded for materials using fewer internal variables. We may remove this restriction in the future to improve memory efficiency*).\n", + "Under the hood, this function tells optimism everything it needs to solve the problem through automatic differentiation. The derivative of `mechanicsFunctions.compute_strain_energy(U, internalVariables)` with respect to the nodal displacements `U` gives the residual function, and the Hessian yields the stiffness matrix (more precisely, the Hessian-vector product is evaluated to give the action of the stiffness matrix on a vector). All of these derivatives are taken automatically by optimism in the solver, so we don't need to worry about them at the application level.\n", "\n", "All of the minimization solvers in optimism require the objective function to have a signature like this:\n", "```python\n", "energy_function(Uu, p)\n", "```\n", - "where `Uu` is the array of unknowns, and `p` is called the parameter set. The parameter set essentially holds all of the information that is needed to specify the problem but *isn't* the set of unknowns. These are things like values of the boundary conditions and the internal variables, as well as some other categories. Parameter sets are constructed by calling the `Params` function in the `Objective` module. This is to help organize them in certain categories that the solver needs to be aware of, such as which ones have derivatives taken with respect to inside the solver. We need to cast our energy function in this form. Let's write it like this and work out what the intervening steps must be:\n", + "where `Uu` is the array of unknowns, and `p` is called the parameter set. The parameter set essentially holds all of the information that is needed to specify the problem but *isn't* the set of unknowns. These are things like values of the boundary conditions and the internal variables, as well as some other categories. We need to cast our energy function in this form. Let's write it like this and work out what the intervening steps must be:\n", "\n", "```python\n", "def energy_function(Uu, p):\n", @@ -604,15 +614,15 @@ " internalVariables = p.state_data\n", " return mechanicsFunctions.compute_strain_energy(U, internalVariables)\n", "```\n", - "On the first line, we use the parameter set to pass from the array of unknowns to the full array of nodal displacements. This means we need to fill in the values of the inhomogeneous boundary conditions. Next, we pull out the internal variables from the parameter set. Finally, we use the canned `compute_strain_energy(...)` function with these variables in order to compute the total energy.\n", + "On the first line, we use the parameter set to extend the array of unknown displacements to the full array of nodal displacements. This means we need to fill in the values of the inhomogeneous boundary conditions, which is what we'll do when we implement `create_field(...)`. Next, we pull out the internal variables from the parameter set. Finally, we use the canned `compute_strain_energy(...)` function with these variables in order to compute the total energy.\n", "\n", "The inhomogeneous boundary condition part is handled like so:" ] }, { "cell_type": "code", - "execution_count": 34, - "id": "217248f1", + "execution_count": null, + "id": "379d4116", "metadata": {}, "outputs": [], "source": [ @@ -627,73 +637,80 @@ }, { "cell_type": "markdown", - "id": "c578d78c", + "id": "559b6254", "metadata": {}, "source": [ - "We will store the applied displacement in the first slot of the parameter set. In line 4 above we extract it. Then we make an array of the same size as the nodal displacements, set it to zero, and replace the values in the DOFs on the top edge with the value of the applied displacement.\n", + "We will store the applied displacement in the first slot of the parameter set, `p[0]`. In line 4 above we extract it. Then we make an array of the same size as the nodal displacements, set it to zero, and replace the values in the DOFs on the top edge with the value of the applied displacement.\n", "\n", - "Now we can write the `create_field` function shown above in the proposed `energy_function`:" + "Now we can write the `create_field(...)` function shown above in the proposed `energy_function(...)`:" ] }, { "cell_type": "code", - "execution_count": 36, - "id": "2755c6db", + "execution_count": 35, + "id": "3a23d5ce", "metadata": {}, "outputs": [], "source": [ "# helper function to go from unknowns to full DoF array\n", "def create_field(Uu, p):\n", - " return dofManager.create_field(Uu, get_ubcs(p))\n", + " Uebc = get_ubcs(p)\n", + " return dofManager.create_field(Uu, Uebc)\n", " \n" ] }, { "cell_type": "markdown", - "id": "34511c5d", + "id": "fab84ba9", "metadata": {}, "source": [ - "The pieces are finally in place to define the energy function for our problem:" + "The `create_field(...)` method in the `dofManager` takes in the unknown displacements `Uu` and an array of the same size with the values of the essential boundary conditions `Uebc` and packages them up together to create an array of all DOF values. The pieces are finally in place to define the energy function for our problem:" ] }, { "cell_type": "code", "execution_count": 37, - "id": "be482478", + "id": "6ebd34c5", "metadata": {}, "outputs": [], "source": [ "# write the energy to minimize\n", "def energy_function(Uu, p):\n", " U = create_field(Uu, p)\n", - " internalVariables = p.state_data\n", - " return mechanicsFunctions.compute_strain_energy(U, internalVariables)\n", - "\n", - "\n" + " internalVariables = p[1]\n", + " return mechanicsFunctions.compute_strain_energy(U, internalVariables)\n" ] }, { "cell_type": "markdown", - "id": "5fab41f4", + "id": "7d271df0", + "metadata": {}, + "source": [ + "The slot `p[1]` is always reserved for the internal variable field." + ] + }, + { + "cell_type": "markdown", + "id": "2815942c", "metadata": {}, "source": [ "## Set up the optimization solver\n", - "We have an objective function - `energy_function`, which we will hand to a routine that will find the unknowns displacements that minimize it. In this section, we specficy which optimization solver we want to use, and tell it how to work. \n", + "We have an objective function - `energy_function(...)`, which we will hand to an optimization routine that will find the unknowns displacements. In this section, we specficy which optimization solver we want to use, and tell it how to work. \n", "\n", - "We will use the Steighaug trust region method. This method uses linear conjugate gradient iterations as part of its algorithm, which in turn need to be preconditioned in order to be effective. Currently, the only available preconditioner in optimism is a Cholesky factorization on the stiffness matrix. We need to intruct the solver how to assemble the stiffness matrix like this:" + "We will use the Steighaug trust region method. This method uses linear conjugate gradient iterations as part of its algorithm, which in turn need to be preconditioned in order to be effective. Currently, the only available preconditioner in optimism is a Cholesky factorization of the stiffness matrix. We need to intruct the solver how to assemble the stiffness matrix like this:" ] }, { "cell_type": "code", "execution_count": 38, - "id": "cd2bb5dd", + "id": "396fcac4", "metadata": {}, "outputs": [], "source": [ "# Tell objective how to assemble preconditioner matrix\n", "def assemble_sparse_preconditioner_matrix(Uu, p):\n", " U = create_field(Uu, p)\n", - " internalVariables = p.state_data\n", + " internalVariables = p[1]\n", " elementStiffnesses = mechanicsFunctions.compute_element_stiffnesses(U, internalVariables)\n", " return SparseMatrixAssembler.assemble_sparse_stiffness_matrix(\n", " elementStiffnesses, mesh.conns, dofManager)\n" @@ -701,65 +718,63 @@ }, { "cell_type": "markdown", - "id": "2b5da4b6", + "id": "50fbde23", + "metadata": {}, + "source": [ + "We see once again that `mechanicsFunctions` provides a helper function. In this case, the function `compute_element_stiffnesses(...)` takes the same inputs as the energy function, but instead of returning the total energy, it returns an array containing the stiffness matrix for each element. The elemental stiffness matrices are the Hessians of the total energy in each element, and automatic differentiation is again used to perform this calculation. The `assemble_sparse_stiffness_matrix(...)` function takes these elemental stiffness matrices and contructs the system's stiffness matrix using a sparse matrix data structure (from SciPy). We tell the solver how to use this capability by building something called a `PrecondStrategy`: " + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "id": "7aeca02f", + "metadata": {}, + "outputs": [], + "source": [ + "precondStrategy = Objective.PrecondStrategy(assemble_sparse_preconditioner_matrix)\n" + ] + }, + { + "cell_type": "markdown", + "id": "6e4b621c", "metadata": {}, "source": [ - "We see once again that `mechanicsFunctions` provides a helper function. In this case, the function `compute_element_stiffnesses` takes the same inputs as the energy function, but instead of returning the total energy, it returns an array containing the stiffness matrix for each element. The elemental stiffness matrices are the Hessians of the total energy in each element, and automatic differentiation is again used to perform this calculation. The `assemble_sparse_stiffness_matrix(...)` function takes these elemental stiffness matrices and contructs the system's stiffness matrix using a sparse matrix data structure (from SciPy). We tell the solver how to use this capability by building something called a `PrecondStrategy`: " + "Finally, we can specify custom settings for the solver is we wish." ] }, { "cell_type": "code", "execution_count": 40, - "id": "b1bc05ed", + "id": "01556d5f", "metadata": {}, "outputs": [], "source": [ - "precondStrategy = Objective.PrecondStrategy(assemble_sparse_preconditioner_matrix)\n", - "\n", "# solver settings\n", "solverSettings = EquationSolver.get_settings(use_preconditioned_inner_product_for_cg=True)\n" ] }, { "cell_type": "markdown", - "id": "d4b55505", + "id": "dbf08b3f", "metadata": {}, "source": [ - "Finally, we can specify custom settings for the solver is we wish." + "## Solve!" ] }, { "cell_type": "markdown", - "id": "3e10766a", + "id": "031e405b", "metadata": {}, "source": [ - "## Solve!" + "Parameter sets are constructed by calling the `Params` function in the `Objective` module. This is to help organize them in certain categories that the solver needs to be aware of, such as which ones have derivatives taken with respect to inside the solver." ] }, { "cell_type": "code", - "execution_count": 30, - "id": "28b893dd", + "execution_count": null, + "id": "c18895d2", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Assembling preconditioner 0\n", - "min, max diagonal stiffness = 0.9999999999999998 1.0000000000000002\n", - "Factorizing preconditioner\n", - "num warm start cg iters = 1\n", - "Assembling preconditioner 0\n", - "min, max diagonal stiffness = 0.9875365254464693 1.0110599625088337\n", - "Factorizing preconditioner\n", - "\n", - "Initial objective, residual = 0.04170630570636782 1.5580526760827706e-05\n", - "obj= -1.47133479e-10 , model obj= -1.47130159e-10 , res= 5.10352716e-10 , model res= 7.43900588e-21 , cg iters= 1 , tr size= 2.0 , interior , accepted= True\n", - "\n" - ] - } - ], + "outputs": [], "source": [ "# initialize unknown displacements to zero\n", "Uu = dofManager.get_unknown_values(np.zeros(fieldShape))\n", @@ -783,7 +798,7 @@ }, { "cell_type": "markdown", - "id": "2b1968f3", + "id": "782f1ac5", "metadata": {}, "source": [ "## Post-process\n", @@ -792,8 +807,8 @@ }, { "cell_type": "code", - "execution_count": 31, - "id": "e14a722a", + "execution_count": 41, + "id": "f5ac29dc", "metadata": {}, "outputs": [], "source": [ @@ -801,22 +816,23 @@ "writer = VTKWriter.VTKWriter(mesh, baseFileName='uniaxial_two_material')\n", "\n", "U = create_field(Uu, p)\n", - "writer.add_nodal_field(name='displacement', nodalData=U, fieldType=VTKWriter.VTKFieldType.VECTORS)\n", + "writer.add_nodal_field(name='displacement', nodalData=U, \n", + " fieldType=VTKWriter.VTKFieldType.VECTORS)\n", "\n" ] }, { "cell_type": "markdown", - "id": "da5415dd", + "id": "be97ee7c", "metadata": {}, "source": [ - "Let's also show an example of plotting quadrature field data. A commonly needed output is the stress field. We make use of another function in `mechanicsFunctions` to help us. However, before we write out any quadrature field data, we should update the internal variables. Currently, `state` still holds the initial conditions of the internal variables. That is, the solver finds the equilibrium displacement field, but doesn't change `state` in place. This is again due to Jax's functional programming paradigm. The following function call returns the updated internal variables using the new displacement field and the old internal variables as inputs:" + "Let's also show an example of plotting quadrature field data. A commonly needed output is the stress field. We make use of another function in `mechanicsFunctions` to help us. However, before we write out any quadrature field data, we should update the internal variables. Currently, `state` still holds the initial values of the internal variables (before the load step was taken). That is, the solver finds the equilibrium displacement field, but doesn't change `state` in place. This is again due to Jax's functional programming paradigm. The following function call returns the updated internal variables using the new equilibrium displacement field and the old internal variables as inputs:" ] }, { "cell_type": "code", "execution_count": null, - "id": "e46298c0", + "id": "d2926f79", "metadata": {}, "outputs": [], "source": [ @@ -827,7 +843,7 @@ }, { "cell_type": "markdown", - "id": "d23fc4a2", + "id": "57ea016b", "metadata": {}, "source": [ "We make use of another of the `mechanicsFunctions` member functions to get the stress field, using the updated internal variables. Note that we don't pay the cost of iteratively solving the history-dependent material model response again. This function expects that the internal variables are already updated when it is called, and simply uses theses values in the evaluation of the material model energy density functions." @@ -836,7 +852,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d3b3e743", + "id": "dc41d135", "metadata": {}, "outputs": [], "source": [ @@ -846,18 +862,18 @@ }, { "cell_type": "markdown", - "id": "b4221fcc", + "id": "b2cae332", "metadata": {}, "source": [ - "As the left-hand side of this assignment implies, the scalar energy density field (at every quadrature point) is returned in addition to the stresses. In fact, only the scalar-valued energy density function for each material is implemented in each material model. The stress field is the derivative of the energy density function with respect to the displacement gradient, and optimism uses Jax's automatic differentiation to compute this derivative. When computing this derivative, evaluating the energy density at the same time is essentially free, so the above function coputes them both.\n", + "As the left-hand side of this assignment implies, the scalar energy density field (at every quadrature point) is returned in addition to the stresses. In fact, only the scalar-valued energy density function for each material is implemented in each material model. The stress field is the derivative of the energy density function with respect to the displacement gradient, and optimism uses Jax's automatic differentiation to compute this derivative. When computing this derivative, evaluating the energy density at the same time is essentially free, so the above function computes them both.\n", "\n", - "The `stresses` array contains the stress tensor for every quadrature point of every element. There is no interpolation associated with a quadrature field, so in order to visualize the stresses, this must be remedied, either on the physics solver side or in the visualization software. One simple way to accomplish this is to compute element averages of the field and plot it as a cell-based field (instead of a nodal field). Optimism provides a method to do this:" + "The `stresses` array contains the stress tensor at every quadrature point of every element. There is no interpolation associated with a quadrature field, so in order to visualize the stresses, one must be created, either on the application side or in the visualization software, if it has this capability. One simple way to accomplish this is to compute element-wise averages of the quadrature field and plot it as a cell-based field (instead of a nodal field). Optimism provides a method to do this called `FunctionSpace.project_quadrature_field_to_element_field(...)`:" ] }, { "cell_type": "code", "execution_count": null, - "id": "e0d3a7bb", + "id": "55774221", "metadata": {}, "outputs": [], "source": [ @@ -872,22 +888,30 @@ }, { "cell_type": "markdown", - "id": "2558d9ca", + "id": "1f3b4a7e", "metadata": {}, "source": [ - "The above code can be wrapped in a loop that will compute the response at multiple time steps. Note that before moving to the next time step, you should update the internal variables in the parameter set:" + "This concludes the tutorial. As an enhancement, the load stepping part above can be wrapped in a loop that will compute the response over multiple time steps. Note that before moving to the next time step, you should update the internal variables in the parameter set:" ] }, { "cell_type": "code", - "execution_count": 33, - "id": "7e0acf41", + "execution_count": 42, + "id": "6fcb6668", "metadata": {}, "outputs": [], "source": [ "# update the state variables in the new equilibrium configuration\n", "p = Objective.param_index_update(p, 1, state)" ] + }, + { + "cell_type": "markdown", + "id": "782062d0", + "metadata": {}, + "source": [ + "You can look at other example simulations to see how to do this." + ] } ], "metadata": {