From 11b88a7af030c62a8672f3093a5dc03b3a414ad6 Mon Sep 17 00:00:00 2001 From: Brad Aagaard <baagaard@usgs.gov> Date: Mon, 28 Aug 2023 11:46:48 -0600 Subject: [PATCH 1/7] FIX: Fix missing arguments and typos in KinSrcTimeHistory Python interface. --- pylith/faults/KinSrcTimeHistory.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pylith/faults/KinSrcTimeHistory.py b/pylith/faults/KinSrcTimeHistory.py index 4a00c7e449..72af5c776f 100644 --- a/pylith/faults/KinSrcTimeHistory.py +++ b/pylith/faults/KinSrcTimeHistory.py @@ -61,12 +61,12 @@ def __init__(self, name="kinsrctimehistory"): KinSrc.__init__(self, name) return - def preinitialize(self): + def preinitialize(self, problem): """Do pre-initialization setup. """ - KinSrc.preinitialize(self) + KinSrc.preinitialize(self, problem) - ModuleKinSrc.setTimeHistory(self, self.dbTimeHistory) + ModuleKinSrc.setTimeHistoryDB(self, self.dbTimeHistory) return # PRIVATE METHODS //////////////////////////////////////////////////// From 5a3bbec5a3ba1a410aa884b4ef78b3d9b50cf680 Mon Sep 17 00:00:00 2001 From: Brad Aagaard <baagaard@usgs.gov> Date: Mon, 28 Aug 2023 11:49:00 -0600 Subject: [PATCH 2/7] Update full-scale test to use KinSrcTimeHistory. --- tests/fullscale/linearelasticity/faults-2d/Makefile.am | 3 ++- .../linearelasticity/faults-2d/shearnoslip.cfg | 6 ++++++ .../linearelasticity/faults-2d/zeroslipfn.timedb | 10 ++++++++++ 3 files changed, 18 insertions(+), 1 deletion(-) create mode 100644 tests/fullscale/linearelasticity/faults-2d/zeroslipfn.timedb diff --git a/tests/fullscale/linearelasticity/faults-2d/Makefile.am b/tests/fullscale/linearelasticity/faults-2d/Makefile.am index 2fe40b9b85..ff2a499995 100644 --- a/tests/fullscale/linearelasticity/faults-2d/Makefile.am +++ b/tests/fullscale/linearelasticity/faults-2d/Makefile.am @@ -49,7 +49,8 @@ dist_noinst_DATA = \ threeblocks_ic_tri.cfg \ shearnoslip.cfg \ shearnoslip_quad.cfg \ - shearnoslip_tri.cfg + shearnoslip_tri.cfg \ + zeroslipfn.timedb noinst_TMP = diff --git a/tests/fullscale/linearelasticity/faults-2d/shearnoslip.cfg b/tests/fullscale/linearelasticity/faults-2d/shearnoslip.cfg index 2a2ba9fb3f..426e726b0c 100644 --- a/tests/fullscale/linearelasticity/faults-2d/shearnoslip.cfg +++ b/tests/fullscale/linearelasticity/faults-2d/shearnoslip.cfg @@ -20,12 +20,18 @@ features = [ # ---------------------------------------------------------------------- # faults # ---------------------------------------------------------------------- +[pylithapp.problem.interfaces.fault_xmid.eq_ruptures] +rupture = pylith.faults.KinSrcTimeHistory + [pylithapp.problem.interfaces.fault_xmid.eq_ruptures.rupture] db_auxiliary_field = spatialdata.spatialdb.UniformDB db_auxiliary_field.description = Fault rupture auxiliary field spatial database db_auxiliary_field.values = [initiation_time, final_slip_left_lateral, final_slip_opening] db_auxiliary_field.data = [0.0*s, 0.0*m, 0.0*m] +time_history.description = Slip time function time history +time_history.filename = zeroslipfn.timedb + [pylithapp.problem.interfaces.fault_xneg.eq_ruptures.rupture] db_auxiliary_field = spatialdata.spatialdb.UniformDB db_auxiliary_field.description = Fault rupture auxiliary field spatial database diff --git a/tests/fullscale/linearelasticity/faults-2d/zeroslipfn.timedb b/tests/fullscale/linearelasticity/faults-2d/zeroslipfn.timedb new file mode 100644 index 0000000000..e13154a424 --- /dev/null +++ b/tests/fullscale/linearelasticity/faults-2d/zeroslipfn.timedb @@ -0,0 +1,10 @@ +#TIME HISTORY ascii +TimeHistory { + num-points = 5 + time-units = year +} + 0.0 0.0 + 1.0 0.0 + 4.0 0.0 + 8.0 0.0 +10.0 0.0 From 0434e9c088bcb5490ffee706f0fe8e734c24a209 Mon Sep 17 00:00:00 2001 From: Brad Aagaard <baagaard@usgs.gov> Date: Fri, 1 Sep 2023 10:13:17 -0600 Subject: [PATCH 3/7] FIX: Use public VecView() instead of private VecView_Seq() and VecView_MPI() in HDF5 writers. --- libsrc/pylith/meshio/DataWriter.cc | 29 +++++++++++++++++++++++ libsrc/pylith/meshio/DataWriter.hh | 12 ++++++++++ libsrc/pylith/meshio/DataWriterHDF5.cc | 25 ++----------------- libsrc/pylith/meshio/DataWriterHDF5Ext.cc | 25 ++----------------- 4 files changed, 45 insertions(+), 46 deletions(-) diff --git a/libsrc/pylith/meshio/DataWriter.cc b/libsrc/pylith/meshio/DataWriter.cc index 199898dc48..545d0bc37c 100644 --- a/libsrc/pylith/meshio/DataWriter.cc +++ b/libsrc/pylith/meshio/DataWriter.cc @@ -248,4 +248,33 @@ pylith::meshio::DataWriter::getCoordsGlobalVec(PetscVec* coordsGlobalVec, } // getCoordsGlobalVec +// ---------------------------------------------------------------------- +void +pylith::meshio::DataWriter::_writeVec(PetscVec vector, + PetscViewer viewer) { + // From plexhdf5.c DMPlexGlobalVectorView_HDF5_Internal + + /* To save vec in where we want, we create a new Vec (temp) with */ + /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */ + PetscVec temp; + const PetscScalar *array; + const char* vecName; + PetscLayout map; + PetscErrorCode err = PETSC_SUCCESS; + + err = VecCreate(PetscObjectComm((PetscObject)vector), &temp);PYLITH_CHECK_ERROR(err); + err = PetscObjectGetName((PetscObject)vector, &vecName);PYLITH_CHECK_ERROR(err); + err = PetscObjectSetName((PetscObject)temp, vecName);PYLITH_CHECK_ERROR(err); + err = VecGetLayout(vector, &map);PYLITH_CHECK_ERROR(err); + err = VecSetLayout(temp, map);PYLITH_CHECK_ERROR(err); + err = VecSetUp(temp);PYLITH_CHECK_ERROR(err); + err = VecGetArrayRead(vector, &array);PYLITH_CHECK_ERROR(err); + err = VecPlaceArray(temp, array);PYLITH_CHECK_ERROR(err); + err = VecView(temp, viewer);PYLITH_CHECK_ERROR(err); + err = VecResetArray(temp);PYLITH_CHECK_ERROR(err); + err = VecRestoreArrayRead(vector, &array);PYLITH_CHECK_ERROR(err); + err = VecDestroy(&temp);PYLITH_CHECK_ERROR(err); +} // _writeVec + + // End of file diff --git a/libsrc/pylith/meshio/DataWriter.hh b/libsrc/pylith/meshio/DataWriter.hh index 3dc8469c14..a8c95aa6a5 100644 --- a/libsrc/pylith/meshio/DataWriter.hh +++ b/libsrc/pylith/meshio/DataWriter.hh @@ -146,6 +146,18 @@ protected: void getCoordsGlobalVec(PetscVec* coordinatesVec, const pylith::topology::Mesh& mesh); + /** Write PETSc Vec using viewer. + * + * This method allows us to circumvent the usual place a vector is written and write it to + * a desired location. See DMPlexGlobalVectorView_HDF5_Internal() in plexhdf5.c in PETSc. + * + * @param[in] vector PETSc Vec to write. + * @param[inout] viewer PETSc viewer to use to write vector. + */ + static + void _writeVec(PetscVec vector, + PetscViewer viewer); + // NOT IMPLEMENTED ////////////////////////////////////////////////////// private: diff --git a/libsrc/pylith/meshio/DataWriterHDF5.cc b/libsrc/pylith/meshio/DataWriterHDF5.cc index 08c92dc811..03210aaf3a 100644 --- a/libsrc/pylith/meshio/DataWriterHDF5.cc +++ b/libsrc/pylith/meshio/DataWriterHDF5.cc @@ -40,15 +40,6 @@ #include <sstream> // USES std::ostringstream #include <stdexcept> // USES std::runtime_error -extern "C" { - extern PetscErrorCode VecView_Seq(Vec, - PetscViewer); - - extern PetscErrorCode VecView_MPI(Vec, - PetscViewer); - -} - #if H5_VERS_MAJOR == 1 && H5_VERS_MINOR >= 8 #define PYLITH_HDF5_USE_API_18 #endif @@ -207,13 +198,7 @@ pylith::meshio::DataWriterHDF5::writeVertexField(const PylithScalar t, err = PetscViewerHDF5SetTimestep(_viewer, istep);PYLITH_CHECK_ERROR(err); PetscVec vector = subfield.getVector();assert(vector); - PetscBool isseq; - err = PetscObjectTypeCompare((PetscObject) vector, VECSEQ, &isseq);PYLITH_CHECK_ERROR(err); - if (isseq) { - err = VecView_Seq(vector, _viewer);PYLITH_CHECK_ERROR(err); - } else { - err = VecView_MPI(vector, _viewer);PYLITH_CHECK_ERROR(err); - } + DataWriter::_writeVec(vector, _viewer);PYLITH_CHECK_ERROR(err); err = PetscViewerHDF5PopTimestepping(_viewer);PYLITH_CHECK_ERROR(err); err = PetscViewerHDF5PopGroup(_viewer);PYLITH_CHECK_ERROR(err); @@ -276,13 +261,7 @@ pylith::meshio::DataWriterHDF5::writeCellField(const PylithScalar t, err = PetscViewerHDF5SetTimestep(_viewer, istep);PYLITH_CHECK_ERROR(err); PetscVec vector = subfield.getVector();assert(vector); - PetscBool isseq; - err = PetscObjectTypeCompare((PetscObject) vector, VECSEQ, &isseq);PYLITH_CHECK_ERROR(err); - if (isseq) { - err = VecView_Seq(vector, _viewer);PYLITH_CHECK_ERROR(err); - } else { - err = VecView_MPI(vector, _viewer);PYLITH_CHECK_ERROR(err); - } // if/else + DataWriter::_writeVec(vector, _viewer); err = PetscViewerHDF5PopTimestepping(_viewer);PYLITH_CHECK_ERROR(err); err = PetscViewerHDF5PopGroup(_viewer);PYLITH_CHECK_ERROR(err); diff --git a/libsrc/pylith/meshio/DataWriterHDF5Ext.cc b/libsrc/pylith/meshio/DataWriterHDF5Ext.cc index 6fb83543ab..b230417f5f 100644 --- a/libsrc/pylith/meshio/DataWriterHDF5Ext.cc +++ b/libsrc/pylith/meshio/DataWriterHDF5Ext.cc @@ -37,15 +37,6 @@ #include <sstream> // USES std::ostringstream #include <stdexcept> // USES std::runtime_error -extern "C" { - extern PetscErrorCode VecView_Seq(Vec, - PetscViewer); - - extern PetscErrorCode VecView_MPI(Vec, - PetscViewer); - -} - // ---------------------------------------------------------------------- // Constructor pylith::meshio::DataWriterHDF5Ext::DataWriterHDF5Ext(void) : @@ -195,13 +186,7 @@ pylith::meshio::DataWriterHDF5Ext::writeVertexField(const PylithScalar t, assert(binaryViewer); PetscVec vector = subfield.getVector();assert(vector); - PetscBool isseq; - err = PetscObjectTypeCompare((PetscObject) vector, VECSEQ, &isseq);PYLITH_CHECK_ERROR(err); - if (isseq) { - err = VecView_Seq(vector, binaryViewer);PYLITH_CHECK_ERROR(err); - } else { - err = VecView_MPI(vector, binaryViewer);PYLITH_CHECK_ERROR(err); - } // if/else + DataWriter::_writeVec(vector, binaryViewer); ExternalDataset& datasetInfo = _datasets[name]; ++datasetInfo.numTimeSteps; @@ -339,13 +324,7 @@ pylith::meshio::DataWriterHDF5Ext::writeCellField(const PylithScalar t, assert(binaryViewer); PetscVec vector = subfield.getVector();assert(vector); - PetscBool isseq; - err = PetscObjectTypeCompare((PetscObject) vector, VECSEQ, &isseq);PYLITH_CHECK_ERROR(err); - if (isseq) { - err = VecView_Seq(vector, binaryViewer);PYLITH_CHECK_ERROR(err); - } else { - err = VecView_MPI(vector, binaryViewer);PYLITH_CHECK_ERROR(err); - } // if/else + DataWriter::_writeVec(vector, binaryViewer); ExternalDataset& datasetInfo = _datasets[name]; ++datasetInfo.numTimeSteps; From 493f0b2f702144ec98a3284d782ee0b289a5c39e Mon Sep 17 00:00:00 2001 From: Brad Aagaard <baagaard@usgs.gov> Date: Fri, 1 Sep 2023 10:13:47 -0600 Subject: [PATCH 4/7] FIX: Allow zero residual in MMS test with rigid block motion. --- tests/mmstests/linearelasticity/faults-2d/ThreeBlocksStatic.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/mmstests/linearelasticity/faults-2d/ThreeBlocksStatic.cc b/tests/mmstests/linearelasticity/faults-2d/ThreeBlocksStatic.cc index f1351501ef..34fa9dd9fd 100644 --- a/tests/mmstests/linearelasticity/faults-2d/ThreeBlocksStatic.cc +++ b/tests/mmstests/linearelasticity/faults-2d/ThreeBlocksStatic.cc @@ -213,6 +213,7 @@ class pylith::_ThreeBlocksStatic { data->journalName = "ThreeBlocksStatic"; data->isJacobianLinear = true; + data->allowZeroResidual = true; data->meshFilename = ":UNKNOWN:"; // Set in child class. From 4abba7e8e41967718df53f15f509914e14186b2f Mon Sep 17 00:00:00 2001 From: Brad Aagaard <baagaard@usgs.gov> Date: Fri, 1 Sep 2023 10:14:01 -0600 Subject: [PATCH 5/7] Update PETSc to 3.19.5. --- configure.ac | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/configure.ac b/configure.ac index 05370cac6d..802637d784 100644 --- a/configure.ac +++ b/configure.ac @@ -109,7 +109,7 @@ AC_CHECK_HEADER([mpi.h], [], [AC_MSG_ERROR([header 'mpi.h' not found])]) dnl PETSC AC_LANG(C) -CIT_PATH_PETSC([3.19.4]) +CIT_PATH_PETSC([3.19.5]) CIT_HEADER_PETSC CIT_CHECK_LIB_PETSC From 9df1d3362b091c890c2c6fc0692fb54fbf0e6b44 Mon Sep 17 00:00:00 2001 From: Brad Aagaard <baagaard@usgs.gov> Date: Fri, 1 Sep 2023 10:12:11 -0600 Subject: [PATCH 6/7] FIX: Fix memory leak associated with buried edges. --- libsrc/pylith/faults/FaultCohesive.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/libsrc/pylith/faults/FaultCohesive.cc b/libsrc/pylith/faults/FaultCohesive.cc index 81000c572c..fefdd41489 100644 --- a/libsrc/pylith/faults/FaultCohesive.cc +++ b/libsrc/pylith/faults/FaultCohesive.cc @@ -375,6 +375,10 @@ pylith::faults::FaultCohesive::createConstraints(const pylith::topology::Field& } // if } // for } // for + if (pointIS) { + err = ISRestoreIndices(pointIS, &points);PYLITH_CHECK_ERROR(err); + err = ISDestroy(&pointIS);PYLITH_CHECK_ERROR(err); + } // if std::vector<pylith::feassemble::Constraint*> constraintArray; pylith::feassemble::ConstraintSimple *constraint = new pylith::feassemble::ConstraintSimple(this);assert(constraint); From 60ce6eaa3c4ccaf16094cd8ae59428d1cc043c28 Mon Sep 17 00:00:00 2001 From: Brad Aagaard <baagaard@usgs.gov> Date: Fri, 1 Sep 2023 15:11:44 -0600 Subject: [PATCH 7/7] Add 3D full-scale test with buried fault edges. --- configure.ac | 1 + tests/fullscale/linearelasticity/Makefile.am | 1 + .../faults-3d-buried/Makefile.am | 46 +++++++ .../faults-3d-buried/TestUniformSlip.py | 93 ++++++++++++++ .../faults-3d-buried/generate_gmsh.py | 103 +++++++++++++++ .../faults-3d-buried/mesh_tet.msh | Bin 0 -> 100773 bytes .../faults-3d-buried/meshes.py | 42 ++++++ .../faults-3d-buried/pylithapp.cfg | 121 ++++++++++++++++++ .../faults-3d-buried/test_pylith.py | 51 ++++++++ .../faults-3d-buried/uniformslip.cfg | 108 ++++++++++++++++ .../faults-3d-buried/uniformslip_soln.py | 96 ++++++++++++++ 11 files changed, 662 insertions(+) create mode 100644 tests/fullscale/linearelasticity/faults-3d-buried/Makefile.am create mode 100644 tests/fullscale/linearelasticity/faults-3d-buried/TestUniformSlip.py create mode 100755 tests/fullscale/linearelasticity/faults-3d-buried/generate_gmsh.py create mode 100644 tests/fullscale/linearelasticity/faults-3d-buried/mesh_tet.msh create mode 100644 tests/fullscale/linearelasticity/faults-3d-buried/meshes.py create mode 100644 tests/fullscale/linearelasticity/faults-3d-buried/pylithapp.cfg create mode 100755 tests/fullscale/linearelasticity/faults-3d-buried/test_pylith.py create mode 100644 tests/fullscale/linearelasticity/faults-3d-buried/uniformslip.cfg create mode 100644 tests/fullscale/linearelasticity/faults-3d-buried/uniformslip_soln.py diff --git a/configure.ac b/configure.ac index 802637d784..fc8a351c7d 100644 --- a/configure.ac +++ b/configure.ac @@ -250,6 +250,7 @@ AC_CONFIG_FILES([Makefile tests/fullscale/linearelasticity/nofaults-3d/Makefile tests/fullscale/linearelasticity/faults-2d/Makefile tests/fullscale/linearelasticity/faults-3d/Makefile + tests/fullscale/linearelasticity/faults-3d-buried/Makefile tests/fullscale/linearelasticity/greensfns-2d/Makefile tests/fullscale/incompressibleelasticity/Makefile tests/fullscale/incompressibleelasticity/nofaults-2d/Makefile diff --git a/tests/fullscale/linearelasticity/Makefile.am b/tests/fullscale/linearelasticity/Makefile.am index e739ab9874..cc988d7ba7 100644 --- a/tests/fullscale/linearelasticity/Makefile.am +++ b/tests/fullscale/linearelasticity/Makefile.am @@ -20,6 +20,7 @@ SUBDIRS = \ nofaults-3d \ faults-2d \ faults-3d \ + faults-3d-buried \ greensfns-2d diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/Makefile.am b/tests/fullscale/linearelasticity/faults-3d-buried/Makefile.am new file mode 100644 index 0000000000..629f948d38 --- /dev/null +++ b/tests/fullscale/linearelasticity/faults-3d-buried/Makefile.am @@ -0,0 +1,46 @@ +# -*- Makefile -*- +# +# ---------------------------------------------------------------------- +# +# Brad T. Aagaard, U.S. Geological Survey +# Charles A. Williams, GNS Science +# Matthew G. Knepley, University at Buffalo +# +# This code was developed as part of the Computational Infrastructure +# for Geodynamics (http://geodynamics.org). +# +# Copyright (c) 2010-2022 University of California, Davis +# +# See LICENSE.md for license information. +# +# ---------------------------------------------------------------------- + +TESTS = test_pylith.py + +dist_check_SCRIPTS = test_pylith.py + +dist_noinst_PYTHON = \ + generate_gmsh.py \ + meshes.py \ + TestUniformSlip.py \ + uniformslip_soln.py + +dist_noinst_DATA = \ + mesh_tet.msh \ + pylithapp.cfg \ + uniformslip.cfg + + +noinst_TMP = + + +export_datadir = $(abs_builddir) +include $(top_srcdir)/tests/data.am + +clean-local: clean-local-tmp clean-data +.PHONY: clean-local-tmp +clean-local-tmp: + $(RM) $(RM_FLAGS) -r output __pycache__ + + +# End of file diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/TestUniformSlip.py b/tests/fullscale/linearelasticity/faults-3d-buried/TestUniformSlip.py new file mode 100644 index 0000000000..3ea7bab78e --- /dev/null +++ b/tests/fullscale/linearelasticity/faults-3d-buried/TestUniformSlip.py @@ -0,0 +1,93 @@ +#!/usr/bin/env nemesis +# +# ---------------------------------------------------------------------- +# +# Brad T. Aagaard, U.S. Geological Survey +# Charles A. Williams, GNS Science +# Matthew G. Knepley, University at Buffalo +# +# This code was developed as part of the Computational Infrastructure +# for Geodynamics (http://geodynamics.org). +# +# Copyright (c) 2010-2022 University of California, Davis +# +# See LICENSE.md for license information. +# +# ---------------------------------------------------------------------- + +import unittest + +from pylith.testing.FullTestApp import (FullTestCase, Check) + +import meshes +import uniformslip_soln + + +# ------------------------------------------------------------------------------------------------- +class TestCase(FullTestCase): + + def setUp(self): + defaults = { + "filename": "output/{name}-{mesh_entity}.h5", + "exact_soln": uniformslip_soln.AnalyticalSoln(), + "mesh": self.mesh, + } + self.checks = [ + Check( + mesh_entities=["domain"], + vertex_fields=[], + defaults=defaults, + ), + Check( + mesh_entities=["mat_elastic"], + filename="output/{name}-{mesh_entity}_info.h5", + cell_fields = ["density", "bulk_modulus", "shear_modulus"], + defaults=defaults, + ), + Check( + mesh_entities=["bc_xneg", "bc_xpos", "bc_yneg", "bc_ypos", "bc_zneg"], + filename="output/{name}-{mesh_entity}_info.h5", + vertex_fields=["initial_amplitude"], + defaults=defaults, + ), + Check( + mesh_entities=["fault"], + vertex_fields=["slip"], + defaults=defaults, + ), + ] + + def run_pylith(self, testName, args, nprocs=1): + FullTestCase.run_pylith(self, testName, args, nprocs=nprocs) + + +# ------------------------------------------------------------------------------------------------- +class TestTetGmsh(TestCase): + + def setUp(self): + self.name = "uniformslip" + self.mesh = meshes.TetGmsh() + super().setUp() + + TestCase.run_pylith(self, self.name, ["uniformslip.cfg"], nprocs=1) + return + + +# ------------------------------------------------------------------------------------------------- +def test_cases(): + return [ + TestTetGmsh, + ] + + +# ------------------------------------------------------------------------------------------------- +if __name__ == '__main__': + FullTestCase.parse_args() + + suite = unittest.TestSuite() + for test in test_cases(): + suite.addTest(unittest.makeSuite(test)) + unittest.TextTestRunner(verbosity=2).run(suite) + + +# End of file diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/generate_gmsh.py b/tests/fullscale/linearelasticity/faults-3d-buried/generate_gmsh.py new file mode 100755 index 0000000000..2df4ec47c1 --- /dev/null +++ b/tests/fullscale/linearelasticity/faults-3d-buried/generate_gmsh.py @@ -0,0 +1,103 @@ +#!/usr/bin/env nemesis + +import numpy +import gmsh +from pylith.meshio.gmsh_utils import (VertexGroup, MaterialGroup, GenerateMesh) + +class App(GenerateMesh): + """ + Block is DOMAIN_X by DOMAIN_Y x DOMAN_Z with discretization size DX. + + -30km <= x <= +30km + -30km <= y <= +30km + -30km <= z <= 0 + """ + DOMAIN_X = DOMAIN_Y = 40.0e+3 + DOMAIN_Z = 20.0e+3 + FAULT_LENGTH = 20.0e+3 + FAULT_WIDTH = 10.0e+3 + DX = 5.0e+3 + + def __init__(self): + self.cell_choices = { + "required": False, + "default": "tet", + "choices": ["tet"], + } + self.filename = "mesh_tet.msh" + + def create_geometry(self): + """Create geometry. + """ + gmsh.initialize() + self.box = gmsh.model.occ.add_box(-0.5*self.DOMAIN_X, -0.5*self.DOMAIN_Y, -self.DOMAIN_Z, + self.DOMAIN_X, self.DOMAIN_Y, self.DOMAIN_Z) + + self.s_xneg = 1 + self.s_xpos = 2 + self.s_yneg = 3 + self.s_ypos = 4 + self.s_zneg = 5 + self.s_zpos = 6 + + p_faulttrace_yneg = gmsh.model.occ.add_point(0.0, -0.5*self.FAULT_LENGTH, 0.0) + p_faulttrace_ypos = gmsh.model.occ.add_point(0.0, +0.5*self.FAULT_LENGTH, 0.0) + p_faultbot_yneg = gmsh.model.occ.add_point(0.0, -0.5*self.FAULT_LENGTH, -self.FAULT_WIDTH) + p_faultbot_ypos = gmsh.model.occ.add_point(0.0, +0.5*self.FAULT_LENGTH, -self.FAULT_WIDTH) + + self.c_faulttrace = gmsh.model.occ.add_line(p_faulttrace_yneg, p_faulttrace_ypos) + self.c_fault_yneg = gmsh.model.occ.add_line(p_faulttrace_yneg, p_faultbot_yneg) + self.c_faultbot = gmsh.model.occ.add_line(p_faultbot_yneg, p_faultbot_ypos) + self.c_fault_ypos = gmsh.model.occ.add_line(p_faulttrace_ypos, p_faultbot_ypos) + + loop_fault = gmsh.model.occ.add_curve_loop([-self.c_faulttrace, self.c_fault_yneg, self.c_faultbot, -self.c_fault_ypos]) + self.s_fault = gmsh.model.occ.add_plane_surface([loop_fault]) + + gmsh.model.occ.synchronize() + gmsh.model.mesh.embed(2, [self.s_fault], 3, self.box) + gmsh.model.mesh.embed(1, [self.c_faulttrace], 2, self.s_zpos) + + def mark(self): + """Mark geometry for materials, boundary conditions, faults, etc. + """ + materials = ( + MaterialGroup(tag=1, entities=[self.box]), + ) + for material in materials: + material.create_physical_group() + + vertex_groups = ( + VertexGroup(name="boundary_xneg", tag=10, dim=2, entities=[self.s_xneg]), + VertexGroup(name="boundary_xpos", tag=11, dim=2, entities=[self.s_xpos]), + VertexGroup(name="boundary_yneg", tag=12, dim=2, entities=[self.s_yneg]), + VertexGroup(name="boundary_ypos", tag=13, dim=2, entities=[self.s_ypos]), + VertexGroup(name="boundary_zneg", tag=14, dim=2, entities=[self.s_zneg]), + VertexGroup(name="boundary_zpos", tag=15, dim=2, entities=[self.s_zpos]), + VertexGroup(name="fault", tag=20, dim=2, entities=[self.s_fault]), + VertexGroup(name="fault_edges", tag=21, dim=1, entities=[self.c_fault_yneg, self.c_faultbot, self.c_fault_ypos]), + ) + for group in vertex_groups: + group.create_physical_group() + + def generate_mesh(self, cell): + """Generate the mesh. Should also include optimizing the mesh quality. + """ + gmsh.option.setNumber("Mesh.MeshSizeMin", self.DX) + gmsh.option.setNumber("Mesh.MeshSizeMax", self.DX) + if cell == "hex": + gmsh.model.mesh.set_transfinite_automatic(recombine=True) + #nnodes = 3 + int(self.DOMAIN_Y/self.DX + 0.5) + #gmsh.model.mesh.set_transfinite_curve(self.l_xneg, nnodes) + #gmsh.model.mesh.set_transfinite_curve(self.l_split0, nnodes) + #gmsh.model.mesh.set_transfinite_surface(self.s_xmid, cornerTags=self.xmid_corners) + #gmsh.model.mesh.set_recombine(3, self.s_xmid) + + gmsh.model.mesh.generate(3) + gmsh.model.mesh.optimize("Laplace2D") + + +if __name__ == "__main__": + App().main() + + +# End of file diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/mesh_tet.msh b/tests/fullscale/linearelasticity/faults-3d-buried/mesh_tet.msh new file mode 100644 index 0000000000000000000000000000000000000000..fdfd2d31bbf99e0a029eedd522224efd421b3698 GIT binary patch literal 100773 zcma(43A~N<*Z+?n8B>&^NQ2BmnUWG+N<tCIkO*ZgN>L(Y9wJndB$=lqnagx>%yY9e z2~A3=k7)MW*E!GUKDylJ`~UC9LoaK+UTeLl^&a-F&wUqf(4j}?I*&ZotzFODRj#gB zq+*e3xwE9Cq~tDs>%;9&<BB)x+`C7Y2ikRQ*sfcL9=R)3$z7pH#R^4=-v7ws54Ugk zSnoDZKHQ;G(a2fxbZ7TRdZalko#yPVvhwLpl~qo2_ECB5=}wfD60LV^_juQyQ3I6{ zt)>}mI<)W9A<<@HBxm+6;XJc<3Fn!;OE}N$UBY>0?-I^4dsnf@X}wDr@!nPXzj~K& zp4q#E^UU5QoM-kf;XJc<3Fn!;t8{wrDkqk_Xkzs{Jl3UM*YaK3-%zn=+S>njy%W{Y zvrEq|iIqP$Q5@0d{B%gk2}VT|eZ-Uz?}up)F=6CLb0mr@fzh8RSM=vR{+yM-w2v{T zVP>&;69Lhmtm%;{r?ngCipS<j5YZpo&89G;y&%n%urrIb-Rv2$$rxS6o|AAyf8yhn zPxSe}v8S~gkF{Or{j}Kh9Eo<!YixYHxe}VRPe!`|c5=M&Slc~2Ju>5;BT_>R*C2nw z5&bEU_!F)3e?B8p1MUAyAGQ0RkNFemI-Q<BYQvgL#@WXIBPOcTwhtx#k(lc^P91;h z)T!a`hhkdZeSK<LpVqE+go%h#5yo@XeRK7sUq9}%)?Dxb{{Q5R>d~B0|NcqjjONIk z<OpK6{SOk-l$2OMqqR3@G=EXf2+a}4ct5qx5sf{n1?>p`PUMXCLA-yl+=00g(*N`` zTJvZP%^B72Ux}O%#&bHK<aRF3n9nHZ>GK)&?~g>zXy2GqjQ7*r3F#UAJbgZ+eH7)4 zFq2%#wJ|5U(6L7SJbgZ+G5?;(8SS%Jzsw~j=QEx&n#<F3Ms59>$Qi9&JZF49%@wcD zDCg<(8P#XBMp6IGnaLaknT<V~%hUTA)nzpHsI7R;ct6b<?`M?r^nOPBCt9P(Z;tSz zgpoKKw9SD(SxeCvqfd0+xL(n^oJi!0`fa@VMDqn@62jr>=i{-DGtLjs#?y2DmXJpM z{{NNJKAPJc@to1R98cto?rC!p6U6&z&V)4O|JBc^exjWJ|D2t!Pn0t{lbp*ahv#@A zRoYn;t<k@Uzei)qq=rKCi3#a{@<!t|Z&a7jo{aW-RQIQ!p~Ra{Ozx#b$ba%iwP@ax ziCob*qA^5mM1J#S((fR%{zmm^-e~QkK1OXtwHlpO=FJ3?dEKwhZ&a7&jn3t$UZc5) z+KA4Ucs^^xeu~zxd867iZ?sm?97Oq|dWh!PJl2+Io_LL!H>z9nMspI?>o19%zb5`f zns|R>c5)sP(*LYm6dUy~`Xo}N)s}akcT3tsH~F$H-tdw3i3e`+^q7=+^w5rKN3^~; zriEraBEvMBf4EVN(f&VU$PpDq^f5*|&6&|okY@7_C({`1|1*XhQNcwYW3<zp8IQy; z&E_9Yt})vGXAJH5Xf;MV78B4!8brk#eT>mgb7u4{zz!qyk!p-~n#O<lu*8lEI{Fx+ zo$lqIJ}fb#Lnrzeqn+;MpSA-=n>_j$qn+;MpSF{TiF9gR|7fE`pJ=m3pJ;O?KN)oq zrupvJm5y<LvR|<;Y9p4UwP}62Mt<STM?)8Y7lP47jy}<YJo;Re4(b22ZTljp-@csT zUlQs{L-qf0CJg-}mIVKylK#~H+keNdq(7Qt=U2R#WCg)OU>YQI^X1L4#r=-$41fAO zyFZ=hOK77ocqw=pcsY0lSOiQ1+c>i?VLG4w7|EC~*PeMcUNpCM5}D_uZ64-1?suLO z{%D>(U!3PFX{RVy%yuKp*yCZEJ7a7?nw@MbIkuQTeV#p^oM+GZ5@|_NQc8lQz|vqD z@G3A3?2~<w)9;#};Xk9E110~#WS%qXIbHAhV32y}IpL4y*|9p$Wzm-d%Y$iP4$-k` zC-IDV4$|Y(=h^ewwy&m*YrqO%MX(ZB8LR@Pfo*2=CH#-yc}~^<dogi_p%a5N+kc#? zPv_6@pHa{0I_Ei|kLKC?#W7z?J5|B!z-r+2;0<7PFb&QeYnbLgbFBE!7;CyeEhzDi zj4`PX(x@~1XN)ynf7)1$i^lBz=$LOzBU4gpfHlFJzyPcT-VCNee9Ut3F(>`8G3yr{ zThfl5S4n?t%=+VFPWYoSd%rv8ThP}AZw2dsw}ExRdSDu4j5$cN?Q_PMWB%lr|Lac+ zPW;2YlFXAF^MCc89dUien3MWw%)U<?bNz%b`rMulDR+Pkz=q(RU?VUM;$xPJk2&d& zjak3w7?O5u%t?Q8odaEI7$0-O8I9TZmFH$-g6;yFfd2u{oc}ogAI52Km8ebMa?#(T z&)w<p-+w0ODgCVV-Df|WqHYE@|8I1<Wk1ao*TntyGvSZ=<~z`Swur|^cCyBnr2WvP z);G&N#IywO1zUlw!TZ2AU|X;qct7|6*dFWvb_6?toxv{PgWyA8SFjuSF!%`A9efmg z4D12+1Rn>V0DFOHkl3Gnuc=>S#rLbjh3(!Odrp?OLj6TX-d3l>k=KHjo97)`FmX}R zzv#ISo7OxuBm8OjkOr04>`waI-}P1MmsQsUmkexN`?k7!lKvTAK30CwxLskE4b>*4 zjNOv-=MUbP{Y$;A!GR~U&l=kI&7{9;?uz>>7TXn6xcRaMCCY3~`b$6X@2v$J&j_CR zqyF&v+cMfsxnxkD8|SwQcQ*f{$J@&`C-s~EJ-mF|m?hz~NoPOu$uAlHB2Sjy5<WCD z98&+D9Jk!II;pQv?&=9Oi>(POm49~KfZQ4VfBKrSb83}M4e#%J&ZB$#X7vBqclR{h z{paqmc;nn-X5OBW=lHJCf4+F(n&9vI*WEd;NyfNt=s06dg|||J(jQeX+kS4CY<KeO zxvC9%XlXF}&!uImG|T8`{)dPCxM=&<aNP5UhFn)`VNyT7#xG?i=G+=q9en??ZuK(y zSz*MnHSbQ}9o%(&*$)5wnbDt(4}UqKTDGmhJ73T3bMq$|`tC<xYSW?2(%{F@r5C)? zCPSYReMZa(7quJq#@}DhNw!<N)%l$+PF)>-eD095yDiJ;|KaC<%GD~d$66KaeBr0H zwk7q~7GCt!m*4LR#VBob|5r=*2c=)QDc_)p3zGV!#X9Yok!x>o_S!F>y0zxsq~G`| zG5yTq-NwFm*XD5fio^F_f7`mG{#f~MtK_fyYEZY*xLPMlhiQA>?~m^JClg(XKE2Z+ zr4RTN*ca>v_6MH^2Y>^?LEvC;2sjiR1`Y>DfFr?Yz)|36@L6ySI2IfSJ_kMzz5tE~ zUj!$BX^^<5Mt8sD=e8HD33D_mT%%;x<em>Q+4Dh*lwJ$VAD$7meB`VKJJ%=oUXaP2 z558OfV1dqUwuW7r7dxD9M}|I=Js)JPKJIv@c6-9re{Sts>eE%B{mkSn2s%G?sNJzf zGlIu{`|_>4+cWZ~JpXXxJd0X|+1_v5Xv*@<$vl@{IAdnX*O!C~D)-Md;jaw;+>eV~ z`^A$p!)oikx~^Nj<Q@p7{x*NbeNESdGhTVI$<3c;<iFs%X8+`TJT(k|XtJSrx8zw6 zWU}XjO!j<`?VJ*O`X%=KzV<yww8_}>!HnXSt4;iMYdG=carcZWx-i-QO!jt=$=(h| zef@i>`b~C+Z)f@N-2uxp`WZBTuXwcQ=MJh`VQR*n4^pDfnR`Asl55KB3BRlk>t~(* z;|D8}dp`K!taB#bo7nSJwj5t~C}YnD?o}~8@B2X}dp^iy&nNT5_Iz;mdA}67>Fdqm z)JwWn9CTB1&j<h4p7$N%o_{H=Gbt$(!I!~F;49!{@KtaM_!>ABoCZz@XMi)oS>SAN z4mcN_2hImy2N!^EfD6Gl!A0O(;9_tIxD;Fl)?AnUg}Vx`3DeBP`8#;S-@EG_-W_&t z+O6<2BW48pYvo*gS-CmZ`i!K0=ixht6gjdeEO_bVIqRQVozw?Q78a`i$D*+Ci_LqN zJMCV~c5GP7ih~w~htto1q~5dOi&3Laty{J$nSbhl0`pVmZ3#bkv}3C*t@kALJ%{Bk zTkOcHuzQ(Pzf|v)(VszU`u+6l^4($kBDwoaymLu1PwBVktXY42bGU!RKgX-wnVR&^ zty^_Z-Ktx{zbB?tpRqK#$HU)C_Bpuum4w}QK#n&WR7kcv|GHfZKKx*FxUc4fQOyS| zN&4&TI`7|nCpU-3%9VcaWdDqL+d1j^g+K3E6@KtXzgb)BuSx3v9iQ{dZ@$?bwklA) ze5tOfN&meg_dl^CadwPo`qupWk8Vx+gIyO){;Al)aN@olKWA^UJL$jn<lKj*C-(5s z=C9RzFi*yNXu#%wAOG%-J>k8Ty41fpCF5Qnb28=Nr<bh>Tg>@mR{07U^LG85soR(I zOTQ0O?#QSI@A<iHuFBo;i;VTHdExhmk0stIxl$6J)9Pd@?}FGKPu}x!zxR9v?)@lT z^bSzx*`U8T_x#nteJ6H&v^mhfx6{(`@7Gut_N?(!mVchv66)W1*6OmK_evjEwlrr- zihF)J_6l$%_%`?sxC&eit^wDA>%b6P4{iWc!HwW1a5K0C+zM_3w}U&to!~BTH@FAf z3%(1!2fhz}0LJ!wkjb78GTHM%CVM`}WX}hg?D-&*Js)JU=Yvf4e2~eW4>H;FK_+`X z$Yjq4ne6!>lRY0~vgd<L_I!}Zo)0qF^FbziKFDOx2bt{oAd@{GWU}XjO!j<`$(~Q% z)9&T8@=LrgGuiV=eQeJMne6!>lRY1tzUQ-3^EZPZf**k&gZscw!2RH-;Ai0H;1}SR z-~sR}@E~{yJPdvfegl3Beg}RJ9sz#<e*}+$KY>4k$G~5}U%~%@$HCvg6W~ekckmDJ zPw+4BZ}1=RU+@%|_G47y?^&{>JI(^Lg4w|AU=A=Rm<v1`O#4rJ6RFO@J{Qab<^|6K z&j<5?7l0Rn`N0C<Mc|}7`_}wC>CAV{fv%n2*)e%{P~@s5wJ(1v<DTrkb;JePn$Her zE@)78@NF6Q@ONuZb$H_Fj9~H6Wwl0EHQ)Q;qKzXyn)=V$pmxbaeQv%pLx0u3R~{;w zxc^(6-Fwxzo7MzrvzEB`?t4A|A7?-MPFT9c{dv~Cm+=m|`{a^+1^ev|Dtw+j@0J}I z`i)heY*PEn-NDv%y&nDFc^UVArEXJ>F1dGi*tqhSyPm0+;eUD9>>ZQFZVhi7^TX%Y zJija1Zu_^g3?4IYcW~^`g2Tm&XWakQCjC03TiMM)$JupSecvnNJ=HSnlvO!PM)&Xe z9cOfX;U@cE`qwe77Co>kXnWR_rg;Zt%-f<JQ=i+>cvbLSo}MLkC%&7hZ*=(k9X*fC z3{n<N?NX`Q`C-!E=*>-sR?iF$6uRiSGW9dwSyQu&>G1YvFNe7v&DT2bV;Onc_O4Vo z*PfTdy-!pdc5-w^yI)kFGIVL{6=BMOMo*1Dk|mkH?X6Q@zw58%;ebQ)t1SE_BY%Y( zU%8~f^((@1H<oy!?%<5^p5LoigMAOJ3MN0YcKBblGv=r5tSieOs5mLCnSJe;Qq3~_ z6VLtXn%UpI6HGn0_WsZJWW1YpPQ7@~qs2CcQ=|XVV|YgY&)<1p_g5dB87wN)^NL}I zHU!SYN56IM^>NP?;qOz<JHJ=|)MWo(n~^&0)bR~rg|@jXEt!{ir$(Pd{TzSr)eE0o zo*K^WHtgmtr8DYe%B-n<%f7iRoPK_-#lP*{7^vUdY1{rC8zzUthF{%a`vVzuvXy>z zPs^JzpZO<T_gRS-6Lqrk+C>`|Wc0s#ll)6kf8G>Ke0$J!=Vo7;%%9`kgL@wgHU$^d z+EwG8%Qq$czdt|v`K9OV31+@q^WFkeGTyn3Mz6Z+>n&@-ebb8+f2&!>e(7HJP`)QV zof(W+Hs*pAFJ`PukpnfFJ+uAg@Uv(>Z_gNSuSt{2uW!C59DUp54hJvHSeMv-Pxfcr zb<_K8{bNQ@`r4|uyqqIrUyLYHX2X!NGlNEVJksZy*E8mSWwku(9(-k0&~shI>i0Ct zn77d#SFJA6adWV`e2<C~$7bjsD_Z)zLlw6K?T_~D)bxUk^Wn?iZn|>YyL-aE^HL`c zJwIa}dRFZ)Y{T&_;gX#nRCx2R4E^(${dVe_aa)7#KlE<)#fgmhnec7i57rb&4Tqd+ zvZBbp(Y{Q467`V(lRJ8h7?$`Cl2_b5ty;@1q5iOX$05(o8xxLOaQu$zwr1$l`5EEd z+8<4tdBfhMzH7DC|2$D}O|WHu!!>np$?&%w&}HF~Vk?8%1AhAHpY<7eyz@od<!%qF z#NPYf`SR)mbfW&7KJK--LYr;r|8Z*a-xIta*1(^x%h=b|Pxq%xO@G2C<|q5nOXg&& zvn^b;dRm^ol{SRtzp+@(R|`F~Ej+Nj#^?2Vqz3xWO8g~9nvs%H6TCP}`d<oyg}_U| z!r-OgW#HxD6<`tYO0Xzc3@i?o084_Uz|vqD@G7t@SPm=?UJYIYRsbu4mB7ki74TZH zDtH}O4ZI$_0jv(*2-X1oFD-Av4!~OA&EPFyZSYpG4tN__7pw=?2X6=O02_b}!8^f5 zU}NwuunBlK_~hdY4ix)&cNi95-15%SoV|$>I8Y;N(^uE-4qCmj@!C_3Idc;JIfHio za$w5ta72@tWp_-N5u{g3LjPN{<x3at4mv;m;T3Z~o}Todm~hLb*G*Uzj?de5&N+wJ z<B6cY3)j52s6is{tRs&PUB;P^@OQ1XqQ`mTR|Ug9Yx{7{4H^C+l~>PiTWC$Nbkn#l z-)zdr<M%FphfTJd?YuiTby&VDm|thjl037wB>kt(9@g-R#P7Oo$878H&~2QtiFWG_ zezyM7hMU7ik2WgScm3j|e%sut?F+xLIZVm1ze&~2_k{^*+NV{uDWg9L)`V^6W?6pp z>X!rDDR=igx_8lM>k?;CmLkorni=Zf`+U9)y$Y@hyLSAk-r8fUlKK&w8oz(l?8JC# z<r!4y4eq8y(2*@qf8B1*=HP?(e=k{m6?a6!-|FX+caJT;C74j>$;X;_CMW#ISFQ~U zH(wj%_+oJPrAIQ>E7$58tv7dC8*J@y`?8mw&*;ySGnf9eZR75s!0}Zl`+HX=@^n8w ze^TDPyMu|BzR^2%Q$|1Ae>HaMkApS`$LDA3JF!^C{Ist&_?rsvZ4NK#yzzw!4O4?e zinPy9?_M$Q(>;kd;ht@UemFKg#GR7kk%@zk-L*CB`SYF6UH*N>`VOnN`~LYarH1wH zzx3<F*|sG0Eq>fT^XGX>!bXh;Jw4@u%}IY+zQk`e|G5Jb&fiCTUh2RLi8Ft|sJes8 zW#mU6PECI!2WfWVADzby-`BR!jG$W6*6-CSpHWu_(%v*nf*YbUy99SZ6x8<qDz7z5 z4W3vW+<Pn~qpoh;HsPsTzS<mCX?f2>1u88LBWd&*bYa*2&EMV}47lOk2?v`$5T-d3 z{}^(6{v#K)TN}*n^G}r<t4<B{r__0T=fu2QgSowP|ND4Nej`kfqhBhWt?0SC!<_RT zD?jdyjQQD6<;n7ozPKf<yXCiRbMu|n{}P{{{n+>|L5DlneO{wL#ykYgD`)%mqr^UH z)N^;u)7CfF#P>(_U%Mv!@#?c$AA5LJXn&TT`1x4LJd1<prrlC$VD`;P|DOY!z1_6S zTj7uI-g51I4KwCv#NbCtb?7%YJpZdYy+&rsI8Ua&@%)fNugnX|we9xAq-<*g^Upka z*5(S87YD(V0*4E4>6OesvDTz&!RHHuz3J<a^dG+U{HNc`x-e|9?87>Zp2*lQrBCkP zFlo^Iux;~8F6h~MPBMSnKYtvvEytX&^BdRS*Ks6o*F@vq<zm8>_VF%HTld7j)5D!o zHT~{R{yRPWvAaIdf8vYf%S(K<F6_K?*<(Scj5_)I*@|a3tsDl0-hFTQC0#St@%+C& z?OU_+!XV7C?Ank2d_9?`?FTysH`<r@&G??XdJKCh49(vn-5=mj>Z3b7CnE3mrdiS* z&A{eh3-BJWB}jNm{_>mZ&6qJW9Jz8$mS>(?7V59oVR@m8U!ECst^7y#B4gf4`aAZ> z^~IWM?}Ytx7o7Fp+y!xeu>Q{;Z%kOVGU?yh<FB1>cU~3Le5-7)Zms4d{Ud(7vC5@g zR)xP0$dkWGyH#P<l=S~TBmpURM*m2B?xnA-z}Db>U>mS4*bclOd;n|@b^tqqoxsju z7w|#wA+Rgh4SX1U1ndqz3O)w*0DFRugHM3Hz$d}pU?1=)urJsT><>N-4gd#&gTTSy z5O63s3>*%Q07rt)fTO_C;IrTua4a|ud=7jbd;uH}z6eeLUjiqBFN2f7SHQ{OtKbyy zHE=384V(_n0B3@;z}es&a4t9x{Q088GYUSF8g5_n%eyt(&ksLtc%s~<ds4&htuJX^ z?O<xqrbpe6J6t(Ec<iZjzrT6<hOq3&8l5J0New=He%+nrrmhWBcHKRDbf?sCaccd0 ztItji*YwR+|FPFogVkB?xHiku)Nn$R?I(Y(wm!Hj->36eZrc=0nEX`7=FQiK*<Sf) zLh6ace*;je=J#X2-yB|bXXA?>e`aI&zuW&^)MWV9FxNMG@9WoeesFSpYPlr`SBEW| zRv(f3;<aJcmh-BXYd9x(=Bb;Xo%Hy&aDD!r4cqQr5xn;P(N8b=az@x~LdDOHmfjYY zDL-q>yqYV6K39CXwfMAcVWZDZgw<YqC;T#X=rxb$+aBcGbMyEb7tRU~q}$=Y7f$-7 z|Boxe#|KvLcw_xAESXw;{g&z5!p(yY?P~V>wxGwxKHnT&wJm(7b(w~#=WY-Go^tz& ze`jn9-u|uq9lOeJ3$tF4`-^(H-VPTQd%fw|Uv3K?duvp&&Z%pH&Y#t4f6X_E{!E(M zds@G3!J>->75lvOobapbZ+NEjODn^MPyIIX&NA<WAMM|tyJNT^oV@+w+xO>wJDk`4 zubYZ=UlJD1wzKg?IkyKp^L{_@;-TBa`Xk%Un>2rWaOHJB)ZKN~>R|Bm`?h>~<=f#6 z#ov3d<;u5%*WP@q$KKSpgUe_B`sjvJ+rseEinZP^u{ii?->aP`SK1cLuTbjS%Hj5K zbK|#<blJ2$+*)+XGvQs!gCd`<z3$5f+rxgTMfRn9xHvdIt8jr2%We-#pImi#Q=i0m zmwg%D`t!<g<%_dEo&CbHV1LDT)|UKfRoJCXwGq$$wk$08d-<Y|b=w~FZIkn(i<Ykp zfBtl2u@?rv9jqK)yWqyguZOi?8TtCcixcA-Is1~Yuh<bDn?3dUi_dx|JaOv<9X>vH zM|jtv0-x>pe0%u&z-QM_uD?30-Ky3jtLv={?;5wN)@AqY2s`Jycv6jzw*@sH`hISX z9y`LG&lmdc@DKCD!$%jDDLZaQ82t0YMSW^)4@+HntjgXI+r#pOYTkWW>m9+k2El{l z2doTdwi@=>;3jW{@7}PWf5}qQ!(SV;x$LIzc7$uMec`~fXKx8!KYs9|Y58{q55BYP z!V_h8gwGYb<;eDFD}%Qt&0qUP{}n;Die<81@#~84o8phAc0RN{9C6|JLS+i}3)|NI z`(Vq0JHi3K-!N%=(H%kaMg!}uY_TJFV_x2G9+|KroV@ne)!Ba95w@sw^5e@3?+mYh zVbk@uT(ms==-aXtD=u3RZpi!jrN!<X5tjJ<?`vB(T^^Rnk$-tmY;KTtDJ1?~-`{_I z+yD7?JHxSEC!JI5=z^g7=n5ZK7`!9Nm+X6O>-M+K3m@pZa^E96c7(<iF8_FgWp6A9 zr|cg&_xtZx2QA9IeYDKGJHn~k-dcBH+RiiZwKqo}nLlQIIHb!ZSC5^&Gu%Hj=g(I+ z*%^*}Z$pjXlWk$omr54dllPTyN4ff6+_+>(SnTm9zG?N#=HR<`cAr)7x*5UF*Pni= zSEXIyh8>sgTldh;aP5TZ_s-}R28&u28e3`RvhdOoS$BUwDGYZv`?1u9RU3lK$N%=u z$aTBIyBoiIs(g`|LH<u3JE!r88Nu++#|~<fXusfre%}uryfc_Ew#y%X<y;q*DF5X> z^9Sq<U!T;X=%8X@@Iv$3OTBpGhA{O+k4iPJ-xW6KI(F;0AHuL!^^2}9|Le~1(g%-R z|8eVOVYSgu)#!Q0<X~8qPXA2po%nCrW)?o0_m&OegGW2{?EdG@pzzKY+bwMQTCk$l zhW%IEv?@IB%0oAFI{0$%%|B(nd#hv^jNJV9<eKO03WKorqhHJ$6)f&Ew_DbBGsEGf z2G?r4YF9X==70yX&Ric(*<Ij&W4CMwh8~{Srsayo;l#Y3<{h|wXK-rg7Y9~+zBbHq z&HQ%36Kle5$DS+K`bgp(UgVZdS2Ya7r@qbGs%5stVduGB7IgW2efYzUgPnKHc`Lkp zR{wKHzPT&x+xU;Nv+ubpyy2qu73=2R9o)Mv@8G`KcZP5LlQQLox7UZ&uDamY;h!!E z-%8mtcSqMH;o6V#e*Hq8rD3C{x4$%h?~?G1ho)7pdOQs8x#*$i&MC1m_^iUPfj=jH zldRJ9-8uCuFAf`h@Z0b^f7%rmJktJyrwY#x%WPeJ{NsZg!lL{07yshV^}#n6wXOSP z`60o5SM|D~-yMs?QR`0p@x|CBLDOtgAK#aMcCcyiD{q~enHsKtFyDf_xnB-nPC2@B zRQX-OD{Y@SQm+1n@QI<5H=h04hVbe(>n=L-{;uHr!Sf#ZuvOxmOL?tSlVuabMkm&W zUH5JXzo_uw@)_^%3Y(pG_QQp?EDkEw8NaS*%Ia{|2iczpKiLq@Z`&jIXa276e@&a5 ztao@rP_NC8#~KaU818#-(BjtFW`;Mv`pu&^{W?AP{Yb$s1<LITv;TSGk*fLD2ZztU z_ooZ;?+We>CjKz%+a+P$-wIqcf7be-(g$yxXjpAmxV_;EM=mb6DJ=V3&B1MkZV3Av zZMEQs+Pi}scVrzF-nKJ5R^_M5ODtFtmg-ny+^bnqgBL5T|8C6Y4dL0<2aUNg_wJy@ zlw6lSdT?{%UR%?l{W}}OhQ)Gr&vjsZ(D3-A)}Ped8Sb0c?T;(J-5sX=j|BOj2Hx7W zL$?kO_ssYo39ikRwrx|mlT)-KZ{#}<)IJ@fzekVc(>0MDeZ=VR(WCQpO=L$OG5UM9 z#GlhOksW=+=<nGRe@@p#cJvXWzef-C(=|Ta#Yhu9>QC4B$QL6`G%crVn5Hy;#7CdA z6Ms(EM0WHMqrc}){5f3{+0jRg{(es4&*_@Tjy_`a_j40}PS-?s^bw=A@{e<tCry#D z7TS67$L%P8^pV#-Ki!?ye)LysCBkT(BXoZ9CH@|L)EIj~x;yzFU#q(izqY#k__fs) zz^|?DBK+FwF2=8|t{{GGYcc9~^!Esz$3kee&3OraZF3gJuWkIL__a@OGs=DXJYI&@ zePmmg<JY#WEAVUERuTN#=DiZXwt0)<*EVl4{MtpZi-X$EQ8bolU!FdXCD6JT?L$fY z+V-Ioer@|u8o##tKWZz=9ijDD2CW+Vdli0d`&JgewtXvyU)#Qw$FFVQuEwuzKd!;A zZ9gjD*Y+%m#uEKKLg%(3TE}RgE8*9+&z13O+vh6ywNGy|%6<AgUW+!G-{|A`s-~Ms z+won8U)%9j!>{dFuE($KSZ=_t?O3Yg*LEy7;@5VLqOnBh^66u$f!00d7;EC!c8oXS z*LI8ner?BC3%|Bwycxf?W4r~wwqvY~U)wc^#uD8lr;q7Yw2?3RI3IP=&7|$vZ^N(c z*z4lgcI@@=YX{i%L2bu=JAQ4)eg}SS_h;1aa}s|dbnP0Twf>!_hWNFer#ta$J133s zYda^6@oPILcj4D|PMYA?c24fbuWkRMv7DRubNZN?qRk3Af6ZdH^Vb}|w)58lzqWIA z4}NXuswIAH=jvYk+Rjxg{Mx6V5qT0^G$z-jHQH#OMIYz&zH~EbJFji<Ydf!P@oPsl z9(|&=@`CQ!b`g>I*rsc6f4V#QKSGOYD*8k@&I9fD0}+w<sBvxDr@K?MU7HU0wOyNz z__d=ok3LbG=Y#gUQ$!>_j?H!KobFE1Hb)c}<%!TSc0n7BC;CKu^m#lTQXT~LM}3Gs z@(-n($+vp6p3&#AbV%t6UIN;8`EK|x0<Gm9*bm42g|Oux!M+%DJ$hnykNJyW%Rh>3 zzoXiSKGvcAl#l!EhjGrM>u5X1sgrjewT;*RM0$75c=Nh0y<)t6<DZ1j2VMZWw|mEU z{pRTti@Olpd4DRF&;0Uz;f>3W?HST9#_Kn}KfG}Tu&w>4WBJs{4}dq$I`cdk7~}OD zKL}pG@vi$|_=2Ey=(#imzqY(}F%;gqHjn#eSS-JJ<cG)N_1nJ@@P$G5q5F0uer<W{ z{u%g7V|hH+M&Z}C9rKNbzbux=yW&~=+O}i9G4Pki^0dVs8}pk-ejNN2pn1HLp2M## zZ~XJ{ML_d-ZoLrmn@4_pEMC9s|02Bihx^a^p8&64o#)z1@a|W2o*NTm@y5x&4DY@( zzwJ+g*RRfV=apEz`{rT%lViMo<6nhu3%Y*ou&039^7rF^4b+ymex|~^E^V>B_ol^o z{oV`H;axZL$j^ZHeD_|r{%6ALSLgiCip6^m$j^@P=5>G0fp^~T!}bj@H)hLwPRxVX zmiPRa5AV5SUiGiX@~d-x7R2J6NBKA4^&9X0SP0)1v<_USH)DS5z<OO2<MmsAZ^4_- z{PK(8jWfUdXbHT2>(2FF8sqgFzYJc#@vif7c=hV6j}`Flf9vOQ{3}8I>WqII-uiKU z-KX!s>vx{ySHbI7=lZUOSFg_bJuAWc=WeuXz+&KdFgv(5X3H1HzYfd+J_lY2hTzqp zy!!R{uL8XTa$;`)wdKp>PX)_@-VwR5H-g&oRq$^D?eA#pqTuEjZ~xV8iRt@eJD*!) zym9Kb!D~Ce#&3srzDHu00C&WA<JIkq>7T?d3GRyV#;MyK(~rR}1?~aW%bS01Oy3*3 zH27|eH(uR)G5s*?GT{55dU^AI5YzX=w%$IB@y4tB2wvOzGyY@va^TaT^|&wQxBlcm zfw!KY!8YIim@RKRpT^?tpSsWBuL1jlSA(C&{MNVp7w{FpA)xubjM?(Ga{#^~*atM< zS24fw@(1B7fdfGE9g5lVwsRQ1GWZl|zOQ5c0Wtm?Y|o2P*yj5-X3N{ocd>ZaU)}fc zgF)-x{d5Gs_Sl&I2W<DPZ(Q|1#^UA8a}?hFWBc~!r&xOvV)~zB_OsaTqhm4Ncy+(P z4*&;(uHUaQ|DYKEKWzQhkMYOhU5^pi1Hj*6{5WiNC$L>FdGnowx80%G*8lI=`bWq3 zKd|+?|BU|=-uWDk?RoSUwtn}${NM2UJs;%%fwz9VKfI6rh1a%z)SrU4|MKenk<0VS ze%RkES<(||JHP7Bg4dRJJ+i`UyT8<DgC9UYUDtuw+2OUVd+RAjjJJN(<%Bn{y!mp$ zt9Rd-|7>_|>%@A>9pjBxcMiO^^Y40`3vZmf`aJN~m;2K7$P2G6?|9FH*LHo>pAT=G z^P~Q(#5!b6{E6^x)+8S>K=R*XzXg8*#EsbM<S)dL2R_OjeHNudN`83#*W#Bi0ACD$ zc5n%J5xoAJWBkSNIq+WzE(Qz2SI1T-UkE-gyyJKk`x1C<dHY`&{yg;Nc^&)GSbp=! zUlxnk@BCg4e<k^fg7d&DV!VFGR|MWT=Xnb8SH^h#<|zuV-+0%j7`*FRmAv!8;xS&o z^H&1ixa){}2`m}o^_!;@ym8ftdkriNZ=5>$GVsP-PuxWCsu-`|__FZE-9X%IupGQ` z>g3DA8)uzPCjROeuiyA<;Po5t`d5IjK|c5KE7%p`wdJk9O7Jz&o97Me%CY?Bk*@-O z6MFN!jD0P<w!HCG;RE#Mc@z7(Sbp=!SA(yG-aM1AuZP!`H~t2A&mHs3CB8bme)GuR z7>n2M{;2_PeOgyj$zL<Z>z@XHQ;gT|{s`dBYaPxYz81W;{7m?p;kD(hzgys)_vyqf z0BghRSLgcP3O@_IyzzBn`PI2Tx4}E_@(amVH<n+W@%7+c7x|^c*N^2lul()s*1h9& ze0RX>SGSD(4d6ZZ<=uY`WBJt?e<!@_;W;pe{EgtX<vrgU!)wcXe%u8=jeOo;)>jjF z{p#%B-SFzwdHy$rcU<zuH-mS5yocOB&EfT{bA4LCJ3e*t_rOnOeC`wXLreH+@akO8 zd*QX^t>0F${N|N!9m{VX$A2HZaqdIwyiJVPZ~nIM`i-}a+rjI%j;xpaW4wOjAAr~I zJUhSb;nk~?KPxe>*%E&uygTtH`g9;BYr2`T7c2yJgs>)pvrw7=8~JwG(|h!)wdm zl1`^|fiI8!9+(q+5MEop5&T2&xv<{_%Yt3u>&1BW-QWwtuLd2@!|?j;zx*Tc&ewM0 z&DTAaU*3F=#^RlSb&tV2o~`6B2lj~Zj$d6*c=c<DbG{#s@y4lp0$$tsH@+9V>$#2m zwZSK2yz{57H@y1g#N7(^iSfp%dkS9tTH@+}ec{#1tM3P|eg$#2f&F8=@#>z2SHFq4 zy5InK_44Wm!mD3Job@^=#v8A0Fub<wZ~PGW+sXS5=zbmwuiyHU9|nI1aXUft4UgrQ zH{S^O2E?rd%{LNWzwz?Vz&9jrH)y_5vHbGp8x4Oaac_g>dlp{5@$zHfJ#V%UZ@#gy z{PN}-7mIiQsCy3H`dmkT>+Ja$AHu79A;!Bu)QyKX-v;8X#~0!Cr^fgR@cOMk<6nYz zzIPD69-IiT?fR;J8Ga*r<K-v8Yr8(`Ux9bNcacYaaxA~`>R*L-J>Dl?eo8FAyzRUO zZ@oLd_2iiful)gj_0!-z@7;g$)8Vy^S3d*Zb=gdQ&x4up`aR#}XTj_D{E(jwAClMl z^8T6wuWkR-&xJS6^I!cuc*kSB{Cs%V$Nl4d^g6t@>!W@FyyKHs{|3DK$$j8{SqSfV z)-w<4--Oq1{mL(b*Os@Px8Tk1_|3N%UcGhgeq0jc&981LytZ{~{4#iL>%@9m9^;Ky zw*p?<`E|Zm!W$>A{;b4&W>5TyaC+iT^m&_@$QOMsNr#lY_}>Aq#jj3&6?_r+Y~XqL zSA+U*iScXTug0GPJRkpBP=ECpzYabZ{wiQT{2}NatxkSDyyGrR{1@Pc7_Z;)r@|X| zDRKXT8{v&pC%*~axXXz95ZoN&^&7th-nh$&I|XiqH%^`WHhAN%AZ|anJ;v)deh0j9 zMTpCee<x_1I{97j+QwTyyWy`PuXXY{_8xd`dFy8{d<FF8Ie`6cEWdf=--EA+-aH>+ zzYni1Z~O=FmC&2#B=(1~{N|DW2);6U^L&Q=F}$|C@%!N27v}ks_)lW_%_F}*7O&s^ z^C^5)@_A0=!2cPjE${yM9R51u&GRYt7x3Eh=KB)98hZ1bga1G*zj@@pg1?@4^X$Vu z2(K+~{2}-o(3|I6{D))t%_ILce0Ac@lMVkjptij6-^SwAS^wX`d+yezoiD-fW4wO% z-w}A@ZYAz7@CSI~)XDz{Z(JSXJ_e7*c>Tuz1aI7J#HHZ>88l9v{4sdr>JtA6_)Cn} zZ~U+D#(9tA#Q#6gICb*J;kAu-|NI8;dG0xV0sa%9{!{cjKmL=TdUf)@!+WkAApR@x zj~K82ApD;(UjIq>zu=w!<NCqB;q|M_jsKrmyz?&qFTDHaH}Yn|e+tyE&i$G8*h$3b zH|{Lrvhbjh*YEy13;qw<b)8+$tnm8P<spAIc*mtqK0CbYx}SC~Bt8ece(%YP;B&&O zS0|qf{!{Y)0y<x3!|PXfHu-bM;*E2D&x!G#%f_Dz@4aH({7yge#CZMAW8N6A|4-u2 zgFk@n{hgKmpAWAs?>U(dUR&Pt^8)xU$nQQaK>iEi9k)8ipFb9_&ik(byn5r!e-XUr ztLI`a+P^r)TOXd^1>yCZ$NDG~i#M<QCGcM{e&^ZwEDV1TUfm(=OX0QUe}umbUR&P$ zm&fv(=P>bC#CZMIQ4x6S?jUjUSHc@-{;$y$h1Zt<8@?F4e#aqSJeJ=)-w<CS#_RtH zzGRHoZ~sfd8)yD+i7ySWE&m;S8F+1Z`+rp|zj^*4zAU``vLDCb%fV~Q8($vYy!PXJ z;;)YJ`j5b06XW%pzXH5*=Kq2CityU<N8u~MYs))7mEmpIc77(lO00d`IRSqyyn5s1 ztH$!H`ycVw#d!VRx7A|2e&^?Uc>Cozy$^1P@%n$I{pvAZzx}@v-uZHWcn{Zr*OvDl zt_iO#@A+^Oy#2EO-ctd*w!H1tg4b_9<!^?!zO7^L?_1#AkLtXqYQt;G8-FXj>*o4; zpVxuc@4WiHyA58wI{CWr_S<#$eN_+M_lY|B`muO*_WyQx*WGpSKDz_neyj7n+8`Eh zocCWtc=g81-wAJBJCELDjo`KAz26$c>vtaI?}GQ<wH@;{fp=czoxi)`9lzsneVW4i zp0Q4RUp9;J`u+aU9Nu?{`JJB@@V+k{pY7iR@BVUM`5tN+<Mk)s$BA<wN8(Q;o5`7Q zFELrL@8B=zW4D68F~+NJ4Sx}QLx`-{_rYt+*N1NdUmUvu@!7H4!fVUl4BrkuC-$Ah zUxj@?e4`kz{sH*L@V63wHg<b>?Yi*lJHS`Qe{C$jV=UhM>N~}F<ExOTGraaKvHV@& z_1i!D`yhN1##4v*+}IDnYu^s9zH5xX4!$bb4PLu;EdRr?_-gRxc_hZ$KXu*V&1?Vc z@1yY6YfbX!V|<Urc<WnT4|w$#6Mq5mJ!8Ca>K=zzUxD}wiGKoKy}bHf@ahW@pP%?A zW4!U|dc&)~lK29|_kmY0ul^}`^;Zz@e(oFNjaSzXUfca`e1G^7^rsy8-S1Ds>vw<4 z4}dR8d_m&PHxORGdE^Jdmm<Ct@#Y&G%P()fA@HS%zl3=64TVqoWBf4qGQ^jMH{bAB zetGkafcHErLcIA#!s|C){+U?3=cD>j@HOaPMe=*zkA~Op{*iwcz9#XP5^ugS@cPXo zKNkKb;ywS&H!hZ6-h9u&2gG|$n(ujd{l?3`0AGvv%jm!P#>4A3kNk`9HxnNaZ@vk! z{PN~|DHd=2tD6X4pZqtG|90$`W4!xM-6VMRg^9m|_*Y`Qaq1?+tM{C5K>Vxl>gCl> zfmiRj-;ns%;Ej`4KNVj6<@Enf;-|%U<JC=vS6_p8-vcw?)yu1&37_=G;%CA89`OFX zi~h}q*RDsr`Z=+9dEW_h;k}RR6K}qGvHW+ztDhg^?cZ(aUWeB<-uw$<@wRV2-+;G0 z<Gsfh#(3{<^SlY~{dfcY@}651<D0;%dkfxsS>Ak$;XT(&6Yu%H1YUoc7{3%=zxRvr z%iz6VixTfWxje=fgIBjA#yf6xE8)ElsuS;c-iG(wzZzcsJMj9AmtO_%{`H;Zeq0Uj zxzLFG>es~L<;}Ae-t)Q=@!kXLV!ZESbzzM6{8qOf-t*LVu;=s!c+b<uw68uj7B6p} zjqsi;)rj}J*aWZt`WU|%Uccvv@mt`nmurZ>0^AC(-*u4R2Cv`qUw%7$5%S#)S}!|d zy#2GkJK?pBH{Y&UyzdM1?1s1A9f#xH18=<-r+@1A!s{;?<KKnXZ+_$7i{&pty!qab z@t!N{K7e;$%G=)$;njOCn*SqsZO=XX`!T$6^6K}&+wRTu-}XL%*KZ#A{qVNu{&b!` zjq#q(wQ1)wcx~g&_c^@%HeUV<c+V@(AJ5w_;kB=#UC-MCF}^Iky02oq`PChSSMUBd z{~>s7&k4`@!!h1?bzj5VUfo##zKOLj@A>jAyteu6?|1Ny>qhe1-uLkOy~pK`!0Wet z`5)jNkMB>%|0BHTh5O(8<S4wh{jk43#rRtEU)|3!-u&v0!CSA^pXc2#@ZJ}`11i$a zukiZ){v!WBc<&SUk@=3t^83Cp-*527$s2znmj71zrT!$m_lM*5y!}1K*CAfrAMo1F zpYeafdw$qI+y5&Te^+dM{)TtG<gKTF;O&>Z{rMN(aeF@5endyV#YO0E2B(Pejoh3* z`J0yihv~KD@8->&g&!OJjp!cYeWRWQ-vYh7`mFHgYf1b$=(E9V%QuD34&NNR1@Y%% z=YZFiw?8@IeS<e8-ucJ{uix>@pAGMLn-Oom+_C)f<~t`A@BFAc7v6Rqzw?y`-uW<Z zGsd4c#v8BhJb3lC>wKPXJmZpApU-%B=hyfP;PpE`=l{Z3e*5Qm^T+Zh<7qFVqrYZL z{F%=W*9C}q9L-R$05}Z12;JbAe`+FE^cjQyVsve>&D${+R}ftr{DZ*0;7qVkOxF$n z^)dY===7`e+%a!qbiLBMozk5=?ZAHErRbi-cHEB3dAST-HFUP)cpS&&=yuWWBk1Gf zzapmFj9#5_MbO3jV?SMwE79Ffd)Crp*e`%Z(Ph@&Q|OAJn?k?(gQLOD=!(bsXI)tj z&R+?1UlHGd_y@pe!IJ1&W5;z5;xC2Hcy+Dd)s;p!kakDJ@{YqU6KmIa^Y=t|6}rOM zJIUJvyEj-Co&7P7dCgZ2-LK?v-rPqM$X6cSRCI&!Ter^h)#$7v>(u>XonC{^JK$mR z#p5cV8$r8W&^eFJU&WYCejst`Dxuqle<Jv3Y&@0GS<k-zt!wMK3cC3GyDqPh_gZw$ ztMeA$S5?vZesJF6=gW2Is*>Mx!2V6aUk#n-hWT6<*W-G0qnHo(q4PbLd^g1Mna6hR zM|E_E$!Gn!E~ChIBRc2VJomtRe$+rWoqSV4>&|niCOXIAzpHi}j_W3Lp3k1gp3kq) zet^z(vhE*(x7}LkoL}z+>$yMrn`1it>RgXo&`o5V@%kJ~zS`)%LU#yszCD+2Md$oE zj;ZJ!M;&zLF|YYNCvQXNc|4Cip4Xk}U)`9l6MpsT>Y>|59`}dio=RMObgAgPr@Zg1 zzuVF2SMR-J+#Tq;kk5NMUKb6}#p~C3d6oV(MCbW@6MoNw-q?4dyA|8Fm38$p{zm9L zC*5z>i{ovKPQP{jKI5{^?uxZ*dlPBTcAKEH?!6Z}knaifccZI>t>3)XRa12SHzN<D zkMGN7=uXnk08pLt-yGd*@a~gy$ZP$yK<E59Kelfj+=I?>JRKX4<7kPl3cUNsx4--H zUUc`P^L+8YdGcOuh3;N-eL&*|V7Er+`Q*7+9A4dh=-y;r?4SL$e{Ik?52MKA{94a# z(K&C{k#*&H)efEO<@$X@KRk!-N2k9Jym9ItK=%~9`P5kt?a{gZuKNt)9B&77lkm4^ zyaDJr)DfNQ;W{|)&T}Vp_RD_TKj*PCx@xrJ{owuXe(r+q1pTv4tQ+g)L3CBgw-UYQ zz!dxsq4Qj`|Gm(A&UHm6@A!(5XE@jmUHm>jL_Wv+Fgo9B1K<b7bdR9(e73$kpRKFz zv38T~&@b!uQFQOnPW-%m4gF*2+&7LhKK>r)_LHv_=)CnJUr%(tr<>uQ47!gVNB0$Z z4uXe3>*on{)=M4ov_NOQ^g`$FP?@#=B)V4QasJ#tzKeUKbAMch-u>Zz>4UBv{jol+ zt5=A7D%P&;xz4_i`l36VxB}o%(Dz+GbjF*{ae0pQN7vnX1--wXAn()YrosCzvOUlH z0qEKi*Aeu5@SYop?pFAEp!d4>>L7G4)2`>U=e6h6V07)!xsRP^=Y0q|_l@VW^Yj$` z7>aHm`Mf9N=jJeU?hp4xQSx~23`h4b`PzdIfzI~`bpAGFyyvRt%}8|4r|aaryB^P= z8$dhmn}OKY_b7DP>327H_p|$VG`eZ<u7AAUXVKZOr^y@FjX~%59Ix}}c*ml1-CQr% z&+~5_I_qj8?OI>f)pO{4KiH0Ss_uDoj@P=M8XNx$=sb5Ehu?eW(2w!x1`+S?3$Dv} zbT6VigziW1E70#16VPqN?|Z{Ko=d)$(Akdd4kXTVZ6Z3?$@+3$T+f$dI@h-~`G(T& zBy`S;-}|e=d;h)?(>ZVIoxjQG+#l9Wd>vm!=Y2Vd_N*K0WeU3ZJ$MuP-r#HKmN7o( zXIgB$Q_<BX&i*@o$2AR|?`(PREBiYgo#*u<<aHdb^9*#ptJ|ZCzhh^jdm=WD`1h1q z=-mIF2j1(R2eZ-ne(|0eKz`4cIp~}(=g;}{yqb&dZu;-NcRjZCJr7+W{GJyrX?G<4 z`RK03c7Dw}8UO3(0&M56F?MHg0lKBMqfUMnx;N0-p6xV)_nos4o$KlPS{K&Ao9LPl z@4l>r?Y>-u?g9Mn_tw~+8*ib@jlV1S9QYErIMzSwuob-HTY^r%^Xs_OEk$RYm`|Pc zybPWF8bW*a*M2TXcOAU@$n(YdTY>H?)@?L;=f`<jiS7$@-pk(8o@;NT^IY_uV_h5f z4!Y0C=l37Kr#TO+(CtI_8|e45sq}X>x_a2&Zv(Ksx7NgT^37wqwdg#zJJRm!;M3^W zp>w|6*RHqsM~LoK;>~Lw=Wji_P3XMWt$*ua13LF%Kz}X+-G8a*+#l`-_lNCoMCb2m zoyg;Ub{}s-=RPx^`JMO8=ms+m=O=!@Z$X!b_Wa+9#ed(~itavi_k*5qexKTg&U4#( zYm44_-j1#u{*m<4ar)laf$jzPcfowbxgI;wc@Nkh?}O>+ccCke?f1r_*p6p6I`1RL z<M$EgeGfY8$#ELz`tC(%e&^5m^j-EYy7;<V53}g!d+6N%o(sO?o%i?Ajbh&1*W<9~ zlJ5g_j?=v6aos*d=l2)KGn_oHfFGfoOFNT5_lfKLF}m;JeJ?pb6Ues@o$bkc&iOs+ z6LeF^=YERM=YDkl4sr;+_tPla{S=*bYQGD@Tc@9)vmN)p<2LSdbiQA#(~HPwy?%l2 zPx|FK9bd;U(fNJD{Sbe5A3%3I?O1Omv3>u2h0gt;&UsgN5S{zYem_cor<3;(I?omB zFMhrpMmL!FfuQG>-wVD*=lVN;@!#3LK^I@&gXA4S|Gq`H2EXq(=Wi~$@6de|>zDVY z`MyVY8*$CRyTHld5p?c1&kfHJ_uCKX{9VC!v-Rcg3_qgt94J6L@%?)g-HVLFb7eHP zb@~%J-vid2b?7?&jLvtq?<@In^y3)1)wJXH7whIVbibgpU-rlTy5D|9=e~45K1qIc z|3l~edTz(h=i}&%H(x3GW&Qnz?p4}xy<W%mygz}?emlSY(K)V@=sd^lpZB}>`S0jF zpM2j~cfQ;IK<7K%cV`R6;d}Q_bidOd*UNqC`}Z$&)`9ni^<doJ=;HUI_2<3t54!)+ zuKPB=pZ`VI51n;YgmJh}PoeYNU5MWO;eF}<^Kk304La-3dd!k7op9eeU-A8W7P>X$ zc?Ps!bLnSRbYH`}euuEVXR@L5`_QArJC14C+0j|g*4r}jjX{?Ko#*{M`0In7>p9Um zFWbq}6?9&5q1%Xm5P4Q&Pr*JL-C_KWJ3ijr=;G(qzr<NL=b$?Wzx%@PHG^sQTy&ez z`Ck77+j}t&I?q$*cNRL&)4b?>-@9McPomxPV(t2UAvM;I^U+;Gob@~mo#$jebQ|&S z2AyZ;{Q`9M+uxxbm*cn)opojZtuOaSesoXMp6_k<x${>5-IM4%AID&OPF{r0eeON& zK7WRJxES4F;s?g^I4=d!c}`j%?gP)wLg*eOkN%pp=YF{aT}yPnZyv@TiM}wpY52XD zZP(wKFGc4*Yk}T#*Zp=Gx-8`1htBnwivMzS-fzyg-;=!Gu0SX6z2!QnD}t^A?Ro!s ze>t8j(Ru#6kNkf49Q`Yb?g#SuZg*cgj$-Kiz0>{SdvzxHilcL%o99V*_jw6)u3Ki` zLrS7^zm+GC`_27V3Z3IP8XJ$}C>?9pdMQNS;b0kb1Igp}ZR@)Ox~tH6f7&14QQjYA z(Y;DO&-rQC)?+zz?vMC6`7-U6M^^^hdhvX-Uam&Bgm!$lOu<%v4Z5qz=Q!NAzIQ92 za~`Zu=f`#{qMJxN{&yiWu&swm=&q&RyWy)~TQ8N-t%R=%-vjjCs)DW+I)7L8`;6~| zYthBe3F~1xd8?vZLB16H@%!mIbk2kIF&e$|QVpH^(R(13{4b%u9$ow#ssrzP<_2`_ ziT7N12;1|uI=al>OOE$Obnc($Y0q=PeN+RT?~M;=uQz(X2h>FO3%u*@xiAd<P3XMm zou54L-dh1W?<?!CIy&pM7CP6{d#elmaJ)C8b6>s=@BZ|jxdmNz+VOk;McCF=ZFHVj z_A}nUThYbi?Wgx@9dw>|MQE=d{dAw+hOQHS&oj?g@4>q0ypJ+FKkA`tM!o^GYaMtF z)kilMo$vQa*q$%9qg#XD@8!O?Jtyx#cP{ODe#|D$`?UeO__^DPIP0JxI^WMF;cJ8b z4ss{Dp|sNp-tT?x-$v+M*Z6zMb{nIM*Wdl*^>?nj(0xX`{tn@|Mxbwkt|+#7U1!(# zZgj`Umm9tJz2{3)bWPxWSNYDe?wg@=zj~f{KTROMIXcf%`xoDDEzo%%ji;a9OV-Oh z=<bDY55~V&v_z->Lv+R^`-{$Z?1l8(@2&2;R_F$hrw=#~^jvR^?l$uHz4AtE@5}qp z{Y1R)G0!#Er42gk+jG))l;>_+bp9@2d{^3aecPe)d&Fh%?n}Qv+>dSp{qr6g6wCJj zy8G~3$DUJj@V7^2U0Mgl;D>`9&^e#plYYN>7F|bl@pW?j+((_zonjpOLEizcZ)bGQ zgXh&k^v+8cbe-sj=Ysw8Ja`b@mFOI|=dk<oA#{Jzj=!6_&%NKeqH9XN=Ah@V=WaK2 z?n8O!+j)K%T_4(YzWx2v@jQa=Z*+&i2WWpF*d1Mbovhnw=pMCQY~K-0$!k47hR*Y# z0ea5`&x0Q5yzhMfc;36;dZIf)JJwBp@>(yCqw~G%x$Qk|z9-OmpL-AYCckmL(0NXJ zPsYdlWUO8D+OGL}qjO)nuiclvcl)68o^DRR+@F45eG1(h@Ru?l{@ySKyDz%7=-f~7 zebf(~_l)O|?>76}AKfVOJq=p#-uq9ZvyR;V*4ZrjHvpaIqUWS>W6=#n=Q(M=Jb&%) zAavG&_2s&_9)r;(e~(}s@%wT}tlc5>$Me~ChobXd_?>+5b72^|_<iO*XFU%`SDyBY zfYz(y8iDQ{_<gkd8|ZwFL}&j!=dCNx`)APU_kOVM)Qv({nRfkNdpGU+J#RER*U9&f z`}|enpG8*`+jGzJ_67W7(0R_9*ZDU8SahC~o{tCUpXcN_bnZ9zx#OQfyU(HXKCwQn zThGbo(b>)r+Or?F`vSTL$?yGB5POn-boR&l!8-B1_98m>cX#63(tc8p&ikktx}Kox z@e(@M$$Q58!+D;F&h>ZQt&1__eHq>7^yhs1eX+flC!zEG)fe53;AHR>boUVNdpEhC z)uZ#@lMF!b?>W}dtLVxQ=X!L;c0Hz`^Ic%Q7DH#<zJ_i!?fJWX8|>cbr=l~?b#k3w zK{pLuTjE_`-{ZzjN9Q>(5&gO3v2JIe^L_8{b>_1_GtoIO)~EC1yv#y3nf4ZguCwbq z8=dn%5ZwfD9Q~Yw&Uy7cR04iDI2T>Ko^OM<uIHiqn|6K#eHRTxKOdd@&3X15^xS+M z-E-t|oQ~V^E<m>)T@TvXiS6$yZ=kyy-t*pbzCZqj=++TG1;2UC_a?ee(E0mT3VQF) zMd&WU-y3}?Y}ey0bpDR?9lAc~d=D;0=e)SToge3A2|CB&ed>4|$5M2jYi(&aeqSy_ z=e;Vg-u<y0T{Ze$3#<cr-mE~koA$1Q9|X<>SEBQKhvS<^obQ~s(cMa%-vhiaJwM+; z=kII&@5aw2&iCjlblyMiYwsuDO{>whA>MbI^)wg%8g!ns&V%nL&)v1?tPksD4CA$4 z)}izJY8BeK478six(o4l2VHO1Z9Te=X=g3C9JGEmpu3nnH9_wQ-*>6#yw5x*nv!oC zaU0RK#{UG^5quimCUo)h$#ZKWy3ObwA>Mw*`?m#MapJsZN@DvS+lsC%e$Qja<GO7_ zXZ<{h-a4C#y&c_Q`se%V8|+c&cA(3T?e{X*(RJI2&bqhWD-&nk??Pvt+b{b&4E=6& z`i*xV_<ed0y07WC>umo$ulAzLY<^#$-FMM>UzpGQ=6eraY2tj}xh}59`{-K5bRDqW zM<1Z`KJvWvec=1=Lv&?{xBYm#AEC4V)`|UimH3a*`F{C;dE5^Ae%XhvHF3U|JO{_) z{{-FJ#9xlz^VRcoKf1Z-8ldy{1oMB2&T~Y+{czl$q4Rg@_sP2z^c?&gT?^WCf3-k2 z2mcr7{5_~8yzi)y=)OefdRhne&v`k3ZZ>)1@3eX3`wE@kbKIAQi5r3LAiCSJ>w_gh z$8`vu@4n&i50KAy=V5fdo2(o6fpz~iI>+bwk0!6<{RW-uZXI~PS_j{v%TK)LO8j2< z4xM#U4*f{RVO@NWt|h$ZhWD`N#u0R$PlM543;KQL2Xxl6<B7iyenfXA@t!k<u&tM) z=-Sfn<U1U@ANrrrc^>=w!$0U|p!*qJX>9ASAh!GC7`iO@{oTOt58eyEpz~elJJ)xP z>-{S_>t!H$tsCp*f9O7<-v>eWi}idQ-Q(nO9(?Cnr@x_VMjp>ge?N0xPM{l(&i8dc zY|n#}=saJnYwKwi{rMf8?bxpSYzVqP&^5*`1DenE_!FJ;DR29>`xm-)<gvZ`u)Vkb zMmL`ExliKz;~#WQh_3|J1byfHi_UejU-5l@3Z3s9-(8+Bo;PX#PlywJ^Io=|t#8*Y zOLhi8e{zHVUilLCS?F%Wt_Avz?Teiio$pECk$=%lS9ICXx&H^_e;J+UL3VVv!h62l zjP3g`2fFh3F9GeJ^_UY~Hu%Z37yllW3*9RCoao=h_Pu{Ly1w{72dy*LFE={p|1i3r zz=z3y4m$5Y?-B1o>+)Q5zUzAtZ=F7aod?~swBtPZJAv<^yy(2wiW6`B``$heT^sWG zUh#LB>BOCn&VA{Ae-Yk&nGaphnBQ}CGX4wD`JPfA*IkJ2IQ`2_ymjZfksn=kc;_LW zuK+r~zj)qy9@>wK(5<Fjzh}EHvuWpIbk@}-^q+v%RY7!Fh<g+q2zq}MLT4T2fp`D9 z?=L~;?+ZWBp1(WvA+9hw>&W`D?p)tX(Rt6DB%k+<_snH6o#)$E=-nTeqbp56N`c<1 zuFDnZ{C(B>bUm!oBIs7o&L%Jw^mm#o(RmJe{#;7?o<l{^IS<Z{?+ovSV(9Lto$jFb ztoLehboSdivc6uW{SxS^Vf${r4to$-5?u<m>+d){$4kZfXTRe2era_1$mhA|IcS}h zLFa!ja6kGZjQbhvtI$ovKNGaB{avLjx^Lk<FXQj9a_F2F=fQcINxS9I-Hh%1;eG7= zaW%U2wCnGP?_j&1uR(Vye&;Rz-&IsV*B)JF`ikhh5BwhDe)e9dgwB2I_jJF{dmdCq zSDEqpyMXoWJzE8x?=|lO?}xd>UyIK5_4f$Z*Y&Q7&VAvz=lL<6eAl6Cj$IFI3c7!* zp$qUg0OvBFzTd7#XM5g9_QQ3%0o|v>KS@64Zy0uUbnf$e@ec>xzc-@0hkVI@KZ5PM ztOh#!WB=^0^HLL?_nY^T>)}0k6S}8qXDjG_bln1U72rLm4`BQ6&1#_=f?s`HcQd*+ z%!B>vf^ECEpnHydevcWA?YpQpI?q|_)jIXuy%pVG#Cwk~$M#;WgU)*NUhux3LqBgr zH-&iXDgHZ1U3C85eHnW9dt2;!=;YP=eOg_8bOVU<zO0Dt`rM9Azj}Gs?GAK*5dRtY zC+NB~KsO2hEbvXxcUePpE%Ez(!uO~3e<!+i=z1`(>#^N$jnECk@4avTz4sfV+l<b7 z_V*;~`mR|2{GMNd_U&I2bnZ9zvG=t5?QV3wyYi67`Lx}p=&UR2sX6_fL0mI*xv;(W ztXu1}x%r6qeb52h`?v+Vi^+Eb=(*x}??Lw@@xGG=V~;1UCAwkw-M@Y>@VvPf-9mKU z%ht&(^0h+e`oDnR^UXSFjcz^p)LTDOh`SG+?`+?r-g7UbYlF`5Iqos!cl>S9IiJqk zz2tcXY=^Es_B~+yeR)4R_qTa3K;I910G;=X_w9`Me4_Ile+J(A_dM-@&VJRQox8xk zU`KSmTUL^<Fue2L37z%qJb#bgdhU$wZsM;8i-7&XF6gX#?=|b+x_=N|U)u4$@%ysp z!9(a=56^Me!}mc~bdJM)>AS@CyP@+Q>O{Zi&=2eRVRX)ub?JR$JwJl3HhHRoEy!nl zcXa2$Z-d_hI!}+H8;9R>djz(1@EE$5<n!Ol_?~rt^g!o3(f5tt^QIEt6P@?+i?lZs zbe<kZcLj0x!S}}YTzvwa?-$>#$?>4?h0b%$^V{>ybL~lV{(D-_0mn0ee7(_~OTPx; z_j`})*aw~Gp7rbd*Li*lo#XI(*<|wiy{s>~W3=o2VY`0s?uRaZPI!Jej{fLE@>y@I zvE4^cqdS1#{o{Ue9}PfvIsNu~(skG~!GY+SVBZW@0<Et>=x)JZ1@!xc_rqXx@%nT9 z9Pbcxs~Cs(a8B}i{|=4iGmqz%`G%nzN<QBm&WHElaCBK{XCUbJ&}XqnptGJ{!teT8 zS0m9apk2@T`1$_~I)5)7iM~JBnf#;Bc@I91-+sD(N28lY{4#Jl===Fubmht8dG0yx zei(z!eedtj&WGpZSahC0u7~}#f8)@7N`GDgotI(s=Q(r};hn#k*kjQ>kIr+$daVin z1o#5Fi?B1hm&T*JoOXP#@5gQhzKG8Eh3}3B$?v;q0y@vf9>jZYdS1SSZUuY_d8{kz zY$7_(d(WZf=wAe1M(6i~ugH4<wEicV55MQG{~p0}_Z4*h?lOvYJkLGXCZqE`H4&ZX zh4;~`=&Yv)&{<d0u&1CKO8>@!?qm1wYv{5Q*ArbguowQR=vrVm1MdNSS4>0q5AlxA z_k{1E>FC-J@BQKa@&1^B&ig37uAUn+(LF%B&cAi+`*aq%rnGAvG$YRYXf`_Q%63!H zSyyw=`5y3G@tkr0&P8_(?OMM#VZQ**L+5$#eI382=cDsI<@X}rX|u`qI=W5t+wt6n z?S5H+?koJ(|3T~y<a-01=RtS;{+?y~3(@(#Cw~5T{=JFL^C5nIw59z;=qh0k0G|R+ zf^VUVzl*Pd_kLfD&i#51y!B(<EJ5eF<2lrjal4L7(Rq%0etVv~f0v>2e7ul6o|De! za&(>>o)<sSPxsY|Si8O(d`FBW|4MXU67T!@Jb3r{+vvVU=Xqs6+`sRjTZ?WTxE!=D zSE2JgbpSrT9;?y$Zg3vmPwt~N==|Q@k9@v6JwMi>`ySr!F+XGbJ!c&{=go7~`(P6D z6RIQrA+Rv`FSs6^_loz8>*0C60iEY$ZrbsjbbV9Nm4vSk76(0dH=^^M-51`v@m$-4 z&iV8Air)0!`*$<CtmN?><Nom;+=9;c*Nf!w{pCHg72PWMcR|0O`##-<t`ofb?n?UM z_oMCTPN3@te*riF+=0&ev@TzR_nh3BNxKEeZ~ME@`EKz(9gEKK?nc)fzxS8x;ri}D zw-ntxa1Q7>xfk6G{DbMg?_~G+yXdOJKM0lv{oeH+I`@V9!TsVse;=Le;d{(=aXmgj zH=F*skL~|F+W!!p^BT`*J0GFjLq6Xn@$>3qbS2Q0qn(1Fzh~`3=YI5j_WqxS{u6Xx zk<a?K-hF57N9R5)L3}ZA0Qf07^SF=V`94GEIblCN2R%1FN0*QMIl({aulwi=bl$s} zt;?6_JfE%S`0vUGV*Rr}_HQu#_zK;}^vC<bedoLBAi6fh`MZYw_x*PWoptW_ovGxv z-NWc^C7<8B>tXwz{u<p5;@5+o556bALFav7ojgE3>*ZT?Q^?aFEDDYWzeDHy+26N( zKYLz(kIwtRd5q8B5p@5O-*fLnY<~~=0bM}98-ry)&*>l0x!>B5$9?Df;3&G`@cypm zJ>vcQ6FSd<fwX76IDbE*b6?wDcJ#J;44wPT{S{x|U(i(|-a4&}JqG+0T`uh7j3@q% z{vWzq(D_bqpZZQXj&3q}+~?Vdn+N`e&i^h^z30AhC(vympWiz+V*5MFNp$h=AAh5J zo_2pnXPw9E$-4goo#%zWBY2Luo_}VN&;9N`{|lY#@A(<uM}MQ+PQUkpzPEj6|AWr^ z^byAQ9B8}$qWg_Jet$}yALvh^bDw)IyAJj@B}aN9zgNf)CBJc5&|O1*?@9YV68~A~ zeE<7>D*innE4q7$^Ig*i+wT?G&{?0>rQ>qHWk+Y-`Oa`W)@crOo|8RjZ$0C+UUH)I zdyf6GpY|^ox_BLpCGRZoY;=vVoyVK7UjTEX^ZVs2`sMh|dk#9^iTR25y!738E;{$k zA@Y3<x{vapdk(+nd%Q04qRT>i?whl)t*i6UdB6F-n@WD`^?Y>Yh%W)}yW4k8zW=XI zvkLRF>f$I$HxAw1-Q6{GGjumYhcr?%bPO;<4IM)Xf*=wK7J>-|Dhf6tSlD2E7vEZ5 zoaf={_v~}liT{bc_jzaF^gCJL&D-7Q2)HuT8{f)eZ$h)e>4)??`XSFjHn>~N%k^4D zKfBQEaK<BlPdn<3M>*hJfBOmh&k1Lqy^ww;qWjQXaQ=SRKRNZD``mER{0{H!2<)Bk zdEmSku7h#iar46Mp#1@K5&guV`QTb&yASS*d3%01{h#AGzT*~v`<U@Y((fPGi?IvB zU84LfYTWgGR0wW6^`6TFY=1{#IQ@5f#tHp+5xASdc*bMLEefX}(0{rw`h#L{E$Pp{ zf4q#HjCmA?)8A&L-0xR@e<=Z1nEGt=GaSuLc}ci|l=nvUE1Tg;!F8v;JLR5R&ueM8 zBslZRr?I^sW#ER<kN(y;q2DbFr~mQX8;^|(<=}=ej=y6Nw)tdvIQ@ok#&fFQr~qeP z70o=u{Hh|{0>+6)=b?U&s063~HqVWqU)Q@b+)mnWL>IuB|5bt8N4@b!|DZpp3K!-P z*J<ZHs0NppaZXWR16w~)9quIcu8Vn@=e`D9DY!iJQyz_>ye6D^S98jJH#kl$xS7=J z7yO?3JLR?ETo?1wE!4X{b>M1Kz8B5F{5D|Mh4VXv=fu3k^HLA42km_ic~254uMcPb z?R%mZ^`7SjaQdOs%twEu|7ZwjKkldZYCHWlg4;wt56e3m!&RnTKJ+l{Jbz8#{QhFx z8U$~=YzpUo>A&4C{dhAt*Twble!4Er;f67vOlV*9Ald@XcUCO((!XrLZV6W(dk$I~ z-Hk@VMbf|i`3$yxx)t1L%5R`~(SOj^aD%a@GR}Bx{can$*_3yHi^F!lZQ;xZ%%8%% zsU2K;`l(4j+t38GJ={EOzo*Q{9!mQTaK*9fpsu&;+!3xW^`%hr1oMMVaJ3l6c{)Gu zWoNkaaPOkV|84Z&1<w4(d(#izc-$4v{MLA*KhU3cgX>Rwzc)Cpc}I6R@6lZ5Qw8-N zMZrZ;u78Qd_WbvNGtY6Ho6JMM+Y_!R?H{%;z2MxBlJM?}{;oHi^XUfH5_P_P;EFJx zt*Gw{{cvA6<B9uZ9N9;|{owQij^ldh5BkFyr~hW$bo8TN8UUA#`d(<54-JHip+1aX z(<vVW=ehLWY^Pp7JQ(g3%JmcajhW1I2%LGl-x0#Rb0}O++C4_UJ<$e~4};SW45wWG z;CUVnH;eY}*CK4&jezq$d_cYN(|b7*&U<El+mG?x?@@5x595pRrat|ThBF_M+ep3l zC>pLO<8(%izsBP+aK_VyaA8~>3wN7w?w~W^<I!<&<|T{aR-+@~#=~`^d@$wn&^UAg z+?#N|Q=i55{W>w2pY!?#ZWsMcg43V6Z~9gDZ8DtsR+xu9px-HQh3LOD<uRynYARf~ z@7{C$)#GsbHT|f5RX;Zk&U2;T*01}nogR$4lzG~Z<IV^;=Or&U6K){CZxGrWHBQfh zGcI4H{v2wYo(*@3a{b_I*jt(J960yW^KlE#{hbS^KO0BAc||P!Jpt#tJArb4x8u%( z`<H%vr(}cooiZP83EX;g4mtvU0i5T+@$>`w`-O1&?Q8U-U$@;Nxa{<Q4=w_G3K|2K z0o%Cl`8U4D!Zm?2UbezEe_0F{MY(awcxZfD0_Q#Oo|w0J4wk~%zB|7|f8~2&8C(m- zZH(&Q^lQuEDp21Q^}RI)dj(t##`AlQ@hb*yC7garzoox2{;q;+K)Wy=jbfbDaISN6 z_^Rk=bPZgr?NI${5?maddHY4UaGlq}xz7F$*V}zt2Ny#>uABbNeOwPW70!I!_-h^? z50@3r`<a{hS3(oua$u)H&Fl1A8{mvreorb0uV35<XMD5V{MdW93GNZv-({Ze?+VH{ z!+C%0PycIt-vYOiex{>-7x6ntBAol>eQ+OrZ*7J9kbcaUQZuhe%D2Jkho(`kKUxRB z9c~%*#_Nx<tD-yL^y8i<{kikq3HJu=!uag(*aheL@Vw}sJukcAT!#+y>pI!r9=Npp zuI%Uy_(j-z;fi3pU!EKHZ6BO|b|Jj})_9l%XFg*53ggy(xH|OH1@(98PY=Kuj~vhO zjSC0i?(_THzi{6U!FiwiGVUf+e{~qnc;NXD<HD10e=_bLsN;^L-y?9wY4b4SwQ>4r zFhA$DKA7KAaGMw}KlP2V8(|*{#{G)=Td3n7hx0snz6#Kf=jv&=&9o2ma`WsHa6M?> z6YY)~k50lR(%%l$c;h*JCgA$Q`CaEA_9;01aBar(_joU#g$w<z{?`2HG+ZzGF)kP< zJV$5X+$Z<V^X0iZ3->(pHNScRdn@CdgG-{^eQu2H{yqoiJ1`gZ=Eu?4&%^b>9*GV> z$20DExFOi9Q1hTY*cagRyU8feh>k;Ffb(9PXV#(Kd;TKa>&)i@s{fb-{}P;eweuec zKO634xSZII7y6l3;9Sqr@a0hR(~EH4o7;@%efIvo3U`V6^Qimbd*KpXAI8~&dM^^O zFT?f2HqINL?e8@>^OYI!p4W}oufv%aHm1BD>iBQKxvxL+JHq{a6Rsxpz89}!n=ifv zXPh!l8aKU{Z^Nah{Q~AS4fQ-;f%ANLZt}vlrTta7(bz*#{fGYL9k?m9_Z;=cjz`~x z8;<RLFb~>K`Fn7lTla4*_3q#MaQadGs^2f_Gv8})5wv%nzH6M{2XHHBw-}8_&96R$ z(+@eH^KhQ4kKo!<?)L`sm;%h_I-K*LOS!+xxbQKY_vbhIy@UEYK7sRl&LQfzqn?W! zaIUZUi}~kH=5Z6QKkcK@5vbg!a8Z<Zr@e7(BKBu+4X~@A{$79o=WyPiuj%)DbPesk zfU}>fl;=ZTzc1m=(S9jv9y$m9E4Uier=r}r8iV~cocB$CF^~G~=q)(c%k}g7xa<9m z9PP}jjf2L;+i>P{Tj9Mg=5^n~`8~z^>)&s^2Y2A~2N(F=uc7vP7w%ii|3=Lx_2b{c zxh}n!pX;MP{T?oodh>xY*yaO2z?EY>_sxBLmGU3qy3yWoqp%(Kr(j&`ZRb4h!A+$9 zMD!@C|N9xPEaUilT4LMZFK|sL_ndleeV_jd*Ay;4S`9V7_zms>^YMLg4cqhaJDlg( zd57ou54fjj-<I*hd-*4v>;4M-8PxCX_u<@k<LhAhjYA*6={Lqwz6|xA{{^>=akito zP~+U+aK<0wk>`Cg{r&^j3ET7L`SZU03-=x4{Dppl`a6=P{(nJQ`g5IKcm2jAaK7^% z_P%`-?gZoL7mXA8wd8Q-55K}+M_rc`aM_t(gzEt3`lN*O9DA-^7teDlxV^OBg4RKe z=c(a*7x+$yfHTib12>Fv&*4FAe}7sy{j+(Jzt?q62RDuOekWOly^i0X9&S73_O}Du z{xZOMem$R_>qPp?2zQ$Heiyocy$voCoaf#1Z@>1J8P0dkGxYN!nm~V9;4)C&1D#L3 zapEyJ_tX6?1?RbsfV)pWsj2t9?mf>6XFl{Cyz_N_+2D*Hx8RH;#*gf9em8iOaqq&L zFXn)A{arWr$9tI*t||R^-}S5d(_C=It*(r>fOh83x#9dB#y|b<QusV@J?Y2uAD;WX zaQY3;@lUk#{N;nw|LBigC-aW{aACQ;c}W4d_RO~-S`alaDG1k!^46&N-4ob_;Jmlq z<0<?O?{8taSjx@2jO$7CR|GCO_9yV4qRS~S3a4Lm{4j16gX>2-@1gmPc~fz?<yZ z=R7=rCE$#+_Ur!JUrD$-^k>{KkDtMOO2KV}t4?`C>?3GtIQ@n3(euBHc4gqaw;$3z zyuW4PT>n_O*{J(dE*N(r<Cx$3d&|S=hql7GPI49CjDz|!<Dq`5BAo9h-%~Z|w*XoR z&hlE6r{j0%4=cm%rhFf29P?dT1<rN$p1a<zb5%IM4}M91U!(3zHMqW%C!*$;<_p!~ zno;lfkek?h`CT>Su=84udai52)u!Baaou)X4)+=3c-}w7_P*7IdzSKxXb05(>cEYq zd<43L@phth;o4B|@6aE;PI)~z_p=7=-Ea4|KHQVEbD#HM>&F|wWuhPR`U7z8M?<*L zl)GQHlWPPwmHuMU-sonuF<fr!qNsV1?V7+%VBE*y`eVC(O@nbA&we)0t{I%~&S?4x z|No^qTr=9a&aS)b+yZVCoZm<KU~fTN!d+y%cToL`=Oq#@7xn&5zyCSkR&d?u-}qx5 z;5)B1oc`K-J%M)Sr)}VhP+l0Vh<aYz!s#c>b29|}wS)70mzLk7KR2&v4`&{0+)qV+ z#{CX(v*>3Lx)AkVc7*d@dLP4k-U-h4r1Nr~&aX3^{@r-5KQuq>0;m7e@7-d4`oFGl z-gED_d9wGs8(a(8Cu2R#dlRs`!{xx%Kc>X?zC^*9?=_-d<LM^s9&in4-xeK3z4=Q| zILqb2Jf|0&=f(54kN!M=z2RJE?^gym<7FQ>&z<kB6wK59`oc}7z4vGgw(HUlt}pF< z4|$(`2lR(?eKW&{>pTF?d*eMS2=6@_2-lYR=<iQq>u(0Zx$nL!JQtpW!El`zX9(I9 zU2i$uHO8?Y`?bHJa2w$cp?*(YK!3yFVyKTrSD<6yhQmc;n;$O3c7I2}8TU`pkMZBQ zKN8OPzLtLa1^tbJ)1Mfx!uU2CZYkrHL+7G?=ZS{PNV)OScxXN|CKxv@^SVL1Mc8BE zmQlVMbzI*8<KWC&TF^cny!T~1oad+=^`0m5#|dx`sCT_R->&mSxVPahqxvh~N0Z>n zQtx_Q!Cs3_hO5u~JD`P7<KYyzR+Q^^!+19pPXFV(C-g&)!#zrW#=|fUPJ{cK`TCyy z7TfbW9j-L}`Q0EFwtiy<oN>i?*c{F{GZW5za9=$C?#C=R<5olFW&F~g&xSK!^W2B$ zZw}ll`Z2Dq$JQUth4a1RJ@ejqf1iNUujprtU;3AMaE>?M^+lJnKJ(#x@B5y=3pWpL z0o>!14@L)~=A8@SD$}kC8vdPR5nN8%8%N&4?vKX68He^#z7O3`|FLl1pLeMb@6lqo zZq)lOHm~!0&JsB9YZ}J!+<0D=!kwo-&(rhRa?9ZSJ~*BFc+|XhIh?=WygZDnE8x60 zx9R5=x{LX)gu6)jIrJ*3KUf8)KPX1~)~NaKYB=-6Pw1xtngxChTpR3WXgWCkX&jvX z;797y*$=uFZa?*F(8B03bRAq@Y~x56Ki0$bgY$io1KWHe9`0f1#C1-9GjB@4yfdJl z(+zMnDK~ByM~ojE;amsrxA)z9z6mahd1pZNqxy}_aI>k8LzkkSgDr6Sv+C3zK#j+V zaK;Pchxf?yvK4Lz?Gu9jVpxZ5aK>rlSvB~_;I_l1r(FN!ecl1L1I}~ldGXwOPItn& z-{CoEPk+1M@-VLJqyP7O?S?xCcM*LDb)EOXeNOow=$GhX#@!3&d)sy0K)vtKeQ+%) zuZNaJeJ3Zu=}*kB3o)*Kct6}&%15GpU-7%Z0XWZ{=g@bC=kFk#{-rqc($Dz24#DaF zj0ZUx|15eKPJhyq{`5EcmnY%OUyMuj=*N8K2;4Brr=p`!-(^SPN-$nWRR3nYehN;1 zrhogI{?^d$7@X%*|5g}oIeHw<yd*8<zCXHRKP`uSAN5{&k50fDAI(2K@2zQf67Dr@ z|K9Qv_Evt!GjQfP<~P&e%yUk`9i)6SdIEKwpM^8eEo(ihe>e^2d9T2H>Y|?eGjIhc zH$Se8y&F9X_n%*5JpGIJ?HrtbpgLSC`q5uL2N!-9a37D*?s+)BW0;5NPmMF@;k<_v zXzzXTUS5D3PQASGZVB^z0q!zfZd5<$xp)z-6Ya{NrBU<om*Bj&zT=zH?<MqQxGUIa z&{k+_`h5k?^XU2X+%|-}2<P`TzYlm%y_c`Tb!Hs>WhDK$@0Z}bXWp+Uxb^5|IQ^)< zw-Ej7PhW%6--Le1d+<7(`J?yE`{#Xo1J3tWa^~Z^%Y5TaxSh<$d**m6Y4;YK{?B=4 zrk(fo?O=Wvs6UTxqrWR~uE$O4!*g&I?q}-1MSnq;(%(C9Vf^*IH5&e1IOC9cknu^s z`yQP6jQ7|$zKC((hpU6#1nq(vhp)kf=gRYBe)IucCE6E7?LQX#L%1#cu48Cj)coZm zxD>RL54r1bo?p+q=iBrAF`W0%c;P+uUVdVK%*WrSze;32H{cpzd+v?L`um%3Q|K=| zf1b-v;ml(@&~F*q#i5_U6{0@mZ2vjjMcVmY+kDNq@CBUr#{1wr{T*Mzc@E5XS~CyV z=_@$%zuy_(Jkb2_Yq(L=yFRXm<KBYH!8qQhG}zvwZ{W=9j0?l*-?(ra?jGfL(BIH_ z`u!Fzn)1tN7)S2F`ThPIxW7>2-CejA^lx654$i#pJGi2h&q24M$I<WM_F(Tq5233V z=Lb0V{|)NHbMPaa^D^Ibf4k7`C%7TldC`36Kj=NU-a&nS?Dx>0;oirtgjPmppufQF z#`gF6`~AJY!ntm)<2u^8e!sz8rQG{<4!b=1JDm5-{Hi1M7tuf9KE<wx`g;~){|R>; zJ1zTUeieg#AI^K`y$k*B1Gob8=lcGOJpla+&ihz}@+VN&`)@eM=|+7o=4GDv4_q<o zm!R$7%Ax<leMkLI=&$H@xMXSmU(lEGOlU?l5q$*CbMPeP`_Mh~_b8lk(Dn!5jDyMH zCQ@#^bslml;EXfIm*$LT9+eXAF#VaI=wHoCQo-GYdsx1c8qW3dJErUD`lW%3U_STg zH}q?1;XI$S;66bspy}XhVb?}KKsTW2;e40)F8AHyJ0%0$3fje^)6vz;J0o29{&T&J z=b7Mg&^|lb4=qT4nc)UgJ{}!}`h7PG+yu%;pd-<p^!FHC9OdWG^5{}D0<HyiOY|8u z8U1C2^Bj5KJHdJ0vcZ{OWTD>tV-$9FxU1Out&Z3i(HwC4#fRPFoNzU1r{Ai9oq*<o zi^X;w=3*!Ddvn9}r#z?a(X(hCIPZmdpzneG^p_W|3GMVhN!X3x^1+p$yg2ILU8Z5@ zhjX7^p!_*Bf&L1>{Y<%eeVCUNgsTfz2i=1Dy9&W|r#uSv-Mou&3&VA$K9X|Ji+;NZ zob}$P@E#R~D?@*!(L(e)8@m`>Gwi0Q@1IN9#o>z6-gkxX0N)iQ;2zNKB-$KxoRV<8 zD9?swMRPMwDY#OUmqhbWub(Mxf0U<(GtT*)wG5p1&-HaZUFWiJ1?k5){1J8^v>crG z)_XjKdhc&}xHFWuLVce^Vpo7$j$Ifng6>Bv!aa?>4>fM?VI3;Lng5yZb*A3DwKAOh z?Y^C--@Wiv;M!Ae{@DTB?^RXd9M5-z^KkrXaJ~<Gzh|O9-!0YQ22!ry5B+@&xSntk z=woP3en(BXB+8GY`WgLBEx4Ddk3>I2TcWk$O3;5r)bFi+)<jdtD({g82`E?g)2 z@t&8#HeS|)+f0AHSL<QdSK&POVLbBOH-HQGD||0BgiB8UlhC(Ozk4-;OHH|P%QzNC zd1JVv*nZ!56ni4t1de4$=6!pLdhcCRxRjI|55s+K2A2x%7`g`ayG3)j)0DSH8#0kL z*e&2DV`rv51@>ELOSrYzS5S(R`5hyYf@Rd-N2{S#(N=KA<#LplMPEf*!!^S8J8xs` zO=ugq8`yc#Jm_e&E!=Qy^R;Z)h57yM;7VhAj!R*0LfgZ+|EVZ<9*eL$z!^`Cuf|os z=X8YY%6xq9cn`cso#2+hnKu>3_TAGNZVu)8<)YYI(JpW^u|G!LuPWGG;W}gceb@V^ zpX>%_oY9|!d180C^YnKCeHY!v`bELTQtr7~jhzJ715SUKM15n_Jgg^NY068X<>@B| zyBFL>?CR7vz;1~4hBK~oq}+G<MeIIso+tgJ@!b9G3wN1zxzSu`G};fY5Vr9@CG#9Z zd4ISJl-HqrFIpQN0Qa!>g!gnH+$q|%L))S?7-tY%b;`XT2jtMfaNDqJqUq_U7WR;! zKkuD>SARSdPCsj$rn_YN-C=O+=wJWqeKNidhl_)|iW(<;ca4CXM7{AQ2exr$B%Jft zPx(70!jFP84tcH|*LjSF^PD)Y^LU&7qT#&%{$79oakw#X`t{J?w}KlB=X>3_w1#@) z(Kxsgl$S(TYO^up;VNNAP@feYhfaVqo|-?nzwY})INxucf8RU4-zLH3ql3eY<2^So znG9!qECgrVoQypMt{V1Ql;Z#5<x~m^Q}6m6!*(4Xhx7Mkq2Ax`d7cL6e44^JU*q|7 zIM28H?s?FE&VcjWdX79_<`pyHUS~Y>i@Z!&zc35Ve6BNmC)B)dHk|9=y3~d9T+M;= zet0j+z{R0+;f`W^&oW|9M4y0jKKkXf)EkHA!MWeAZ%#PZc|IJH%>3pGoabTzfbW!= ztebxIP56ayU9erhF#axrYXoPU3jZ!1182S*_UAaU!MMJQ($Rkl=DQfK4R*-c-x4_2 zwK=@&?m91pi=jR{It)#VE(<v8<)^?chr7qTL(cwI1f2cI+uusKne?9#9gB`bR|TB) z@_#a))o}5&`+@S?=x^wnfU{nH5bfjO7GOL70obF_wE<_n{4}_AaEr0OU_L{!2czo) z&U*RJX%`Q-BIs{k&|gBpS#N)L>2CvESK7})E1-?hjc}IB*=`KnCb&h|Cm6@S3npT3 z4mj)OUC%9We&;Jmzs5!5VIrL6a<=n(;Z`{P>?Hd0{XQLgTfq4al<Nq$9nNw&`?tRx za64(|_ob1n#{_g|z*#Rpp7y)oVzG0g<LGB9_U?eQUf#T658OuD8`q8J)#-09oaJ)1 zGY{AYryn-Ioligd>!g6QUfz6pKin$XJKiDYy9YfGaMsIvuMfg~M7uGl?*Z5Q5S-<5 zwoA!;4#Vxh?n-|{u-Btc2AuWs>1lTa&U`+Ba{bD6?4xj&%h}HN!BcQYXy<;q9<I-^ zfU{oScz7Jndp4bZ!#wqAIqLu8SP#GRo`8E5&Ny$JF#ex}vs})0b7}t!oZkt0FrINF z3hq?ESud}jdloJO?LyA}P6wR*$lKo;xMlR`yCM8tXW=ZDv)z5#pM%p6<>mJ_Lnoon b!C5Y6JN;U^V$H`6X&yav%+RruC#Cxzl7>kA literal 0 HcmV?d00001 diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/meshes.py b/tests/fullscale/linearelasticity/faults-3d-buried/meshes.py new file mode 100644 index 0000000000..46cb7ee9b8 --- /dev/null +++ b/tests/fullscale/linearelasticity/faults-3d-buried/meshes.py @@ -0,0 +1,42 @@ +# +# ---------------------------------------------------------------------- +# +# Brad T. Aagaard, U.S. Geological Survey +# Charles A. Williams, GNS Science +# Matthew G. Knepley, University at Buffalo +# +# This code was developed as part of the Computational Infrastructure +# for Geodynamics (http://geodynamics.org). +# +# Copyright (c) 2010-2022 University of California, Davis +# +# See LICENSE.md for license information. +# +# ---------------------------------------------------------------------- + +from pylith.testing.FullTestApp import MeshEntity + + +class TetGmsh(object): + """Mesh information for tet mesh using Gmsh. + """ + ENTITIES = { + "domain": MeshEntity(ncells=1450, ncorners=4, nvertices=425+9), + + # Materials + "mat_elastic": MeshEntity(ncells=1450, ncorners=4, nvertices=425+9), + + # Faults + "fault": MeshEntity(ncells=22, ncorners=3, nvertices=18), + + # Boundaries + "bc_xneg": MeshEntity(ncells=84, ncorners=3, nvertices=55), + "bc_xpos": MeshEntity(ncells=84, ncorners=3, nvertices=55), + "bc_yneg": MeshEntity(ncells=84, ncorners=3, nvertices=55), + "bc_ypos": MeshEntity(ncells=84, ncorners=3, nvertices=55), + "bc_zneg": MeshEntity(ncells=162, ncorners=3, nvertices=98), + "bc_zpos": MeshEntity(ncells=164, ncorners=3, nvertices=94), + } + + +# End of file diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/pylithapp.cfg b/tests/fullscale/linearelasticity/faults-3d-buried/pylithapp.cfg new file mode 100644 index 0000000000..40cd56965f --- /dev/null +++ b/tests/fullscale/linearelasticity/faults-3d-buried/pylithapp.cfg @@ -0,0 +1,121 @@ +[pylithapp.metadata] +keywords = [full-scale test, 3D, box, fault, buried edges, field split conditioner, schur complement] +features = [ + Field split preconditioner, + Schur preconditioner, + pylith.meshio.MeshIOPetsc, + pylith.problems.SolnDispLagrange, + pylith.problems.TimeDependent, + pylith.meshio.DataWriterHDF5, + pylith.meshio.OutputSolnBoundary, + spatialdata.spatialdb.UniformDB + ] + +[pylithapp] +# dump_parameters.filename = output/SIMULATION-parameters.json +# problem.progress_monitor.filename = output/SIMULATION-progress.txt + +[pylithapp.launcher] # WARNING: THIS IS NOT PORTABLE +command = mpiexec -np ${nodes} + +# ---------------------------------------------------------------------- +# journal +# ---------------------------------------------------------------------- +[pylithapp.journal.info] +timedependent = 1 +solution = 1 +petsc = 1 +meshiopetsc = 1 +isotropiclinearelasticity = 1 +dirichlettimedependent = 1 +faultcohesivekin = 1 + +[pylithapp.journal.debug] +#timedependent = 1 +#solution = 1 +#isotropiclinearelasticity = 1 +#dirichlettimedependent = 1 +#constraintspatialdb = 1 +#faultcohesivekin = 1 +#integratorinterface = 1 +#kinsrcstep = 1 +#outputphysics = 1 +#outputsolndomain = 1 + +# ---------------------------------------------------------------------- +# mesh_generator +# ---------------------------------------------------------------------- +[pylithapp.mesh_generator] +reader = pylith.meshio.MeshIOPetsc + +[pylithapp.mesh_generator.reader] +filename = mesh_tet.msh +coordsys.space_dim = 3 + +# ---------------------------------------------------------------------- +# problem +# ---------------------------------------------------------------------- +[pylithapp.problem] +defaults.quadrature_order = 1 + +# Use nonlinear solver to ensure residual and Jacobian are consistent. +solver = nonlinear + +solution = pylith.problems.SolnDispLagrange + +[pylithapp.problem.solution.subfields] +displacement.basis_order = 1 +lagrange_fault.basis_order = 1 + +[pylithapp.problem] +solution_observers = [domain, bc_ypos] +solution_observers.bc_ypos = pylith.meshio.OutputSolnBoundary + +[pylithapp.problem.solution_observers.bc_ypos] +label = boundary_ypos +label_value = 13 + +# ---------------------------------------------------------------------- +# materials +# ---------------------------------------------------------------------- +[pylithapp.problem] +materials = [mat_elastic] + +[pylithapp.problem.materials.mat_elastic] +description = Elastic material +label_value = 1 + +db_auxiliary_field = spatialdata.spatialdb.UniformDB +db_auxiliary_field.description = Elastic properties +db_auxiliary_field.values = [density, vs, vp] +db_auxiliary_field.data = [2500*kg/m**3, 3.0*km/s, 5.2915026*km/s] + + +auxiliary_subfields.density.basis_order = 0 +bulk_rheology.auxiliary_subfields.bulk_modulus.basis_order = 0 +bulk_rheology.auxiliary_subfields.shear_modulus.basis_order = 0 + +derived_subfields.cauchy_stress.basis_order = 0 +derived_subfields.cauchy_strain.basis_order = 0 + + +# ---------------------------------------------------------------------- +# PETSc +# ---------------------------------------------------------------------- +[pylithapp.problem.petsc_defaults] +solver = True +testing = True +monitors = True + +[pylithapp.petsc] +ksp_rtol = 1.0e-15 +ksp_atol = 1.0e-15 +ksp_max_it = 100 +ksp_gmres_restart = 50 + +snes_atol = 1.0e-12 + +snes_max_it = 1 + + +# End of file diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/test_pylith.py b/tests/fullscale/linearelasticity/faults-3d-buried/test_pylith.py new file mode 100755 index 0000000000..73621bbc4a --- /dev/null +++ b/tests/fullscale/linearelasticity/faults-3d-buried/test_pylith.py @@ -0,0 +1,51 @@ +#!/usr/bin/env nemesis +# +# ====================================================================== +# +# Brad T. Aagaard, U.S. Geological Survey +# Charles A. Williams, GNS Science +# Matthew G. Knepley, University at Buffalo +# +# This code was developed as part of the Computational Infrastructure +# for Geodynamics (http://geodynamics.org). +# +# Copyright (c) 2010-2022 University of California, Davis +# +# See LICENSE.md for license information. +# +# ====================================================================== + +from pylith.testing.FullTestApp import TestDriver, FullTestCase + +import unittest + + +class TestApp(TestDriver): + """Driver application for full-scale tests. + """ + + def __init__(self): + """Constructor. + """ + TestDriver.__init__(self) + return + + def _suite(self): + """Create test suite. + """ + suite = unittest.TestSuite() + + import TestUniformSlip + for test in TestUniformSlip.test_cases(): + suite.addTest(unittest.makeSuite(test)) + + return suite + + +# ---------------------------------------------------------------------- +if __name__ == '__main__': + FullTestCase.parse_args() + TestApp().main() + + +# End of file diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/uniformslip.cfg b/tests/fullscale/linearelasticity/faults-3d-buried/uniformslip.cfg new file mode 100644 index 0000000000..64c9ffca5f --- /dev/null +++ b/tests/fullscale/linearelasticity/faults-3d-buried/uniformslip.cfg @@ -0,0 +1,108 @@ +[pylithapp.metadata] +# There is no analytical solution for this problem. We just check ranges of values in the output. +# +description = Prescribed uniform left-lateral fault slip with Dirichlet boundary conditions. +authors = [Brad Aagaard] +keywords = [fault, prescribed slip, uniform slip, tetrahedral cells] +version = 1.0.0 +pylith_version = [>=3.0, <4.0] +base = [pylithapp.cfg] +arguments = [uniformslip.cfg] + +features = [ + Static simulation, + Field split preconditioner, + pylith.faults.FaultCohesiveKin, + pylith.faults.KinSrcStep, + pylith.materials.Elasticity, + pylith.materials.IsotropicLinearElasticity, + pylith.bc.DirichletTimeDependent, + spatialdata.spatialdb.UniformDB + ] + +# ---------------------------------------------------------------------- +# output +# ---------------------------------------------------------------------- +[pylithapp] +dump_parameters.filename = output/uniformslip-parameters.json +problem.progress_monitor.filename = output/uniformslip-progress.txt + +problem.defaults.name = uniformslip + +# ---------------------------------------------------------------------- +# solution +# ---------------------------------------------------------------------- +[pylithapp.problem.solution.subfields.displacement] +basis_order = 1 + +# ---------------------------------------------------------------------- +# fault +# ---------------------------------------------------------------------- +[pylithapp.problem] +interfaces = [fault] + +[pylithapp.problem.interfaces.fault] +label = fault +label_value = 20 + +edge = fault_edges +edge_value = 21 + +observers.observer.data_fields = [slip, lagrange_multiplier_fault] + +[pylithapp.problem.interfaces.fault.eq_ruptures.rupture] +db_auxiliary_field = spatialdata.spatialdb.UniformDB +db_auxiliary_field.description = Fault rupture auxiliary field spatial database +db_auxiliary_field.values = [initiation_time, final_slip_left_lateral, final_slip_reverse, final_slip_opening] +db_auxiliary_field.data = [0.0*s, +2.0*m, 0.0*m, 0.0*m] + + +# ---------------------------------------------------------------------- +# boundary conditions +# ---------------------------------------------------------------------- +[pylithapp.problem] +bc = [bc_xneg, bc_xpos, bc_yneg, bc_ypos, bc_zneg] +bc.bc_xneg = pylith.bc.DirichletTimeDependent +bc.bc_xpos = pylith.bc.DirichletTimeDependent +bc.bc_yneg = pylith.bc.DirichletTimeDependent +bc.bc_ypos = pylith.bc.DirichletTimeDependent +bc.bc_zneg = pylith.bc.DirichletTimeDependent + +[pylithapp.problem.bc.bc_xneg] +label = boundary_xneg +label_value = 10 +constrained_dof = [0] +db_auxiliary_field = pylith.bc.ZeroDB +db_auxiliary_field.description = Dirichlet BC -x edge + +[pylithapp.problem.bc.bc_xpos] +label = boundary_xpos +label_value = 11 +constrained_dof = [0] +db_auxiliary_field = pylith.bc.ZeroDB +db_auxiliary_field.description = Dirichlet BC +x edge + +[pylithapp.problem.bc.bc_yneg] +label = boundary_yneg +label_value = 12 +constrained_dof = [1] +db_auxiliary_field = pylith.bc.ZeroDB +db_auxiliary_field.description = Dirichlet BC -y edge + +[pylithapp.problem.bc.bc_ypos] +label = boundary_ypos +label_value = 13 +constrained_dof = [1] +db_auxiliary_field = pylith.bc.ZeroDB +db_auxiliary_field.description = Dirichlet BC +y edge + +[pylithapp.problem.bc.bc_zneg] +label = boundary_zneg +label_value = 14 +constrained_dof = [2] +db_auxiliary_field = pylith.bc.ZeroDB +db_auxiliary_field.description = Dirichlet BC -z edge + + + +# End of file diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/uniformslip_soln.py b/tests/fullscale/linearelasticity/faults-3d-buried/uniformslip_soln.py new file mode 100644 index 0000000000..1a9d30c93e --- /dev/null +++ b/tests/fullscale/linearelasticity/faults-3d-buried/uniformslip_soln.py @@ -0,0 +1,96 @@ +# ---------------------------------------------------------------------- +# +# Brad T. Aagaard, U.S. Geological Survey +# Charles A. Williams, GNS Science +# Matthew G. Knepley, University at Buffalo +# +# This code was developed as part of the Computational Infrastructure +# for Geodynamics (http://geodynamics.org). +# +# Copyright (c) 2010-2022 University of California, Davis +# +# See LICENSE.md for license information. +# +# ---------------------------------------------------------------------- +# +# @brief Uniform slip on a buried fault (no analytical solution). +# +# We just check info fields. + +import numpy + +# Physical properties +p_density = 2500.0 +p_vs = 3000.0 +p_vp = 5291.502622129181 + +p_mu = p_density * p_vs**2 +p_lambda = p_density * p_vp**2 - 2 * p_mu + +# ---------------------------------------------------------------------- +class AnalyticalSoln(object): + """Analytical solution to axial/shear displacement problem. + """ + SPACE_DIM = 3 + TENSOR_SIZE = 6 + + def __init__(self): + self.fields = { + "density": self.density, + "shear_modulus": self.shear_modulus, + "bulk_modulus": self.bulk_modulus, + "initial_amplitude": self.zero_displacement, + "slip": self.slip, + } + return + + def getField(self, name, mesh_entity, pts): + return self.fields[name](pts) + + def getMask(self, name, mesh_entity, pts): + mask = None + if name == "displacement": + fdist = numpy.dot(pts-fault_center, fault_normal) + index = numpy.argmin(numpy.abs(numpy.min(fdist))) + mask = numpy.abs(fdist) < 1.0 + return mask + + def zero_displacement(self, locs): + """Compute displacement field at locations. + """ + (npts, dim) = locs.shape + + disp = numpy.zeros((1, npts, self.SPACE_DIM), dtype=numpy.float64) + return disp + + def density(self, locs): + """Compute density field at locations. + """ + (npts, dim) = locs.shape + density = p_density * numpy.ones((1, npts, 1), dtype=numpy.float64) + return density + + def shear_modulus(self, locs): + """Compute shear modulus field at locations. + """ + (npts, dim) = locs.shape + shear_modulus = p_mu * numpy.ones((1, npts, 1), dtype=numpy.float64) + return shear_modulus + + def bulk_modulus(self, locs): + """Compute bulk modulus field at locations. + """ + (npts, dim) = locs.shape + bulk_modulus = (p_lambda + 2.0 / 3.0 * p_mu) * numpy.ones((1, npts, 1), dtype=numpy.float64) + return bulk_modulus + + def slip(self, locs): + """Compute slip field at locations. + """ + (npts, dim) = locs.shape + slip = numpy.zeros((1, npts, 3), dtype=numpy.float64) + slip[:,:,1] = +2.0 + return slip + + +# End of file