From e5e98c717ba1d8b499243fc5f39e07664fc4d2e0 Mon Sep 17 00:00:00 2001 From: Tim Brooks <41971846+timryanb@users.noreply.github.com> Date: Sat, 3 Jun 2023 06:12:38 -0400 Subject: [PATCH] Adding buckling eigenvalue solver (#213) * Adding geometric stiffness to shell elements * Adding TACSBuckling API to cython wrapper * Adding BucklingProblem class * Adding buckling integration test * clang format fixes * loosening up shell test tol * Adding buckling solver to mphys wrapper * Minor mphys refactor * Clang format fixes * Reducing mesh size for mphys buckling test * Increasing testflo CI timelimit * Updating sigma value in eigenvalue test * Adding rtol setting to ModalProblem setup * Adding buckling example and benchmark * Reducing N_PROCS to 1 for mphys buckling test --- .github/workflows/unit_tests.yml | 2 +- docs/source/pytacs/buckling.rst | 19 + docs/source/pytacs/problems.rst | 3 +- examples/plate/benchmark/benchmark_buckle.py | 44 + examples/plate/plate_buckle.py | 50 + examples/plate/ss_plate_buckle.bdf | 2736 +++++++++++++++++ src/TACSBuckling.cpp | 221 +- src/TACSBuckling.h | 7 +- src/elements/shell/TACSShellElement.h | 69 +- tacs/TACS.pxd | 9 +- tacs/TACS.pyx | 127 +- tacs/mphys/mphys_tacs.py | 152 +- tacs/problems/__init__.py | 9 +- tacs/problems/base.py | 5 + tacs/problems/buckling.py | 972 ++++++ tacs/problems/modal.py | 10 +- tacs/pytacs.py | 40 + .../shell_tests/test_shell_element.py | 2 +- .../input_files/ss_plate.bdf | 200 ++ .../test_mphys_struct_buckling.py | 149 + .../test_shell_plate_buckling.py | 75 + .../test_shell_plate_quad.py | 2 +- 22 files changed, 4863 insertions(+), 40 deletions(-) create mode 100644 docs/source/pytacs/buckling.rst create mode 100644 examples/plate/benchmark/benchmark_buckle.py create mode 100755 examples/plate/plate_buckle.py create mode 100644 examples/plate/ss_plate_buckle.bdf create mode 100644 tacs/problems/buckling.py create mode 100644 tests/integration_tests/input_files/ss_plate.bdf create mode 100644 tests/integration_tests/test_mphys_struct_buckling.py create mode 100644 tests/integration_tests/test_shell_plate_buckling.py diff --git a/.github/workflows/unit_tests.yml b/.github/workflows/unit_tests.yml index 560518f8d..67961cb1a 100644 --- a/.github/workflows/unit_tests.yml +++ b/.github/workflows/unit_tests.yml @@ -133,7 +133,7 @@ jobs: if [[ ${{ matrix.NAME }} == 'Real' ]]; then source ${TACS_DIR}/extern/ESP122/ESPenv.sh fi - testflo --timeout 120 .; + testflo --timeout 240 .; - name: Build docs run: | cd $TACS_DIR/docs; diff --git a/docs/source/pytacs/buckling.rst b/docs/source/pytacs/buckling.rst new file mode 100644 index 000000000..a99fe0a83 --- /dev/null +++ b/docs/source/pytacs/buckling.rst @@ -0,0 +1,19 @@ +BucklingProblem +------------ +.. automodule:: tacs.problems.buckling + +Options +^^^^^^^ +Options can be set for :class:`~tacs.problems.BucklingProblem` at time of creation for the class in the +:meth:`pyTACS.createBucklingProblem ` method or using the +:meth:`BucklingProblem.setOption ` method. Current option values for a class +instance can be printed out using the :meth:`BucklingProblem.printOption ` method. +The following options, their default values and descriptions are listed below: + +.. program-output:: python -c "from tacs.problems import BucklingProblem; BucklingProblem.printDefaultOptions()" + +API Reference +^^^^^^^^^^^^^ +.. autoclass:: tacs.problems.BucklingProblem + :members: + :inherited-members: \ No newline at end of file diff --git a/docs/source/pytacs/problems.rst b/docs/source/pytacs/problems.rst index 2d644702b..645d1ec32 100644 --- a/docs/source/pytacs/problems.rst +++ b/docs/source/pytacs/problems.rst @@ -6,4 +6,5 @@ Problem classes static transient - modal \ No newline at end of file + modal + buckling \ No newline at end of file diff --git a/examples/plate/benchmark/benchmark_buckle.py b/examples/plate/benchmark/benchmark_buckle.py new file mode 100644 index 000000000..84c52d720 --- /dev/null +++ b/examples/plate/benchmark/benchmark_buckle.py @@ -0,0 +1,44 @@ +""" +This script is used to regression test the example against historical values. +""" + +import numpy as np +import unittest +import sys +import os + +# Set the path to the example script we're testing +example_path = os.path.join(os.path.dirname(__file__), "..") +sys.path.append(example_path) + +# Reference values for eval functions +FUNC_REF = { + "buckle_eigsb.0": 20.082400947865022, + "buckle_eigsb.1": 43.01440496545738, + "buckle_eigsb.2": 80.31233170890538, + "buckle_eigsb.3": 84.69751385025641, + "buckle_eigsb.4": 86.82024838978569, +} + + +class ExampleBenchmark(unittest.TestCase): + N_PROCS = 8 # this is how many MPI processes to use for this TestCase. + + def setUp(self): + # Import the example to automatically run the script + import plate_buckle + + self.example = plate_buckle + + def benchmark_funcs(self): + """ + Test the example eval functions against reference values + """ + func_dict = self.example.funcs + + # Test functions values against historical values + for func_name in func_dict: + with self.subTest(function=func_name): + np.testing.assert_allclose( + func_dict[func_name], FUNC_REF[func_name], rtol=1e-6, atol=1e-6 + ) diff --git a/examples/plate/plate_buckle.py b/examples/plate/plate_buckle.py new file mode 100755 index 000000000..a326b8ab9 --- /dev/null +++ b/examples/plate/plate_buckle.py @@ -0,0 +1,50 @@ +""" +This problem performs a buckling analysis on a 4m x 3m flat plate. +The perimeter of the plate is pinned and loaded in compression (20kN/m) on its horizontal edges. +We use TACS buckling eigenvalue solver through the pyTACS BucklingProblem interface. +""" +# ============================================================================== +# Standard Python modules +# ============================================================================== +from __future__ import print_function +import os + +# ============================================================================== +# External Python modules +# ============================================================================== +from pprint import pprint + +import numpy as np +from mpi4py import MPI + +# ============================================================================== +# Extension modules +# ============================================================================== +from tacs import pyTACS + +comm = MPI.COMM_WORLD + +# Instantiate FEAAssembler +bdfFile = os.path.join(os.path.dirname(__file__), "ss_plate_buckle.bdf") +FEAAssembler = pyTACS(bdfFile, comm=comm) +# Set up constitutive objects and elements +FEAAssembler.initialize() + +# Setup static problem +bucklingProb = FEAAssembler.createBucklingProblem(name="buckle", sigma=10.0, numEigs=5) +bucklingProb.setOption("printLevel", 2) +# Add Loads +bucklingProb.addLoadFromBDF(loadID=1) + +# solve and evaluate functions/sensitivities +funcs = {} +funcsSens = {} +bucklingProb.solve() +bucklingProb.evalFunctions(funcs) +bucklingProb.evalFunctionsSens(funcsSens) +bucklingProb.writeSolution(outputDir=os.path.dirname(__file__)) + + +if comm.rank == 0: + pprint(funcs) + pprint(funcsSens) diff --git a/examples/plate/ss_plate_buckle.bdf b/examples/plate/ss_plate_buckle.bdf new file mode 100644 index 000000000..40288c746 --- /dev/null +++ b/examples/plate/ss_plate_buckle.bdf @@ -0,0 +1,2736 @@ +INIT MASTER(S) +NASTRAN SYSTEM(442)=-1,SYSTEM(319)=1 +ID FEMAP,FEMAP +SOL SEBUCKL +CEND + TITLE = Simcenter Nastran Buckling Analysis Set + ECHO = NONE + DISPLACEMENT(PLOT) = ALL + SPCFORCE(PLOT) = ALL + STRESS(PLOT,CORNER) = ALL + SPC = 1 +SUBCASE 1 + LOAD = 1 +SUBCASE 2 + LABEL = BUCKLING CASE + METHOD = 1 +BEGIN BULK +$ *************************************************************************** +$ Written by : Femap +$ Version : 2023.1.0 +$ Translator : Simcenter Nastran +$ From Model : +$ Date : Fri May 19 13:59:11 2023 +$ Output To : C:\Users\trbrooks\AppData\Local\Temp\2\ +$ *************************************************************************** +$ +PARAM,PRGPST,NO +PARAM,POST,-1 +PARAM,OGEOM,NO +PARAM,AUTOSPC,YES +PARAM,GRDPNT,0 +EIGRL 1 3000. 1 0 MAX +CORD2C 1 0 0. 0. 0. 0. 0. 1.+FEMAPC1 ++FEMAPC1 1. 0. 1. +CORD2S 2 0 0. 0. 0. 0. 0. 1.+FEMAPC2 ++FEMAPC2 1. 0. 1. +$ Femap Load Set 1 : Untitled +FORCE 1 71 0 1. 0. -1000. 0. +FORCE 1 72 0 1. 0. -2000. 0. +FORCE 1 73 0 1. 0. -2000. 0. +FORCE 1 74 0 1. 0. -2000. 0. +FORCE 1 75 0 1. 0. -2000. 0. +FORCE 1 76 0 1. 0. -2000. 0. +FORCE 1 77 0 1. 0. -2000. 0. +FORCE 1 78 0 1. 0. -2000. 0. +FORCE 1 79 0 1. 0. -2000. 0. +FORCE 1 80 0 1. 0. -2000. 0. +FORCE 1 81 0 1. 0. -2000. 0. +FORCE 1 82 0 1. 0. -2000. 0. +FORCE 1 83 0 1. 0. -2000. 0. +FORCE 1 84 0 1. 0. -2000. 0. +FORCE 1 85 0 1. 0. -2000. 0. +FORCE 1 86 0 1. 0. -2000. 0. +FORCE 1 87 0 1. 0. -2000. 0. +FORCE 1 88 0 1. 0. -2000. 0. +FORCE 1 89 0 1. 0. -2000. 0. +FORCE 1 90 0 1. 0. -2000. 0. +FORCE 1 91 0 1. 0. -2000. 0. +FORCE 1 92 0 1. 0. -2000. 0. +FORCE 1 93 0 1. 0. -2000. 0. +FORCE 1 94 0 1. 0. -2000. 0. +FORCE 1 95 0 1. 0. -2000. 0. +FORCE 1 96 0 1. 0. -2000. 0. +FORCE 1 97 0 1. 0. -2000. 0. +FORCE 1 98 0 1. 0. -2000. 0. +FORCE 1 99 0 1. 0. -2000. 0. +FORCE 1 100 0 1. 0. -2000. 0. +FORCE 1 101 0 1. 0. -2000. 0. +FORCE 1 102 0 1. 0. -2000. 0. +FORCE 1 103 0 1. 0. -2000. 0. +FORCE 1 104 0 1. 0. -2000. 0. +FORCE 1 105 0 1. 0. -2000. 0. +FORCE 1 106 0 1. 0. -2000. 0. +FORCE 1 107 0 1. 0. -2000. 0. +FORCE 1 108 0 1. 0. -2000. 0. +FORCE 1 109 0 1. 0. -2000. 0. +FORCE 1 110 0 1. 0. -2000. 0. +FORCE 1 111 0 1. 0. -1000. 0. +FORCE 1 1 0 1. 0. 1000. 0. +FORCE 1 2 0 1. 0. 2000. 0. +FORCE 1 3 0 1. 0. 2000. 0. +FORCE 1 4 0 1. 0. 2000. 0. +FORCE 1 5 0 1. 0. 2000. 0. +FORCE 1 6 0 1. 0. 2000. 0. +FORCE 1 7 0 1. 0. 2000. 0. +FORCE 1 8 0 1. 0. 2000. 0. +FORCE 1 9 0 1. 0. 2000. 0. +FORCE 1 10 0 1. 0. 2000. 0. +FORCE 1 11 0 1. 0. 2000. 0. +FORCE 1 12 0 1. 0. 2000. 0. +FORCE 1 13 0 1. 0. 2000. 0. +FORCE 1 14 0 1. 0. 2000. 0. +FORCE 1 15 0 1. 0. 2000. 0. +FORCE 1 16 0 1. 0. 2000. 0. +FORCE 1 17 0 1. 0. 2000. 0. +FORCE 1 18 0 1. 0. 2000. 0. +FORCE 1 19 0 1. 0. 2000. 0. +FORCE 1 20 0 1. 0. 2000. 0. +FORCE 1 21 0 1. 0. 2000. 0. +FORCE 1 22 0 1. 0. 2000. 0. +FORCE 1 23 0 1. 0. 2000. 0. +FORCE 1 24 0 1. 0. 2000. 0. +FORCE 1 25 0 1. 0. 2000. 0. +FORCE 1 26 0 1. 0. 2000. 0. +FORCE 1 27 0 1. 0. 2000. 0. +FORCE 1 28 0 1. 0. 2000. 0. +FORCE 1 29 0 1. 0. 2000. 0. +FORCE 1 30 0 1. 0. 2000. 0. +FORCE 1 31 0 1. 0. 2000. 0. +FORCE 1 32 0 1. 0. 2000. 0. +FORCE 1 33 0 1. 0. 2000. 0. +FORCE 1 34 0 1. 0. 2000. 0. +FORCE 1 35 0 1. 0. 2000. 0. +FORCE 1 36 0 1. 0. 2000. 0. +FORCE 1 37 0 1. 0. 2000. 0. +FORCE 1 38 0 1. 0. 2000. 0. +FORCE 1 39 0 1. 0. 2000. 0. +FORCE 1 40 0 1. 0. 2000. 0. +FORCE 1 41 0 1. 0. 1000. 0. +$ Femap Constraint Set 1 : Untitled +SPC1 1 123456 1 +SPC1 1 236 2 +SPC1 1 236 3 +SPC1 1 236 4 +SPC1 1 236 5 +SPC1 1 236 6 +SPC1 1 236 7 +SPC1 1 236 8 +SPC1 1 236 9 +SPC1 1 236 10 +SPC1 1 236 11 +SPC1 1 236 12 +SPC1 1 236 13 +SPC1 1 236 14 +SPC1 1 236 15 +SPC1 1 236 16 +SPC1 1 236 17 +SPC1 1 236 18 +SPC1 1 236 19 +SPC1 1 236 20 +SPC1 1 236 21 +SPC1 1 236 22 +SPC1 1 236 23 +SPC1 1 236 24 +SPC1 1 236 25 +SPC1 1 236 26 +SPC1 1 236 27 +SPC1 1 236 28 +SPC1 1 236 29 +SPC1 1 236 30 +SPC1 1 236 31 +SPC1 1 236 32 +SPC1 1 236 33 +SPC1 1 236 34 +SPC1 1 236 35 +SPC1 1 236 36 +SPC1 1 236 37 +SPC1 1 236 38 +SPC1 1 236 39 +SPC1 1 236 40 +SPC1 1 236 41 +SPC1 1 36 42 +SPC1 1 36 43 +SPC1 1 36 44 +SPC1 1 36 45 +SPC1 1 36 46 +SPC1 1 36 47 +SPC1 1 36 48 +SPC1 1 36 49 +SPC1 1 36 50 +SPC1 1 36 51 +SPC1 1 36 52 +SPC1 1 36 53 +SPC1 1 36 54 +SPC1 1 36 55 +SPC1 1 36 56 +SPC1 1 36 57 +SPC1 1 36 58 +SPC1 1 36 59 +SPC1 1 36 60 +SPC1 1 36 61 +SPC1 1 36 62 +SPC1 1 36 63 +SPC1 1 36 64 +SPC1 1 36 65 +SPC1 1 36 66 +SPC1 1 36 67 +SPC1 1 36 68 +SPC1 1 36 69 +SPC1 1 36 70 +SPC1 1 36 71 +SPC1 1 36 72 +SPC1 1 36 73 +SPC1 1 36 74 +SPC1 1 36 75 +SPC1 1 36 76 +SPC1 1 36 77 +SPC1 1 36 78 +SPC1 1 36 79 +SPC1 1 36 80 +SPC1 1 36 81 +SPC1 1 36 82 +SPC1 1 36 83 +SPC1 1 36 84 +SPC1 1 36 85 +SPC1 1 36 86 +SPC1 1 36 87 +SPC1 1 36 88 +SPC1 1 36 89 +SPC1 1 36 90 +SPC1 1 36 91 +SPC1 1 36 92 +SPC1 1 36 93 +SPC1 1 36 94 +SPC1 1 36 95 +SPC1 1 36 96 +SPC1 1 36 97 +SPC1 1 36 98 +SPC1 1 36 99 +SPC1 1 36 100 +SPC1 1 36 101 +SPC1 1 36 102 +SPC1 1 36 103 +SPC1 1 36 104 +SPC1 1 36 105 +SPC1 1 36 106 +SPC1 1 36 107 +SPC1 1 36 108 +SPC1 1 36 109 +SPC1 1 36 110 +SPC1 1 36 111 +SPC1 1 36 112 +SPC1 1 36 113 +SPC1 1 36 114 +SPC1 1 36 115 +SPC1 1 36 116 +SPC1 1 36 117 +SPC1 1 36 118 +SPC1 1 36 119 +SPC1 1 36 120 +SPC1 1 36 121 +SPC1 1 36 122 +SPC1 1 36 123 +SPC1 1 36 124 +SPC1 1 36 125 +SPC1 1 36 126 +SPC1 1 36 127 +SPC1 1 36 128 +SPC1 1 36 129 +SPC1 1 36 130 +SPC1 1 36 131 +SPC1 1 36 132 +SPC1 1 36 133 +SPC1 1 36 134 +SPC1 1 36 135 +SPC1 1 36 136 +SPC1 1 36 137 +SPC1 1 36 138 +SPC1 1 36 139 +SPC1 1 36 140 +$ Femap Property 1 : PLATE Property +PSHELL 1 1 .02 1 1 0. +$ Femap Material 1 : ISOTROPIC Material +MAT1 1 2.05+11 .3 0. 0. 0. +GRID 1 0 0. 0. 0. 0 +GRID 2 0 .1 0. 0. 0 +GRID 3 0 .2 0. 0. 0 +GRID 4 0 .3 0. 0. 0 +GRID 5 0 .4 0. 0. 0 +GRID 6 0 .5 0. 0. 0 +GRID 7 0 .6 0. 0. 0 +GRID 8 0 .7 0. 0. 0 +GRID 9 0 .8 0. 0. 0 +GRID 10 0 .9 0. 0. 0 +GRID 11 0 1. 0. 0. 0 +GRID 12 0 1.1 0. 0. 0 +GRID 13 0 1.2 0. 0. 0 +GRID 14 0 1.3 0. 0. 0 +GRID 15 0 1.4 0. 0. 0 +GRID 16 0 1.5 0. 0. 0 +GRID 17 0 1.6 0. 0. 0 +GRID 18 0 1.7 0. 0. 0 +GRID 19 0 1.8 0. 0. 0 +GRID 20 0 1.9 0. 0. 0 +GRID 21 0 2. 0. 0. 0 +GRID 22 0 2.1 0. 0. 0 +GRID 23 0 2.2 0. 0. 0 +GRID 24 0 2.3 0. 0. 0 +GRID 25 0 2.4 0. 0. 0 +GRID 26 0 2.5 0. 0. 0 +GRID 27 0 2.6 0. 0. 0 +GRID 28 0 2.7 0. 0. 0 +GRID 29 0 2.8 0. 0. 0 +GRID 30 0 2.9 0. 0. 0 +GRID 31 0 3. 0. 0. 0 +GRID 32 0 3.1 0. 0. 0 +GRID 33 0 3.2 0. 0. 0 +GRID 34 0 3.3 0. 0. 0 +GRID 35 0 3.4 0. 0. 0 +GRID 36 0 3.5 0. 0. 0 +GRID 37 0 3.6 0. 0. 0 +GRID 38 0 3.7 0. 0. 0 +GRID 39 0 3.8 0. 0. 0 +GRID 40 0 3.9 0. 0. 0 +GRID 41 0 4. 0. 0. 0 +GRID 42 0 4. .1 0. 0 +GRID 43 0 4. .2 0. 0 +GRID 44 0 4. .3 0. 0 +GRID 45 0 4. .4 0. 0 +GRID 46 0 4. .5 0. 0 +GRID 47 0 4. .6 0. 0 +GRID 48 0 4. .7 0. 0 +GRID 49 0 4. .8 0. 0 +GRID 50 0 4. .9 0. 0 +GRID 51 0 4. 1. 0. 0 +GRID 52 0 4. 1.1 0. 0 +GRID 53 0 4. 1.2 0. 0 +GRID 54 0 4. 1.3 0. 0 +GRID 55 0 4. 1.4 0. 0 +GRID 56 0 4. 1.5 0. 0 +GRID 57 0 4. 1.6 0. 0 +GRID 58 0 4. 1.7 0. 0 +GRID 59 0 4. 1.8 0. 0 +GRID 60 0 4. 1.9 0. 0 +GRID 61 0 4. 2. 0. 0 +GRID 62 0 4. 2.1 0. 0 +GRID 63 0 4. 2.2 0. 0 +GRID 64 0 4. 2.3 0. 0 +GRID 65 0 4. 2.4 0. 0 +GRID 66 0 4. 2.5 0. 0 +GRID 67 0 4. 2.6 0. 0 +GRID 68 0 4. 2.7 0. 0 +GRID 69 0 4. 2.8 0. 0 +GRID 70 0 4. 2.9 0. 0 +GRID 71 0 4. 3. 0. 0 +GRID 72 0 3.9 3. 0. 0 +GRID 73 0 3.8 3. 0. 0 +GRID 74 0 3.7 3. 0. 0 +GRID 75 0 3.6 3. 0. 0 +GRID 76 0 3.5 3. 0. 0 +GRID 77 0 3.4 3. 0. 0 +GRID 78 0 3.3 3. 0. 0 +GRID 79 0 3.2 3. 0. 0 +GRID 80 0 3.1 3. 0. 0 +GRID 81 0 3. 3. 0. 0 +GRID 82 0 2.9 3. 0. 0 +GRID 83 0 2.8 3. 0. 0 +GRID 84 0 2.7 3. 0. 0 +GRID 85 0 2.6 3. 0. 0 +GRID 86 0 2.5 3. 0. 0 +GRID 87 0 2.4 3. 0. 0 +GRID 88 0 2.3 3. 0. 0 +GRID 89 0 2.2 3. 0. 0 +GRID 90 0 2.1 3. 0. 0 +GRID 91 0 2. 3. 0. 0 +GRID 92 0 1.9 3. 0. 0 +GRID 93 0 1.8 3. 0. 0 +GRID 94 0 1.7 3. 0. 0 +GRID 95 0 1.6 3. 0. 0 +GRID 96 0 1.5 3. 0. 0 +GRID 97 0 1.4 3. 0. 0 +GRID 98 0 1.3 3. 0. 0 +GRID 99 0 1.2 3. 0. 0 +GRID 100 0 1.1 3. 0. 0 +GRID 101 0 1. 3. 0. 0 +GRID 102 0 .9 3. 0. 0 +GRID 103 0 .8 3. 0. 0 +GRID 104 0 .7 3. 0. 0 +GRID 105 0 .6 3. 0. 0 +GRID 106 0 .5 3. 0. 0 +GRID 107 0 .4 3. 0. 0 +GRID 108 0 .3 3. 0. 0 +GRID 109 0 .2 3. 0. 0 +GRID 110 0 .1 3. 0. 0 +GRID 111 0 0. 3. 0. 0 +GRID 112 0 0. 2.9 0. 0 +GRID 113 0 0. 2.8 0. 0 +GRID 114 0 0. 2.7 0. 0 +GRID 115 0 0. 2.6 0. 0 +GRID 116 0 0. 2.5 0. 0 +GRID 117 0 0. 2.4 0. 0 +GRID 118 0 0. 2.3 0. 0 +GRID 119 0 0. 2.2 0. 0 +GRID 120 0 0. 2.1 0. 0 +GRID 121 0 0. 2. 0. 0 +GRID 122 0 0. 1.9 0. 0 +GRID 123 0 0. 1.8 0. 0 +GRID 124 0 0. 1.7 0. 0 +GRID 125 0 0. 1.6 0. 0 +GRID 126 0 0. 1.5 0. 0 +GRID 127 0 0. 1.4 0. 0 +GRID 128 0 0. 1.3 0. 0 +GRID 129 0 0. 1.2 0. 0 +GRID 130 0 0. 1.1 0. 0 +GRID 131 0 0. 1. 0. 0 +GRID 132 0 0. .9 0. 0 +GRID 133 0 0. .8 0. 0 +GRID 134 0 0. .7 0. 0 +GRID 135 0 0. .6 0. 0 +GRID 136 0 0. .5 0. 0 +GRID 137 0 0. .4 0. 0 +GRID 138 0 0. .3 0. 0 +GRID 139 0 0. .2 0. 0 +GRID 140 0 0. .1 0. 0 +GRID 141 0 .1 .1 0. 0 +GRID 142 0 .2 .1 0. 0 +GRID 143 0 .3 .1 0. 0 +GRID 144 0 .4 .1 0. 0 +GRID 145 0 .5 .1 0. 0 +GRID 146 0 .6 .1 0. 0 +GRID 147 0 .7 .1 0. 0 +GRID 148 0 .8 .1 0. 0 +GRID 149 0 .9 .1 0. 0 +GRID 150 0 1. .1 0. 0 +GRID 151 0 1.1 .1 0. 0 +GRID 152 0 1.2 .1 0. 0 +GRID 153 0 1.3 .1 0. 0 +GRID 154 0 1.4 .1 0. 0 +GRID 155 0 1.5 .1 0. 0 +GRID 156 0 1.6 .1 0. 0 +GRID 157 0 1.7 .1 0. 0 +GRID 158 0 1.8 .1 0. 0 +GRID 159 0 1.9 .1 0. 0 +GRID 160 0 2. .1 0. 0 +GRID 161 0 2.1 .1 0. 0 +GRID 162 0 2.2 .1 0. 0 +GRID 163 0 2.3 .1 0. 0 +GRID 164 0 2.4 .1 0. 0 +GRID 165 0 2.5 .1 0. 0 +GRID 166 0 2.6 .1 0. 0 +GRID 167 0 2.7 .1 0. 0 +GRID 168 0 2.8 .1 0. 0 +GRID 169 0 2.9 .1 0. 0 +GRID 170 0 3. .1 0. 0 +GRID 171 0 3.1 .1 0. 0 +GRID 172 0 3.2 .1 0. 0 +GRID 173 0 3.3 .1 0. 0 +GRID 174 0 3.4 .1 0. 0 +GRID 175 0 3.5 .1 0. 0 +GRID 176 0 3.6 .1 0. 0 +GRID 177 0 3.7 .1 0. 0 +GRID 178 0 3.8 .1 0. 0 +GRID 179 0 3.9 .1 0. 0 +GRID 180 0 .1 .2 0. 0 +GRID 181 0 .2 .2 0. 0 +GRID 182 0 .3 .2 0. 0 +GRID 183 0 .4 .2 0. 0 +GRID 184 0 .5 .2 0. 0 +GRID 185 0 .6 .2 0. 0 +GRID 186 0 .7 .2 0. 0 +GRID 187 0 .8 .2 0. 0 +GRID 188 0 .9 .2 0. 0 +GRID 189 0 1. .2 0. 0 +GRID 190 0 1.1 .2 0. 0 +GRID 191 0 1.2 .2 0. 0 +GRID 192 0 1.3 .2 0. 0 +GRID 193 0 1.4 .2 0. 0 +GRID 194 0 1.5 .2 0. 0 +GRID 195 0 1.6 .2 0. 0 +GRID 196 0 1.7 .2 0. 0 +GRID 197 0 1.8 .2 0. 0 +GRID 198 0 1.9 .2 0. 0 +GRID 199 0 2. .2 0. 0 +GRID 200 0 2.1 .2 0. 0 +GRID 201 0 2.2 .2 0. 0 +GRID 202 0 2.3 .2 0. 0 +GRID 203 0 2.4 .2 0. 0 +GRID 204 0 2.5 .2 0. 0 +GRID 205 0 2.6 .2 0. 0 +GRID 206 0 2.7 .2 0. 0 +GRID 207 0 2.8 .2 0. 0 +GRID 208 0 2.9 .2 0. 0 +GRID 209 0 3. .2 0. 0 +GRID 210 0 3.1 .2 0. 0 +GRID 211 0 3.2 .2 0. 0 +GRID 212 0 3.3 .2 0. 0 +GRID 213 0 3.4 .2 0. 0 +GRID 214 0 3.5 .2 0. 0 +GRID 215 0 3.6 .2 0. 0 +GRID 216 0 3.7 .2 0. 0 +GRID 217 0 3.8 .2 0. 0 +GRID 218 0 3.9 .2 0. 0 +GRID 219 0 .1 .3 0. 0 +GRID 220 0 .2 .3 0. 0 +GRID 221 0 .3 .3 0. 0 +GRID 222 0 .4 .3 0. 0 +GRID 223 0 .5 .3 0. 0 +GRID 224 0 .6 .3 0. 0 +GRID 225 0 .7 .3 0. 0 +GRID 226 0 .8 .3 0. 0 +GRID 227 0 .9 .3 0. 0 +GRID 228 0 1. .3 0. 0 +GRID 229 0 1.1 .3 0. 0 +GRID 230 0 1.2 .3 0. 0 +GRID 231 0 1.3 .3 0. 0 +GRID 232 0 1.4 .3 0. 0 +GRID 233 0 1.5 .3 0. 0 +GRID 234 0 1.6 .3 0. 0 +GRID 235 0 1.7 .3 0. 0 +GRID 236 0 1.8 .3 0. 0 +GRID 237 0 1.9 .3 0. 0 +GRID 238 0 2. .3 0. 0 +GRID 239 0 2.1 .3 0. 0 +GRID 240 0 2.2 .3 0. 0 +GRID 241 0 2.3 .3 0. 0 +GRID 242 0 2.4 .3 0. 0 +GRID 243 0 2.5 .3 0. 0 +GRID 244 0 2.6 .3 0. 0 +GRID 245 0 2.7 .3 0. 0 +GRID 246 0 2.8 .3 0. 0 +GRID 247 0 2.9 .3 0. 0 +GRID 248 0 3. .3 0. 0 +GRID 249 0 3.1 .3 0. 0 +GRID 250 0 3.2 .3 0. 0 +GRID 251 0 3.3 .3 0. 0 +GRID 252 0 3.4 .3 0. 0 +GRID 253 0 3.5 .3 0. 0 +GRID 254 0 3.6 .3 0. 0 +GRID 255 0 3.7 .3 0. 0 +GRID 256 0 3.8 .3 0. 0 +GRID 257 0 3.9 .3 0. 0 +GRID 258 0 .1 .4 0. 0 +GRID 259 0 .2 .4 0. 0 +GRID 260 0 .3 .4 0. 0 +GRID 261 0 .4 .4 0. 0 +GRID 262 0 .5 .4 0. 0 +GRID 263 0 .6 .4 0. 0 +GRID 264 0 .7 .4 0. 0 +GRID 265 0 .8 .4 0. 0 +GRID 266 0 .9 .4 0. 0 +GRID 267 0 1. .4 0. 0 +GRID 268 0 1.1 .4 0. 0 +GRID 269 0 1.2 .4 0. 0 +GRID 270 0 1.3 .4 0. 0 +GRID 271 0 1.4 .4 0. 0 +GRID 272 0 1.5 .4 0. 0 +GRID 273 0 1.6 .4 0. 0 +GRID 274 0 1.7 .4 0. 0 +GRID 275 0 1.8 .4 0. 0 +GRID 276 0 1.9 .4 0. 0 +GRID 277 0 2. .4 0. 0 +GRID 278 0 2.1 .4 0. 0 +GRID 279 0 2.2 .4 0. 0 +GRID 280 0 2.3 .4 0. 0 +GRID 281 0 2.4 .4 0. 0 +GRID 282 0 2.5 .4 0. 0 +GRID 283 0 2.6 .4 0. 0 +GRID 284 0 2.7 .4 0. 0 +GRID 285 0 2.8 .4 0. 0 +GRID 286 0 2.9 .4 0. 0 +GRID 287 0 3. .4 0. 0 +GRID 288 0 3.1 .4 0. 0 +GRID 289 0 3.2 .4 0. 0 +GRID 290 0 3.3 .4 0. 0 +GRID 291 0 3.4 .4 0. 0 +GRID 292 0 3.5 .4 0. 0 +GRID 293 0 3.6 .4 0. 0 +GRID 294 0 3.7 .4 0. 0 +GRID 295 0 3.8 .4 0. 0 +GRID 296 0 3.9 .4 0. 0 +GRID 297 0 .1 .5 0. 0 +GRID 298 0 .2 .5 0. 0 +GRID 299 0 .3 .5 0. 0 +GRID 300 0 .4 .5 0. 0 +GRID 301 0 .5 .5 0. 0 +GRID 302 0 .6 .5 0. 0 +GRID 303 0 .7 .5 0. 0 +GRID 304 0 .8 .5 0. 0 +GRID 305 0 .9 .5 0. 0 +GRID 306 0 1. .5 0. 0 +GRID 307 0 1.1 .5 0. 0 +GRID 308 0 1.2 .5 0. 0 +GRID 309 0 1.3 .5 0. 0 +GRID 310 0 1.4 .5 0. 0 +GRID 311 0 1.5 .5 0. 0 +GRID 312 0 1.6 .5 0. 0 +GRID 313 0 1.7 .5 0. 0 +GRID 314 0 1.8 .5 0. 0 +GRID 315 0 1.9 .5 0. 0 +GRID 316 0 2. .5 0. 0 +GRID 317 0 2.1 .5 0. 0 +GRID 318 0 2.2 .5 0. 0 +GRID 319 0 2.3 .5 0. 0 +GRID 320 0 2.4 .5 0. 0 +GRID 321 0 2.5 .5 0. 0 +GRID 322 0 2.6 .5 0. 0 +GRID 323 0 2.7 .5 0. 0 +GRID 324 0 2.8 .5 0. 0 +GRID 325 0 2.9 .5 0. 0 +GRID 326 0 3. .5 0. 0 +GRID 327 0 3.1 .5 0. 0 +GRID 328 0 3.2 .5 0. 0 +GRID 329 0 3.3 .5 0. 0 +GRID 330 0 3.4 .5 0. 0 +GRID 331 0 3.5 .5 0. 0 +GRID 332 0 3.6 .5 0. 0 +GRID 333 0 3.7 .5 0. 0 +GRID 334 0 3.8 .5 0. 0 +GRID 335 0 3.9 .5 0. 0 +GRID 336 0 .1 .6 0. 0 +GRID 337 0 .2 .6 0. 0 +GRID 338 0 .3 .6 0. 0 +GRID 339 0 .4 .6 0. 0 +GRID 340 0 .5 .6 0. 0 +GRID 341 0 .6 .6 0. 0 +GRID 342 0 .7 .6 0. 0 +GRID 343 0 .8 .6 0. 0 +GRID 344 0 .9 .6 0. 0 +GRID 345 0 1. .6 0. 0 +GRID 346 0 1.1 .6 0. 0 +GRID 347 0 1.2 .6 0. 0 +GRID 348 0 1.3 .6 0. 0 +GRID 349 0 1.4 .6 0. 0 +GRID 350 0 1.5 .6 0. 0 +GRID 351 0 1.6 .6 0. 0 +GRID 352 0 1.7 .6 0. 0 +GRID 353 0 1.8 .6 0. 0 +GRID 354 0 1.9 .6 0. 0 +GRID 355 0 2. .6 0. 0 +GRID 356 0 2.1 .6 0. 0 +GRID 357 0 2.2 .6 0. 0 +GRID 358 0 2.3 .6 0. 0 +GRID 359 0 2.4 .6 0. 0 +GRID 360 0 2.5 .6 0. 0 +GRID 361 0 2.6 .6 0. 0 +GRID 362 0 2.7 .6 0. 0 +GRID 363 0 2.8 .6 0. 0 +GRID 364 0 2.9 .6 0. 0 +GRID 365 0 3. .6 0. 0 +GRID 366 0 3.1 .6 0. 0 +GRID 367 0 3.2 .6 0. 0 +GRID 368 0 3.3 .6 0. 0 +GRID 369 0 3.4 .6 0. 0 +GRID 370 0 3.5 .6 0. 0 +GRID 371 0 3.6 .6 0. 0 +GRID 372 0 3.7 .6 0. 0 +GRID 373 0 3.8 .6 0. 0 +GRID 374 0 3.9 .6 0. 0 +GRID 375 0 .1 .7 0. 0 +GRID 376 0 .2 .7 0. 0 +GRID 377 0 .3 .7 0. 0 +GRID 378 0 .4 .7 0. 0 +GRID 379 0 .5 .7 0. 0 +GRID 380 0 .6 .7 0. 0 +GRID 381 0 .7 .7 0. 0 +GRID 382 0 .8 .7 0. 0 +GRID 383 0 .9 .7 0. 0 +GRID 384 0 1. .7 0. 0 +GRID 385 0 1.1 .7 0. 0 +GRID 386 0 1.2 .7 0. 0 +GRID 387 0 1.3 .7 0. 0 +GRID 388 0 1.4 .7 0. 0 +GRID 389 0 1.5 .7 0. 0 +GRID 390 0 1.6 .7 0. 0 +GRID 391 0 1.7 .7 0. 0 +GRID 392 0 1.8 .7 0. 0 +GRID 393 0 1.9 .7 0. 0 +GRID 394 0 2. .7 0. 0 +GRID 395 0 2.1 .7 0. 0 +GRID 396 0 2.2 .7 0. 0 +GRID 397 0 2.3 .7 0. 0 +GRID 398 0 2.4 .7 0. 0 +GRID 399 0 2.5 .7 0. 0 +GRID 400 0 2.6 .7 0. 0 +GRID 401 0 2.7 .7 0. 0 +GRID 402 0 2.8 .7 0. 0 +GRID 403 0 2.9 .7 0. 0 +GRID 404 0 3. .7 0. 0 +GRID 405 0 3.1 .7 0. 0 +GRID 406 0 3.2 .7 0. 0 +GRID 407 0 3.3 .7 0. 0 +GRID 408 0 3.4 .7 0. 0 +GRID 409 0 3.5 .7 0. 0 +GRID 410 0 3.6 .7 0. 0 +GRID 411 0 3.7 .7 0. 0 +GRID 412 0 3.8 .7 0. 0 +GRID 413 0 3.9 .7 0. 0 +GRID 414 0 .1 .8 0. 0 +GRID 415 0 .2 .8 0. 0 +GRID 416 0 .3 .8 0. 0 +GRID 417 0 .4 .8 0. 0 +GRID 418 0 .5 .8 0. 0 +GRID 419 0 .6 .8 0. 0 +GRID 420 0 .7 .8 0. 0 +GRID 421 0 .8 .8 0. 0 +GRID 422 0 .9 .8 0. 0 +GRID 423 0 1. .8 0. 0 +GRID 424 0 1.1 .8 0. 0 +GRID 425 0 1.2 .8 0. 0 +GRID 426 0 1.3 .8 0. 0 +GRID 427 0 1.4 .8 0. 0 +GRID 428 0 1.5 .8 0. 0 +GRID 429 0 1.6 .8 0. 0 +GRID 430 0 1.7 .8 0. 0 +GRID 431 0 1.8 .8 0. 0 +GRID 432 0 1.9 .8 0. 0 +GRID 433 0 2. .8 0. 0 +GRID 434 0 2.1 .8 0. 0 +GRID 435 0 2.2 .8 0. 0 +GRID 436 0 2.3 .8 0. 0 +GRID 437 0 2.4 .8 0. 0 +GRID 438 0 2.5 .8 0. 0 +GRID 439 0 2.6 .8 0. 0 +GRID 440 0 2.7 .8 0. 0 +GRID 441 0 2.8 .8 0. 0 +GRID 442 0 2.9 .8 0. 0 +GRID 443 0 3. .8 0. 0 +GRID 444 0 3.1 .8 0. 0 +GRID 445 0 3.2 .8 0. 0 +GRID 446 0 3.3 .8 0. 0 +GRID 447 0 3.4 .8 0. 0 +GRID 448 0 3.5 .8 0. 0 +GRID 449 0 3.6 .8 0. 0 +GRID 450 0 3.7 .8 0. 0 +GRID 451 0 3.8 .8 0. 0 +GRID 452 0 3.9 .8 0. 0 +GRID 453 0 .1 .9 0. 0 +GRID 454 0 .2 .9 0. 0 +GRID 455 0 .3 .9 0. 0 +GRID 456 0 .4 .9 0. 0 +GRID 457 0 .5 .9 0. 0 +GRID 458 0 .6 .9 0. 0 +GRID 459 0 .7 .9 0. 0 +GRID 460 0 .8 .9 0. 0 +GRID 461 0 .9 .9 0. 0 +GRID 462 0 1. .9 0. 0 +GRID 463 0 1.1 .9 0. 0 +GRID 464 0 1.2 .9 0. 0 +GRID 465 0 1.3 .9 0. 0 +GRID 466 0 1.4 .9 0. 0 +GRID 467 0 1.5 .9 0. 0 +GRID 468 0 1.6 .9 0. 0 +GRID 469 0 1.7 .9 0. 0 +GRID 470 0 1.8 .9 0. 0 +GRID 471 0 1.9 .9 0. 0 +GRID 472 0 2. .9 0. 0 +GRID 473 0 2.1 .9 0. 0 +GRID 474 0 2.2 .9 0. 0 +GRID 475 0 2.3 .9 0. 0 +GRID 476 0 2.4 .9 0. 0 +GRID 477 0 2.5 .9 0. 0 +GRID 478 0 2.6 .9 0. 0 +GRID 479 0 2.7 .9 0. 0 +GRID 480 0 2.8 .9 0. 0 +GRID 481 0 2.9 .9 0. 0 +GRID 482 0 3. .9 0. 0 +GRID 483 0 3.1 .9 0. 0 +GRID 484 0 3.2 .9 0. 0 +GRID 485 0 3.3 .9 0. 0 +GRID 486 0 3.4 .9 0. 0 +GRID 487 0 3.5 .9 0. 0 +GRID 488 0 3.6 .9 0. 0 +GRID 489 0 3.7 .9 0. 0 +GRID 490 0 3.8 .9 0. 0 +GRID 491 0 3.9 .9 0. 0 +GRID 492 0 .1 1. 0. 0 +GRID 493 0 .2 1. 0. 0 +GRID 494 0 .3 1. 0. 0 +GRID 495 0 .4 1. 0. 0 +GRID 496 0 .5 1. 0. 0 +GRID 497 0 .6 1. 0. 0 +GRID 498 0 .7 1. 0. 0 +GRID 499 0 .8 1. 0. 0 +GRID 500 0 .9 1. 0. 0 +GRID 501 0 1. 1. 0. 0 +GRID 502 0 1.1 1. 0. 0 +GRID 503 0 1.2 1. 0. 0 +GRID 504 0 1.3 1. 0. 0 +GRID 505 0 1.4 1. 0. 0 +GRID 506 0 1.5 1. 0. 0 +GRID 507 0 1.6 1. 0. 0 +GRID 508 0 1.7 1. 0. 0 +GRID 509 0 1.8 1. 0. 0 +GRID 510 0 1.9 1. 0. 0 +GRID 511 0 2. 1. 0. 0 +GRID 512 0 2.1 1. 0. 0 +GRID 513 0 2.2 1. 0. 0 +GRID 514 0 2.3 1. 0. 0 +GRID 515 0 2.4 1. 0. 0 +GRID 516 0 2.5 1. 0. 0 +GRID 517 0 2.6 1. 0. 0 +GRID 518 0 2.7 1. 0. 0 +GRID 519 0 2.8 1. 0. 0 +GRID 520 0 2.9 1. 0. 0 +GRID 521 0 3. 1. 0. 0 +GRID 522 0 3.1 1. 0. 0 +GRID 523 0 3.2 1. 0. 0 +GRID 524 0 3.3 1. 0. 0 +GRID 525 0 3.4 1. 0. 0 +GRID 526 0 3.5 1. 0. 0 +GRID 527 0 3.6 1. 0. 0 +GRID 528 0 3.7 1. 0. 0 +GRID 529 0 3.8 1. 0. 0 +GRID 530 0 3.9 1. 0. 0 +GRID 531 0 .1 1.1 0. 0 +GRID 532 0 .2 1.1 0. 0 +GRID 533 0 .3 1.1 0. 0 +GRID 534 0 .4 1.1 0. 0 +GRID 535 0 .5 1.1 0. 0 +GRID 536 0 .6 1.1 0. 0 +GRID 537 0 .7 1.1 0. 0 +GRID 538 0 .8 1.1 0. 0 +GRID 539 0 .9 1.1 0. 0 +GRID 540 0 1. 1.1 0. 0 +GRID 541 0 1.1 1.1 0. 0 +GRID 542 0 1.2 1.1 0. 0 +GRID 543 0 1.3 1.1 0. 0 +GRID 544 0 1.4 1.1 0. 0 +GRID 545 0 1.5 1.1 0. 0 +GRID 546 0 1.6 1.1 0. 0 +GRID 547 0 1.7 1.1 0. 0 +GRID 548 0 1.8 1.1 0. 0 +GRID 549 0 1.9 1.1 0. 0 +GRID 550 0 2. 1.1 0. 0 +GRID 551 0 2.1 1.1 0. 0 +GRID 552 0 2.2 1.1 0. 0 +GRID 553 0 2.3 1.1 0. 0 +GRID 554 0 2.4 1.1 0. 0 +GRID 555 0 2.5 1.1 0. 0 +GRID 556 0 2.6 1.1 0. 0 +GRID 557 0 2.7 1.1 0. 0 +GRID 558 0 2.8 1.1 0. 0 +GRID 559 0 2.9 1.1 0. 0 +GRID 560 0 3. 1.1 0. 0 +GRID 561 0 3.1 1.1 0. 0 +GRID 562 0 3.2 1.1 0. 0 +GRID 563 0 3.3 1.1 0. 0 +GRID 564 0 3.4 1.1 0. 0 +GRID 565 0 3.5 1.1 0. 0 +GRID 566 0 3.6 1.1 0. 0 +GRID 567 0 3.7 1.1 0. 0 +GRID 568 0 3.8 1.1 0. 0 +GRID 569 0 3.9 1.1 0. 0 +GRID 570 0 .1 1.2 0. 0 +GRID 571 0 .2 1.2 0. 0 +GRID 572 0 .3 1.2 0. 0 +GRID 573 0 .4 1.2 0. 0 +GRID 574 0 .5 1.2 0. 0 +GRID 575 0 .6 1.2 0. 0 +GRID 576 0 .7 1.2 0. 0 +GRID 577 0 .8 1.2 0. 0 +GRID 578 0 .9 1.2 0. 0 +GRID 579 0 1. 1.2 0. 0 +GRID 580 0 1.1 1.2 0. 0 +GRID 581 0 1.2 1.2 0. 0 +GRID 582 0 1.3 1.2 0. 0 +GRID 583 0 1.4 1.2 0. 0 +GRID 584 0 1.5 1.2 0. 0 +GRID 585 0 1.6 1.2 0. 0 +GRID 586 0 1.7 1.2 0. 0 +GRID 587 0 1.8 1.2 0. 0 +GRID 588 0 1.9 1.2 0. 0 +GRID 589 0 2. 1.2 0. 0 +GRID 590 0 2.1 1.2 0. 0 +GRID 591 0 2.2 1.2 0. 0 +GRID 592 0 2.3 1.2 0. 0 +GRID 593 0 2.4 1.2 0. 0 +GRID 594 0 2.5 1.2 0. 0 +GRID 595 0 2.6 1.2 0. 0 +GRID 596 0 2.7 1.2 0. 0 +GRID 597 0 2.8 1.2 0. 0 +GRID 598 0 2.9 1.2 0. 0 +GRID 599 0 3. 1.2 0. 0 +GRID 600 0 3.1 1.2 0. 0 +GRID 601 0 3.2 1.2 0. 0 +GRID 602 0 3.3 1.2 0. 0 +GRID 603 0 3.4 1.2 0. 0 +GRID 604 0 3.5 1.2 0. 0 +GRID 605 0 3.6 1.2 0. 0 +GRID 606 0 3.7 1.2 0. 0 +GRID 607 0 3.8 1.2 0. 0 +GRID 608 0 3.9 1.2 0. 0 +GRID 609 0 .1 1.3 0. 0 +GRID 610 0 .2 1.3 0. 0 +GRID 611 0 .3 1.3 0. 0 +GRID 612 0 .4 1.3 0. 0 +GRID 613 0 .5 1.3 0. 0 +GRID 614 0 .6 1.3 0. 0 +GRID 615 0 .7 1.3 0. 0 +GRID 616 0 .8 1.3 0. 0 +GRID 617 0 .9 1.3 0. 0 +GRID 618 0 1. 1.3 0. 0 +GRID 619 0 1.1 1.3 0. 0 +GRID 620 0 1.2 1.3 0. 0 +GRID 621 0 1.3 1.3 0. 0 +GRID 622 0 1.4 1.3 0. 0 +GRID 623 0 1.5 1.3 0. 0 +GRID 624 0 1.6 1.3 0. 0 +GRID 625 0 1.7 1.3 0. 0 +GRID 626 0 1.8 1.3 0. 0 +GRID 627 0 1.9 1.3 0. 0 +GRID 628 0 2. 1.3 0. 0 +GRID 629 0 2.1 1.3 0. 0 +GRID 630 0 2.2 1.3 0. 0 +GRID 631 0 2.3 1.3 0. 0 +GRID 632 0 2.4 1.3 0. 0 +GRID 633 0 2.5 1.3 0. 0 +GRID 634 0 2.6 1.3 0. 0 +GRID 635 0 2.7 1.3 0. 0 +GRID 636 0 2.8 1.3 0. 0 +GRID 637 0 2.9 1.3 0. 0 +GRID 638 0 3. 1.3 0. 0 +GRID 639 0 3.1 1.3 0. 0 +GRID 640 0 3.2 1.3 0. 0 +GRID 641 0 3.3 1.3 0. 0 +GRID 642 0 3.4 1.3 0. 0 +GRID 643 0 3.5 1.3 0. 0 +GRID 644 0 3.6 1.3 0. 0 +GRID 645 0 3.7 1.3 0. 0 +GRID 646 0 3.8 1.3 0. 0 +GRID 647 0 3.9 1.3 0. 0 +GRID 648 0 .1 1.4 0. 0 +GRID 649 0 .2 1.4 0. 0 +GRID 650 0 .3 1.4 0. 0 +GRID 651 0 .4 1.4 0. 0 +GRID 652 0 .5 1.4 0. 0 +GRID 653 0 .6 1.4 0. 0 +GRID 654 0 .7 1.4 0. 0 +GRID 655 0 .8 1.4 0. 0 +GRID 656 0 .9 1.4 0. 0 +GRID 657 0 1. 1.4 0. 0 +GRID 658 0 1.1 1.4 0. 0 +GRID 659 0 1.2 1.4 0. 0 +GRID 660 0 1.3 1.4 0. 0 +GRID 661 0 1.4 1.4 0. 0 +GRID 662 0 1.5 1.4 0. 0 +GRID 663 0 1.6 1.4 0. 0 +GRID 664 0 1.7 1.4 0. 0 +GRID 665 0 1.8 1.4 0. 0 +GRID 666 0 1.9 1.4 0. 0 +GRID 667 0 2. 1.4 0. 0 +GRID 668 0 2.1 1.4 0. 0 +GRID 669 0 2.2 1.4 0. 0 +GRID 670 0 2.3 1.4 0. 0 +GRID 671 0 2.4 1.4 0. 0 +GRID 672 0 2.5 1.4 0. 0 +GRID 673 0 2.6 1.4 0. 0 +GRID 674 0 2.7 1.4 0. 0 +GRID 675 0 2.8 1.4 0. 0 +GRID 676 0 2.9 1.4 0. 0 +GRID 677 0 3. 1.4 0. 0 +GRID 678 0 3.1 1.4 0. 0 +GRID 679 0 3.2 1.4 0. 0 +GRID 680 0 3.3 1.4 0. 0 +GRID 681 0 3.4 1.4 0. 0 +GRID 682 0 3.5 1.4 0. 0 +GRID 683 0 3.6 1.4 0. 0 +GRID 684 0 3.7 1.4 0. 0 +GRID 685 0 3.8 1.4 0. 0 +GRID 686 0 3.9 1.4 0. 0 +GRID 687 0 .1 1.5 0. 0 +GRID 688 0 .2 1.5 0. 0 +GRID 689 0 .3 1.5 0. 0 +GRID 690 0 .4 1.5 0. 0 +GRID 691 0 .5 1.5 0. 0 +GRID 692 0 .6 1.5 0. 0 +GRID 693 0 .7 1.5 0. 0 +GRID 694 0 .8 1.5 0. 0 +GRID 695 0 .9 1.5 0. 0 +GRID 696 0 1. 1.5 0. 0 +GRID 697 0 1.1 1.5 0. 0 +GRID 698 0 1.2 1.5 0. 0 +GRID 699 0 1.3 1.5 0. 0 +GRID 700 0 1.4 1.5 0. 0 +GRID 701 0 1.5 1.5 0. 0 +GRID 702 0 1.6 1.5 0. 0 +GRID 703 0 1.7 1.5 0. 0 +GRID 704 0 1.8 1.5 0. 0 +GRID 705 0 1.9 1.5 0. 0 +GRID 706 0 2. 1.5 0. 0 +GRID 707 0 2.1 1.5 0. 0 +GRID 708 0 2.2 1.5 0. 0 +GRID 709 0 2.3 1.5 0. 0 +GRID 710 0 2.4 1.5 0. 0 +GRID 711 0 2.5 1.5 0. 0 +GRID 712 0 2.6 1.5 0. 0 +GRID 713 0 2.7 1.5 0. 0 +GRID 714 0 2.8 1.5 0. 0 +GRID 715 0 2.9 1.5 0. 0 +GRID 716 0 3. 1.5 0. 0 +GRID 717 0 3.1 1.5 0. 0 +GRID 718 0 3.2 1.5 0. 0 +GRID 719 0 3.3 1.5 0. 0 +GRID 720 0 3.4 1.5 0. 0 +GRID 721 0 3.5 1.5 0. 0 +GRID 722 0 3.6 1.5 0. 0 +GRID 723 0 3.7 1.5 0. 0 +GRID 724 0 3.8 1.5 0. 0 +GRID 725 0 3.9 1.5 0. 0 +GRID 726 0 .1 1.6 0. 0 +GRID 727 0 .2 1.6 0. 0 +GRID 728 0 .3 1.6 0. 0 +GRID 729 0 .4 1.6 0. 0 +GRID 730 0 .5 1.6 0. 0 +GRID 731 0 .6 1.6 0. 0 +GRID 732 0 .7 1.6 0. 0 +GRID 733 0 .8 1.6 0. 0 +GRID 734 0 .9 1.6 0. 0 +GRID 735 0 1. 1.6 0. 0 +GRID 736 0 1.1 1.6 0. 0 +GRID 737 0 1.2 1.6 0. 0 +GRID 738 0 1.3 1.6 0. 0 +GRID 739 0 1.4 1.6 0. 0 +GRID 740 0 1.5 1.6 0. 0 +GRID 741 0 1.6 1.6 0. 0 +GRID 742 0 1.7 1.6 0. 0 +GRID 743 0 1.8 1.6 0. 0 +GRID 744 0 1.9 1.6 0. 0 +GRID 745 0 2. 1.6 0. 0 +GRID 746 0 2.1 1.6 0. 0 +GRID 747 0 2.2 1.6 0. 0 +GRID 748 0 2.3 1.6 0. 0 +GRID 749 0 2.4 1.6 0. 0 +GRID 750 0 2.5 1.6 0. 0 +GRID 751 0 2.6 1.6 0. 0 +GRID 752 0 2.7 1.6 0. 0 +GRID 753 0 2.8 1.6 0. 0 +GRID 754 0 2.9 1.6 0. 0 +GRID 755 0 3. 1.6 0. 0 +GRID 756 0 3.1 1.6 0. 0 +GRID 757 0 3.2 1.6 0. 0 +GRID 758 0 3.3 1.6 0. 0 +GRID 759 0 3.4 1.6 0. 0 +GRID 760 0 3.5 1.6 0. 0 +GRID 761 0 3.6 1.6 0. 0 +GRID 762 0 3.7 1.6 0. 0 +GRID 763 0 3.8 1.6 0. 0 +GRID 764 0 3.9 1.6 0. 0 +GRID 765 0 .1 1.7 0. 0 +GRID 766 0 .2 1.7 0. 0 +GRID 767 0 .3 1.7 0. 0 +GRID 768 0 .4 1.7 0. 0 +GRID 769 0 .5 1.7 0. 0 +GRID 770 0 .6 1.7 0. 0 +GRID 771 0 .7 1.7 0. 0 +GRID 772 0 .8 1.7 0. 0 +GRID 773 0 .9 1.7 0. 0 +GRID 774 0 1. 1.7 0. 0 +GRID 775 0 1.1 1.7 0. 0 +GRID 776 0 1.2 1.7 0. 0 +GRID 777 0 1.3 1.7 0. 0 +GRID 778 0 1.4 1.7 0. 0 +GRID 779 0 1.5 1.7 0. 0 +GRID 780 0 1.6 1.7 0. 0 +GRID 781 0 1.7 1.7 0. 0 +GRID 782 0 1.8 1.7 0. 0 +GRID 783 0 1.9 1.7 0. 0 +GRID 784 0 2. 1.7 0. 0 +GRID 785 0 2.1 1.7 0. 0 +GRID 786 0 2.2 1.7 0. 0 +GRID 787 0 2.3 1.7 0. 0 +GRID 788 0 2.4 1.7 0. 0 +GRID 789 0 2.5 1.7 0. 0 +GRID 790 0 2.6 1.7 0. 0 +GRID 791 0 2.7 1.7 0. 0 +GRID 792 0 2.8 1.7 0. 0 +GRID 793 0 2.9 1.7 0. 0 +GRID 794 0 3. 1.7 0. 0 +GRID 795 0 3.1 1.7 0. 0 +GRID 796 0 3.2 1.7 0. 0 +GRID 797 0 3.3 1.7 0. 0 +GRID 798 0 3.4 1.7 0. 0 +GRID 799 0 3.5 1.7 0. 0 +GRID 800 0 3.6 1.7 0. 0 +GRID 801 0 3.7 1.7 0. 0 +GRID 802 0 3.8 1.7 0. 0 +GRID 803 0 3.9 1.7 0. 0 +GRID 804 0 .1 1.8 0. 0 +GRID 805 0 .2 1.8 0. 0 +GRID 806 0 .3 1.8 0. 0 +GRID 807 0 .4 1.8 0. 0 +GRID 808 0 .5 1.8 0. 0 +GRID 809 0 .6 1.8 0. 0 +GRID 810 0 .7 1.8 0. 0 +GRID 811 0 .8 1.8 0. 0 +GRID 812 0 .9 1.8 0. 0 +GRID 813 0 1. 1.8 0. 0 +GRID 814 0 1.1 1.8 0. 0 +GRID 815 0 1.2 1.8 0. 0 +GRID 816 0 1.3 1.8 0. 0 +GRID 817 0 1.4 1.8 0. 0 +GRID 818 0 1.5 1.8 0. 0 +GRID 819 0 1.6 1.8 0. 0 +GRID 820 0 1.7 1.8 0. 0 +GRID 821 0 1.8 1.8 0. 0 +GRID 822 0 1.9 1.8 0. 0 +GRID 823 0 2. 1.8 0. 0 +GRID 824 0 2.1 1.8 0. 0 +GRID 825 0 2.2 1.8 0. 0 +GRID 826 0 2.3 1.8 0. 0 +GRID 827 0 2.4 1.8 0. 0 +GRID 828 0 2.5 1.8 0. 0 +GRID 829 0 2.6 1.8 0. 0 +GRID 830 0 2.7 1.8 0. 0 +GRID 831 0 2.8 1.8 0. 0 +GRID 832 0 2.9 1.8 0. 0 +GRID 833 0 3. 1.8 0. 0 +GRID 834 0 3.1 1.8 0. 0 +GRID 835 0 3.2 1.8 0. 0 +GRID 836 0 3.3 1.8 0. 0 +GRID 837 0 3.4 1.8 0. 0 +GRID 838 0 3.5 1.8 0. 0 +GRID 839 0 3.6 1.8 0. 0 +GRID 840 0 3.7 1.8 0. 0 +GRID 841 0 3.8 1.8 0. 0 +GRID 842 0 3.9 1.8 0. 0 +GRID 843 0 .1 1.9 0. 0 +GRID 844 0 .2 1.9 0. 0 +GRID 845 0 .3 1.9 0. 0 +GRID 846 0 .4 1.9 0. 0 +GRID 847 0 .5 1.9 0. 0 +GRID 848 0 .6 1.9 0. 0 +GRID 849 0 .7 1.9 0. 0 +GRID 850 0 .8 1.9 0. 0 +GRID 851 0 .9 1.9 0. 0 +GRID 852 0 1. 1.9 0. 0 +GRID 853 0 1.1 1.9 0. 0 +GRID 854 0 1.2 1.9 0. 0 +GRID 855 0 1.3 1.9 0. 0 +GRID 856 0 1.4 1.9 0. 0 +GRID 857 0 1.5 1.9 0. 0 +GRID 858 0 1.6 1.9 0. 0 +GRID 859 0 1.7 1.9 0. 0 +GRID 860 0 1.8 1.9 0. 0 +GRID 861 0 1.9 1.9 0. 0 +GRID 862 0 2. 1.9 0. 0 +GRID 863 0 2.1 1.9 0. 0 +GRID 864 0 2.2 1.9 0. 0 +GRID 865 0 2.3 1.9 0. 0 +GRID 866 0 2.4 1.9 0. 0 +GRID 867 0 2.5 1.9 0. 0 +GRID 868 0 2.6 1.9 0. 0 +GRID 869 0 2.7 1.9 0. 0 +GRID 870 0 2.8 1.9 0. 0 +GRID 871 0 2.9 1.9 0. 0 +GRID 872 0 3. 1.9 0. 0 +GRID 873 0 3.1 1.9 0. 0 +GRID 874 0 3.2 1.9 0. 0 +GRID 875 0 3.3 1.9 0. 0 +GRID 876 0 3.4 1.9 0. 0 +GRID 877 0 3.5 1.9 0. 0 +GRID 878 0 3.6 1.9 0. 0 +GRID 879 0 3.7 1.9 0. 0 +GRID 880 0 3.8 1.9 0. 0 +GRID 881 0 3.9 1.9 0. 0 +GRID 882 0 .1 2. 0. 0 +GRID 883 0 .2 2. 0. 0 +GRID 884 0 .3 2. 0. 0 +GRID 885 0 .4 2. 0. 0 +GRID 886 0 .5 2. 0. 0 +GRID 887 0 .6 2. 0. 0 +GRID 888 0 .7 2. 0. 0 +GRID 889 0 .8 2. 0. 0 +GRID 890 0 .9 2. 0. 0 +GRID 891 0 1. 2. 0. 0 +GRID 892 0 1.1 2. 0. 0 +GRID 893 0 1.2 2. 0. 0 +GRID 894 0 1.3 2. 0. 0 +GRID 895 0 1.4 2. 0. 0 +GRID 896 0 1.5 2. 0. 0 +GRID 897 0 1.6 2. 0. 0 +GRID 898 0 1.7 2. 0. 0 +GRID 899 0 1.8 2. 0. 0 +GRID 900 0 1.9 2. 0. 0 +GRID 901 0 2. 2. 0. 0 +GRID 902 0 2.1 2. 0. 0 +GRID 903 0 2.2 2. 0. 0 +GRID 904 0 2.3 2. 0. 0 +GRID 905 0 2.4 2. 0. 0 +GRID 906 0 2.5 2. 0. 0 +GRID 907 0 2.6 2. 0. 0 +GRID 908 0 2.7 2. 0. 0 +GRID 909 0 2.8 2. 0. 0 +GRID 910 0 2.9 2. 0. 0 +GRID 911 0 3. 2. 0. 0 +GRID 912 0 3.1 2. 0. 0 +GRID 913 0 3.2 2. 0. 0 +GRID 914 0 3.3 2. 0. 0 +GRID 915 0 3.4 2. 0. 0 +GRID 916 0 3.5 2. 0. 0 +GRID 917 0 3.6 2. 0. 0 +GRID 918 0 3.7 2. 0. 0 +GRID 919 0 3.8 2. 0. 0 +GRID 920 0 3.9 2. 0. 0 +GRID 921 0 .1 2.1 0. 0 +GRID 922 0 .2 2.1 0. 0 +GRID 923 0 .3 2.1 0. 0 +GRID 924 0 .4 2.1 0. 0 +GRID 925 0 .5 2.1 0. 0 +GRID 926 0 .6 2.1 0. 0 +GRID 927 0 .7 2.1 0. 0 +GRID 928 0 .8 2.1 0. 0 +GRID 929 0 .9 2.1 0. 0 +GRID 930 0 1. 2.1 0. 0 +GRID 931 0 1.1 2.1 0. 0 +GRID 932 0 1.2 2.1 0. 0 +GRID 933 0 1.3 2.1 0. 0 +GRID 934 0 1.4 2.1 0. 0 +GRID 935 0 1.5 2.1 0. 0 +GRID 936 0 1.6 2.1 0. 0 +GRID 937 0 1.7 2.1 0. 0 +GRID 938 0 1.8 2.1 0. 0 +GRID 939 0 1.9 2.1 0. 0 +GRID 940 0 2. 2.1 0. 0 +GRID 941 0 2.1 2.1 0. 0 +GRID 942 0 2.2 2.1 0. 0 +GRID 943 0 2.3 2.1 0. 0 +GRID 944 0 2.4 2.1 0. 0 +GRID 945 0 2.5 2.1 0. 0 +GRID 946 0 2.6 2.1 0. 0 +GRID 947 0 2.7 2.1 0. 0 +GRID 948 0 2.8 2.1 0. 0 +GRID 949 0 2.9 2.1 0. 0 +GRID 950 0 3. 2.1 0. 0 +GRID 951 0 3.1 2.1 0. 0 +GRID 952 0 3.2 2.1 0. 0 +GRID 953 0 3.3 2.1 0. 0 +GRID 954 0 3.4 2.1 0. 0 +GRID 955 0 3.5 2.1 0. 0 +GRID 956 0 3.6 2.1 0. 0 +GRID 957 0 3.7 2.1 0. 0 +GRID 958 0 3.8 2.1 0. 0 +GRID 959 0 3.9 2.1 0. 0 +GRID 960 0 .1 2.2 0. 0 +GRID 961 0 .2 2.2 0. 0 +GRID 962 0 .3 2.2 0. 0 +GRID 963 0 .4 2.2 0. 0 +GRID 964 0 .5 2.2 0. 0 +GRID 965 0 .6 2.2 0. 0 +GRID 966 0 .7 2.2 0. 0 +GRID 967 0 .8 2.2 0. 0 +GRID 968 0 .9 2.2 0. 0 +GRID 969 0 1. 2.2 0. 0 +GRID 970 0 1.1 2.2 0. 0 +GRID 971 0 1.2 2.2 0. 0 +GRID 972 0 1.3 2.2 0. 0 +GRID 973 0 1.4 2.2 0. 0 +GRID 974 0 1.5 2.2 0. 0 +GRID 975 0 1.6 2.2 0. 0 +GRID 976 0 1.7 2.2 0. 0 +GRID 977 0 1.8 2.2 0. 0 +GRID 978 0 1.9 2.2 0. 0 +GRID 979 0 2. 2.2 0. 0 +GRID 980 0 2.1 2.2 0. 0 +GRID 981 0 2.2 2.2 0. 0 +GRID 982 0 2.3 2.2 0. 0 +GRID 983 0 2.4 2.2 0. 0 +GRID 984 0 2.5 2.2 0. 0 +GRID 985 0 2.6 2.2 0. 0 +GRID 986 0 2.7 2.2 0. 0 +GRID 987 0 2.8 2.2 0. 0 +GRID 988 0 2.9 2.2 0. 0 +GRID 989 0 3. 2.2 0. 0 +GRID 990 0 3.1 2.2 0. 0 +GRID 991 0 3.2 2.2 0. 0 +GRID 992 0 3.3 2.2 0. 0 +GRID 993 0 3.4 2.2 0. 0 +GRID 994 0 3.5 2.2 0. 0 +GRID 995 0 3.6 2.2 0. 0 +GRID 996 0 3.7 2.2 0. 0 +GRID 997 0 3.8 2.2 0. 0 +GRID 998 0 3.9 2.2 0. 0 +GRID 999 0 .1 2.3 0. 0 +GRID 1000 0 .2 2.3 0. 0 +GRID 1001 0 .3 2.3 0. 0 +GRID 1002 0 .4 2.3 0. 0 +GRID 1003 0 .5 2.3 0. 0 +GRID 1004 0 .6 2.3 0. 0 +GRID 1005 0 .7 2.3 0. 0 +GRID 1006 0 .8 2.3 0. 0 +GRID 1007 0 .9 2.3 0. 0 +GRID 1008 0 1. 2.3 0. 0 +GRID 1009 0 1.1 2.3 0. 0 +GRID 1010 0 1.2 2.3 0. 0 +GRID 1011 0 1.3 2.3 0. 0 +GRID 1012 0 1.4 2.3 0. 0 +GRID 1013 0 1.5 2.3 0. 0 +GRID 1014 0 1.6 2.3 0. 0 +GRID 1015 0 1.7 2.3 0. 0 +GRID 1016 0 1.8 2.3 0. 0 +GRID 1017 0 1.9 2.3 0. 0 +GRID 1018 0 2. 2.3 0. 0 +GRID 1019 0 2.1 2.3 0. 0 +GRID 1020 0 2.2 2.3 0. 0 +GRID 1021 0 2.3 2.3 0. 0 +GRID 1022 0 2.4 2.3 0. 0 +GRID 1023 0 2.5 2.3 0. 0 +GRID 1024 0 2.6 2.3 0. 0 +GRID 1025 0 2.7 2.3 0. 0 +GRID 1026 0 2.8 2.3 0. 0 +GRID 1027 0 2.9 2.3 0. 0 +GRID 1028 0 3. 2.3 0. 0 +GRID 1029 0 3.1 2.3 0. 0 +GRID 1030 0 3.2 2.3 0. 0 +GRID 1031 0 3.3 2.3 0. 0 +GRID 1032 0 3.4 2.3 0. 0 +GRID 1033 0 3.5 2.3 0. 0 +GRID 1034 0 3.6 2.3 0. 0 +GRID 1035 0 3.7 2.3 0. 0 +GRID 1036 0 3.8 2.3 0. 0 +GRID 1037 0 3.9 2.3 0. 0 +GRID 1038 0 .1 2.4 0. 0 +GRID 1039 0 .2 2.4 0. 0 +GRID 1040 0 .3 2.4 0. 0 +GRID 1041 0 .4 2.4 0. 0 +GRID 1042 0 .5 2.4 0. 0 +GRID 1043 0 .6 2.4 0. 0 +GRID 1044 0 .7 2.4 0. 0 +GRID 1045 0 .8 2.4 0. 0 +GRID 1046 0 .9 2.4 0. 0 +GRID 1047 0 1. 2.4 0. 0 +GRID 1048 0 1.1 2.4 0. 0 +GRID 1049 0 1.2 2.4 0. 0 +GRID 1050 0 1.3 2.4 0. 0 +GRID 1051 0 1.4 2.4 0. 0 +GRID 1052 0 1.5 2.4 0. 0 +GRID 1053 0 1.6 2.4 0. 0 +GRID 1054 0 1.7 2.4 0. 0 +GRID 1055 0 1.8 2.4 0. 0 +GRID 1056 0 1.9 2.4 0. 0 +GRID 1057 0 2. 2.4 0. 0 +GRID 1058 0 2.1 2.4 0. 0 +GRID 1059 0 2.2 2.4 0. 0 +GRID 1060 0 2.3 2.4 0. 0 +GRID 1061 0 2.4 2.4 0. 0 +GRID 1062 0 2.5 2.4 0. 0 +GRID 1063 0 2.6 2.4 0. 0 +GRID 1064 0 2.7 2.4 0. 0 +GRID 1065 0 2.8 2.4 0. 0 +GRID 1066 0 2.9 2.4 0. 0 +GRID 1067 0 3. 2.4 0. 0 +GRID 1068 0 3.1 2.4 0. 0 +GRID 1069 0 3.2 2.4 0. 0 +GRID 1070 0 3.3 2.4 0. 0 +GRID 1071 0 3.4 2.4 0. 0 +GRID 1072 0 3.5 2.4 0. 0 +GRID 1073 0 3.6 2.4 0. 0 +GRID 1074 0 3.7 2.4 0. 0 +GRID 1075 0 3.8 2.4 0. 0 +GRID 1076 0 3.9 2.4 0. 0 +GRID 1077 0 .1 2.5 0. 0 +GRID 1078 0 .2 2.5 0. 0 +GRID 1079 0 .3 2.5 0. 0 +GRID 1080 0 .4 2.5 0. 0 +GRID 1081 0 .5 2.5 0. 0 +GRID 1082 0 .6 2.5 0. 0 +GRID 1083 0 .7 2.5 0. 0 +GRID 1084 0 .8 2.5 0. 0 +GRID 1085 0 .9 2.5 0. 0 +GRID 1086 0 1. 2.5 0. 0 +GRID 1087 0 1.1 2.5 0. 0 +GRID 1088 0 1.2 2.5 0. 0 +GRID 1089 0 1.3 2.5 0. 0 +GRID 1090 0 1.4 2.5 0. 0 +GRID 1091 0 1.5 2.5 0. 0 +GRID 1092 0 1.6 2.5 0. 0 +GRID 1093 0 1.7 2.5 0. 0 +GRID 1094 0 1.8 2.5 0. 0 +GRID 1095 0 1.9 2.5 0. 0 +GRID 1096 0 2. 2.5 0. 0 +GRID 1097 0 2.1 2.5 0. 0 +GRID 1098 0 2.2 2.5 0. 0 +GRID 1099 0 2.3 2.5 0. 0 +GRID 1100 0 2.4 2.5 0. 0 +GRID 1101 0 2.5 2.5 0. 0 +GRID 1102 0 2.6 2.5 0. 0 +GRID 1103 0 2.7 2.5 0. 0 +GRID 1104 0 2.8 2.5 0. 0 +GRID 1105 0 2.9 2.5 0. 0 +GRID 1106 0 3. 2.5 0. 0 +GRID 1107 0 3.1 2.5 0. 0 +GRID 1108 0 3.2 2.5 0. 0 +GRID 1109 0 3.3 2.5 0. 0 +GRID 1110 0 3.4 2.5 0. 0 +GRID 1111 0 3.5 2.5 0. 0 +GRID 1112 0 3.6 2.5 0. 0 +GRID 1113 0 3.7 2.5 0. 0 +GRID 1114 0 3.8 2.5 0. 0 +GRID 1115 0 3.9 2.5 0. 0 +GRID 1116 0 .1 2.6 0. 0 +GRID 1117 0 .2 2.6 0. 0 +GRID 1118 0 .3 2.6 0. 0 +GRID 1119 0 .4 2.6 0. 0 +GRID 1120 0 .5 2.6 0. 0 +GRID 1121 0 .6 2.6 0. 0 +GRID 1122 0 .7 2.6 0. 0 +GRID 1123 0 .8 2.6 0. 0 +GRID 1124 0 .9 2.6 0. 0 +GRID 1125 0 1. 2.6 0. 0 +GRID 1126 0 1.1 2.6 0. 0 +GRID 1127 0 1.2 2.6 0. 0 +GRID 1128 0 1.3 2.6 0. 0 +GRID 1129 0 1.4 2.6 0. 0 +GRID 1130 0 1.5 2.6 0. 0 +GRID 1131 0 1.6 2.6 0. 0 +GRID 1132 0 1.7 2.6 0. 0 +GRID 1133 0 1.8 2.6 0. 0 +GRID 1134 0 1.9 2.6 0. 0 +GRID 1135 0 2. 2.6 0. 0 +GRID 1136 0 2.1 2.6 0. 0 +GRID 1137 0 2.2 2.6 0. 0 +GRID 1138 0 2.3 2.6 0. 0 +GRID 1139 0 2.4 2.6 0. 0 +GRID 1140 0 2.5 2.6 0. 0 +GRID 1141 0 2.6 2.6 0. 0 +GRID 1142 0 2.7 2.6 0. 0 +GRID 1143 0 2.8 2.6 0. 0 +GRID 1144 0 2.9 2.6 0. 0 +GRID 1145 0 3. 2.6 0. 0 +GRID 1146 0 3.1 2.6 0. 0 +GRID 1147 0 3.2 2.6 0. 0 +GRID 1148 0 3.3 2.6 0. 0 +GRID 1149 0 3.4 2.6 0. 0 +GRID 1150 0 3.5 2.6 0. 0 +GRID 1151 0 3.6 2.6 0. 0 +GRID 1152 0 3.7 2.6 0. 0 +GRID 1153 0 3.8 2.6 0. 0 +GRID 1154 0 3.9 2.6 0. 0 +GRID 1155 0 .1 2.7 0. 0 +GRID 1156 0 .2 2.7 0. 0 +GRID 1157 0 .3 2.7 0. 0 +GRID 1158 0 .4 2.7 0. 0 +GRID 1159 0 .5 2.7 0. 0 +GRID 1160 0 .6 2.7 0. 0 +GRID 1161 0 .7 2.7 0. 0 +GRID 1162 0 .8 2.7 0. 0 +GRID 1163 0 .9 2.7 0. 0 +GRID 1164 0 1. 2.7 0. 0 +GRID 1165 0 1.1 2.7 0. 0 +GRID 1166 0 1.2 2.7 0. 0 +GRID 1167 0 1.3 2.7 0. 0 +GRID 1168 0 1.4 2.7 0. 0 +GRID 1169 0 1.5 2.7 0. 0 +GRID 1170 0 1.6 2.7 0. 0 +GRID 1171 0 1.7 2.7 0. 0 +GRID 1172 0 1.8 2.7 0. 0 +GRID 1173 0 1.9 2.7 0. 0 +GRID 1174 0 2. 2.7 0. 0 +GRID 1175 0 2.1 2.7 0. 0 +GRID 1176 0 2.2 2.7 0. 0 +GRID 1177 0 2.3 2.7 0. 0 +GRID 1178 0 2.4 2.7 0. 0 +GRID 1179 0 2.5 2.7 0. 0 +GRID 1180 0 2.6 2.7 0. 0 +GRID 1181 0 2.7 2.7 0. 0 +GRID 1182 0 2.8 2.7 0. 0 +GRID 1183 0 2.9 2.7 0. 0 +GRID 1184 0 3. 2.7 0. 0 +GRID 1185 0 3.1 2.7 0. 0 +GRID 1186 0 3.2 2.7 0. 0 +GRID 1187 0 3.3 2.7 0. 0 +GRID 1188 0 3.4 2.7 0. 0 +GRID 1189 0 3.5 2.7 0. 0 +GRID 1190 0 3.6 2.7 0. 0 +GRID 1191 0 3.7 2.7 0. 0 +GRID 1192 0 3.8 2.7 0. 0 +GRID 1193 0 3.9 2.7 0. 0 +GRID 1194 0 .1 2.8 0. 0 +GRID 1195 0 .2 2.8 0. 0 +GRID 1196 0 .3 2.8 0. 0 +GRID 1197 0 .4 2.8 0. 0 +GRID 1198 0 .5 2.8 0. 0 +GRID 1199 0 .6 2.8 0. 0 +GRID 1200 0 .7 2.8 0. 0 +GRID 1201 0 .8 2.8 0. 0 +GRID 1202 0 .9 2.8 0. 0 +GRID 1203 0 1. 2.8 0. 0 +GRID 1204 0 1.1 2.8 0. 0 +GRID 1205 0 1.2 2.8 0. 0 +GRID 1206 0 1.3 2.8 0. 0 +GRID 1207 0 1.4 2.8 0. 0 +GRID 1208 0 1.5 2.8 0. 0 +GRID 1209 0 1.6 2.8 0. 0 +GRID 1210 0 1.7 2.8 0. 0 +GRID 1211 0 1.8 2.8 0. 0 +GRID 1212 0 1.9 2.8 0. 0 +GRID 1213 0 2. 2.8 0. 0 +GRID 1214 0 2.1 2.8 0. 0 +GRID 1215 0 2.2 2.8 0. 0 +GRID 1216 0 2.3 2.8 0. 0 +GRID 1217 0 2.4 2.8 0. 0 +GRID 1218 0 2.5 2.8 0. 0 +GRID 1219 0 2.6 2.8 0. 0 +GRID 1220 0 2.7 2.8 0. 0 +GRID 1221 0 2.8 2.8 0. 0 +GRID 1222 0 2.9 2.8 0. 0 +GRID 1223 0 3. 2.8 0. 0 +GRID 1224 0 3.1 2.8 0. 0 +GRID 1225 0 3.2 2.8 0. 0 +GRID 1226 0 3.3 2.8 0. 0 +GRID 1227 0 3.4 2.8 0. 0 +GRID 1228 0 3.5 2.8 0. 0 +GRID 1229 0 3.6 2.8 0. 0 +GRID 1230 0 3.7 2.8 0. 0 +GRID 1231 0 3.8 2.8 0. 0 +GRID 1232 0 3.9 2.8 0. 0 +GRID 1233 0 .1 2.9 0. 0 +GRID 1234 0 .2 2.9 0. 0 +GRID 1235 0 .3 2.9 0. 0 +GRID 1236 0 .4 2.9 0. 0 +GRID 1237 0 .5 2.9 0. 0 +GRID 1238 0 .6 2.9 0. 0 +GRID 1239 0 .7 2.9 0. 0 +GRID 1240 0 .8 2.9 0. 0 +GRID 1241 0 .9 2.9 0. 0 +GRID 1242 0 1. 2.9 0. 0 +GRID 1243 0 1.1 2.9 0. 0 +GRID 1244 0 1.2 2.9 0. 0 +GRID 1245 0 1.3 2.9 0. 0 +GRID 1246 0 1.4 2.9 0. 0 +GRID 1247 0 1.5 2.9 0. 0 +GRID 1248 0 1.6 2.9 0. 0 +GRID 1249 0 1.7 2.9 0. 0 +GRID 1250 0 1.8 2.9 0. 0 +GRID 1251 0 1.9 2.9 0. 0 +GRID 1252 0 2. 2.9 0. 0 +GRID 1253 0 2.1 2.9 0. 0 +GRID 1254 0 2.2 2.9 0. 0 +GRID 1255 0 2.3 2.9 0. 0 +GRID 1256 0 2.4 2.9 0. 0 +GRID 1257 0 2.5 2.9 0. 0 +GRID 1258 0 2.6 2.9 0. 0 +GRID 1259 0 2.7 2.9 0. 0 +GRID 1260 0 2.8 2.9 0. 0 +GRID 1261 0 2.9 2.9 0. 0 +GRID 1262 0 3. 2.9 0. 0 +GRID 1263 0 3.1 2.9 0. 0 +GRID 1264 0 3.2 2.9 0. 0 +GRID 1265 0 3.3 2.9 0. 0 +GRID 1266 0 3.4 2.9 0. 0 +GRID 1267 0 3.5 2.9 0. 0 +GRID 1268 0 3.6 2.9 0. 0 +GRID 1269 0 3.7 2.9 0. 0 +GRID 1270 0 3.8 2.9 0. 0 +GRID 1271 0 3.9 2.9 0. 0 +CQUAD4 1 1 1 2 141 140 +CQUAD4 2 1 2 3 142 141 +CQUAD4 3 1 3 4 143 142 +CQUAD4 4 1 4 5 144 143 +CQUAD4 5 1 5 6 145 144 +CQUAD4 6 1 6 7 146 145 +CQUAD4 7 1 7 8 147 146 +CQUAD4 8 1 8 9 148 147 +CQUAD4 9 1 9 10 149 148 +CQUAD4 10 1 10 11 150 149 +CQUAD4 11 1 11 12 151 150 +CQUAD4 12 1 12 13 152 151 +CQUAD4 13 1 13 14 153 152 +CQUAD4 14 1 14 15 154 153 +CQUAD4 15 1 15 16 155 154 +CQUAD4 16 1 16 17 156 155 +CQUAD4 17 1 17 18 157 156 +CQUAD4 18 1 18 19 158 157 +CQUAD4 19 1 19 20 159 158 +CQUAD4 20 1 20 21 160 159 +CQUAD4 21 1 21 22 161 160 +CQUAD4 22 1 22 23 162 161 +CQUAD4 23 1 23 24 163 162 +CQUAD4 24 1 24 25 164 163 +CQUAD4 25 1 25 26 165 164 +CQUAD4 26 1 26 27 166 165 +CQUAD4 27 1 27 28 167 166 +CQUAD4 28 1 28 29 168 167 +CQUAD4 29 1 29 30 169 168 +CQUAD4 30 1 30 31 170 169 +CQUAD4 31 1 31 32 171 170 +CQUAD4 32 1 32 33 172 171 +CQUAD4 33 1 33 34 173 172 +CQUAD4 34 1 34 35 174 173 +CQUAD4 35 1 35 36 175 174 +CQUAD4 36 1 36 37 176 175 +CQUAD4 37 1 37 38 177 176 +CQUAD4 38 1 38 39 178 177 +CQUAD4 39 1 39 40 179 178 +CQUAD4 40 1 40 41 42 179 +CQUAD4 41 1 140 141 180 139 +CQUAD4 42 1 141 142 181 180 +CQUAD4 43 1 142 143 182 181 +CQUAD4 44 1 143 144 183 182 +CQUAD4 45 1 144 145 184 183 +CQUAD4 46 1 145 146 185 184 +CQUAD4 47 1 146 147 186 185 +CQUAD4 48 1 147 148 187 186 +CQUAD4 49 1 148 149 188 187 +CQUAD4 50 1 149 150 189 188 +CQUAD4 51 1 150 151 190 189 +CQUAD4 52 1 151 152 191 190 +CQUAD4 53 1 152 153 192 191 +CQUAD4 54 1 153 154 193 192 +CQUAD4 55 1 154 155 194 193 +CQUAD4 56 1 155 156 195 194 +CQUAD4 57 1 156 157 196 195 +CQUAD4 58 1 157 158 197 196 +CQUAD4 59 1 158 159 198 197 +CQUAD4 60 1 159 160 199 198 +CQUAD4 61 1 160 161 200 199 +CQUAD4 62 1 161 162 201 200 +CQUAD4 63 1 162 163 202 201 +CQUAD4 64 1 163 164 203 202 +CQUAD4 65 1 164 165 204 203 +CQUAD4 66 1 165 166 205 204 +CQUAD4 67 1 166 167 206 205 +CQUAD4 68 1 167 168 207 206 +CQUAD4 69 1 168 169 208 207 +CQUAD4 70 1 169 170 209 208 +CQUAD4 71 1 170 171 210 209 +CQUAD4 72 1 171 172 211 210 +CQUAD4 73 1 172 173 212 211 +CQUAD4 74 1 173 174 213 212 +CQUAD4 75 1 174 175 214 213 +CQUAD4 76 1 175 176 215 214 +CQUAD4 77 1 176 177 216 215 +CQUAD4 78 1 177 178 217 216 +CQUAD4 79 1 178 179 218 217 +CQUAD4 80 1 179 42 43 218 +CQUAD4 81 1 139 180 219 138 +CQUAD4 82 1 180 181 220 219 +CQUAD4 83 1 181 182 221 220 +CQUAD4 84 1 182 183 222 221 +CQUAD4 85 1 183 184 223 222 +CQUAD4 86 1 184 185 224 223 +CQUAD4 87 1 185 186 225 224 +CQUAD4 88 1 186 187 226 225 +CQUAD4 89 1 187 188 227 226 +CQUAD4 90 1 188 189 228 227 +CQUAD4 91 1 189 190 229 228 +CQUAD4 92 1 190 191 230 229 +CQUAD4 93 1 191 192 231 230 +CQUAD4 94 1 192 193 232 231 +CQUAD4 95 1 193 194 233 232 +CQUAD4 96 1 194 195 234 233 +CQUAD4 97 1 195 196 235 234 +CQUAD4 98 1 196 197 236 235 +CQUAD4 99 1 197 198 237 236 +CQUAD4 100 1 198 199 238 237 +CQUAD4 101 1 199 200 239 238 +CQUAD4 102 1 200 201 240 239 +CQUAD4 103 1 201 202 241 240 +CQUAD4 104 1 202 203 242 241 +CQUAD4 105 1 203 204 243 242 +CQUAD4 106 1 204 205 244 243 +CQUAD4 107 1 205 206 245 244 +CQUAD4 108 1 206 207 246 245 +CQUAD4 109 1 207 208 247 246 +CQUAD4 110 1 208 209 248 247 +CQUAD4 111 1 209 210 249 248 +CQUAD4 112 1 210 211 250 249 +CQUAD4 113 1 211 212 251 250 +CQUAD4 114 1 212 213 252 251 +CQUAD4 115 1 213 214 253 252 +CQUAD4 116 1 214 215 254 253 +CQUAD4 117 1 215 216 255 254 +CQUAD4 118 1 216 217 256 255 +CQUAD4 119 1 217 218 257 256 +CQUAD4 120 1 218 43 44 257 +CQUAD4 121 1 138 219 258 137 +CQUAD4 122 1 219 220 259 258 +CQUAD4 123 1 220 221 260 259 +CQUAD4 124 1 221 222 261 260 +CQUAD4 125 1 222 223 262 261 +CQUAD4 126 1 223 224 263 262 +CQUAD4 127 1 224 225 264 263 +CQUAD4 128 1 225 226 265 264 +CQUAD4 129 1 226 227 266 265 +CQUAD4 130 1 227 228 267 266 +CQUAD4 131 1 228 229 268 267 +CQUAD4 132 1 229 230 269 268 +CQUAD4 133 1 230 231 270 269 +CQUAD4 134 1 231 232 271 270 +CQUAD4 135 1 232 233 272 271 +CQUAD4 136 1 233 234 273 272 +CQUAD4 137 1 234 235 274 273 +CQUAD4 138 1 235 236 275 274 +CQUAD4 139 1 236 237 276 275 +CQUAD4 140 1 237 238 277 276 +CQUAD4 141 1 238 239 278 277 +CQUAD4 142 1 239 240 279 278 +CQUAD4 143 1 240 241 280 279 +CQUAD4 144 1 241 242 281 280 +CQUAD4 145 1 242 243 282 281 +CQUAD4 146 1 243 244 283 282 +CQUAD4 147 1 244 245 284 283 +CQUAD4 148 1 245 246 285 284 +CQUAD4 149 1 246 247 286 285 +CQUAD4 150 1 247 248 287 286 +CQUAD4 151 1 248 249 288 287 +CQUAD4 152 1 249 250 289 288 +CQUAD4 153 1 250 251 290 289 +CQUAD4 154 1 251 252 291 290 +CQUAD4 155 1 252 253 292 291 +CQUAD4 156 1 253 254 293 292 +CQUAD4 157 1 254 255 294 293 +CQUAD4 158 1 255 256 295 294 +CQUAD4 159 1 256 257 296 295 +CQUAD4 160 1 257 44 45 296 +CQUAD4 161 1 137 258 297 136 +CQUAD4 162 1 258 259 298 297 +CQUAD4 163 1 259 260 299 298 +CQUAD4 164 1 260 261 300 299 +CQUAD4 165 1 261 262 301 300 +CQUAD4 166 1 262 263 302 301 +CQUAD4 167 1 263 264 303 302 +CQUAD4 168 1 264 265 304 303 +CQUAD4 169 1 265 266 305 304 +CQUAD4 170 1 266 267 306 305 +CQUAD4 171 1 267 268 307 306 +CQUAD4 172 1 268 269 308 307 +CQUAD4 173 1 269 270 309 308 +CQUAD4 174 1 270 271 310 309 +CQUAD4 175 1 271 272 311 310 +CQUAD4 176 1 272 273 312 311 +CQUAD4 177 1 273 274 313 312 +CQUAD4 178 1 274 275 314 313 +CQUAD4 179 1 275 276 315 314 +CQUAD4 180 1 276 277 316 315 +CQUAD4 181 1 277 278 317 316 +CQUAD4 182 1 278 279 318 317 +CQUAD4 183 1 279 280 319 318 +CQUAD4 184 1 280 281 320 319 +CQUAD4 185 1 281 282 321 320 +CQUAD4 186 1 282 283 322 321 +CQUAD4 187 1 283 284 323 322 +CQUAD4 188 1 284 285 324 323 +CQUAD4 189 1 285 286 325 324 +CQUAD4 190 1 286 287 326 325 +CQUAD4 191 1 287 288 327 326 +CQUAD4 192 1 288 289 328 327 +CQUAD4 193 1 289 290 329 328 +CQUAD4 194 1 290 291 330 329 +CQUAD4 195 1 291 292 331 330 +CQUAD4 196 1 292 293 332 331 +CQUAD4 197 1 293 294 333 332 +CQUAD4 198 1 294 295 334 333 +CQUAD4 199 1 295 296 335 334 +CQUAD4 200 1 296 45 46 335 +CQUAD4 201 1 136 297 336 135 +CQUAD4 202 1 297 298 337 336 +CQUAD4 203 1 298 299 338 337 +CQUAD4 204 1 299 300 339 338 +CQUAD4 205 1 300 301 340 339 +CQUAD4 206 1 301 302 341 340 +CQUAD4 207 1 302 303 342 341 +CQUAD4 208 1 303 304 343 342 +CQUAD4 209 1 304 305 344 343 +CQUAD4 210 1 305 306 345 344 +CQUAD4 211 1 306 307 346 345 +CQUAD4 212 1 307 308 347 346 +CQUAD4 213 1 308 309 348 347 +CQUAD4 214 1 309 310 349 348 +CQUAD4 215 1 310 311 350 349 +CQUAD4 216 1 311 312 351 350 +CQUAD4 217 1 312 313 352 351 +CQUAD4 218 1 313 314 353 352 +CQUAD4 219 1 314 315 354 353 +CQUAD4 220 1 315 316 355 354 +CQUAD4 221 1 316 317 356 355 +CQUAD4 222 1 317 318 357 356 +CQUAD4 223 1 318 319 358 357 +CQUAD4 224 1 319 320 359 358 +CQUAD4 225 1 320 321 360 359 +CQUAD4 226 1 321 322 361 360 +CQUAD4 227 1 322 323 362 361 +CQUAD4 228 1 323 324 363 362 +CQUAD4 229 1 324 325 364 363 +CQUAD4 230 1 325 326 365 364 +CQUAD4 231 1 326 327 366 365 +CQUAD4 232 1 327 328 367 366 +CQUAD4 233 1 328 329 368 367 +CQUAD4 234 1 329 330 369 368 +CQUAD4 235 1 330 331 370 369 +CQUAD4 236 1 331 332 371 370 +CQUAD4 237 1 332 333 372 371 +CQUAD4 238 1 333 334 373 372 +CQUAD4 239 1 334 335 374 373 +CQUAD4 240 1 335 46 47 374 +CQUAD4 241 1 135 336 375 134 +CQUAD4 242 1 336 337 376 375 +CQUAD4 243 1 337 338 377 376 +CQUAD4 244 1 338 339 378 377 +CQUAD4 245 1 339 340 379 378 +CQUAD4 246 1 340 341 380 379 +CQUAD4 247 1 341 342 381 380 +CQUAD4 248 1 342 343 382 381 +CQUAD4 249 1 343 344 383 382 +CQUAD4 250 1 344 345 384 383 +CQUAD4 251 1 345 346 385 384 +CQUAD4 252 1 346 347 386 385 +CQUAD4 253 1 347 348 387 386 +CQUAD4 254 1 348 349 388 387 +CQUAD4 255 1 349 350 389 388 +CQUAD4 256 1 350 351 390 389 +CQUAD4 257 1 351 352 391 390 +CQUAD4 258 1 352 353 392 391 +CQUAD4 259 1 353 354 393 392 +CQUAD4 260 1 354 355 394 393 +CQUAD4 261 1 355 356 395 394 +CQUAD4 262 1 356 357 396 395 +CQUAD4 263 1 357 358 397 396 +CQUAD4 264 1 358 359 398 397 +CQUAD4 265 1 359 360 399 398 +CQUAD4 266 1 360 361 400 399 +CQUAD4 267 1 361 362 401 400 +CQUAD4 268 1 362 363 402 401 +CQUAD4 269 1 363 364 403 402 +CQUAD4 270 1 364 365 404 403 +CQUAD4 271 1 365 366 405 404 +CQUAD4 272 1 366 367 406 405 +CQUAD4 273 1 367 368 407 406 +CQUAD4 274 1 368 369 408 407 +CQUAD4 275 1 369 370 409 408 +CQUAD4 276 1 370 371 410 409 +CQUAD4 277 1 371 372 411 410 +CQUAD4 278 1 372 373 412 411 +CQUAD4 279 1 373 374 413 412 +CQUAD4 280 1 374 47 48 413 +CQUAD4 281 1 134 375 414 133 +CQUAD4 282 1 375 376 415 414 +CQUAD4 283 1 376 377 416 415 +CQUAD4 284 1 377 378 417 416 +CQUAD4 285 1 378 379 418 417 +CQUAD4 286 1 379 380 419 418 +CQUAD4 287 1 380 381 420 419 +CQUAD4 288 1 381 382 421 420 +CQUAD4 289 1 382 383 422 421 +CQUAD4 290 1 383 384 423 422 +CQUAD4 291 1 384 385 424 423 +CQUAD4 292 1 385 386 425 424 +CQUAD4 293 1 386 387 426 425 +CQUAD4 294 1 387 388 427 426 +CQUAD4 295 1 388 389 428 427 +CQUAD4 296 1 389 390 429 428 +CQUAD4 297 1 390 391 430 429 +CQUAD4 298 1 391 392 431 430 +CQUAD4 299 1 392 393 432 431 +CQUAD4 300 1 393 394 433 432 +CQUAD4 301 1 394 395 434 433 +CQUAD4 302 1 395 396 435 434 +CQUAD4 303 1 396 397 436 435 +CQUAD4 304 1 397 398 437 436 +CQUAD4 305 1 398 399 438 437 +CQUAD4 306 1 399 400 439 438 +CQUAD4 307 1 400 401 440 439 +CQUAD4 308 1 401 402 441 440 +CQUAD4 309 1 402 403 442 441 +CQUAD4 310 1 403 404 443 442 +CQUAD4 311 1 404 405 444 443 +CQUAD4 312 1 405 406 445 444 +CQUAD4 313 1 406 407 446 445 +CQUAD4 314 1 407 408 447 446 +CQUAD4 315 1 408 409 448 447 +CQUAD4 316 1 409 410 449 448 +CQUAD4 317 1 410 411 450 449 +CQUAD4 318 1 411 412 451 450 +CQUAD4 319 1 412 413 452 451 +CQUAD4 320 1 413 48 49 452 +CQUAD4 321 1 133 414 453 132 +CQUAD4 322 1 414 415 454 453 +CQUAD4 323 1 415 416 455 454 +CQUAD4 324 1 416 417 456 455 +CQUAD4 325 1 417 418 457 456 +CQUAD4 326 1 418 419 458 457 +CQUAD4 327 1 419 420 459 458 +CQUAD4 328 1 420 421 460 459 +CQUAD4 329 1 421 422 461 460 +CQUAD4 330 1 422 423 462 461 +CQUAD4 331 1 423 424 463 462 +CQUAD4 332 1 424 425 464 463 +CQUAD4 333 1 425 426 465 464 +CQUAD4 334 1 426 427 466 465 +CQUAD4 335 1 427 428 467 466 +CQUAD4 336 1 428 429 468 467 +CQUAD4 337 1 429 430 469 468 +CQUAD4 338 1 430 431 470 469 +CQUAD4 339 1 431 432 471 470 +CQUAD4 340 1 432 433 472 471 +CQUAD4 341 1 433 434 473 472 +CQUAD4 342 1 434 435 474 473 +CQUAD4 343 1 435 436 475 474 +CQUAD4 344 1 436 437 476 475 +CQUAD4 345 1 437 438 477 476 +CQUAD4 346 1 438 439 478 477 +CQUAD4 347 1 439 440 479 478 +CQUAD4 348 1 440 441 480 479 +CQUAD4 349 1 441 442 481 480 +CQUAD4 350 1 442 443 482 481 +CQUAD4 351 1 443 444 483 482 +CQUAD4 352 1 444 445 484 483 +CQUAD4 353 1 445 446 485 484 +CQUAD4 354 1 446 447 486 485 +CQUAD4 355 1 447 448 487 486 +CQUAD4 356 1 448 449 488 487 +CQUAD4 357 1 449 450 489 488 +CQUAD4 358 1 450 451 490 489 +CQUAD4 359 1 451 452 491 490 +CQUAD4 360 1 452 49 50 491 +CQUAD4 361 1 132 453 492 131 +CQUAD4 362 1 453 454 493 492 +CQUAD4 363 1 454 455 494 493 +CQUAD4 364 1 455 456 495 494 +CQUAD4 365 1 456 457 496 495 +CQUAD4 366 1 457 458 497 496 +CQUAD4 367 1 458 459 498 497 +CQUAD4 368 1 459 460 499 498 +CQUAD4 369 1 460 461 500 499 +CQUAD4 370 1 461 462 501 500 +CQUAD4 371 1 462 463 502 501 +CQUAD4 372 1 463 464 503 502 +CQUAD4 373 1 464 465 504 503 +CQUAD4 374 1 465 466 505 504 +CQUAD4 375 1 466 467 506 505 +CQUAD4 376 1 467 468 507 506 +CQUAD4 377 1 468 469 508 507 +CQUAD4 378 1 469 470 509 508 +CQUAD4 379 1 470 471 510 509 +CQUAD4 380 1 471 472 511 510 +CQUAD4 381 1 472 473 512 511 +CQUAD4 382 1 473 474 513 512 +CQUAD4 383 1 474 475 514 513 +CQUAD4 384 1 475 476 515 514 +CQUAD4 385 1 476 477 516 515 +CQUAD4 386 1 477 478 517 516 +CQUAD4 387 1 478 479 518 517 +CQUAD4 388 1 479 480 519 518 +CQUAD4 389 1 480 481 520 519 +CQUAD4 390 1 481 482 521 520 +CQUAD4 391 1 482 483 522 521 +CQUAD4 392 1 483 484 523 522 +CQUAD4 393 1 484 485 524 523 +CQUAD4 394 1 485 486 525 524 +CQUAD4 395 1 486 487 526 525 +CQUAD4 396 1 487 488 527 526 +CQUAD4 397 1 488 489 528 527 +CQUAD4 398 1 489 490 529 528 +CQUAD4 399 1 490 491 530 529 +CQUAD4 400 1 491 50 51 530 +CQUAD4 401 1 131 492 531 130 +CQUAD4 402 1 492 493 532 531 +CQUAD4 403 1 493 494 533 532 +CQUAD4 404 1 494 495 534 533 +CQUAD4 405 1 495 496 535 534 +CQUAD4 406 1 496 497 536 535 +CQUAD4 407 1 497 498 537 536 +CQUAD4 408 1 498 499 538 537 +CQUAD4 409 1 499 500 539 538 +CQUAD4 410 1 500 501 540 539 +CQUAD4 411 1 501 502 541 540 +CQUAD4 412 1 502 503 542 541 +CQUAD4 413 1 503 504 543 542 +CQUAD4 414 1 504 505 544 543 +CQUAD4 415 1 505 506 545 544 +CQUAD4 416 1 506 507 546 545 +CQUAD4 417 1 507 508 547 546 +CQUAD4 418 1 508 509 548 547 +CQUAD4 419 1 509 510 549 548 +CQUAD4 420 1 510 511 550 549 +CQUAD4 421 1 511 512 551 550 +CQUAD4 422 1 512 513 552 551 +CQUAD4 423 1 513 514 553 552 +CQUAD4 424 1 514 515 554 553 +CQUAD4 425 1 515 516 555 554 +CQUAD4 426 1 516 517 556 555 +CQUAD4 427 1 517 518 557 556 +CQUAD4 428 1 518 519 558 557 +CQUAD4 429 1 519 520 559 558 +CQUAD4 430 1 520 521 560 559 +CQUAD4 431 1 521 522 561 560 +CQUAD4 432 1 522 523 562 561 +CQUAD4 433 1 523 524 563 562 +CQUAD4 434 1 524 525 564 563 +CQUAD4 435 1 525 526 565 564 +CQUAD4 436 1 526 527 566 565 +CQUAD4 437 1 527 528 567 566 +CQUAD4 438 1 528 529 568 567 +CQUAD4 439 1 529 530 569 568 +CQUAD4 440 1 530 51 52 569 +CQUAD4 441 1 130 531 570 129 +CQUAD4 442 1 531 532 571 570 +CQUAD4 443 1 532 533 572 571 +CQUAD4 444 1 533 534 573 572 +CQUAD4 445 1 534 535 574 573 +CQUAD4 446 1 535 536 575 574 +CQUAD4 447 1 536 537 576 575 +CQUAD4 448 1 537 538 577 576 +CQUAD4 449 1 538 539 578 577 +CQUAD4 450 1 539 540 579 578 +CQUAD4 451 1 540 541 580 579 +CQUAD4 452 1 541 542 581 580 +CQUAD4 453 1 542 543 582 581 +CQUAD4 454 1 543 544 583 582 +CQUAD4 455 1 544 545 584 583 +CQUAD4 456 1 545 546 585 584 +CQUAD4 457 1 546 547 586 585 +CQUAD4 458 1 547 548 587 586 +CQUAD4 459 1 548 549 588 587 +CQUAD4 460 1 549 550 589 588 +CQUAD4 461 1 550 551 590 589 +CQUAD4 462 1 551 552 591 590 +CQUAD4 463 1 552 553 592 591 +CQUAD4 464 1 553 554 593 592 +CQUAD4 465 1 554 555 594 593 +CQUAD4 466 1 555 556 595 594 +CQUAD4 467 1 556 557 596 595 +CQUAD4 468 1 557 558 597 596 +CQUAD4 469 1 558 559 598 597 +CQUAD4 470 1 559 560 599 598 +CQUAD4 471 1 560 561 600 599 +CQUAD4 472 1 561 562 601 600 +CQUAD4 473 1 562 563 602 601 +CQUAD4 474 1 563 564 603 602 +CQUAD4 475 1 564 565 604 603 +CQUAD4 476 1 565 566 605 604 +CQUAD4 477 1 566 567 606 605 +CQUAD4 478 1 567 568 607 606 +CQUAD4 479 1 568 569 608 607 +CQUAD4 480 1 569 52 53 608 +CQUAD4 481 1 129 570 609 128 +CQUAD4 482 1 570 571 610 609 +CQUAD4 483 1 571 572 611 610 +CQUAD4 484 1 572 573 612 611 +CQUAD4 485 1 573 574 613 612 +CQUAD4 486 1 574 575 614 613 +CQUAD4 487 1 575 576 615 614 +CQUAD4 488 1 576 577 616 615 +CQUAD4 489 1 577 578 617 616 +CQUAD4 490 1 578 579 618 617 +CQUAD4 491 1 579 580 619 618 +CQUAD4 492 1 580 581 620 619 +CQUAD4 493 1 581 582 621 620 +CQUAD4 494 1 582 583 622 621 +CQUAD4 495 1 583 584 623 622 +CQUAD4 496 1 584 585 624 623 +CQUAD4 497 1 585 586 625 624 +CQUAD4 498 1 586 587 626 625 +CQUAD4 499 1 587 588 627 626 +CQUAD4 500 1 588 589 628 627 +CQUAD4 501 1 589 590 629 628 +CQUAD4 502 1 590 591 630 629 +CQUAD4 503 1 591 592 631 630 +CQUAD4 504 1 592 593 632 631 +CQUAD4 505 1 593 594 633 632 +CQUAD4 506 1 594 595 634 633 +CQUAD4 507 1 595 596 635 634 +CQUAD4 508 1 596 597 636 635 +CQUAD4 509 1 597 598 637 636 +CQUAD4 510 1 598 599 638 637 +CQUAD4 511 1 599 600 639 638 +CQUAD4 512 1 600 601 640 639 +CQUAD4 513 1 601 602 641 640 +CQUAD4 514 1 602 603 642 641 +CQUAD4 515 1 603 604 643 642 +CQUAD4 516 1 604 605 644 643 +CQUAD4 517 1 605 606 645 644 +CQUAD4 518 1 606 607 646 645 +CQUAD4 519 1 607 608 647 646 +CQUAD4 520 1 608 53 54 647 +CQUAD4 521 1 128 609 648 127 +CQUAD4 522 1 609 610 649 648 +CQUAD4 523 1 610 611 650 649 +CQUAD4 524 1 611 612 651 650 +CQUAD4 525 1 612 613 652 651 +CQUAD4 526 1 613 614 653 652 +CQUAD4 527 1 614 615 654 653 +CQUAD4 528 1 615 616 655 654 +CQUAD4 529 1 616 617 656 655 +CQUAD4 530 1 617 618 657 656 +CQUAD4 531 1 618 619 658 657 +CQUAD4 532 1 619 620 659 658 +CQUAD4 533 1 620 621 660 659 +CQUAD4 534 1 621 622 661 660 +CQUAD4 535 1 622 623 662 661 +CQUAD4 536 1 623 624 663 662 +CQUAD4 537 1 624 625 664 663 +CQUAD4 538 1 625 626 665 664 +CQUAD4 539 1 626 627 666 665 +CQUAD4 540 1 627 628 667 666 +CQUAD4 541 1 628 629 668 667 +CQUAD4 542 1 629 630 669 668 +CQUAD4 543 1 630 631 670 669 +CQUAD4 544 1 631 632 671 670 +CQUAD4 545 1 632 633 672 671 +CQUAD4 546 1 633 634 673 672 +CQUAD4 547 1 634 635 674 673 +CQUAD4 548 1 635 636 675 674 +CQUAD4 549 1 636 637 676 675 +CQUAD4 550 1 637 638 677 676 +CQUAD4 551 1 638 639 678 677 +CQUAD4 552 1 639 640 679 678 +CQUAD4 553 1 640 641 680 679 +CQUAD4 554 1 641 642 681 680 +CQUAD4 555 1 642 643 682 681 +CQUAD4 556 1 643 644 683 682 +CQUAD4 557 1 644 645 684 683 +CQUAD4 558 1 645 646 685 684 +CQUAD4 559 1 646 647 686 685 +CQUAD4 560 1 647 54 55 686 +CQUAD4 561 1 127 648 687 126 +CQUAD4 562 1 648 649 688 687 +CQUAD4 563 1 649 650 689 688 +CQUAD4 564 1 650 651 690 689 +CQUAD4 565 1 651 652 691 690 +CQUAD4 566 1 652 653 692 691 +CQUAD4 567 1 653 654 693 692 +CQUAD4 568 1 654 655 694 693 +CQUAD4 569 1 655 656 695 694 +CQUAD4 570 1 656 657 696 695 +CQUAD4 571 1 657 658 697 696 +CQUAD4 572 1 658 659 698 697 +CQUAD4 573 1 659 660 699 698 +CQUAD4 574 1 660 661 700 699 +CQUAD4 575 1 661 662 701 700 +CQUAD4 576 1 662 663 702 701 +CQUAD4 577 1 663 664 703 702 +CQUAD4 578 1 664 665 704 703 +CQUAD4 579 1 665 666 705 704 +CQUAD4 580 1 666 667 706 705 +CQUAD4 581 1 667 668 707 706 +CQUAD4 582 1 668 669 708 707 +CQUAD4 583 1 669 670 709 708 +CQUAD4 584 1 670 671 710 709 +CQUAD4 585 1 671 672 711 710 +CQUAD4 586 1 672 673 712 711 +CQUAD4 587 1 673 674 713 712 +CQUAD4 588 1 674 675 714 713 +CQUAD4 589 1 675 676 715 714 +CQUAD4 590 1 676 677 716 715 +CQUAD4 591 1 677 678 717 716 +CQUAD4 592 1 678 679 718 717 +CQUAD4 593 1 679 680 719 718 +CQUAD4 594 1 680 681 720 719 +CQUAD4 595 1 681 682 721 720 +CQUAD4 596 1 682 683 722 721 +CQUAD4 597 1 683 684 723 722 +CQUAD4 598 1 684 685 724 723 +CQUAD4 599 1 685 686 725 724 +CQUAD4 600 1 686 55 56 725 +CQUAD4 601 1 126 687 726 125 +CQUAD4 602 1 687 688 727 726 +CQUAD4 603 1 688 689 728 727 +CQUAD4 604 1 689 690 729 728 +CQUAD4 605 1 690 691 730 729 +CQUAD4 606 1 691 692 731 730 +CQUAD4 607 1 692 693 732 731 +CQUAD4 608 1 693 694 733 732 +CQUAD4 609 1 694 695 734 733 +CQUAD4 610 1 695 696 735 734 +CQUAD4 611 1 696 697 736 735 +CQUAD4 612 1 697 698 737 736 +CQUAD4 613 1 698 699 738 737 +CQUAD4 614 1 699 700 739 738 +CQUAD4 615 1 700 701 740 739 +CQUAD4 616 1 701 702 741 740 +CQUAD4 617 1 702 703 742 741 +CQUAD4 618 1 703 704 743 742 +CQUAD4 619 1 704 705 744 743 +CQUAD4 620 1 705 706 745 744 +CQUAD4 621 1 706 707 746 745 +CQUAD4 622 1 707 708 747 746 +CQUAD4 623 1 708 709 748 747 +CQUAD4 624 1 709 710 749 748 +CQUAD4 625 1 710 711 750 749 +CQUAD4 626 1 711 712 751 750 +CQUAD4 627 1 712 713 752 751 +CQUAD4 628 1 713 714 753 752 +CQUAD4 629 1 714 715 754 753 +CQUAD4 630 1 715 716 755 754 +CQUAD4 631 1 716 717 756 755 +CQUAD4 632 1 717 718 757 756 +CQUAD4 633 1 718 719 758 757 +CQUAD4 634 1 719 720 759 758 +CQUAD4 635 1 720 721 760 759 +CQUAD4 636 1 721 722 761 760 +CQUAD4 637 1 722 723 762 761 +CQUAD4 638 1 723 724 763 762 +CQUAD4 639 1 724 725 764 763 +CQUAD4 640 1 725 56 57 764 +CQUAD4 641 1 125 726 765 124 +CQUAD4 642 1 726 727 766 765 +CQUAD4 643 1 727 728 767 766 +CQUAD4 644 1 728 729 768 767 +CQUAD4 645 1 729 730 769 768 +CQUAD4 646 1 730 731 770 769 +CQUAD4 647 1 731 732 771 770 +CQUAD4 648 1 732 733 772 771 +CQUAD4 649 1 733 734 773 772 +CQUAD4 650 1 734 735 774 773 +CQUAD4 651 1 735 736 775 774 +CQUAD4 652 1 736 737 776 775 +CQUAD4 653 1 737 738 777 776 +CQUAD4 654 1 738 739 778 777 +CQUAD4 655 1 739 740 779 778 +CQUAD4 656 1 740 741 780 779 +CQUAD4 657 1 741 742 781 780 +CQUAD4 658 1 742 743 782 781 +CQUAD4 659 1 743 744 783 782 +CQUAD4 660 1 744 745 784 783 +CQUAD4 661 1 745 746 785 784 +CQUAD4 662 1 746 747 786 785 +CQUAD4 663 1 747 748 787 786 +CQUAD4 664 1 748 749 788 787 +CQUAD4 665 1 749 750 789 788 +CQUAD4 666 1 750 751 790 789 +CQUAD4 667 1 751 752 791 790 +CQUAD4 668 1 752 753 792 791 +CQUAD4 669 1 753 754 793 792 +CQUAD4 670 1 754 755 794 793 +CQUAD4 671 1 755 756 795 794 +CQUAD4 672 1 756 757 796 795 +CQUAD4 673 1 757 758 797 796 +CQUAD4 674 1 758 759 798 797 +CQUAD4 675 1 759 760 799 798 +CQUAD4 676 1 760 761 800 799 +CQUAD4 677 1 761 762 801 800 +CQUAD4 678 1 762 763 802 801 +CQUAD4 679 1 763 764 803 802 +CQUAD4 680 1 764 57 58 803 +CQUAD4 681 1 124 765 804 123 +CQUAD4 682 1 765 766 805 804 +CQUAD4 683 1 766 767 806 805 +CQUAD4 684 1 767 768 807 806 +CQUAD4 685 1 768 769 808 807 +CQUAD4 686 1 769 770 809 808 +CQUAD4 687 1 770 771 810 809 +CQUAD4 688 1 771 772 811 810 +CQUAD4 689 1 772 773 812 811 +CQUAD4 690 1 773 774 813 812 +CQUAD4 691 1 774 775 814 813 +CQUAD4 692 1 775 776 815 814 +CQUAD4 693 1 776 777 816 815 +CQUAD4 694 1 777 778 817 816 +CQUAD4 695 1 778 779 818 817 +CQUAD4 696 1 779 780 819 818 +CQUAD4 697 1 780 781 820 819 +CQUAD4 698 1 781 782 821 820 +CQUAD4 699 1 782 783 822 821 +CQUAD4 700 1 783 784 823 822 +CQUAD4 701 1 784 785 824 823 +CQUAD4 702 1 785 786 825 824 +CQUAD4 703 1 786 787 826 825 +CQUAD4 704 1 787 788 827 826 +CQUAD4 705 1 788 789 828 827 +CQUAD4 706 1 789 790 829 828 +CQUAD4 707 1 790 791 830 829 +CQUAD4 708 1 791 792 831 830 +CQUAD4 709 1 792 793 832 831 +CQUAD4 710 1 793 794 833 832 +CQUAD4 711 1 794 795 834 833 +CQUAD4 712 1 795 796 835 834 +CQUAD4 713 1 796 797 836 835 +CQUAD4 714 1 797 798 837 836 +CQUAD4 715 1 798 799 838 837 +CQUAD4 716 1 799 800 839 838 +CQUAD4 717 1 800 801 840 839 +CQUAD4 718 1 801 802 841 840 +CQUAD4 719 1 802 803 842 841 +CQUAD4 720 1 803 58 59 842 +CQUAD4 721 1 123 804 843 122 +CQUAD4 722 1 804 805 844 843 +CQUAD4 723 1 805 806 845 844 +CQUAD4 724 1 806 807 846 845 +CQUAD4 725 1 807 808 847 846 +CQUAD4 726 1 808 809 848 847 +CQUAD4 727 1 809 810 849 848 +CQUAD4 728 1 810 811 850 849 +CQUAD4 729 1 811 812 851 850 +CQUAD4 730 1 812 813 852 851 +CQUAD4 731 1 813 814 853 852 +CQUAD4 732 1 814 815 854 853 +CQUAD4 733 1 815 816 855 854 +CQUAD4 734 1 816 817 856 855 +CQUAD4 735 1 817 818 857 856 +CQUAD4 736 1 818 819 858 857 +CQUAD4 737 1 819 820 859 858 +CQUAD4 738 1 820 821 860 859 +CQUAD4 739 1 821 822 861 860 +CQUAD4 740 1 822 823 862 861 +CQUAD4 741 1 823 824 863 862 +CQUAD4 742 1 824 825 864 863 +CQUAD4 743 1 825 826 865 864 +CQUAD4 744 1 826 827 866 865 +CQUAD4 745 1 827 828 867 866 +CQUAD4 746 1 828 829 868 867 +CQUAD4 747 1 829 830 869 868 +CQUAD4 748 1 830 831 870 869 +CQUAD4 749 1 831 832 871 870 +CQUAD4 750 1 832 833 872 871 +CQUAD4 751 1 833 834 873 872 +CQUAD4 752 1 834 835 874 873 +CQUAD4 753 1 835 836 875 874 +CQUAD4 754 1 836 837 876 875 +CQUAD4 755 1 837 838 877 876 +CQUAD4 756 1 838 839 878 877 +CQUAD4 757 1 839 840 879 878 +CQUAD4 758 1 840 841 880 879 +CQUAD4 759 1 841 842 881 880 +CQUAD4 760 1 842 59 60 881 +CQUAD4 761 1 122 843 882 121 +CQUAD4 762 1 843 844 883 882 +CQUAD4 763 1 844 845 884 883 +CQUAD4 764 1 845 846 885 884 +CQUAD4 765 1 846 847 886 885 +CQUAD4 766 1 847 848 887 886 +CQUAD4 767 1 848 849 888 887 +CQUAD4 768 1 849 850 889 888 +CQUAD4 769 1 850 851 890 889 +CQUAD4 770 1 851 852 891 890 +CQUAD4 771 1 852 853 892 891 +CQUAD4 772 1 853 854 893 892 +CQUAD4 773 1 854 855 894 893 +CQUAD4 774 1 855 856 895 894 +CQUAD4 775 1 856 857 896 895 +CQUAD4 776 1 857 858 897 896 +CQUAD4 777 1 858 859 898 897 +CQUAD4 778 1 859 860 899 898 +CQUAD4 779 1 860 861 900 899 +CQUAD4 780 1 861 862 901 900 +CQUAD4 781 1 862 863 902 901 +CQUAD4 782 1 863 864 903 902 +CQUAD4 783 1 864 865 904 903 +CQUAD4 784 1 865 866 905 904 +CQUAD4 785 1 866 867 906 905 +CQUAD4 786 1 867 868 907 906 +CQUAD4 787 1 868 869 908 907 +CQUAD4 788 1 869 870 909 908 +CQUAD4 789 1 870 871 910 909 +CQUAD4 790 1 871 872 911 910 +CQUAD4 791 1 872 873 912 911 +CQUAD4 792 1 873 874 913 912 +CQUAD4 793 1 874 875 914 913 +CQUAD4 794 1 875 876 915 914 +CQUAD4 795 1 876 877 916 915 +CQUAD4 796 1 877 878 917 916 +CQUAD4 797 1 878 879 918 917 +CQUAD4 798 1 879 880 919 918 +CQUAD4 799 1 880 881 920 919 +CQUAD4 800 1 881 60 61 920 +CQUAD4 801 1 121 882 921 120 +CQUAD4 802 1 882 883 922 921 +CQUAD4 803 1 883 884 923 922 +CQUAD4 804 1 884 885 924 923 +CQUAD4 805 1 885 886 925 924 +CQUAD4 806 1 886 887 926 925 +CQUAD4 807 1 887 888 927 926 +CQUAD4 808 1 888 889 928 927 +CQUAD4 809 1 889 890 929 928 +CQUAD4 810 1 890 891 930 929 +CQUAD4 811 1 891 892 931 930 +CQUAD4 812 1 892 893 932 931 +CQUAD4 813 1 893 894 933 932 +CQUAD4 814 1 894 895 934 933 +CQUAD4 815 1 895 896 935 934 +CQUAD4 816 1 896 897 936 935 +CQUAD4 817 1 897 898 937 936 +CQUAD4 818 1 898 899 938 937 +CQUAD4 819 1 899 900 939 938 +CQUAD4 820 1 900 901 940 939 +CQUAD4 821 1 901 902 941 940 +CQUAD4 822 1 902 903 942 941 +CQUAD4 823 1 903 904 943 942 +CQUAD4 824 1 904 905 944 943 +CQUAD4 825 1 905 906 945 944 +CQUAD4 826 1 906 907 946 945 +CQUAD4 827 1 907 908 947 946 +CQUAD4 828 1 908 909 948 947 +CQUAD4 829 1 909 910 949 948 +CQUAD4 830 1 910 911 950 949 +CQUAD4 831 1 911 912 951 950 +CQUAD4 832 1 912 913 952 951 +CQUAD4 833 1 913 914 953 952 +CQUAD4 834 1 914 915 954 953 +CQUAD4 835 1 915 916 955 954 +CQUAD4 836 1 916 917 956 955 +CQUAD4 837 1 917 918 957 956 +CQUAD4 838 1 918 919 958 957 +CQUAD4 839 1 919 920 959 958 +CQUAD4 840 1 920 61 62 959 +CQUAD4 841 1 120 921 960 119 +CQUAD4 842 1 921 922 961 960 +CQUAD4 843 1 922 923 962 961 +CQUAD4 844 1 923 924 963 962 +CQUAD4 845 1 924 925 964 963 +CQUAD4 846 1 925 926 965 964 +CQUAD4 847 1 926 927 966 965 +CQUAD4 848 1 927 928 967 966 +CQUAD4 849 1 928 929 968 967 +CQUAD4 850 1 929 930 969 968 +CQUAD4 851 1 930 931 970 969 +CQUAD4 852 1 931 932 971 970 +CQUAD4 853 1 932 933 972 971 +CQUAD4 854 1 933 934 973 972 +CQUAD4 855 1 934 935 974 973 +CQUAD4 856 1 935 936 975 974 +CQUAD4 857 1 936 937 976 975 +CQUAD4 858 1 937 938 977 976 +CQUAD4 859 1 938 939 978 977 +CQUAD4 860 1 939 940 979 978 +CQUAD4 861 1 940 941 980 979 +CQUAD4 862 1 941 942 981 980 +CQUAD4 863 1 942 943 982 981 +CQUAD4 864 1 943 944 983 982 +CQUAD4 865 1 944 945 984 983 +CQUAD4 866 1 945 946 985 984 +CQUAD4 867 1 946 947 986 985 +CQUAD4 868 1 947 948 987 986 +CQUAD4 869 1 948 949 988 987 +CQUAD4 870 1 949 950 989 988 +CQUAD4 871 1 950 951 990 989 +CQUAD4 872 1 951 952 991 990 +CQUAD4 873 1 952 953 992 991 +CQUAD4 874 1 953 954 993 992 +CQUAD4 875 1 954 955 994 993 +CQUAD4 876 1 955 956 995 994 +CQUAD4 877 1 956 957 996 995 +CQUAD4 878 1 957 958 997 996 +CQUAD4 879 1 958 959 998 997 +CQUAD4 880 1 959 62 63 998 +CQUAD4 881 1 119 960 999 118 +CQUAD4 882 1 960 961 1000 999 +CQUAD4 883 1 961 962 1001 1000 +CQUAD4 884 1 962 963 1002 1001 +CQUAD4 885 1 963 964 1003 1002 +CQUAD4 886 1 964 965 1004 1003 +CQUAD4 887 1 965 966 1005 1004 +CQUAD4 888 1 966 967 1006 1005 +CQUAD4 889 1 967 968 1007 1006 +CQUAD4 890 1 968 969 1008 1007 +CQUAD4 891 1 969 970 1009 1008 +CQUAD4 892 1 970 971 1010 1009 +CQUAD4 893 1 971 972 1011 1010 +CQUAD4 894 1 972 973 1012 1011 +CQUAD4 895 1 973 974 1013 1012 +CQUAD4 896 1 974 975 1014 1013 +CQUAD4 897 1 975 976 1015 1014 +CQUAD4 898 1 976 977 1016 1015 +CQUAD4 899 1 977 978 1017 1016 +CQUAD4 900 1 978 979 1018 1017 +CQUAD4 901 1 979 980 1019 1018 +CQUAD4 902 1 980 981 1020 1019 +CQUAD4 903 1 981 982 1021 1020 +CQUAD4 904 1 982 983 1022 1021 +CQUAD4 905 1 983 984 1023 1022 +CQUAD4 906 1 984 985 1024 1023 +CQUAD4 907 1 985 986 1025 1024 +CQUAD4 908 1 986 987 1026 1025 +CQUAD4 909 1 987 988 1027 1026 +CQUAD4 910 1 988 989 1028 1027 +CQUAD4 911 1 989 990 1029 1028 +CQUAD4 912 1 990 991 1030 1029 +CQUAD4 913 1 991 992 1031 1030 +CQUAD4 914 1 992 993 1032 1031 +CQUAD4 915 1 993 994 1033 1032 +CQUAD4 916 1 994 995 1034 1033 +CQUAD4 917 1 995 996 1035 1034 +CQUAD4 918 1 996 997 1036 1035 +CQUAD4 919 1 997 998 1037 1036 +CQUAD4 920 1 998 63 64 1037 +CQUAD4 921 1 118 999 1038 117 +CQUAD4 922 1 999 1000 1039 1038 +CQUAD4 923 1 1000 1001 1040 1039 +CQUAD4 924 1 1001 1002 1041 1040 +CQUAD4 925 1 1002 1003 1042 1041 +CQUAD4 926 1 1003 1004 1043 1042 +CQUAD4 927 1 1004 1005 1044 1043 +CQUAD4 928 1 1005 1006 1045 1044 +CQUAD4 929 1 1006 1007 1046 1045 +CQUAD4 930 1 1007 1008 1047 1046 +CQUAD4 931 1 1008 1009 1048 1047 +CQUAD4 932 1 1009 1010 1049 1048 +CQUAD4 933 1 1010 1011 1050 1049 +CQUAD4 934 1 1011 1012 1051 1050 +CQUAD4 935 1 1012 1013 1052 1051 +CQUAD4 936 1 1013 1014 1053 1052 +CQUAD4 937 1 1014 1015 1054 1053 +CQUAD4 938 1 1015 1016 1055 1054 +CQUAD4 939 1 1016 1017 1056 1055 +CQUAD4 940 1 1017 1018 1057 1056 +CQUAD4 941 1 1018 1019 1058 1057 +CQUAD4 942 1 1019 1020 1059 1058 +CQUAD4 943 1 1020 1021 1060 1059 +CQUAD4 944 1 1021 1022 1061 1060 +CQUAD4 945 1 1022 1023 1062 1061 +CQUAD4 946 1 1023 1024 1063 1062 +CQUAD4 947 1 1024 1025 1064 1063 +CQUAD4 948 1 1025 1026 1065 1064 +CQUAD4 949 1 1026 1027 1066 1065 +CQUAD4 950 1 1027 1028 1067 1066 +CQUAD4 951 1 1028 1029 1068 1067 +CQUAD4 952 1 1029 1030 1069 1068 +CQUAD4 953 1 1030 1031 1070 1069 +CQUAD4 954 1 1031 1032 1071 1070 +CQUAD4 955 1 1032 1033 1072 1071 +CQUAD4 956 1 1033 1034 1073 1072 +CQUAD4 957 1 1034 1035 1074 1073 +CQUAD4 958 1 1035 1036 1075 1074 +CQUAD4 959 1 1036 1037 1076 1075 +CQUAD4 960 1 1037 64 65 1076 +CQUAD4 961 1 117 1038 1077 116 +CQUAD4 962 1 1038 1039 1078 1077 +CQUAD4 963 1 1039 1040 1079 1078 +CQUAD4 964 1 1040 1041 1080 1079 +CQUAD4 965 1 1041 1042 1081 1080 +CQUAD4 966 1 1042 1043 1082 1081 +CQUAD4 967 1 1043 1044 1083 1082 +CQUAD4 968 1 1044 1045 1084 1083 +CQUAD4 969 1 1045 1046 1085 1084 +CQUAD4 970 1 1046 1047 1086 1085 +CQUAD4 971 1 1047 1048 1087 1086 +CQUAD4 972 1 1048 1049 1088 1087 +CQUAD4 973 1 1049 1050 1089 1088 +CQUAD4 974 1 1050 1051 1090 1089 +CQUAD4 975 1 1051 1052 1091 1090 +CQUAD4 976 1 1052 1053 1092 1091 +CQUAD4 977 1 1053 1054 1093 1092 +CQUAD4 978 1 1054 1055 1094 1093 +CQUAD4 979 1 1055 1056 1095 1094 +CQUAD4 980 1 1056 1057 1096 1095 +CQUAD4 981 1 1057 1058 1097 1096 +CQUAD4 982 1 1058 1059 1098 1097 +CQUAD4 983 1 1059 1060 1099 1098 +CQUAD4 984 1 1060 1061 1100 1099 +CQUAD4 985 1 1061 1062 1101 1100 +CQUAD4 986 1 1062 1063 1102 1101 +CQUAD4 987 1 1063 1064 1103 1102 +CQUAD4 988 1 1064 1065 1104 1103 +CQUAD4 989 1 1065 1066 1105 1104 +CQUAD4 990 1 1066 1067 1106 1105 +CQUAD4 991 1 1067 1068 1107 1106 +CQUAD4 992 1 1068 1069 1108 1107 +CQUAD4 993 1 1069 1070 1109 1108 +CQUAD4 994 1 1070 1071 1110 1109 +CQUAD4 995 1 1071 1072 1111 1110 +CQUAD4 996 1 1072 1073 1112 1111 +CQUAD4 997 1 1073 1074 1113 1112 +CQUAD4 998 1 1074 1075 1114 1113 +CQUAD4 999 1 1075 1076 1115 1114 +CQUAD4 1000 1 1076 65 66 1115 +CQUAD4 1001 1 116 1077 1116 115 +CQUAD4 1002 1 1077 1078 1117 1116 +CQUAD4 1003 1 1078 1079 1118 1117 +CQUAD4 1004 1 1079 1080 1119 1118 +CQUAD4 1005 1 1080 1081 1120 1119 +CQUAD4 1006 1 1081 1082 1121 1120 +CQUAD4 1007 1 1082 1083 1122 1121 +CQUAD4 1008 1 1083 1084 1123 1122 +CQUAD4 1009 1 1084 1085 1124 1123 +CQUAD4 1010 1 1085 1086 1125 1124 +CQUAD4 1011 1 1086 1087 1126 1125 +CQUAD4 1012 1 1087 1088 1127 1126 +CQUAD4 1013 1 1088 1089 1128 1127 +CQUAD4 1014 1 1089 1090 1129 1128 +CQUAD4 1015 1 1090 1091 1130 1129 +CQUAD4 1016 1 1091 1092 1131 1130 +CQUAD4 1017 1 1092 1093 1132 1131 +CQUAD4 1018 1 1093 1094 1133 1132 +CQUAD4 1019 1 1094 1095 1134 1133 +CQUAD4 1020 1 1095 1096 1135 1134 +CQUAD4 1021 1 1096 1097 1136 1135 +CQUAD4 1022 1 1097 1098 1137 1136 +CQUAD4 1023 1 1098 1099 1138 1137 +CQUAD4 1024 1 1099 1100 1139 1138 +CQUAD4 1025 1 1100 1101 1140 1139 +CQUAD4 1026 1 1101 1102 1141 1140 +CQUAD4 1027 1 1102 1103 1142 1141 +CQUAD4 1028 1 1103 1104 1143 1142 +CQUAD4 1029 1 1104 1105 1144 1143 +CQUAD4 1030 1 1105 1106 1145 1144 +CQUAD4 1031 1 1106 1107 1146 1145 +CQUAD4 1032 1 1107 1108 1147 1146 +CQUAD4 1033 1 1108 1109 1148 1147 +CQUAD4 1034 1 1109 1110 1149 1148 +CQUAD4 1035 1 1110 1111 1150 1149 +CQUAD4 1036 1 1111 1112 1151 1150 +CQUAD4 1037 1 1112 1113 1152 1151 +CQUAD4 1038 1 1113 1114 1153 1152 +CQUAD4 1039 1 1114 1115 1154 1153 +CQUAD4 1040 1 1115 66 67 1154 +CQUAD4 1041 1 115 1116 1155 114 +CQUAD4 1042 1 1116 1117 1156 1155 +CQUAD4 1043 1 1117 1118 1157 1156 +CQUAD4 1044 1 1118 1119 1158 1157 +CQUAD4 1045 1 1119 1120 1159 1158 +CQUAD4 1046 1 1120 1121 1160 1159 +CQUAD4 1047 1 1121 1122 1161 1160 +CQUAD4 1048 1 1122 1123 1162 1161 +CQUAD4 1049 1 1123 1124 1163 1162 +CQUAD4 1050 1 1124 1125 1164 1163 +CQUAD4 1051 1 1125 1126 1165 1164 +CQUAD4 1052 1 1126 1127 1166 1165 +CQUAD4 1053 1 1127 1128 1167 1166 +CQUAD4 1054 1 1128 1129 1168 1167 +CQUAD4 1055 1 1129 1130 1169 1168 +CQUAD4 1056 1 1130 1131 1170 1169 +CQUAD4 1057 1 1131 1132 1171 1170 +CQUAD4 1058 1 1132 1133 1172 1171 +CQUAD4 1059 1 1133 1134 1173 1172 +CQUAD4 1060 1 1134 1135 1174 1173 +CQUAD4 1061 1 1135 1136 1175 1174 +CQUAD4 1062 1 1136 1137 1176 1175 +CQUAD4 1063 1 1137 1138 1177 1176 +CQUAD4 1064 1 1138 1139 1178 1177 +CQUAD4 1065 1 1139 1140 1179 1178 +CQUAD4 1066 1 1140 1141 1180 1179 +CQUAD4 1067 1 1141 1142 1181 1180 +CQUAD4 1068 1 1142 1143 1182 1181 +CQUAD4 1069 1 1143 1144 1183 1182 +CQUAD4 1070 1 1144 1145 1184 1183 +CQUAD4 1071 1 1145 1146 1185 1184 +CQUAD4 1072 1 1146 1147 1186 1185 +CQUAD4 1073 1 1147 1148 1187 1186 +CQUAD4 1074 1 1148 1149 1188 1187 +CQUAD4 1075 1 1149 1150 1189 1188 +CQUAD4 1076 1 1150 1151 1190 1189 +CQUAD4 1077 1 1151 1152 1191 1190 +CQUAD4 1078 1 1152 1153 1192 1191 +CQUAD4 1079 1 1153 1154 1193 1192 +CQUAD4 1080 1 1154 67 68 1193 +CQUAD4 1081 1 114 1155 1194 113 +CQUAD4 1082 1 1155 1156 1195 1194 +CQUAD4 1083 1 1156 1157 1196 1195 +CQUAD4 1084 1 1157 1158 1197 1196 +CQUAD4 1085 1 1158 1159 1198 1197 +CQUAD4 1086 1 1159 1160 1199 1198 +CQUAD4 1087 1 1160 1161 1200 1199 +CQUAD4 1088 1 1161 1162 1201 1200 +CQUAD4 1089 1 1162 1163 1202 1201 +CQUAD4 1090 1 1163 1164 1203 1202 +CQUAD4 1091 1 1164 1165 1204 1203 +CQUAD4 1092 1 1165 1166 1205 1204 +CQUAD4 1093 1 1166 1167 1206 1205 +CQUAD4 1094 1 1167 1168 1207 1206 +CQUAD4 1095 1 1168 1169 1208 1207 +CQUAD4 1096 1 1169 1170 1209 1208 +CQUAD4 1097 1 1170 1171 1210 1209 +CQUAD4 1098 1 1171 1172 1211 1210 +CQUAD4 1099 1 1172 1173 1212 1211 +CQUAD4 1100 1 1173 1174 1213 1212 +CQUAD4 1101 1 1174 1175 1214 1213 +CQUAD4 1102 1 1175 1176 1215 1214 +CQUAD4 1103 1 1176 1177 1216 1215 +CQUAD4 1104 1 1177 1178 1217 1216 +CQUAD4 1105 1 1178 1179 1218 1217 +CQUAD4 1106 1 1179 1180 1219 1218 +CQUAD4 1107 1 1180 1181 1220 1219 +CQUAD4 1108 1 1181 1182 1221 1220 +CQUAD4 1109 1 1182 1183 1222 1221 +CQUAD4 1110 1 1183 1184 1223 1222 +CQUAD4 1111 1 1184 1185 1224 1223 +CQUAD4 1112 1 1185 1186 1225 1224 +CQUAD4 1113 1 1186 1187 1226 1225 +CQUAD4 1114 1 1187 1188 1227 1226 +CQUAD4 1115 1 1188 1189 1228 1227 +CQUAD4 1116 1 1189 1190 1229 1228 +CQUAD4 1117 1 1190 1191 1230 1229 +CQUAD4 1118 1 1191 1192 1231 1230 +CQUAD4 1119 1 1192 1193 1232 1231 +CQUAD4 1120 1 1193 68 69 1232 +CQUAD4 1121 1 113 1194 1233 112 +CQUAD4 1122 1 1194 1195 1234 1233 +CQUAD4 1123 1 1195 1196 1235 1234 +CQUAD4 1124 1 1196 1197 1236 1235 +CQUAD4 1125 1 1197 1198 1237 1236 +CQUAD4 1126 1 1198 1199 1238 1237 +CQUAD4 1127 1 1199 1200 1239 1238 +CQUAD4 1128 1 1200 1201 1240 1239 +CQUAD4 1129 1 1201 1202 1241 1240 +CQUAD4 1130 1 1202 1203 1242 1241 +CQUAD4 1131 1 1203 1204 1243 1242 +CQUAD4 1132 1 1204 1205 1244 1243 +CQUAD4 1133 1 1205 1206 1245 1244 +CQUAD4 1134 1 1206 1207 1246 1245 +CQUAD4 1135 1 1207 1208 1247 1246 +CQUAD4 1136 1 1208 1209 1248 1247 +CQUAD4 1137 1 1209 1210 1249 1248 +CQUAD4 1138 1 1210 1211 1250 1249 +CQUAD4 1139 1 1211 1212 1251 1250 +CQUAD4 1140 1 1212 1213 1252 1251 +CQUAD4 1141 1 1213 1214 1253 1252 +CQUAD4 1142 1 1214 1215 1254 1253 +CQUAD4 1143 1 1215 1216 1255 1254 +CQUAD4 1144 1 1216 1217 1256 1255 +CQUAD4 1145 1 1217 1218 1257 1256 +CQUAD4 1146 1 1218 1219 1258 1257 +CQUAD4 1147 1 1219 1220 1259 1258 +CQUAD4 1148 1 1220 1221 1260 1259 +CQUAD4 1149 1 1221 1222 1261 1260 +CQUAD4 1150 1 1222 1223 1262 1261 +CQUAD4 1151 1 1223 1224 1263 1262 +CQUAD4 1152 1 1224 1225 1264 1263 +CQUAD4 1153 1 1225 1226 1265 1264 +CQUAD4 1154 1 1226 1227 1266 1265 +CQUAD4 1155 1 1227 1228 1267 1266 +CQUAD4 1156 1 1228 1229 1268 1267 +CQUAD4 1157 1 1229 1230 1269 1268 +CQUAD4 1158 1 1230 1231 1270 1269 +CQUAD4 1159 1 1231 1232 1271 1270 +CQUAD4 1160 1 1232 69 70 1271 +CQUAD4 1161 1 112 1233 110 111 +CQUAD4 1162 1 1233 1234 109 110 +CQUAD4 1163 1 1234 1235 108 109 +CQUAD4 1164 1 1235 1236 107 108 +CQUAD4 1165 1 1236 1237 106 107 +CQUAD4 1166 1 1237 1238 105 106 +CQUAD4 1167 1 1238 1239 104 105 +CQUAD4 1168 1 1239 1240 103 104 +CQUAD4 1169 1 1240 1241 102 103 +CQUAD4 1170 1 1241 1242 101 102 +CQUAD4 1171 1 1242 1243 100 101 +CQUAD4 1172 1 1243 1244 99 100 +CQUAD4 1173 1 1244 1245 98 99 +CQUAD4 1174 1 1245 1246 97 98 +CQUAD4 1175 1 1246 1247 96 97 +CQUAD4 1176 1 1247 1248 95 96 +CQUAD4 1177 1 1248 1249 94 95 +CQUAD4 1178 1 1249 1250 93 94 +CQUAD4 1179 1 1250 1251 92 93 +CQUAD4 1180 1 1251 1252 91 92 +CQUAD4 1181 1 1252 1253 90 91 +CQUAD4 1182 1 1253 1254 89 90 +CQUAD4 1183 1 1254 1255 88 89 +CQUAD4 1184 1 1255 1256 87 88 +CQUAD4 1185 1 1256 1257 86 87 +CQUAD4 1186 1 1257 1258 85 86 +CQUAD4 1187 1 1258 1259 84 85 +CQUAD4 1188 1 1259 1260 83 84 +CQUAD4 1189 1 1260 1261 82 83 +CQUAD4 1190 1 1261 1262 81 82 +CQUAD4 1191 1 1262 1263 80 81 +CQUAD4 1192 1 1263 1264 79 80 +CQUAD4 1193 1 1264 1265 78 79 +CQUAD4 1194 1 1265 1266 77 78 +CQUAD4 1195 1 1266 1267 76 77 +CQUAD4 1196 1 1267 1268 75 76 +CQUAD4 1197 1 1268 1269 74 75 +CQUAD4 1198 1 1269 1270 73 74 +CQUAD4 1199 1 1270 1271 72 73 +CQUAD4 1200 1 1271 70 71 72 +ENDDATA 2d9c2e3b diff --git a/src/TACSBuckling.cpp b/src/TACSBuckling.cpp index 8984233f4..12f5ed55b 100644 --- a/src/TACSBuckling.cpp +++ b/src/TACSBuckling.cpp @@ -193,24 +193,29 @@ void TACSLinearBuckling::setSigma(TacsScalar _sigma) { (K + sigma G)^{-1} K x = lambda/(lambda - sigma) x */ -void TACSLinearBuckling::solve(TACSVec *rhs, KSMPrint *ksm_print) { +void TACSLinearBuckling::solve(TACSVec *rhs, TACSVec *u0, KSMPrint *ksm_print) { // Zero the variables assembler->zeroVariables(); if (mg) { double alpha = 1.0, beta = 0.0, gamma = 0.0; - mg->assembleJacobian(alpha, beta, gamma, res); - mg->factor(); - // Add the right-hand-side due to external forces - if (rhs) { - // res->axpy(1.0, rhs); - assembler->applyBCs(res); - } + if (u0) { + path->copyValues(u0); + } else { + mg->assembleJacobian(alpha, beta, gamma, res); + mg->factor(); + + // Add the right-hand-side due to external forces + if (rhs) { + // res->axpy(1.0, rhs); + assembler->applyBCs(res); + } - // Solve for the load path - solver->solve(res, path); - path->scale(-1.0); + // Solve for the load path + solver->solve(res, path); + path->scale(-1.0); + } assembler->setVariables(path); // Assemble the linear combination of the stiffness matrix @@ -229,18 +234,24 @@ void TACSLinearBuckling::solve(TACSVec *rhs, KSMPrint *ksm_print) { assembler->assembleMatType(TACS_STIFFNESS_MATRIX, kmat); aux_mat->copyValues(kmat); - pc->factor(); - assembler->assembleRes(res); + if (u0) { + path->copyValues(u0); + assembler->applyBCs(path); + } else { + pc->factor(); + assembler->assembleRes(res); + + // If need to add rhs + if (rhs) { + res->axpy(-1.0, rhs); + assembler->applyBCs(res); + } - // If need to add rhs - if (rhs) { - res->axpy(-1.0, rhs); - assembler->applyBCs(res); + // Solve for the load path and set the variables + solver->solve(res, path); + path->scale(-1.0); } - // Solve for the load path and set the variables - solver->solve(res, path); - path->scale(-1.0); assembler->setVariables(path); // Assemble the stiffness and geometric stiffness matrix @@ -403,6 +414,176 @@ void TACSLinearBuckling::evalEigenDVSens(int n, TACSBVec *dfdx) { dfdx->scale(-1.0 / scale); } +/* + The function computes the derivatives of the buckling eigenvalues. + + Compute the derivative of the eignevalues w.r.t. the design + variables. This function must be called after the solve function has + been called. The stiffness matrix and geometric stiffness matrix + cannot be modified from the previous call to solve. + + The original eigenvalue problem is + + K*u + lambda*G*u = 0 + + The derivative of the eigenvalue problem is given as follows: + + d(lambda)/dx = - u^{T}*(dK/dx + lambda*dG/dx)*u/(u^{T}*G*u) + + The difficulty is that the load path is determined by solving an + auxiliary linear system: + + K*path = f + + Since the geometric stiffness matrix is a function of the path, we + must compute the total derivative of the inner product of the + geometric stiffness matrix as follows: + + d(u^{T}*G*u)/dx = [ p(u^{T}*G*u)/px - psi*d(K*path)/dx ] + + where the adjoint variables psi are found by solving the linear + system: + + K*psi = d(u^{T}*G*u)/d(path) +*/ +void TACSLinearBuckling::evalEigenXptSens(int n, TACSBVec *dfdX) { + // Zero the derivative + dfdX->zeroEntries(); + + // Copy over the values of the stiffness matrix, factor + // the stiffness matrix. + aux_mat->copyValues(kmat); + pc->factor(); + + // Get the eigenvalue and eigenvector + TacsScalar error; + TacsScalar eig = extractEigenvector(n, eigvec, &error); + + // Evaluate the partial derivative for the stiffness matrix + assembler->addMatXptSensInnerProduct(1.0, TACS_STIFFNESS_MATRIX, eigvec, + eigvec, dfdX); + + // Evaluate the derivative of the geometric stiffness matrix + assembler->addMatXptSensInnerProduct( + TacsRealPart(eig), TACS_GEOMETRIC_STIFFNESS_MATRIX, eigvec, eigvec, dfdX); + + // Evaluate derivative of the inner product with respect to + // the path variables + assembler->evalMatSVSensInnerProduct(TACS_GEOMETRIC_STIFFNESS_MATRIX, eigvec, + eigvec, res); + + // Solve for the adjoint vector and evaluate the derivative of + // the adjoint-residual inner product + solver->solve(res, update); + assembler->addAdjointResXptSensProducts(-TacsRealPart(eig), 1, &update, + &dfdX); + + // Now compute the inner product: u^{T}*G*u + gmat->mult(eigvec, res); + TacsScalar scale = res->dot(eigvec); + + dfdX->beginSetValues(TACS_ADD_VALUES); + dfdX->endSetValues(TACS_ADD_VALUES); + + // Scale the final result + dfdX->scale(-1.0 / scale); +} + +/* + The function computes the partial derivatives of the buckling eigenvalues. + + Add the contribution from the partial derivative of the eigenvalues w.r.t. the + design variables. This function must be called after the solve function has + been called. The stiffness matrix and geometric stiffness matrix + cannot be modified from the previous call to solve. +*/ +void TACSLinearBuckling::addEigenDVSens(TacsScalar coef, int n, + TACSBVec *dfdx) { + // Get the eigenvalue and eigenvector + TacsScalar error; + TacsScalar eig = extractEigenvector(n, eigvec, &error); + + // Now compute the inner product: u^{T}*G*u + gmat->mult(eigvec, res); + TacsScalar scale = res->dot(eigvec); + + // Evaluate the partial derivative for the stiffness matrix + assembler->addMatDVSensInnerProduct(-coef / scale, TACS_STIFFNESS_MATRIX, + eigvec, eigvec, dfdx); + + // Evaluate the derivative of the geometric stiffness matrix + assembler->addMatDVSensInnerProduct(-eig * coef / scale, + TACS_GEOMETRIC_STIFFNESS_MATRIX, eigvec, + eigvec, dfdx); +} + +/* + The function computes the partial derivatives of the buckling eigenvalues. + + Add the contribution from the partial derivative of the eigenvalues w.r.t. the + nodal coordinates. This function must be called after the solve function has + been called. The stiffness matrix and geometric stiffness matrix + cannot be modified from the previous call to solve. +*/ +void TACSLinearBuckling::addEigenXptSens(TacsScalar coef, int n, + TACSBVec *dfdX) { + // Get the eigenvalue and eigenvector + TacsScalar error; + TacsScalar eig = extractEigenvector(n, eigvec, &error); + + // Now compute the inner product: u^{T}*G*u + gmat->mult(eigvec, res); + TacsScalar scale = res->dot(eigvec); + + // Evaluate the partial derivative for the stiffness matrix + assembler->addMatXptSensInnerProduct(-coef / scale, TACS_STIFFNESS_MATRIX, + eigvec, eigvec, dfdX); + + // Evaluate the derivative of the geometric stiffness matrix + assembler->addMatXptSensInnerProduct(-coef * eig / scale, + TACS_GEOMETRIC_STIFFNESS_MATRIX, eigvec, + eigvec, dfdX); +} + +/* + The function computes the partial derivatives of the buckling eigenvalues. + + Compute the derivative of the eigenvalues w.r.t. the design + variables. This function must be called after the solve function has + been called. The stiffness matrix and geometric stiffness matrix + cannot be modified from the previous call to solve. + + The original eigenvalue problem is + + K*u + lambda*G*u = 0 + + The derivative of the eigenvalue problem is given as follows: + + d(lambda)/du = -lambda*(u^{T}*dG/du*u)/(u^{T}*G*u) +*/ +void TACSLinearBuckling::evalEigenSVSens(int n, TACSBVec *dfdu) { + dfdu->zeroEntries(); + + // Get the eigenvalue and eigenvector + TacsScalar error; + TacsScalar eig = extractEigenvector(n, eigvec, &error); + + // Evaluate derivative of the inner product with respect to + // the path variables + assembler->evalMatSVSensInnerProduct(TACS_GEOMETRIC_STIFFNESS_MATRIX, eigvec, + eigvec, dfdu); + + // Now compute the inner product: u^{T}*G*u + gmat->mult(eigvec, res); + TacsScalar scale = res->dot(eigvec); + + dfdu->beginSetValues(TACS_ADD_VALUES); + dfdu->endSetValues(TACS_ADD_VALUES); + + // Scale the final result + dfdu->scale(-eig / scale); +} + /*! The following code computes the eigenvalues and eigenvectors for the natural frequency eigenproblem: diff --git a/src/TACSBuckling.h b/src/TACSBuckling.h index 6b4abe54d..83de19bda 100644 --- a/src/TACSBuckling.h +++ b/src/TACSBuckling.h @@ -59,8 +59,13 @@ class TACSLinearBuckling : public TACSObject { // Solve the eigenvalue problem // ---------------------------- - void solve(TACSVec *rhs = NULL, KSMPrint *ksm_print = NULL); + void solve(TACSVec *rhs = NULL, TACSVec *u0 = NULL, + KSMPrint *ksm_print = NULL); void evalEigenDVSens(int n, TACSBVec *dfdx); + void evalEigenXptSens(int n, TACSBVec *dfdX); + void evalEigenSVSens(int n, TACSBVec *dfdu); + void addEigenDVSens(TacsScalar coef, int n, TACSBVec *dfdx); + void addEigenXptSens(TacsScalar coef, int n, TACSBVec *dfdX); // Extract the eigenvalue or check the solution // -------------------------------------------- diff --git a/src/elements/shell/TACSShellElement.h b/src/elements/shell/TACSShellElement.h index 3feadc925..5cfbbeb4d 100644 --- a/src/elements/shell/TACSShellElement.h +++ b/src/elements/shell/TACSShellElement.h @@ -11,6 +11,7 @@ #include "TACSShellElementModel.h" #include "TACSShellElementTransform.h" #include "TACSShellInertialForce.h" +#include "TACSShellInplaneElementModel.h" #include "TACSShellPressure.h" #include "TACSShellTraction.h" #include "TACSShellUtilities.h" @@ -42,6 +43,21 @@ class TACSShellElement : public TACSElement { con = _con; con->incref(); + + // For linear models, we'll need to switch to a nonlinear implementation to + // capture geometric effects + if (typeid(model) == typeid(TACSShellLinearModel)) { + nlElem = new TACSShellElement(transform, con); + } else if (typeid(model) == typeid(TACSShellInplaneLinearModel)) { + nlElem = + new TACSShellElement(transform, con); + } + // For nonlinear models we can use the current class instance + else { + nlElem = this; + } } ~TACSShellElement() { @@ -52,6 +68,11 @@ class TACSShellElement : public TACSElement { if (con) { con->decref(); } + + // free nonlinear element pointer + if (nlElem != this) { + delete nlElem; + } } const char *getObjectName() { return "TACSShellElement"; } @@ -180,6 +201,7 @@ class TACSShellElement : public TACSElement { TACSShellTransform *transform; TACSShellConstitutive *con; + TACSElement *nlElem; }; /* @@ -638,20 +660,57 @@ void TACSShellElement::getMatType( memset(mat, 0, vars_per_node * num_nodes * vars_per_node * num_nodes * sizeof(TacsScalar)); - TacsScalar alpha, beta, gamma; + TacsScalar *path; + TacsScalar alpha, beta, gamma, dh, norm; alpha = beta = gamma = 0.0; + // Create dummy residual vector + TacsScalar res[vars_per_node * num_nodes]; + memset(res, 0, vars_per_node * num_nodes * sizeof(TacsScalar)); + dh = 1e-4; // Set alpha or gamma based on if this is a stiffness or mass matrix if (matType == TACS_STIFFNESS_MATRIX) { alpha = 1.0; } else if (matType == TACS_MASS_MATRIX) { gamma = 1.0; } else { // TACS_GEOMETRIC_STIFFNESS_MATRIX - // Not implimented + // Approximate geometric stiffness using directional derivative of + // tangential stiffness projected along path of current state vars + + // compute norm for normalizing path vec + norm = 0.0; + for (int i = 0; i < vars_per_node * num_nodes; i++) { + norm += vars[i] * vars[i]; + } + norm = sqrt(norm); + + if (TacsRealPart(norm) == 0.0) { + norm = 1.0; + } + + // Central difference the tangent stiffness matrix + alpha = 0.5 * norm / dh; + + // fwd step + path = new TacsScalar[vars_per_node * num_nodes]; + for (int i = 0; i < vars_per_node * num_nodes; i++) { + path[i] = dh * vars[i] / norm; + } + + nlElem->addJacobian(elemIndex, time, alpha, beta, gamma, Xpts, path, vars, + vars, res, mat); + + // bwd step + for (int i = 0; i < vars_per_node * num_nodes; i++) { + path[i] = -dh * vars[i] / norm; + } + + nlElem->addJacobian(elemIndex, time, -alpha, beta, gamma, Xpts, path, vars, + vars, res, mat); + + delete[] path; + return; } - // Create dummy residual vector - TacsScalar res[vars_per_node * num_nodes]; - memset(res, 0, vars_per_node * num_nodes * sizeof(TacsScalar)); // Add appropriate Jacobian to matrix addJacobian(elemIndex, time, alpha, beta, gamma, Xpts, vars, vars, vars, res, mat); diff --git a/tacs/TACS.pxd b/tacs/TACS.pxd index 0d6ccbd80..4c5d9f8ab 100644 --- a/tacs/TACS.pxd +++ b/tacs/TACS.pxd @@ -644,7 +644,7 @@ cdef extern from "TACSBuckling.h": TacsScalar extractEigenvalue(int, TacsScalar*) TacsScalar extractEigenvector(int, TACSBVec*, TacsScalar*) void evalEigenDVSens(int, TACSBVec*) - void evalEigenXptSens(int, TACSBVec *) + void evalEigenXptSens(int, TACSBVec*) cdef cppclass TACSLinearBuckling(TACSObject): TACSLinearBuckling( TACSAssembler *, @@ -655,10 +655,15 @@ cdef extern from "TACSBuckling.h": TACSAssembler* getAssembler() TacsScalar getSigma() void setSigma(TacsScalar) - void solve(TACSVec*, KSMPrint*) + void solve(TACSVec*, TACSVec*, KSMPrint*) void evalEigenDVSens(int, TacsScalar, int) TacsScalar extractEigenvalue(int, TacsScalar*) TacsScalar extractEigenvector(int, TACSBVec*, TacsScalar*) + void evalEigenDVSens(int, TACSBVec*) + void evalEigenXptSens(int, TACSBVec*) + void addEigenDVSens(TacsScalar, int, TACSBVec*) + void addEigenXptSens(TacsScalar, int, TACSBVec*) + void evalEigenSVSens(int, TACSBVec*) cdef extern from "TACSMeshLoader.h": cdef cppclass TACSMeshLoader(TACSObject): diff --git a/tacs/TACS.pyx b/tacs/TACS.pyx index 681ae96d5..37d842820 100644 --- a/tacs/TACS.pyx +++ b/tacs/TACS.pyx @@ -3284,8 +3284,9 @@ cdef class BucklingAnalysis: def setSigma(self, TacsScalar sigma): self.ptr.setSigma(sigma) - def solve(self, Vec force=None, print_flag=True, int freq=10): + def solve(self, Vec force=None, Vec path=None, print_flag=True, int freq=10): cdef TACSBVec *f = NULL + cdef TACSBVec *u0 = NULL cdef MPI_Comm comm cdef int rank cdef TACSAssembler *assembler = NULL @@ -3294,13 +3295,16 @@ cdef class BucklingAnalysis: if force is not None: f = force.getBVecPtr() + if path is not None: + u0 = path.getBVecPtr() + if print_flag: assembler = self.ptr.getAssembler() comm = assembler.getMPIComm() MPI_Comm_rank(comm, &rank) ksm_print = new KSMPrintStdout("BucklingAnalysis", rank, freq) - self.ptr.solve(f, ksm_print) + self.ptr.solve(f, u0, ksm_print) return def extractEigenvalue(self, int eig): @@ -3315,6 +3319,125 @@ cdef class BucklingAnalysis: eigval = self.ptr.extractEigenvector(eig, vec.getBVecPtr(), &err) return eigval, err + def evalEigenDVSens(self, int index, Vec dvsens): + """ + Compute the derivative of the eigenvalues w.r.t. the design variables + + The original eigenvalue problem is, + + K*u = lambda*G*u + + The derivative of the eigenvalue problem is given as follows, + + dK/dx*u + K*du/dx = + d(lambda)/dx*M*u + lambda*dM/dx*u + lambda*M*du/dx + + Since G = G^{T} and K = K^{T}, pre-multiplying by u^{T} gives, + + u^{T}*dK/dx*u = d(lambda)/dx*(u^{T}*G*u) + lambda*u^{T}*dG/dx*u + + Rearranging gives, + + (u^{T}*G*u)*d(lambda)/dx = u^{T}*(dK/dx - lambda*dG/dx)*u + + Args: + index (int): The index of the desired eigenvalue + dvsens (Vec): The vector in which the design variable sensitivity will be stored + """ + self.ptr.evalEigenDVSens(index, dvsens.getBVecPtr()) + + def evalEigenXptSens(self, int index, Vec xptsens): + """ + The function computes the derivatives of the buckling eigenvalues. + + Compute the derivative of the eignevalues w.r.t. the design + variables. This function must be called after the solve function has + been called. The stiffness matrix and geometric stiffness matrix + cannot be modified from the previous call to solve. + + The original eigenvalue problem is + + K*u + lambda*G*u = 0 + + The derivative of the eigenvalue problem is given as follows: + + d(lambda)/dx = - u^{T}*(dK/dx + lambda*dG/dx)*u/(u^{T}*G*u) + + The difficulty is that the load path is determined by solving an + auxiliary linear system: + + K*path = f + + Since the geometric stiffness matrix is a function of the path, we + must compute the total derivative of the inner product of the + geometric stiffness matrix as follows: + + d(u^{T}*G*u)/dx = [ p(u^{T}*G*u)/px - psi*d(K*path)/dx ] + + where the adjoint variables psi are found by solving the linear + system: + + K*psi = d(u^{T}*G*u)/d(path) + + Args: + index (int): The index of the desired eigenvalue + xptsens (Vec): The vector in which the nodal sensitivity will be stored + """ + self.ptr.evalEigenXptSens(index, xptsens.getBVecPtr()) + + def addEigenDVSens(self, TacsScalar scale, int index, Vec dvsens): + """ + The function computes the partial derivatives of the buckling eigenvalues. + + Add the contribution from the partial derivative of the eigenvalues w.r.t. the design + variables. This function must be called after the solve function has + been called. The stiffness matrix and geometric stiffness matrix + cannot be modified from the previous call to solve. + + Args: + index (int): The index of the desired eigenvalue + dvsens (Vec): The vector in which the state variable sensitivity will be stored + """ + self.ptr.addEigenDVSens(scale, index, dvsens.getBVecPtr()) + + def addEigenXptSens(self, TacsScalar scale, int index, Vec xptsens): + """ + The function computes the partial derivatives of the buckling eigenvalues. + + Add the contribution from the partial derivative of the eigenvalues w.r.t. the nodal + coordinates. This function must be called after the solve function has + been called. The stiffness matrix and geometric stiffness matrix + cannot be modified from the previous call to solve. + + Args: + index (int): The index of the desired eigenvalue + xptsens (Vec): The vector in which the state variable sensitivity will be stored + """ + self.ptr.addEigenXptSens(scale, index, xptsens.getBVecPtr()) + + def evalEigenSVSens(self, int index, Vec svsens): + """ + The function computes the derivatives of the buckling eigenvalues. + + Compute the derivative of the eigenvalues w.r.t. the design + variables. This function must be called after the solve function has + been called. The stiffness matrix and geometric stiffness matrix + cannot be modified from the previous call to solve. + + The original eigenvalue problem is + + K*u + lambda*G*u = 0 + + The derivative of the eigenvalue problem is given as follows: + + d(lambda)/du = -lambda*(u^{T}*dG/du*u)/(u^{T}*G*u) + + Args: + index (int): The index of the desired eigenvalue + svsens (Vec): The vector in which the state variable sensitivity will be stored + """ + self.ptr.evalEigenSVSens(index, svsens.getBVecPtr()) + # A generic abstract class for all integrators implemented in TACS cdef class Integrator: """ diff --git a/tacs/mphys/mphys_tacs.py b/tacs/mphys/mphys_tacs.py index 926b6b9c3..6d6b9e938 100644 --- a/tacs/mphys/mphys_tacs.py +++ b/tacs/mphys/mphys_tacs.py @@ -7,6 +7,7 @@ from mphys.builder import Builder from mphys import MaskedConverter, UnmaskedConverter, MaskedVariableDescription +import tacs.problems from .. import pyTACS, functions @@ -813,6 +814,7 @@ def initialize(self): self.options.declare("problem_setup", default=None) self.options.declare("write_solution") self.options.declare("constraint_setup", default=None) + self.options.declare("buckling_setup", default=None) def setup(self): self.fea_assembler = self.options["fea_assembler"] @@ -879,7 +881,29 @@ def setup(self): promotes_outputs=["*"], ) - self.mass_funcs.mphys_set_sp(sp) + self.mass_funcs.mphys_set_sp(sp) + + # Setup TACS problem with user-defined output functions + buckling_setup = self.options["buckling_setup"] + if buckling_setup is not None: + new_problem = buckling_setup(scenario_name, self.fea_assembler) + # Check if the user provided back a new problem to overwrite the default + if isinstance(new_problem, tacs.problems.BucklingProblem): + bp = new_problem + + # Add buckling evaluation component for eigenvalue outputs + self.add_subsystem( + "buckling", + TacsBuckling( + fea_assembler=self.fea_assembler, + check_partials=self.check_partials, + conduction=self.conduction, + write_solution=self.write_solution, + ), + promotes_inputs=promotes_inputs + promotes_states, + promotes_outputs=["*"], + ) + self.buckling.mphys_set_bp(bp) # Check if there are any user-defined TACS constraints # Constraints behave similar to "mass" functions (i.e. they don't depend on the solution state) @@ -997,6 +1021,125 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): d_inputs["x_struct0"] += Jxpt.T.dot(d_outputs[out_name]) +class TacsBuckling(om.ExplicitComponent): + """ + Component to compute non-mass TACS functions + """ + + def initialize(self): + self.options.declare("fea_assembler", recordable=False) + self.options.declare("conduction", default=False) + self.options.declare("check_partials") + self.options.declare("write_solution") + + self.fea_assembler = None + + self.check_partials = False + + def setup(self): + self.fea_assembler = self.options["fea_assembler"] + self.check_partials = self.options["check_partials"] + self.write_solution = self.options["write_solution"] + self.conduction = self.options["conduction"] + self.solution_counter = 0 + + if self.conduction: + self.states_name = "T_conduct" + else: + self.states_name = "u_struct" + + # TACS part of setup + local_ndvs = self.fea_assembler.getNumDesignVars() + + # OpenMDAO part of setup + self.add_input( + "tacs_dvs", + distributed=True, + shape=local_ndvs, + desc="tacs design variables", + tags=["mphys_coupling"], + ) + self.add_input( + "x_struct0", + distributed=True, + shape_by_conn=True, + desc="structural node coordinates", + tags=["mphys_coordinates"], + ) + self.add_input( + self.states_name, + distributed=True, + shape_by_conn=True, + desc="structural state vector", + tags=["mphys_coupling"], + ) + + def mphys_set_bp(self, bp): + # this is the external function to set the bp to this component + self.bp = bp + + # Add eval funcs as outputs + for func_name in self.bp.functionList: + func_name = func_name.replace(".", "_") + self.add_output( + func_name, distributed=False, shape=1, tags=["mphys_result"] + ) + + def _update_internal(self, inputs): + self.bp.setDesignVars(inputs["tacs_dvs"]) + self.bp.setNodes(inputs["x_struct0"]) + + def compute(self, inputs, outputs): + self._update_internal(inputs) + + # Solve + self.bp.solve(u0=inputs[self.states_name]) + + # Evaluate functions + funcs = {} + self.bp.evalFunctions(funcs, evalFuncs=outputs.keys()) + for out_name in outputs: + eig_name = out_name.replace("_", ".") + func_key = f"{self.bp.name}_{eig_name}" + self.bp.evalFunctions(funcs, evalFuncs=[eig_name]) + outputs[out_name] = funcs[func_key] + + if self.write_solution: + # write the solution files. + self.bp.writeSolution(number=self.solution_counter) + self.solution_counter += 1 + + def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): + if mode == "fwd": + if not self.check_partials: + raise ValueError("TACS forward mode requested but not implemented") + if mode == "rev": + # always update internal because same tacs object could be used by multiple scenarios + # and we need to load this scenario's state back into TACS before doing derivatives + self._update_internal(inputs) + + for func_name in d_outputs: + d_func = d_outputs[func_name] + eig_name = func_name.replace("_", ".") + mode_i = self.bp.functionList[eig_name] + + if d_func[0] != 0.0: + if "tacs_dvs" in d_inputs: + self.bp.addDVSens( + [mode_i], [d_inputs["tacs_dvs"]], scale=d_func + ) + + if "x_struct0" in d_inputs: + self.bp.addXptSens( + [mode_i], [d_inputs["x_struct0"]], scale=d_func + ) + + if self.states_name in d_inputs: + sv_sens = np.zeros_like(d_inputs[self.states_name]) + self.bp.evalSVSens([mode_i], [sv_sens]) + d_inputs[self.states_name][:] += sv_sens * d_func + + class TacsBuilder(Builder): def __init__( self, @@ -1046,6 +1189,12 @@ def initialize(self, comm): else: self.constraint_setup = None + # Load optional user-defined callback function for setting up buckling problem + if "buckling_setup" in pytacs_options: + self.buckling_setup = pytacs_options.pop("buckling_setup") + else: + self.buckling_setup = None + # Create pytacs instance self.fea_assembler = pyTACS(bdf_file, options=pytacs_options, comm=comm) self.comm = comm @@ -1087,6 +1236,7 @@ def get_post_coupling_subsystem(self, scenario_name=None): scenario_name=scenario_name, problem_setup=self.problem_setup, constraint_setup=self.constraint_setup, + buckling_setup=self.buckling_setup, ) def get_ndof(self): diff --git a/tacs/problems/__init__.py b/tacs/problems/__init__.py index 5f9c0e6dd..5c672615d 100644 --- a/tacs/problems/__init__.py +++ b/tacs/problems/__init__.py @@ -1,6 +1,13 @@ from .static import StaticProblem from .transient import TransientProblem from .modal import ModalProblem +from .buckling import BucklingProblem from .base import TACSProblem -__all__ = ["StaticProblem", "TransientProblem", "ModalProblem", "TACSProblem"] +__all__ = [ + "StaticProblem", + "TransientProblem", + "ModalProblem", + "BucklingProblem", + "TACSProblem", +] diff --git a/tacs/problems/base.py b/tacs/problems/base.py index c92c8a790..79231c06b 100644 --- a/tacs/problems/base.py +++ b/tacs/problems/base.py @@ -734,6 +734,11 @@ def _addLoadFromBDF(self, FVec, auxElems, loadID, setScale=1.0): setScale : float Factor to scale the BDF loads by before adding to problem. """ + # Make sure bdf is cross-referenced + if self.bdfInfo.is_xrefed is False: + self.bdfInfo.cross_reference() + self.bdfInfo.is_xrefed = True + vpn = self.assembler.getVarsPerNode() # Get loads and scalers for this load case ID loadSet, loadScale, _ = self.bdfInfo.get_reduced_loads(loadID) diff --git a/tacs/problems/buckling.py b/tacs/problems/buckling.py new file mode 100644 index 000000000..887b83a40 --- /dev/null +++ b/tacs/problems/buckling.py @@ -0,0 +1,972 @@ +""" +The main purpose of this class is to represent all relevant +information for a linearized buckling analysis. + +.. note:: This class should be created using the + :meth:`pyTACS.createBucklingProblem ` method. +""" + +# ============================================================================= +# Imports +# ============================================================================= +import os +import time + +import numpy as np + +import tacs.TACS +from .base import TACSProblem + + +class BucklingProblem(TACSProblem): + # Default Option List + defaultOptions = { + "outputDir": [str, "./", "Output directory for F5 file writer."], + # Solution Options + "L2Convergence": [ + float, + 1e-12, + "Absolute convergence tolerance for Eigenvalue solver based on l2 norm of residual.", + ], + "L2ConvergenceRel": [ + float, + 1e-12, + "Relative convergence tolerance for Eigenvalue solver based on l2 norm of residual.", + ], + "RBEStiffnessScaleFactor": [ + float, + 1e3, + "Constraint matrix scaling factor used in RBE Lagrange multiplier stiffness matrix.", + ], + "RBEArtificialStiffness": [ + float, + 1e-3, + "Artificial constant added to diagonals of RBE Lagrange multiplier stiffness matrix \n" + "\t to stabilize preconditioner.", + ], + "subSpaceSize": [ + int, + 10, + "Subspace size for Krylov solver used by Eigenvalue solver.", + ], + "nRestarts": [ + int, + 15, + "Max number of resets for Krylov solver used by Eigenvalue solver.", + ], + # Output Options + "writeSolution": [bool, True, "Flag for suppressing all f5 file writing."], + "numberSolutions": [ + bool, + True, + "Flag for attaching solution counter index to f5 files.", + ], + "printTiming": [ + bool, + False, + "Flag for printing out timing information for class procedures.", + ], + "printLevel": [ + int, + 0, + "Print level for Eigenvalue solver.\n" + "\t Accepts:\n" + "\t\t 0 : No printing.\n" + "\t\t 1 : Print major iterations.\n" + "\t\t > 1 : Print major + minor iterations.", + ], + } + + def __init__( + self, + name, + sigma, + numEigs, + assembler, + comm, + outputViewer=None, + meshLoader=None, + options=None, + ): + """ + NOTE: This class should not be initialized directly by the user. + Use pyTACS.createBucklingProblem instead. + + Parameters + ---------- + name : str + Name of this tacs problem + + sigma : float + Guess for the lowest eigenvalue. This corresponds to the lowest buckling load factor. + + numEigs : int + Number of eigenvalues to solve for + + assembler : TACS.Assembler + Cython object responsible for creating and setting tacs objects used to solve problem + + comm : mpi4py.MPI.Intracomm + The comm object on which to create the pyTACS object. + + outputViewer : TACS.TACSToFH5 + Cython object used to write out f5 files that can be converted and used for postprocessing. + + meshLoader : pymeshloader.pyMeshLoader + pyMeshLoader object used to create the assembler. + + options : dict + Dictionary holding problem-specific option parameters (case-insensitive). + """ + + # Problem name + self.name = name + + # Default setup for common problem class objects, sets up comm and options + TACSProblem.__init__(self, assembler, comm, options, outputViewer, meshLoader) + + # Set time eigenvalue parameters + self.sigma = sigma + self.numEigs = numEigs + + # String name used in evalFunctions + self.valName = "eigsb" + self._initializeFunctionList() + + # Create problem-specific variables + self._createVariables() + + def _createVariables(self): + """ + Internal to create the objects required by TACS Integrator + """ + + self.callCounter = -1 + + # Buckling load state + self.u0 = self.assembler.createVec() + # Load vector + self.F = self.assembler.createVec() + self.F_array = self.F.getArray() + # RHS vector + self.rhs = self.assembler.createVec() + # Auxiliary element object for applying tractions/pressure + self.auxElems = tacs.TACS.AuxElements() + + self.aux = self.assembler.createSchurMat() + self.G = self.assembler.createSchurMat() + self.K = self.assembler.createSchurMat() + + self.pc = tacs.TACS.Pc(self.K) + + # Set artificial stiffness factors in rbe class + c1 = self.getOption("RBEStiffnessScaleFactor") + c2 = self.getOption("RBEArtificialStiffness") + tacs.elements.RBE2.setScalingParameters(c1, c2) + tacs.elements.RBE3.setScalingParameters(c1, c2) + + # Assemble and factor the stiffness/Jacobian matrix. Factor the + # Jacobian and solve the linear system for the displacements + self.assembler.assembleMatType(tacs.TACS.STIFFNESS_MATRIX, self.K) + self.assembler.assembleMatType(tacs.TACS.GEOMETRIC_STIFFNESS_MATRIX, self.G) + + subspace = self.getOption("subSpaceSize") + restarts = self.getOption("nRestarts") + atol = self.getOption("L2Convergence") + rtol = self.getOption("L2ConvergenceRel") + self.gmres = tacs.TACS.KSM(self.aux, self.pc, subspace, restarts) + self.gmres.setTolerances(rtol, atol) + + # Create the buckling analysis object + self.buckleSolver = tacs.TACS.BucklingAnalysis( + self.assembler, + self.sigma, + self.G, + self.K, + self.gmres, + num_eigs=self.numEigs, + eig_tol=rtol, + ) + + def _initializeFunctionList(self): + """ + Create FunctionList dict which maps eigenvalue strings + to mode indices used in evalFunctions method. + """ + self.functionList = {} + for mode_i in range(self.numEigs): + self.functionList[f"{self.valName}.{mode_i}"] = mode_i + + def setOption(self, name, value): + """ + Set a solver option value. The name is not case sensitive. + + Parameters + ---------- + name : str + Name of option to modify + + value : depends on option + New option value to set + """ + # Default setOption for common problem class objects + TACSProblem.setOption(self, name, value) + + # No need to reset solver for output options + if name.lower() in [ + "writesolution", + "printtiming", + "numbersolutions", + "outputdir", + ]: + pass + # Reset solver for all other option changes + else: + self._createVariables() + + def setValName(self, valName): + """ + Set a name for the eigenvalues. Only needs + to be changed if more than 1 pytacs object is used in an + optimization + + Parameters + ---------- + valName : str + Name of the eigenvalue output used in evalFunctions(). + """ + self.valName = valName + self._initializeFunctionList() + + def getNumEigs(self): + """ + Get the number of eigenvalues requested from solver for this problem. + + Returns + ---------- + numEigs : int + Number of eigenvalues. + """ + return self.numEigs + + def addFunction(self, funcName, funcHandle, compIDs=None, **kwargs): + """ + NOT SUPPORTED FOR THIS PROBLEM + """ + self._TACSWarning("addFunction method not supported for this class.") + + def evalFunctions(self, funcs, evalFuncs=None, ignoreMissing=False): + """ + Evaluate eigenvalues for problem. The functions corresponding to + the integers in evalFuncs are evaluated and updated into + the provided dictionary. + + Parameters + ---------- + funcs : dict + Dictionary into which the functions are saved. + evalFuncs : iterable object containing strings. + If not none, use these functions to evaluate. + ignoreMissing : bool + Flag to supress checking for a valid function. Please use + this option with caution. + + Examples + -------- + >>> funcs = {} + >>> bucklingProblem.solve() + >>> bucklingProblem.evalFunctions(funcs, 'eigsm.0') + >>> funcs + >>> # Result will look like (if bucklingProblem has name of 'c1'): + >>> # {'c1_eigsm.0':12354.10} + """ + # Make sure assembler variables are up to date + self._updateAssemblerVars() + + # Check if user specified which eigvals to output + # Otherwise, output them all + if evalFuncs is None: + evalFuncs = self.functionList + else: + userFuncs = sorted(list(evalFuncs)) + evalFuncs = {} + for func in userFuncs: + if func in self.functionList: + evalFuncs[func] = self.functionList[func] + + if not ignoreMissing: + for f in evalFuncs: + if f not in self.functionList: + raise self._TACSError( + f"Supplied function '{f}' has not been added " + "using addFunction()." + ) + + # Loop through each requested eigenvalue + for funcName in evalFuncs: + mode_i = evalFuncs[funcName] + key = f"{self.name}_{funcName}" + funcs[key], _ = self.getVariables(mode_i) + + def evalFunctionsSens(self, funcsSens, evalFuncs=None): + """ + This is the main routine for returning useful (sensitivity) + information from problem. The derivatives of the functions + corresponding to the strings in evalFuncs are evaluated and + updated into the provided dictionary. The derivitives with + respect to all design variables and node locations are computed. + + Parameters + ---------- + funcsSens : dict + Dictionary into which the derivatives are saved. + evalFuncs : iterable object containing strings + The functions the user wants returned + + Examples + -------- + >>> funcsSens = {} + >>> bucklingProblem.evalFunctionsSens(funcsSens, 'eigsm.0') + >>> funcsSens + >>> # Result will look like (if bucklingProblem has name of 'c1'): + >>> # {'c1_eigsm.0':{'struct':[1.234, ..., 7.89], 'Xpts':[3.14, ..., 1.59]}} + """ + # Make sure assembler variables are up to date + self._updateAssemblerVars() + + # Check if user specified which eigvals to output + # Otherwise, output them all + if evalFuncs is None: + evalFuncs = self.functionList + else: + userFuncs = sorted(list(evalFuncs)) + evalFuncs = {} + for func in userFuncs: + if func in self.functionList: + evalFuncs[func] = self.functionList[func] + + dvSens = self.assembler.createDesignVec() + xptSens = self.assembler.createNodeVec() + + indices = [evalFuncs[funcName] for funcName in evalFuncs] + dvSensList = [self.assembler.createDesignVec() for funcName in evalFuncs] + xptSensList = [self.assembler.createNodeVec() for funcName in evalFuncs] + svSensList = [self.assembler.createVec() for funcName in evalFuncs] + adjointList = [self.assembler.createVec() for funcName in evalFuncs] + + self.addDVSens(indices, dvSensList, scale=1.0) + self.addXptSens(indices, xptSensList, scale=1.0) + self.evalSVSens(indices, svSensList) + + self.aux.copyValues(self.K) + self.pc.factor() + + # Loop through each requested eigenvalue + for i, funcName in enumerate(evalFuncs): + rhs = svSensList[i] + adjoint = adjointList[i] + self.gmres.solve(rhs, adjoint) + + # Evaluate adjoint contribution to nodal sens + xptSens = xptSensList[i] + self.assembler.addAdjointResXptSensProducts([adjoint], [xptSens], -1.0) + xptSens.beginSetValues() + xptSens.endSetValues() + # Evaluate adjoint contribution to dv sens + dvSens = dvSensList[i] + self.assembler.addAdjointResProducts([adjoint], [dvSens], -1.0) + dvSens.beginSetValues() + dvSens.endSetValues() + + key = f"{self.name}_{funcName}" + funcsSens[key] = {} + funcsSens[key][self.varName] = dvSens.getArray().copy() + funcsSens[key][self.coordName] = xptSens.getArray().copy() + + def addLoadToComponents(self, compIDs, F, averageLoad=False): + """ + This method is used to add a *FIXED TOTAL LOAD* on one or more + components, defined by COMPIDs. The purpose of this routine is to add loads that + remain fixed throughout an optimization. An example would be an engine load. + This routine determines all the unique nodes in the FE model that are part of the + requested components, then takes the total 'force' by F and divides by the + number of nodes. This average load is then applied to the nodes. + + Parameters + ---------- + + compIDs : list[int] or int + The components with added loads. Use pyTACS selectCompIDs method + to determine this. + + F : numpy.ndarray 1d or 2d length (varsPerNodes) or (numCompIDs, varsPerNodes) + Vector(s) of 'force' to apply to each component. If only one force vector is provided, + force will be copied uniformly across all components. + + averageLoad : bool + Flag to determine whether load should be split evenly across all components (True) + or copied and applied individually to each component (False). Defaults to False. + + Notes + ----- + + The units of the entries of the 'force' vector F are not + necessarily physical forces and their interpretation depends + on the physics problem being solved and the dofs included + in the model. + + A couple of examples of force vector components for common problems are listed below: + + In Heat Conduction with varsPerNode = 1 + F = [Qdot] # heat rate + In Elasticity with varsPerNode = 3, + F = [fx, fy, fz] # forces + In Elasticity with varsPerNode = 6, + F = [fx, fy, fz, mx, my, mz] # forces + moments + In Thermoelasticity with varsPerNode = 4, + F = [fx, fy, fz, Qdot] # forces + heat rate + In Thermoelasticity with varsPerNode = 7, + F = [fx, fy, fz, mx, my, mz, Qdot] # forces + moments + heat rate + """ + self._addLoadToComponents(self.F, compIDs, F, averageLoad) + + def addLoadToNodes(self, nodeIDs, F, nastranOrdering=False): + """ + This method is used to add a fixed point load of F to the + selected node IDs. + + Parameters + ---------- + + nodeIDs : list[int] + The nodes IDs with added loads. + + F : Numpy 1d or 2d array length (varsPerNodes) or (numNodeIDs, varsPerNodes) + Array of force vectors, one for each node. If only one force vector is provided, + force will be copied uniformly across all nodes. + + nastranOrdering : bool + Flag signaling whether nodeIDs are in TACS (default) + or NASTRAN (grid IDs in bdf file) ordering + + Notes + ----- + + The units of the entries of the 'force' vector F are not + necessarily physical forces and their interpretation depends + on the physics problem being solved and the dofs included + in the model. + + A couple of examples of force vector components for common problems are listed below: + + In Heat Conduction with varsPerNode = 1 + F = [Qdot] # heat rate + In Elasticity with varsPerNode = 3, + F = [fx, fy, fz] # forces + In Elasticity with varsPerNode = 6, + F = [fx, fy, fz, mx, my, mz] # forces + moments + In Thermoelasticity with varsPerNode = 4, + F = [fx, fy, fz, Qdot] # forces + heat rate + In Thermoelasticity with varsPerNode = 7, + F = [fx, fy, fz, mx, my, mz, Qdot] # forces + moments + heat rate + """ + + self._addLoadToNodes(self.F, nodeIDs, F, nastranOrdering) + + def addLoadToRHS(self, Fapplied): + """ + This method is used to add a *FIXED TOTAL LOAD* directly to the + right hand side vector given the equation below: + + K*u = f + + Where: + - K : Stiffness matrix for problem + - u : State variables for problem + - f : Right-hand side vector to add loads to + + Parameters + ---------- + + Fapplied : numpy.ndarray or tacs.TACS.Vec + Distributed array containing loads to applied to RHS of the problem. + + """ + self._addLoadToRHS(self.F, Fapplied) + + def addTractionToComponents(self, compIDs, tractions, faceIndex=0): + """ + This method is used to add a *FIXED TOTAL TRACTION* on one or more + components, defined by COMPIDs. The purpose of this routine is + to add loads that remain fixed throughout an optimization. + + Parameters + ---------- + + compIDs : list[int] or int + The components with added loads. Use pyTACS selectCompIDs method + to determine this. + + tractions : Numpy array length 1 or compIDs + Array of traction vectors for each component + + faceIndex : int + Indicates which face (side) of element to apply traction to. + Note: not required for certain elements (i.e. shells) + """ + self._addTractionToComponents(self.auxElems, compIDs, tractions, faceIndex) + + def addTractionToElements( + self, elemIDs, tractions, faceIndex=0, nastranOrdering=False + ): + """ + This method is used to add a fixed traction to the + selected element IDs. Tractions can be specified on an + element by element basis (if tractions is a 2d array) or + set to a uniform value (if tractions is a 1d array) + + Parameters + ---------- + + elemIDs : list[int] + The global element ID numbers for which to apply the traction. + + tractions : Numpy 1d or 2d array length varsPerNodes or (elemIDs, varsPerNodes) + Array of traction vectors for each element + + faceIndex : int + Indicates which face (side) of element to apply traction to. + Note: not required for certain elements (i.e. shells) + + nastranOrdering : bool + Flag signaling whether elemIDs are in TACS (default) + or NASTRAN ordering + """ + + self._addTractionToElements( + self.auxElems, elemIDs, tractions, faceIndex, nastranOrdering + ) + + def addPressureToComponents(self, compIDs, pressures, faceIndex=0): + """ + This method is used to add a *FIXED TOTAL PRESSURE* on one or more + components, defined by COMPIds. The purpose of this routine is + to add loads that remain fixed throughout an optimization. An example + would be a fuel load. + + Parameters + ---------- + + compIDs : list[int] or int + The components with added loads. Use pyTACS selectCompIDs method + to determine this. + + pressures : Numpy array length 1 or compIDs + Array of pressure values for each component + + faceIndex : int + Indicates which face (side) of element to apply pressure to. + Note: not required for certain elements (i.e. shells) + """ + self._addPressureToComponents(self.auxElems, compIDs, pressures, faceIndex) + + def addPressureToElements( + self, elemIDs, pressures, faceIndex=0, nastranOrdering=False + ): + """ + This method is used to add a fixed presure to the + selected element IDs. Pressures can be specified on an + element by element basis (if pressures is an array) or + set to a uniform value (if pressures is a scalar) + + Parameters + ---------- + + elemIDs : list[int] + The global element ID numbers for which to apply the pressure. + + pressures : Numpy array length 1 or elemIDs + Array of pressure values for each element + + faceIndex : int + Indicates which face (side) of element to apply pressure to. + Note: not required for certain elements (i.e. shells) + + nastranOrdering : bool + Flag signaling whether elemIDs are in TACS (default) + or NASTRAN ordering + """ + + self._addPressureToElements( + self.auxElems, elemIDs, pressures, faceIndex, nastranOrdering + ) + + def addInertialLoad(self, inertiaVector): + """ + This method is used to add a fixed inertial load due to + a uniform acceleration over the entire model. + This is most commonly used to model gravity loads on a model. + + Parameters + ---------- + inertiaVector : numpy.ndarray + Acceleration vector used to define inertial load. + """ + self._addInertialLoad(self.auxElems, inertiaVector) + + def addCentrifugalLoad(self, omegaVector, rotCenter, firstOrder=False): + """ + This method is used to add a fixed centrifugal load due to a + uniform rotational velocity over the entire model. + This is most commonly used to model rotors, rolling aircraft, etc. + + Parameters + ---------- + + omegaVector : numpy.ndarray + Rotational velocity vector (rad/s) used to define centrifugal load. + + rotCenter : numpy.ndarray + Location of center of rotation used to define centrifugal load. + + firstOrder : bool, optional + Whether to use first order approximation for centrifugal load, + which computes the force in the displaced position. By default False + """ + self._addCentrifugalLoad(self.auxElems, omegaVector, rotCenter, firstOrder) + + def addLoadFromBDF(self, loadID, scale=1.0): + """ + This method is used to add a fixed load set defined in the BDF file to the problem. + Currently, only supports LOAD, FORCE, MOMENT, GRAV, RFORCE, PLOAD2, and PLOAD4. + + Parameters + ---------- + + loadID : int + Load identification number of load set in BDF file user wishes to add to problem. + + scale : float + Factor to scale the BDF loads by before adding to problem. + """ + self._addLoadFromBDF(self.F, self.auxElems, loadID, scale) + + def zeroLoads(self): + """ + Zero all applied loads + """ + self.F.zeroEntries() + self.auxElems = tacs.TACS.AuxElements() + + ####### Buckling solver methods ######## + + def _updateAssemblerVars(self): + """ + Make sure that the assembler is using + the input variables associated with this problem + """ + + self.assembler.setDesignVars(self.x) + self.assembler.setNodes(self.Xpts) + # Set state variables + self.assembler.setVariables(self.u0) + # Zero any time derivative terms + self.assembler.zeroDotVariables() + self.assembler.zeroDDotVariables() + # Make sure previous auxiliary loads are removed + self.assembler.setAuxElements(self.auxElems) + # Set artificial stiffness factors in rbe class + c1 = self.getOption("RBEStiffnessScaleFactor") + c2 = self.getOption("RBEArtificialStiffness") + tacs.elements.RBE2.setScalingParameters(c1, c2) + tacs.elements.RBE3.setScalingParameters(c1, c2) + + def solve(self, Fext=None, u0=None): + """ + Solve the eigenvalue problem. The + forces must already be set. + + Parameters + ---------- + Fext : tacs.TACS.Vec or numpy.ndarray, optional + Distributed array containing additional loads (ex. aerodynamic forces for aerostructural coupling) + to applied to RHS of the static problem. + + u0 : tacs.TACS.Vec or numpy.ndarray, optional + Distributed array containing additional loads (ex. aerodynamic forces for aerostructural coupling) + to applied to RHS of the static problem. Alternate to Fext. + """ + startTime = time.time() + + self.callCounter += 1 + + setupProblemTime = time.time() + + # Set problem vars to assembler + self._updateAssemblerVars() + + # states were not prescribed, pass in forces + if u0 is None: + # Sum the forces from the loads not handled by TACS + self.rhs.copyValues(self.F) # Fixed loads + + # Add external loads, if specified + if Fext is not None: + if isinstance(Fext, tacs.TACS.Vec): + self.rhs.axpy(1.0, Fext) + elif isinstance(Fext, np.ndarray): + rhsArray = self.rhs.getArray() + rhsArray[:] = rhsArray[:] + Fext[:] + + force = self.rhs + path = None + + # states already prescribed, we don't need forces + else: + # Convert to bvec + if isinstance(u0, np.ndarray): + path = self._arrayToVec(u0) + else: + path = u0 + force = None + + initSolveTime = time.time() + + # Solve the buckling analysis problem + self.buckleSolver.solve( + force=force, path=path, print_flag=self.getOption("printLevel") + ) + + # Save state vars + self.assembler.getVariables(self.u0) + + solveTime = time.time() + + # If timing was was requested print it, if the solution is nonlinear + # print this information automatically if prinititerations was requested. + if self.getOption("printTiming"): + self._pp("+--------------------------------------------------+") + self._pp("|") + self._pp("| TACS Solve Times:") + self._pp("|") + self._pp( + "| %-30s: %10.3f sec" + % ("TACS Setup Time", setupProblemTime - startTime) + ) + self._pp( + "| %-30s: %10.3f sec" + % ("TACS Solve Init Time", initSolveTime - setupProblemTime) + ) + self._pp( + "| %-30s: %10.3f sec" % ("TACS Solve Time", solveTime - initSolveTime) + ) + self._pp("|") + self._pp( + "| %-30s: %10.3f sec" + % ("TACS Total Solution Time", solveTime - startTime) + ) + self._pp("+--------------------------------------------------+") + + return + + def getVariables(self, index, states=None): + """ + Return the current state values for one mode of the current problem + + Parameters + ---------- + index : int + Mode index to return solution for. + + states : tacs.TACS.Vec or numpy.ndarray or None + Place eigenvector for mode into this array (optional). + + Returns + -------- + eigVal: float + Eigenvalue for mode corresponds to buckling load factor + + states : numpy.ndarray + Eigenvector for mode + """ + eigVal, err = self.buckleSolver.extractEigenvalue(index) + eigVector = self.assembler.createVec() + self.buckleSolver.extractEigenvector(index, eigVector) + # Inplace assignment if vectors were provided + if isinstance(states, tacs.TACS.Vec): + states.copyValues(eigVector) + elif isinstance(states, np.ndarray): + states[:] = eigVector.getArray() + return eigVal, eigVector.getArray() + + def addXptSens(self, indices, xptSensList, scale=1.0): + """ + Add partial sensitivity contribution due to nodal coordinates for eigenvalue + + Parameters + ---------- + indices : list[int] + Mode indices to return sensitivity for. + + xptSensList : list[tacs.TACS.Vec] or list[numpy.ndarray] + List of sensitivity vectors to add partial sensitivity to + + scale : float + Scalar to multiply partial sensitivity by. Defaults to 1.0 + """ + # Set problem vars to assembler + self._updateAssemblerVars() + + for index, xptSens in zip(indices, xptSensList): + # Create a tacs BVec copy for the operation if the output is a numpy array + if isinstance(xptSens, np.ndarray): + xptSensBVec = self._arrayToNodeVec(xptSens) + # Otherwise the input is already a BVec and we can do the operation in place + else: + xptSensBVec = xptSens + + self.buckleSolver.addEigenXptSens(scale, index, xptSensBVec) + + # Finalize sensitivity arrays across all procs + xptSensBVec.beginSetValues() + xptSensBVec.endSetValues() + + # Update from the BVec values, if the input was a numpy array + if isinstance(xptSens, np.ndarray): + # Copy values to numpy array + xptSens[:] = xptSensBVec.getArray() + + def addDVSens(self, indices, dvSensList, scale=1.0): + """ + Add partial sensitivity contribution due to design variables for eigenvalue + + Parameters + ---------- + indices : list[int] + Mode indices to return sensitivity for. + + dvSensList : list[tacs.TACS.Vec] or list[numpy.ndarray] + List of sensitivity vectors to add partial sensitivity to + + scale : float + Scalar to multiply partial sensitivity by. Defaults to 1.0 + """ + # Set problem vars to assembler + self._updateAssemblerVars() + + for index, dvSens in zip(indices, dvSensList): + # Create a tacs BVec copy for the operation if the output is a numpy array + if isinstance(dvSens, np.ndarray): + dvSensBVec = self._arrayToDesignVec(dvSens) + # Otherwise the input is already a BVec and we can do the operation in place + else: + dvSensBVec = dvSens + + self.buckleSolver.addEigenDVSens(scale, index, dvSensBVec) + + # Finalize sensitivity arrays across all procs + dvSensBVec.beginSetValues() + dvSensBVec.endSetValues() + + # Update from the BVec values, if the input was a numpy array + if isinstance(dvSens, np.ndarray): + # Copy values to numpy array + dvSens[:] = dvSensBVec.getArray() + + def evalSVSens(self, indices, svSensList): + """ + Add partial sensitivity contribution due to state variables for eigenvalue + + Parameters + ---------- + indices : list[int] + Mode indices to return sensitivity for. + + svSensList : list[tacs.TACS.Vec] or list[numpy.ndarray] + List of sensitivity vectors to add partial sensitivity to + """ + # Set problem vars to assembler + self._updateAssemblerVars() + + for index, svSens in zip(indices, svSensList): + # Create a tacs BVec copy for the operation if the output is a numpy array + if isinstance(svSens, np.ndarray): + svSensBVec = self._arrayToVec(svSens) + # Otherwise the input is already a BVec and we can do the operation in place + else: + svSensBVec = svSens + + self.buckleSolver.evalEigenSVSens(index, svSensBVec) + + # Finalize sensitivity arrays across all procs + svSensBVec.beginSetValues() + svSensBVec.endSetValues() + + # Update from the BVec values, if the input was a numpy array + if isinstance(svSens, np.ndarray): + # Copy values to numpy array + svSens[:] = svSensBVec.getArray() + + def writeSolution(self, outputDir=None, baseName=None, number=None, indices=None): + """ + This is a generic shell function that writes the output + file(s). The intent is that the user or calling program can + call this function and pyTACS writes all the files that the + user has defined. It is recommended that this function is used + along with the associated logical flags in the options to + determine the desired writing procedure + + Parameters + ---------- + outputDir : str or None + Use the supplied output directory + baseName : str or None + Use this supplied string for the base filename. Typically + only used from an external solver. + number : int or None + Use the user supplied number to index solution. Again, only + typically used from an external solver + indices : int or list[int] or None + Mode index or indices to get state variables for. + If None, returns a solution for all modes. + Defaults to None. + """ + # Make sure assembler variables are up-to-date + self._updateAssemblerVars() + + # Check input + if outputDir is None: + outputDir = self.getOption("outputDir") + + if baseName is None: + baseName = self.name + + # If we are numbering solution, it saving the sequence of + # calls, add the call number + if number is not None: + # We need number based on the provided number: + baseName = baseName + "_%3.3d" % number + else: + # if number is none, i.e. standalone, but we need to + # number solutions, use internal counter + if self.getOption("numberSolutions"): + baseName = baseName + "_%3.3d" % self.callCounter + + # Unless the writeSolution option is off write actual file: + if self.getOption("writeSolution"): + # If indices is None, output all modes + if indices is None: + indices = np.arange(self.numEigs) + + # Write out each specified mode + indices = np.atleast_1d(indices) + vec = self.assembler.createVec() + for index in indices: + # Extract eigenvector + eig, _ = self.getVariables(index, states=vec) + # Set eigen mode in assembler + self.assembler.setVariables(vec) + # Write out mode shape as f5 file + modeName = baseName + "_%3.3d" % index + fileName = os.path.join(outputDir, modeName) + ".f5" + self.outputViewer.writeToFile(fileName) diff --git a/tacs/problems/modal.py b/tacs/problems/modal.py index a2541d9c3..a8bb77a2a 100644 --- a/tacs/problems/modal.py +++ b/tacs/problems/modal.py @@ -164,7 +164,8 @@ def _createVariables(self): restarts = self.getOption("nRestarts") self.gmres = tacs.TACS.KSM(self.K, self.pc, subspace, restarts) - eigTol = self.getOption("L2Convergence") + atol = self.getOption("L2Convergence") + rtol = self.getOption("L2ConvergenceRel") # Create the frequency analysis object self.freqSolver = tacs.TACS.FrequencyAnalysis( @@ -174,7 +175,9 @@ def _createVariables(self): self.K, self.gmres, num_eigs=self.numEigs, - eig_tol=eigTol, + eig_tol=atol, + eig_atol=atol, + eig_rtol=rtol, ) def _initializeFunctionList(self): @@ -369,8 +372,7 @@ def _updateAssemblerVars(self): def solve(self): """ - Solve the time integrated transient problem. The - forces must already be set. + Solve the eigenvalue problem. """ startTime = time.time() diff --git a/tacs/pytacs.py b/tacs/pytacs.py index ab84cf46a..0f8103125 100755 --- a/tacs/pytacs.py +++ b/tacs/pytacs.py @@ -1478,6 +1478,46 @@ def createModalProblem(self, name, sigma, numEigs, options=None): problem.setNodes(self.Xpts0) return problem + @postinitialize_method + def createBucklingProblem(self, name, sigma, numEigs, options=None): + """ + Create a new BucklingProblem for performing linearized buckling analysis. + This problem can be used to identify the buckling load factors and mode + shapes of the model through eigenvalue analysis. + + Parameters + ---------- + name : str + Name to assign problem. + sigma : float + Guess for the lowest eigenvalue. + This corresponds to the lowest expected buckling load factor. + numEigs : int + Number of eigenvalues to solve for. + options : dict + Problem-specific options to pass to ModalProblem instance (case-insensitive). + Defaults to None. + + Returns + ------- + problem : tacs.problems.BucklingProblem + BucklingProblem object used for performing buckling eigenvalue analysis. + """ + problem = tacs.problems.buckling.BucklingProblem( + name, + sigma, + numEigs, + self.assembler, + self.comm, + self.outputViewer, + self.meshLoader, + options, + ) + # Set with original design vars and coordinates, in case they have changed + problem.setDesignVars(self.x0) + problem.setNodes(self.Xpts0) + return problem + @postinitialize_method def createTACSProbsFromBDF(self): """ diff --git a/tests/element_tests/shell_tests/test_shell_element.py b/tests/element_tests/shell_tests/test_shell_element.py index 23e6cdd24..626b747ae 100644 --- a/tests/element_tests/shell_tests/test_shell_element.py +++ b/tests/element_tests/shell_tests/test_shell_element.py @@ -14,7 +14,7 @@ def setUp(self): # fd/cs step size if TACS.dtype is complex: self.dh = 1e-50 - self.rtol = 1e-10 + self.rtol = 1e-9 else: self.dh = 1e-5 self.rtol = 1e-2 diff --git a/tests/integration_tests/input_files/ss_plate.bdf b/tests/integration_tests/input_files/ss_plate.bdf new file mode 100644 index 000000000..28651a2a2 --- /dev/null +++ b/tests/integration_tests/input_files/ss_plate.bdf @@ -0,0 +1,200 @@ +INIT MASTER(S) +NASTRAN SYSTEM(442)=-1,SYSTEM(319)=1 +ID FEMAP,FEMAP +SOL SEBUCKL +CEND + TITLE = Simcenter Nastran Buckling Analysis Set + ECHO = NONE + DISPLACEMENT(PLOT) = ALL + SPCFORCE(PLOT) = ALL + STRESS(PLOT,CORNER) = ALL + SPC = 1 +SUBCASE 1 + LOAD = 1 +SUBCASE 2 + LABEL = BUCKLING CASE + METHOD = 1 +BEGIN BULK +$ *************************************************************************** +$ Written by : Femap +$ Version : 2023.1.0 +$ Translator : Simcenter Nastran +$ From Model : +$ Date : Fri May 19 14:25:00 2023 +$ Output To : C:\Users\trbrooks\AppData\Local\Temp\2\ +$ *************************************************************************** +$ +PARAM,PRGPST,NO +PARAM,POST,-1 +PARAM,OGEOM,NO +PARAM,AUTOSPC,YES +PARAM,GRDPNT,0 +EIGRL 1 3000. 1 0 MAX +CORD2C 1 0 0. 0. 0. 0. 0. 1.+FEMAPC1 ++FEMAPC1 1. 0. 1. +CORD2S 2 0 0. 0. 0. 0. 0. 1.+FEMAPC2 ++FEMAPC2 1. 0. 1. +$ Femap Load Set 1 : Untitled +FORCE 1 15 0 1. 0. -5000. 0. +FORCE 1 16 0 1. 0. -10000. 0. +FORCE 1 17 0 1. 0. -10000. 0. +FORCE 1 18 0 1. 0. -10000. 0. +FORCE 1 19 0 1. 0. -10000. 0. +FORCE 1 20 0 1. 0. -10000. 0. +FORCE 1 21 0 1. 0. -10000. 0. +FORCE 1 22 0 1. 0. -10000. 0. +FORCE 1 23 0 1. 0. -5000. 0. +FORCE 1 1 0 1. 0. 5000. 0. +FORCE 1 2 0 1. 0. 10000. 0. +FORCE 1 3 0 1. 0. 10000. 0. +FORCE 1 4 0 1. 0. 10000. 0. +FORCE 1 5 0 1. 0. 10000. 0. +FORCE 1 6 0 1. 0. 10000. 0. +FORCE 1 7 0 1. 0. 10000. 0. +FORCE 1 8 0 1. 0. 10000. 0. +FORCE 1 9 0 1. 0. 5000. 0. +$ Femap Constraint Set 1 : Untitled +SPC1 1 123456 1 +SPC1 1 236 2 +SPC1 1 236 3 +SPC1 1 236 4 +SPC1 1 236 5 +SPC1 1 236 6 +SPC1 1 236 7 +SPC1 1 236 8 +SPC1 1 236 9 +SPC1 1 36 10 +SPC1 1 36 11 +SPC1 1 36 12 +SPC1 1 36 13 +SPC1 1 36 14 +SPC1 1 36 15 +SPC1 1 36 16 +SPC1 1 36 17 +SPC1 1 36 18 +SPC1 1 36 19 +SPC1 1 36 20 +SPC1 1 36 21 +SPC1 1 36 22 +SPC1 1 36 23 +SPC1 1 36 24 +SPC1 1 36 25 +SPC1 1 36 26 +SPC1 1 36 27 +SPC1 1 36 28 +$ Femap Property 1 : PLATE Property +PSHELL 1 1 .02 1 1 0. +$ Femap Material 1 : ISOTROPIC Material +MAT1 1 2.05+11 .3 0. 0. 0. +GRID 1 0 0. 0. 0. 0 +GRID 2 0 .5 0. 0. 0 +GRID 3 0 1. 0. 0. 0 +GRID 4 0 1.5 0. 0. 0 +GRID 5 0 2. 0. 0. 0 +GRID 6 0 2.5 0. 0. 0 +GRID 7 0 3. 0. 0. 0 +GRID 8 0 3.5 0. 0. 0 +GRID 9 0 4. 0. 0. 0 +GRID 10 0 4. .5 0. 0 +GRID 11 0 4. 1. 0. 0 +GRID 12 0 4. 1.5 0. 0 +GRID 13 0 4. 2. 0. 0 +GRID 14 0 4. 2.5 0. 0 +GRID 15 0 4. 3. 0. 0 +GRID 16 0 3.5 3. 0. 0 +GRID 17 0 3. 3. 0. 0 +GRID 18 0 2.5 3. 0. 0 +GRID 19 0 2. 3. 0. 0 +GRID 20 0 1.5 3. 0. 0 +GRID 21 0 1. 3. 0. 0 +GRID 22 0 .5 3. 0. 0 +GRID 23 0 0. 3. 0. 0 +GRID 24 0 0. 2.5 0. 0 +GRID 25 0 0. 2. 0. 0 +GRID 26 0 0. 1.5 0. 0 +GRID 27 0 0. 1. 0. 0 +GRID 28 0 0. .5 0. 0 +GRID 29 0 .5 .5 0. 0 +GRID 30 0 1. .5 0. 0 +GRID 31 0 1.5 .5 0. 0 +GRID 32 0 2. .5 0. 0 +GRID 33 0 2.5 .5 0. 0 +GRID 34 0 3. .5 0. 0 +GRID 35 0 3.5 .5 0. 0 +GRID 36 0 .5 1. 0. 0 +GRID 37 0 1. 1. 0. 0 +GRID 38 0 1.5 1. 0. 0 +GRID 39 0 2. 1. 0. 0 +GRID 40 0 2.5 1. 0. 0 +GRID 41 0 3. 1. 0. 0 +GRID 42 0 3.5 1. 0. 0 +GRID 43 0 .5 1.5 0. 0 +GRID 44 0 1. 1.5 0. 0 +GRID 45 0 1.5 1.5 0. 0 +GRID 46 0 2. 1.5 0. 0 +GRID 47 0 2.5 1.5 0. 0 +GRID 48 0 3. 1.5 0. 0 +GRID 49 0 3.5 1.5 0. 0 +GRID 50 0 .5 2. 0. 0 +GRID 51 0 1. 2. 0. 0 +GRID 52 0 1.5 2. 0. 0 +GRID 53 0 2. 2. 0. 0 +GRID 54 0 2.5 2. 0. 0 +GRID 55 0 3. 2. 0. 0 +GRID 56 0 3.5 2. 0. 0 +GRID 57 0 .5 2.5 0. 0 +GRID 58 0 1. 2.5 0. 0 +GRID 59 0 1.5 2.5 0. 0 +GRID 60 0 2. 2.5 0. 0 +GRID 61 0 2.5 2.5 0. 0 +GRID 62 0 3. 2.5 0. 0 +GRID 63 0 3.5 2.5 0. 0 +CQUAD4 1 1 1 2 29 28 +CQUAD4 2 1 2 3 30 29 +CQUAD4 3 1 3 4 31 30 +CQUAD4 4 1 4 5 32 31 +CQUAD4 5 1 5 6 33 32 +CQUAD4 6 1 6 7 34 33 +CQUAD4 7 1 7 8 35 34 +CQUAD4 8 1 8 9 10 35 +CQUAD4 9 1 28 29 36 27 +CQUAD4 10 1 29 30 37 36 +CQUAD4 11 1 30 31 38 37 +CQUAD4 12 1 31 32 39 38 +CQUAD4 13 1 32 33 40 39 +CQUAD4 14 1 33 34 41 40 +CQUAD4 15 1 34 35 42 41 +CQUAD4 16 1 35 10 11 42 +CQUAD4 17 1 27 36 43 26 +CQUAD4 18 1 36 37 44 43 +CQUAD4 19 1 37 38 45 44 +CQUAD4 20 1 38 39 46 45 +CQUAD4 21 1 39 40 47 46 +CQUAD4 22 1 40 41 48 47 +CQUAD4 23 1 41 42 49 48 +CQUAD4 24 1 42 11 12 49 +CQUAD4 25 1 26 43 50 25 +CQUAD4 26 1 43 44 51 50 +CQUAD4 27 1 44 45 52 51 +CQUAD4 28 1 45 46 53 52 +CQUAD4 29 1 46 47 54 53 +CQUAD4 30 1 47 48 55 54 +CQUAD4 31 1 48 49 56 55 +CQUAD4 32 1 49 12 13 56 +CQUAD4 33 1 25 50 57 24 +CQUAD4 34 1 50 51 58 57 +CQUAD4 35 1 51 52 59 58 +CQUAD4 36 1 52 53 60 59 +CQUAD4 37 1 53 54 61 60 +CQUAD4 38 1 54 55 62 61 +CQUAD4 39 1 55 56 63 62 +CQUAD4 40 1 56 13 14 63 +CQUAD4 41 1 24 57 22 23 +CQUAD4 42 1 57 58 21 22 +CQUAD4 43 1 58 59 20 21 +CQUAD4 44 1 59 60 19 20 +CQUAD4 45 1 60 61 18 19 +CQUAD4 46 1 61 62 17 18 +CQUAD4 47 1 62 63 16 17 +CQUAD4 48 1 63 14 15 16 +ENDDATA diff --git a/tests/integration_tests/test_mphys_struct_buckling.py b/tests/integration_tests/test_mphys_struct_buckling.py new file mode 100644 index 000000000..77212df38 --- /dev/null +++ b/tests/integration_tests/test_mphys_struct_buckling.py @@ -0,0 +1,149 @@ +import os +import unittest + +import numpy as np +import openmdao.api as om +from mphys.multipoint import Multipoint +from mphys.scenario_structural import ScenarioStructural + +import tacs.mphys +from openmdao_analysis_base_test import OpenMDAOTestCase +from tacs import elements, constitutive, TACS + +""" +This is a simple 1m by 2m plate made up of four quad shell elements. +The plate is structurally loaded under a compression load and a unit force, +"f_struct", is applied on on every node. The mass and KSFailure of the plate +are evaluated as outputs and have their partial and total sensitivities checked. +""" + +base_dir = os.path.dirname(os.path.abspath(__file__)) +bdf_file = os.path.join(base_dir, "./input_files/debug_plate.bdf") + +# Historical reference values for function outputs +FUNC_REFS = { + "analysis.eigsb_0": -1.08864541, + "analysis.eigsb_1": 1.08938765, +} + +# Inputs to check total sensitivities wrt +wrt = ["mesh.fea_mesh.x_struct0", "dv_struct", "f_struct"] + + +class ProblemTest(OpenMDAOTestCase.OpenMDAOTest): + N_PROCS = 1 # this is how many MPI processes to use for this TestCase. + + def setup_problem(self, dtype): + """ + Setup openmdao problem object we will be testing. + """ + + # Overwrite default tolerances + if dtype == complex: + self.rtol = 1e-7 + self.dh = 1e-50 + else: + self.rtol = 1e-1 + self.dh = 1e-7 + + # Callback function used to setup TACS element objects and DVs + def element_callback( + dv_num, comp_id, comp_descript, elem_descripts, special_dvs, **kwargs + ): + rho = 2780.0 # density, kg/m^3 + E = 73.1e9 # elastic modulus, Pa + nu = 0.33 # poisson's ratio + ys = 324.0e6 # yield stress, Pa + thickness = 0.012 + min_thickness = 0.002 + max_thickness = 0.05 + + # Setup (isotropic) property and constitutive objects + prop = constitutive.MaterialProperties(rho=rho, E=E, nu=nu, ys=ys) + # Set one thickness dv for every component + con = constitutive.IsoShellConstitutive( + prop, t=thickness, tNum=dv_num, tlb=min_thickness, tub=max_thickness + ) + + # For each element type in this component, + # pass back the appropriate tacs element object + transform = None + elem = elements.Quad4Shell(transform, con) + + return elem + + def problem_setup(scenario_name, fea_assembler, problem): + """ + Helper function to add fixed forces and eval functions + to structural problems used in tacs builder + """ + # Set convergence to be tight for test + problem.setOption("L2Convergence", 1e-20) + problem.setOption("L2ConvergenceRel", 1e-20) + + def buckling_setup(scenario_name, fea_assembler): + """ + Helper function to add fixed forces and eval functions + to structural problems used in tacs builder + """ + problem = fea_assembler.createBucklingProblem( + "buckling", sigma=1e0, numEigs=2 + ) + problem.setOption("L2Convergence", 1e-20) + problem.setOption("L2ConvergenceRel", 1e-20) + return problem + + class Top(Multipoint): + def setup(self): + tacs_options = { + "element_callback": element_callback, + "problem_setup": problem_setup, + "buckling_setup": buckling_setup, + "mesh_file": bdf_file, + } + + tacs_builder = tacs.mphys.TacsBuilder( + tacs_options, + check_partials=True, + coupled=True, + write_solution=False, + ) + tacs_builder.initialize(self.comm) + ndv_struct = tacs_builder.get_ndv() + + dvs = self.add_subsystem("dvs", om.IndepVarComp(), promotes=["*"]) + dvs.add_output("dv_struct", np.array(ndv_struct * [0.01])) + + f_size = tacs_builder.get_ndof() * tacs_builder.get_number_of_nodes() + forces = self.add_subsystem("forces", om.IndepVarComp(), promotes=["*"]) + f = np.zeros(f_size) + f[1::3] = -1e0 + forces.add_output("f_struct", val=np.ones(f_size), distributed=True) + + self.add_subsystem("mesh", tacs_builder.get_mesh_coordinate_subsystem()) + self.mphys_add_scenario( + "analysis", ScenarioStructural(struct_builder=tacs_builder) + ) + self.connect("mesh.x_struct0", "analysis.x_struct0") + self.connect("dv_struct", "analysis.dv_struct") + self.connect("f_struct", "analysis.f_struct") + + prob = om.Problem() + prob.model = Top() + + return prob + + def setup_funcs(self): + """ + Create a dict of functions to be tested and a list of inputs + to test their sensitivities with respect to. + """ + return FUNC_REFS, wrt + + # This test is very difficult to verify with FD, so we only run it w/ CS + @unittest.skipIf(TACS.dtype != complex, "Skipping test in real mode.") + def test_partials(self): + """ + Test partial sensitivities using fd/cs + """ + return OpenMDAOTestCase.OpenMDAOTest.test_partials(self) diff --git a/tests/integration_tests/test_shell_plate_buckling.py b/tests/integration_tests/test_shell_plate_buckling.py new file mode 100644 index 000000000..df759a9da --- /dev/null +++ b/tests/integration_tests/test_shell_plate_buckling.py @@ -0,0 +1,75 @@ +import os + +from pytacs_analysis_base_test import PyTACSTestCase +from tacs import pytacs, elements, constitutive + +"""" +The nominal case is a 4m x 3m flat plate under a buckling analysis. The +perimeter of the plate is pinned and loaded in compression on its horizontal edges. +This tests the eigenvalues and eigenvalue sensitivities +""" + +base_dir = os.path.dirname(os.path.abspath(__file__)) +bdf_file = os.path.join(base_dir, "./input_files/ss_plate.bdf") + + +class ProblemTest(PyTACSTestCase.PyTACSTest): + N_PROCS = 2 # this is how many MPI processes to use for this TestCase. + + FUNC_REFS = { + "buckling_eigsb.0": 20.18364503170502, + "buckling_eigsb.1": 48.93915410646456, + "buckling_eigsb.2": 82.32683532479277, + "buckling_eigsb.3": 86.45944604770253, + "buckling_eigsb.4": 126.77071579688763, + } + + def setup_tacs_problems(self, comm): + """ + Setup pytacs object for problems we will be testing. + """ + + # Overwrite default check values + if self.dtype == complex: + self.rtol = 1e-8 + self.atol = 1e-8 + self.dh = 1e-50 + else: + self.rtol = 2e-1 + self.atol = 1e-4 + self.dh = 1e-5 + + # Instantiate FEA Assembler + fea_assembler = pytacs.pyTACS(bdf_file, comm) + + def elem_call_back( + dv_num, comp_id, comp_descript, elem_descripts, global_dvs, **kwargs + ): + # Material properties + rho = 2500.0 # density kg/m^3 + E = 205e9 # Young's modulus (Pa) + nu = 0.3 # Poisson's ratio + ys = 464.0e6 # yield stress + + # Plate geometry + tplate = 0.020 # 20 mm + + # Set up property model + prop = constitutive.MaterialProperties(rho=rho, E=E, nu=nu, ys=ys) + # Set up constitutive model + con = constitutive.IsoShellConstitutive(prop, t=tplate, tNum=dv_num) + transform = None + # Set up element + elem = elements.Quad4Shell(transform, con) + scale = [100.0] + return elem, scale + + # Set up constitutive objects and elements + fea_assembler.initialize(elem_call_back) + + buckle_prob = fea_assembler.createBucklingProblem("buckling", 10.0, 10) + buckle_prob.setOption("L2Convergence", 1e-20) + buckle_prob.setOption("L2ConvergenceRel", 1e-20) + buckle_prob.addLoadFromBDF(loadID=1) + + return [buckle_prob], fea_assembler diff --git a/tests/integration_tests/test_shell_plate_quad.py b/tests/integration_tests/test_shell_plate_quad.py index e501a6fee..f94047a32 100644 --- a/tests/integration_tests/test_shell_plate_quad.py +++ b/tests/integration_tests/test_shell_plate_quad.py @@ -191,7 +191,7 @@ def elem_call_back( ) # Add modal problem - mp = fea_assembler.createModalProblem("modal", 7e4, 10) + mp = fea_assembler.createModalProblem("modal", 2e5, 10) mp.setOption("printLevel", 2) tacs_probs.append(mp)