diff --git a/Thermo_electro_elasticity/CMakeLists.txt b/Thermo_electro_elasticity/CMakeLists.txt new file mode 100755 index 00000000..fc49ead4 --- /dev/null +++ b/Thermo_electro_elasticity/CMakeLists.txt @@ -0,0 +1,72 @@ +## +# CMake script for the step-8 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "thermo-electro-elastic_pump") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) + +FIND_PACKAGE(deal.II 9.1.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() + +# +# Are all dependencies fulfilled? +# +IF(NOT DEAL_II_WITH_MPI OR NOT DEAL_II_WITH_P4EST OR NOT DEAL_II_WITH_TRILINOS) + MESSAGE(FATAL_ERROR " +Error! The deal.II library found at ${DEAL_II_PATH} was not configured with + DEAL_II_WITH_MPI = ON + DEAL_II_WITH_P4EST = ON + DEAL_II_WITH_TRILINOS = ON +One or all of these are OFF in your installation but are required for this tutorial step." + ) +ENDIF() + +ADD_CUSTOM_TARGET(clean_output + COMMAND ${CMAKE_COMMAND} -E remove *.vtu *.pvtu *.pvd *.vtk + COMMENT "clean_output invoked" +) + +ADD_CUSTOM_TARGET(run_mpi_2 + COMMAND mpirun -np 2 ${TARGET} + COMMENT "run_mpi_2 invoked" +) + +ADD_CUSTOM_TARGET(run_mpi_4 + COMMAND mpirun -np 4 ${TARGET} + COMMENT "run_mpi_4 invoked" +) + +ADD_CUSTOM_TARGET(run_mpi_8 + COMMAND mpirun -np 8 ${TARGET} + COMMENT "run_mpi_8 invoked" +) diff --git a/Thermo_electro_elasticity/Pump_um_coarse.inp b/Thermo_electro_elasticity/Pump_um_coarse.inp new file mode 100644 index 00000000..b8a7237f --- /dev/null +++ b/Thermo_electro_elasticity/Pump_um_coarse.inp @@ -0,0 +1,1891 @@ +*Heading +** Job name: Pump_um_coarse Model name: Model-1 +** Generated by: Abaqus/CAE 6.13-3 +*Preprint, echo=NO, model=NO, history=NO, contact=NO +** +** PARTS +** +*Part, name=Part-5 +*Node + 1, 490., 0., 70. + 2, 490., 0., 50. + 3, 540., 0., 50. + 4, 540., 0., 70. + 5, 490., 50., 0. + 6, 490., 70., 0. + 7, 540., 70., 0. + 8, 540., 50., 0. + 9, 490., 0., 80. + 10, 490., 99.4987411, 0. + 11, 490., 99.4987411, 80. + 12, 469.57428, 99.4987411, 0. + 13, 469.57428, 99.4987411, 80. + 14, 480., 0., 80. + 15, 480., 0., 50. + 16, 477.388733, 50., 0. + 17, 0., 0., 80. + 18, 0., 480., 80. + 19, 0., 500., 80. + 20, 0., 0., 100. + 21, 0., 500., 100. + 22, 490., 99.4987411, 100. + 23, 490., 0., 100. + 24, 0., 500., 120. + 25, 0., 0., 120. + 26, 490., 0., 120. + 27, 490., 99.4987411, 120. + 28, 0., 480., 0. + 29, 0., 500., 0. + 30, 515., 0., 50. + 31, 515., 0., 70. + 32, 540., 46.1939774, 19.1341724 + 33, 540., 35.3553391, 35.3553391 + 34, 540., 19.1341724, 46.1939774 + 35, 490., 19.9819832, 45.8336143 + 36, 490., 36.1765594, 34.5145836 + 37, 490., 46.45504, 18.4913311 + 38, 490., 29.1102028, 63.6600037 + 39, 490., 51.6898003, 47.2034378 + 40, 490., 65.3111115, 25.1884727 + 41, 540., 64.6715698, 26.7878399 + 42, 540., 49.4974747, 49.4974747 + 43, 540., 26.7878399, 64.6715698 + 44, 515., 70., 0. + 45, 515., 50., 0. + 46, 490., 69.590332, 80. + 47, 490., 38.4070473, 80. + 48, 490., 0., 74.4686127 + 49, 490., 83.4945602, 0. + 50, 490., 99.4987411, 27.2618198 + 51, 490., 99.4987411, 54.0209007 + 52, 469.57428, 99.4987411, 27.2618198 + 53, 469.57428, 99.4987411, 54.0209007 + 54, 475.309845, 66.9369583, 80. + 55, 478.651611, 35.9530525, 80. + 56, 480., 0., 69.8448181 + 57, 480., 0., 60. + 58, 479.564056, 20.4530792, 45.6253395 + 59, 478.604462, 36.5755119, 34.0915222 + 60, 477.737518, 46.5497169, 18.2516804 + 61, 475.175751, 67.8822479, 0. + 62, 473.179718, 80.6284409, 0. + 63, 0., 30.1678772, 80. + 64, 0., 60.3421745, 80. + 65, 0., 90.5150452, 80. + 66, 0., 120.679443, 80. + 67, 0., 150.828003, 80. + 68, 0., 180.952621, 80. + 69, 0., 211.044479, 80. + 70, 0., 241.093918, 80. + 71, 0., 271.090668, 80. + 72, 0., 301.005829, 80. + 73, 0., 330.824829, 80. + 74, 0., 360.565765, 80. + 75, 0., 390.298248, 80. + 76, 0., 419.836212, 80. + 77, 0., 449.562958, 80. + 78, 28.7764206, 479.136627, 80. + 79, 57.5026703, 476.543213, 80. + 80, 86.0282516, 472.227844, 80. + 81, 114.222107, 466.21167, 80. + 82, 141.967514, 458.525055, 80. + 83, 169.157501, 449.205688, 80. + 84, 195.692413, 438.297241, 80. + 85, 221.478134, 425.849091, 80. + 86, 246.425385, 411.91568, 80. + 87, 270.449432, 396.556549, 80. + 88, 293.408508, 379.883484, 80. + 89, 315.142822, 362.05661, 80. + 90, 335.586884, 343.192993, 80. + 91, 354.918976, 323.160187, 80. + 92, 373.17923, 301.889496, 80. + 93, 390.288452, 279.418915, 80. + 94, 406.115387, 255.871658, 80. + 95, 420.563751, 231.357162, 80. + 96, 433.562225, 205.970383, 80. + 97, 445.04422, 179.821167, 80. + 98, 454.941071, 153.064117, 80. + 99, 463.192505, 125.907555, 80. + 100, 449.595734, 0., 80. + 101, 419.83667, 0., 80. + 102, 390.395508, 0., 80. + 103, 360.947052, 0., 80. + 104, 331.330811, 0., 80. + 105, 301.546997, 0., 80. + 106, 271.632629, 0., 80. + 107, 241.62236, 0., 80. + 108, 211.541504, 0., 80. + 109, 181.4077, 0., 80. + 110, 151.233871, 0., 80. + 111, 121.030128, 0., 80. + 112, 90.8028488, 0., 80. + 113, 60.5537109, 0., 80. + 114, 30.2841644, 0., 80. + 115, 482.865601, 129.772141, 80. + 116, 474.158966, 158.660889, 80. + 117, 463.843109, 186.680389, 80. + 118, 451.91861, 213.938263, 80. + 119, 438.411163, 240.407242, 80. + 120, 423.361603, 266.016815, 80. + 121, 406.820221, 290.684204, 80. + 122, 388.844299, 314.32486, 80. + 123, 369.496796, 336.856232, 80. + 124, 348.845856, 358.199066, 80. + 125, 326.964264, 378.278168, 80. + 126, 303.929199, 397.022705, 80. + 127, 279.821899, 414.366638, 80. + 128, 254.727432, 430.248688, 80. + 129, 228.734299, 444.612885, 80. + 130, 201.934189, 457.408539, 80. + 131, 174.421707, 468.590515, 80. + 132, 146.293869, 478.119354, 80. + 133, 117.649979, 485.961395, 80. + 134, 88.5909805, 492.08905, 80. + 135, 59.2194176, 496.480682, 80. + 136, 29.6389694, 499.120758, 80. + 137, 0., 30.1678772, 100. + 138, 0., 60.3421745, 100. + 139, 0., 90.5150452, 100. + 140, 0., 120.679443, 100. + 141, 0., 150.828003, 100. + 142, 0., 180.952621, 100. + 143, 0., 211.044479, 100. + 144, 0., 241.093918, 100. + 145, 0., 271.090668, 100. + 146, 0., 301.005829, 100. + 147, 0., 330.824829, 100. + 148, 0., 360.565765, 100. + 149, 0., 390.298248, 100. + 150, 0., 419.836212, 100. + 151, 0., 449.562958, 100. + 152, 0., 480., 100. + 153, 29.6389694, 499.120758, 100. + 154, 59.2194176, 496.480682, 100. + 155, 88.5909805, 492.08905, 100. + 156, 117.649979, 485.961395, 100. + 157, 146.293869, 478.119354, 100. + 158, 174.421707, 468.590515, 100. + 159, 201.934189, 457.408539, 100. + 160, 228.734299, 444.612885, 100. + 161, 254.727432, 430.248688, 100. + 162, 279.821899, 414.366638, 100. + 163, 303.929199, 397.022705, 100. + 164, 326.964264, 378.278168, 100. + 165, 348.845856, 358.199066, 100. + 166, 369.496796, 336.856232, 100. + 167, 388.844299, 314.32486, 100. + 168, 406.820221, 290.684204, 100. + 169, 423.361603, 266.016815, 100. + 170, 438.411163, 240.407242, 100. + 171, 451.91861, 213.938263, 100. + 172, 463.843109, 186.680389, 100. + 173, 474.158966, 158.660889, 100. + 174, 482.865601, 129.772141, 100. + 175, 490., 69.590332, 100. + 176, 490., 38.4070473, 100. + 177, 480., 0., 100. + 178, 449.595734, 0., 100. + 179, 419.83667, 0., 100. + 180, 390.395508, 0., 100. + 181, 360.947052, 0., 100. + 182, 331.330811, 0., 100. + 183, 301.546997, 0., 100. + 184, 271.632629, 0., 100. + 185, 241.62236, 0., 100. + 186, 211.541504, 0., 100. + 187, 181.4077, 0., 100. + 188, 151.233871, 0., 100. + 189, 121.030128, 0., 100. + 190, 90.8028488, 0., 100. + 191, 60.5537109, 0., 100. + 192, 30.2841644, 0., 100. + 193, 0., 480., 120. + 194, 0., 449.562958, 120. + 195, 0., 419.836212, 120. + 196, 0., 390.298248, 120. + 197, 0., 360.565765, 120. + 198, 0., 330.824829, 120. + 199, 0., 301.005829, 120. + 200, 0., 271.090668, 120. + 201, 0., 241.093918, 120. + 202, 0., 211.044479, 120. + 203, 0., 180.952621, 120. + 204, 0., 150.828003, 120. + 205, 0., 120.679443, 120. + 206, 0., 90.5150452, 120. + 207, 0., 60.3421745, 120. + 208, 0., 30.1678772, 120. + 209, 30.2841644, 0., 120. + 210, 60.5537109, 0., 120. + 211, 90.8028488, 0., 120. + 212, 121.030128, 0., 120. + 213, 151.233871, 0., 120. + 214, 181.4077, 0., 120. + 215, 211.541504, 0., 120. + 216, 241.62236, 0., 120. + 217, 271.632629, 0., 120. + 218, 301.546997, 0., 120. + 219, 331.330811, 0., 120. + 220, 360.947052, 0., 120. + 221, 390.395508, 0., 120. + 222, 419.83667, 0., 120. + 223, 449.595734, 0., 120. + 224, 480., 0., 120. + 225, 490., 38.4070473, 120. + 226, 490., 69.590332, 120. + 227, 482.865601, 129.772141, 120. + 228, 474.158966, 158.660889, 120. + 229, 463.843109, 186.680389, 120. + 230, 451.91861, 213.938263, 120. + 231, 438.411163, 240.407242, 120. + 232, 423.361603, 266.016815, 120. + 233, 406.820221, 290.684204, 120. + 234, 388.844299, 314.32486, 120. + 235, 369.496796, 336.856232, 120. + 236, 348.845856, 358.199066, 120. + 237, 326.964264, 378.278168, 120. + 238, 303.929199, 397.022705, 120. + 239, 279.821899, 414.366638, 120. + 240, 254.727432, 430.248688, 120. + 241, 228.734299, 444.612885, 120. + 242, 201.934189, 457.408539, 120. + 243, 174.421707, 468.590515, 120. + 244, 146.293869, 478.119354, 120. + 245, 117.649979, 485.961395, 120. + 246, 88.5909805, 492.08905, 120. + 247, 59.2194176, 496.480682, 120. + 248, 29.6389694, 499.120758, 120. + 249, 0., 500., 26.666666 + 250, 0., 500., 53.3333321 + 251, 0., 480., 26.666666 + 252, 0., 480., 53.3333321 + 253, 463.192505, 125.907555, 0. + 254, 454.941071, 153.064117, 0. + 255, 445.04422, 179.821167, 0. + 256, 433.562225, 205.970383, 0. + 257, 420.563751, 231.357162, 0. + 258, 406.115387, 255.871658, 0. + 259, 390.288452, 279.418915, 0. + 260, 373.17923, 301.889496, 0. + 261, 354.918976, 323.160187, 0. + 262, 335.586884, 343.192993, 0. + 263, 315.142822, 362.05661, 0. + 264, 293.408508, 379.883484, 0. + 265, 270.449432, 396.556549, 0. + 266, 246.425385, 411.91568, 0. + 267, 221.478134, 425.849091, 0. + 268, 195.692413, 438.297241, 0. + 269, 169.157501, 449.205688, 0. + 270, 141.967514, 458.525055, 0. + 271, 114.222107, 466.21167, 0. + 272, 86.0282516, 472.227844, 0. + 273, 57.5026703, 476.543213, 0. + 274, 28.7764206, 479.136627, 0. + 275, 29.6389694, 499.120758, 0. + 276, 59.2194176, 496.480682, 0. + 277, 88.5909805, 492.08905, 0. + 278, 117.649979, 485.961395, 0. + 279, 146.293869, 478.119354, 0. + 280, 174.421707, 468.590515, 0. + 281, 201.934189, 457.408539, 0. + 282, 228.734299, 444.612885, 0. + 283, 254.727432, 430.248688, 0. + 284, 279.821899, 414.366638, 0. + 285, 303.929199, 397.022705, 0. + 286, 326.964264, 378.278168, 0. + 287, 348.845856, 358.199066, 0. + 288, 369.496796, 336.856232, 0. + 289, 388.844299, 314.32486, 0. + 290, 406.820221, 290.684204, 0. + 291, 423.361603, 266.016815, 0. + 292, 438.411163, 240.407242, 0. + 293, 451.91861, 213.938263, 0. + 294, 463.843109, 186.680389, 0. + 295, 474.158966, 158.660889, 0. + 296, 482.865601, 129.772141, 0. + 297, 515., 19.1341724, 46.1939774 + 298, 515., 35.3553391, 35.3553391 + 299, 515., 46.1939774, 19.1341724 + 300, 515., 26.7878399, 64.6715698 + 301, 515., 49.4974747, 49.4974747 + 302, 515., 64.6715698, 26.7878399 + 303, 490., 37.6345978, 69.6669769 + 304, 490., 70.8865662, 56.2451134 + 305, 490., 81.0312805, 28.8407288 + 306, 479.093292, 29.4890614, 59.4982147 + 307, 477.230713, 51.4863052, 46.319046 + 308, 475.686493, 64.2055435, 24.6511402 + 309, 478.602783, 36.5973854, 68.2451859 + 310, 475.064911, 68.6537094, 56.6212616 + 311, 473.530762, 78.5405807, 29.0001736 + 312, 308.408386, 310.652832, 80. + 313, 26.3173389, 389.829346, 80. + 314, 52.6818466, 388.653961, 80. + 315, 78.9932709, 386.550964, 80. + 316, 105.167198, 383.335907, 80. + 317, 131.121841, 378.891907, 80. + 318, 156.775497, 373.167786, 80. + 319, 182.044037, 366.15274, 80. + 320, 206.836548, 357.865356, 80. + 321, 231.052521, 348.379974, 80. + 322, 254.458435, 337.954651, 80. + 323, 276.378052, 327.271606, 80. + 324, 295.395996, 317.642548, 80. + 325, 315.614349, 297.691711, 80. + 326, 325.730713, 279.077209, 80. + 327, 336.847534, 257.60675, 80. + 328, 347.654297, 234.562912, 80. + 329, 357.458038, 210.555557, 80. + 330, 365.999237, 185.864929, 80. + 331, 373.196716, 160.638474, 80. + 332, 379.032166, 134.960037, 80. + 333, 383.522583, 108.856888, 80. + 334, 386.735046, 82.2811966, 80. + 335, 388.803833, 55.1683884, 80. + 336, 389.893951, 27.6878395, 80. + 337, 324.677673, 329.254822, 80. + 338, 314.936401, 318.413727, 80. + 339, 360.740326, 27.1858025, 80. + 340, 331.282654, 27.3117256, 80. + 341, 301.589142, 27.7472534, 80. + 342, 271.725739, 28.2847176, 80. + 343, 241.741211, 28.8117523, 80. + 344, 211.668472, 29.2747231, 80. + 345, 181.529816, 29.650486, 80. + 346, 151.3414, 29.9301357, 80. + 347, 121.11602, 30.1112804, 80. + 348, 90.8629074, 30.1979103, 80. + 349, 60.5878525, 30.2102776, 80. + 350, 30.2969837, 30.1900826, 80. + 351, 360.044556, 54.3420868, 80. + 352, 330.861023, 54.6320419, 80. + 353, 301.350372, 55.4945908, 80. + 354, 271.607513, 56.5540733, 80. + 355, 241.700684, 57.5973282, 80. + 356, 211.674927, 58.51791, 80. + 357, 181.560577, 59.2677231, 80. + 358, 151.379471, 59.8273239, 80. + 359, 121.14904, 60.1912689, 80. + 360, 90.8833008, 60.369194, 80. + 361, 60.5941391, 60.4040298, 80. + 362, 30.2944794, 60.3757477, 80. + 363, 358.537354, 81.3123169, 80. + 364, 329.752258, 81.8564682, 80. + 365, 300.526794, 83.1806488, 80. + 366, 270.991333, 84.7731705, 80. + 367, 241.23761, 86.3363419, 80. + 368, 211.32576, 87.7171021, 80. + 369, 181.296173, 88.8440933, 80. + 370, 151.178055, 89.6874771, 80. + 371, 120.994743, 90.2388306, 80. + 372, 90.7659607, 90.5144196, 80. + 373, 60.5101814, 90.5805054, 80. + 374, 30.247448, 90.5526047, 80. + 375, 356.004425, 107.98999, 80. + 376, 327.733093, 108.90448, 80. + 377, 298.899048, 110.745125, 80. + 378, 269.668915, 112.897171, 80. + 379, 240.159576, 114.994896, 80. + 380, 210.447906, 116.84613, 80. + 381, 180.585526, 118.359322, 80. + 382, 150.610001, 119.495255, 80. + 383, 120.551254, 120.242989, 80. + 384, 90.4351273, 120.626007, 80. + 385, 60.2864685, 120.733536, 80. + 386, 30.1319237, 120.714272, 80. + 387, 352.309296, 134.344604, 80. + 388, 324.663239, 135.721466, 80. + 389, 296.328308, 138.132324, 80. + 390, 267.507446, 140.876312, 80. + 391, 238.343216, 143.530441, 80. + 392, 208.929749, 145.869308, 80. + 393, 179.330536, 147.783951, 80. + 394, 149.591949, 149.22702, 80. + 395, 119.7509, 150.185593, 80. + 396, 89.8394928, 150.690521, 80. + 397, 59.8885803, 150.853058, 80. + 398, 29.9307117, 150.852463, 80. + 399, 347.375854, 160.354553, 80. + 400, 320.461884, 162.258362, 80. + 401, 292.733185, 165.28627, 80. + 402, 264.428223, 168.65657, 80. + 403, 235.714813, 171.894257, 80. + 404, 206.704163, 174.744217, 80. + 405, 177.471802, 177.082062, 80. + 406, 148.072998, 178.853287, 80. + 407, 118.551918, 180.043442, 80. + 408, 88.9468689, 180.690521, 80. + 409, 59.2944489, 180.926239, 80. + 410, 29.6324768, 180.957291, 80. + 411, 341.169617, 185.971741, 80. + 412, 315.086395, 188.455673, 80. + 413, 288.068298, 192.144928, 80. + 414, 260.386414, 196.178635, 80. + 415, 232.231873, 200.032379, 80. + 416, 203.732239, 203.423309, 80. + 417, 174.974686, 206.212952, 80. + 418, 146.023346, 208.340302, 80. + 419, 116.929596, 209.789658, 80. + 420, 87.7379913, 210.605667, 80. + 421, 58.4906311, 210.938232, 80. + 422, 29.2301559, 211.017517, 80. + 423, 333.686584, 211.116821, 80. + 424, 308.520233, 214.237442, 80. + 425, 282.311462, 218.636719, 80. + 426, 255.358063, 223.376648, 80. + 427, 227.870941, 227.885544, 80. + 428, 199.992142, 231.854523, 80. + 429, 171.819626, 235.131927, 80. + 430, 143.426025, 237.650848, 80. + 431, 114.869835, 239.394516, 80. + 432, 86.201767, 240.413391, 80. + 433, 57.4692993, 240.872528, 80. + 434, 28.7195873, 241.020782, 80. + 435, 324.99762, 235.660843, 80. + 436, 300.833466, 239.486359, 80. + 437, 275.531067, 244.653351, 80. + 438, 249.407806, 250.15097, 80. + 439, 222.69162, 255.363373, 80. + 440, 195.537811, 259.957214, 80. + 441, 168.054626, 263.768616, 80. + 442, 140.322937, 266.725098, 80. + 443, 112.408203, 268.808777, 80. + 444, 84.3672791, 270.074707, 80. + 445, 56.2529373, 270.699554, 80. + 446, 28.1167679, 270.945374, 80. + 447, 315.477386, 259.325012, 80. + 448, 292.469055, 263.937317, 80. + 449, 268.179993, 269.951599, 80. + 450, 242.966888, 276.282837, 80. + 451, 217.090057, 282.270782, 80. + 452, 190.724731, 287.559174, 80. + 453, 163.991791, 291.973114, 80. + 454, 136.981018, 295.435486, 80. + 455, 109.76503, 297.927429, 80. + 456, 82.4069519, 299.507507, 80. + 457, 54.9653358, 300.361877, 80. + 458, 27.4883575, 300.754333, 80. + 459, 306.185425, 281.512512, 80. + 460, 284.595001, 287.069275, 80. + 461, 261.398132, 294.091431, 80. + 462, 237.087585, 301.402466, 80. + 463, 212.013489, 308.294739, 80. + 464, 186.389832, 314.394531, 80. + 465, 160.35466, 319.522034, 80. + 466, 134.007492, 323.599274, 80. + 467, 107.428268, 326.608826, 80. + 468, 80.6871338, 328.616272, 80. + 469, 53.8459587, 329.808197, 80. + 470, 26.9473057, 330.424438, 80. + 471, 298.903351, 301.300476, 80. + 472, 278.958435, 308.21759, 80. + 473, 256.76825, 316.583801, 80. + 474, 233.173798, 325.140198, 80. + 475, 208.694977, 333.149658, 80. + 476, 183.605652, 340.239746, 80. + 477, 158.05777, 346.242981, 80. + 478, 132.157104, 351.091614, 80. + 479, 105.989746, 354.778778, 80. + 480, 79.6325531, 357.374084, 80. + 481, 53.1538162, 359.036652, 80. + 482, 26.6046143, 359.96405, 80. + 483, 26.7750816, 419.111877, 80. + 484, 27.7743015, 448.630798, 80. + 485, 53.5985489, 417.434937, 80. + 486, 55.5699883, 446.432983, 80. + 487, 80.3450623, 414.615295, 80. + 488, 83.2367935, 442.821136, 80. + 489, 106.907867, 410.521088, 80. + 490, 110.643425, 437.725067, 80. + 491, 133.190231, 405.075562, 80. + 492, 137.674316, 431.118195, 80. + 493, 159.100616, 398.245819, 80. + 494, 164.224823, 423.002411, 80. + 495, 184.549286, 390.027802, 80. + 496, 190.197495, 413.398285, 80. + 497, 209.445496, 380.43866, 80. + 498, 215.499146, 402.339966, 80. + 499, 233.697021, 369.529541, 80. + 500, 240.038361, 389.877319, 80. + 501, 257.088409, 357.470459, 80. + 502, 263.649811, 376.122467, 80. + 503, 279.085999, 344.675598, 80. + 504, 285.982422, 361.318726, 80. + 505, 298.797424, 331.722321, 80. + 506, 306.562988, 345.727051, 80. + 507, 448.509613, 30.5660229, 80. + 508, 418.962402, 28.8342876, 80. + 509, 445.937592, 60.9904785, 80. + 510, 417.099548, 57.4955788, 80. + 511, 441.769073, 90.8010025, 80. + 512, 414.03241, 85.6092377, 80. + 513, 436.391846, 118.202293, 80. + 514, 409.730835, 112.684601, 80. + 515, 429.53244, 144.8508, 80. + 516, 404.071167, 139.068329, 80. + 517, 421.133728, 171.011398, 80. + 518, 396.989044, 164.923126, 80. + 519, 411.21405, 196.618896, 80. + 520, 388.472626, 190.246719, 80. + 521, 399.81778, 221.551819, 80. + 522, 378.544495, 214.959793, 80. + 523, 387.000488, 245.688309, 80. + 524, 367.267426, 238.929092, 80. + 525, 372.842682, 268.875824, 80. + 526, 354.831757, 261.911407, 80. + 527, 357.588928, 290.834473, 80. + 528, 341.694458, 283.422516, 80. + 529, 341.607422, 311.109253, 80. + 530, 328.463165, 302.659302, 80. + 531, 469.57428, 99.4987411, 100. + 532, 475.348907, 66.6589813, 100. + 533, 478.709412, 35.1753426, 100. + 534, 28.7764206, 479.136627, 100. + 535, 57.5026703, 476.543213, 100. + 536, 86.0282516, 472.227844, 100. + 537, 114.222107, 466.21167, 100. + 538, 141.967514, 458.525055, 100. + 539, 169.157501, 449.205688, 100. + 540, 195.692413, 438.297241, 100. + 541, 221.478134, 425.849091, 100. + 542, 246.425385, 411.91568, 100. + 543, 270.449432, 396.556549, 100. + 544, 293.408508, 379.883484, 100. + 545, 315.142822, 362.05661, 100. + 546, 335.586884, 343.192993, 100. + 547, 354.918976, 323.160187, 100. + 548, 373.17923, 301.889496, 100. + 549, 390.288452, 279.418915, 100. + 550, 406.115387, 255.871658, 100. + 551, 420.563751, 231.357162, 100. + 552, 433.562225, 205.970383, 100. + 553, 445.04422, 179.821167, 100. + 554, 454.941071, 153.064117, 100. + 555, 463.192505, 125.907555, 100. + 556, 308.408386, 310.652832, 100. + 557, 26.3173389, 389.829346, 100. + 558, 52.6818466, 388.653961, 100. + 559, 78.9932709, 386.550964, 100. + 560, 105.167198, 383.335907, 100. + 561, 131.121841, 378.891907, 100. + 562, 156.775497, 373.167786, 100. + 563, 182.044037, 366.15274, 100. + 564, 206.836548, 357.865356, 100. + 565, 231.052521, 348.379974, 100. + 566, 254.458435, 337.954651, 100. + 567, 276.378052, 327.271606, 100. + 568, 295.395996, 317.642548, 100. + 569, 315.614349, 297.691711, 100. + 570, 325.730713, 279.077209, 100. + 571, 336.847534, 257.60675, 100. + 572, 347.654297, 234.562912, 100. + 573, 357.458038, 210.555557, 100. + 574, 365.999237, 185.864929, 100. + 575, 373.196716, 160.638474, 100. + 576, 379.032166, 134.960037, 100. + 577, 383.522583, 108.856888, 100. + 578, 386.735046, 82.2811966, 100. + 579, 388.803833, 55.1683884, 100. + 580, 389.893951, 27.6878395, 100. + 581, 324.677673, 329.254822, 100. + 582, 314.936401, 318.413727, 100. + 583, 360.740326, 27.1858025, 100. + 584, 331.282654, 27.3117256, 100. + 585, 301.589142, 27.7472534, 100. + 586, 271.725739, 28.2847176, 100. + 587, 241.741211, 28.8117523, 100. + 588, 211.668472, 29.2747231, 100. + 589, 181.529816, 29.650486, 100. + 590, 151.3414, 29.9301357, 100. + 591, 121.11602, 30.1112804, 100. + 592, 90.8629074, 30.1979103, 100. + 593, 60.5878525, 30.2102776, 100. + 594, 30.2969837, 30.1900826, 100. + 595, 360.044556, 54.3420868, 100. + 596, 330.861023, 54.6320419, 100. + 597, 301.350372, 55.4945908, 100. + 598, 271.607513, 56.5540733, 100. + 599, 241.700684, 57.5973282, 100. + 600, 211.674927, 58.51791, 100. + 601, 181.560577, 59.2677231, 100. + 602, 151.379471, 59.8273239, 100. + 603, 121.14904, 60.1912689, 100. + 604, 90.8833008, 60.369194, 100. + 605, 60.5941391, 60.4040298, 100. + 606, 30.2944794, 60.3757477, 100. + 607, 358.537354, 81.3123169, 100. + 608, 329.752258, 81.8564682, 100. + 609, 300.526794, 83.1806488, 100. + 610, 270.991333, 84.7731705, 100. + 611, 241.23761, 86.3363419, 100. + 612, 211.32576, 87.7171021, 100. + 613, 181.296173, 88.8440933, 100. + 614, 151.178055, 89.6874771, 100. + 615, 120.994743, 90.2388306, 100. + 616, 90.7659607, 90.5144196, 100. + 617, 60.5101814, 90.5805054, 100. + 618, 30.247448, 90.5526047, 100. + 619, 356.004425, 107.98999, 100. + 620, 327.733093, 108.90448, 100. + 621, 298.899048, 110.745125, 100. + 622, 269.668915, 112.897171, 100. + 623, 240.159576, 114.994896, 100. + 624, 210.447906, 116.84613, 100. + 625, 180.585526, 118.359322, 100. + 626, 150.610001, 119.495255, 100. + 627, 120.551254, 120.242989, 100. + 628, 90.4351273, 120.626007, 100. + 629, 60.2864685, 120.733536, 100. + 630, 30.1319237, 120.714272, 100. + 631, 352.309296, 134.344604, 100. + 632, 324.663239, 135.721466, 100. + 633, 296.328308, 138.132324, 100. + 634, 267.507446, 140.876312, 100. + 635, 238.343216, 143.530441, 100. + 636, 208.929749, 145.869308, 100. + 637, 179.330536, 147.783951, 100. + 638, 149.591949, 149.22702, 100. + 639, 119.7509, 150.185593, 100. + 640, 89.8394928, 150.690521, 100. + 641, 59.8885803, 150.853058, 100. + 642, 29.9307117, 150.852463, 100. + 643, 347.375854, 160.354553, 100. + 644, 320.461884, 162.258362, 100. + 645, 292.733185, 165.28627, 100. + 646, 264.428223, 168.65657, 100. + 647, 235.714813, 171.894257, 100. + 648, 206.704163, 174.744217, 100. + 649, 177.471802, 177.082062, 100. + 650, 148.072998, 178.853287, 100. + 651, 118.551918, 180.043442, 100. + 652, 88.9468689, 180.690521, 100. + 653, 59.2944489, 180.926239, 100. + 654, 29.6324768, 180.957291, 100. + 655, 341.169617, 185.971741, 100. + 656, 315.086395, 188.455673, 100. + 657, 288.068298, 192.144928, 100. + 658, 260.386414, 196.178635, 100. + 659, 232.231873, 200.032379, 100. + 660, 203.732239, 203.423309, 100. + 661, 174.974686, 206.212952, 100. + 662, 146.023346, 208.340302, 100. + 663, 116.929596, 209.789658, 100. + 664, 87.7379913, 210.605667, 100. + 665, 58.4906311, 210.938232, 100. + 666, 29.2301559, 211.017517, 100. + 667, 333.686584, 211.116821, 100. + 668, 308.520233, 214.237442, 100. + 669, 282.311462, 218.636719, 100. + 670, 255.358063, 223.376648, 100. + 671, 227.870941, 227.885544, 100. + 672, 199.992142, 231.854523, 100. + 673, 171.819626, 235.131927, 100. + 674, 143.426025, 237.650848, 100. + 675, 114.869835, 239.394516, 100. + 676, 86.201767, 240.413391, 100. + 677, 57.4692993, 240.872528, 100. + 678, 28.7195873, 241.020782, 100. + 679, 324.99762, 235.660843, 100. + 680, 300.833466, 239.486359, 100. + 681, 275.531067, 244.653351, 100. + 682, 249.407806, 250.15097, 100. + 683, 222.69162, 255.363373, 100. + 684, 195.537811, 259.957214, 100. + 685, 168.054626, 263.768616, 100. + 686, 140.322937, 266.725098, 100. + 687, 112.408203, 268.808777, 100. + 688, 84.3672791, 270.074707, 100. + 689, 56.2529373, 270.699554, 100. + 690, 28.1167679, 270.945374, 100. + 691, 315.477386, 259.325012, 100. + 692, 292.469055, 263.937317, 100. + 693, 268.179993, 269.951599, 100. + 694, 242.966888, 276.282837, 100. + 695, 217.090057, 282.270782, 100. + 696, 190.724731, 287.559174, 100. + 697, 163.991791, 291.973114, 100. + 698, 136.981018, 295.435486, 100. + 699, 109.76503, 297.927429, 100. + 700, 82.4069519, 299.507507, 100. + 701, 54.9653358, 300.361877, 100. + 702, 27.4883575, 300.754333, 100. + 703, 306.185425, 281.512512, 100. + 704, 284.595001, 287.069275, 100. + 705, 261.398132, 294.091431, 100. + 706, 237.087585, 301.402466, 100. + 707, 212.013489, 308.294739, 100. + 708, 186.389832, 314.394531, 100. + 709, 160.35466, 319.522034, 100. + 710, 134.007492, 323.599274, 100. + 711, 107.428268, 326.608826, 100. + 712, 80.6871338, 328.616272, 100. + 713, 53.8459587, 329.808197, 100. + 714, 26.9473057, 330.424438, 100. + 715, 298.903351, 301.300476, 100. + 716, 278.958435, 308.21759, 100. + 717, 256.76825, 316.583801, 100. + 718, 233.173798, 325.140198, 100. + 719, 208.694977, 333.149658, 100. + 720, 183.605652, 340.239746, 100. + 721, 158.05777, 346.242981, 100. + 722, 132.157104, 351.091614, 100. + 723, 105.989746, 354.778778, 100. + 724, 79.6325531, 357.374084, 100. + 725, 53.1538162, 359.036652, 100. + 726, 26.6046143, 359.96405, 100. + 727, 26.7750816, 419.111877, 100. + 728, 27.7743015, 448.630798, 100. + 729, 53.5985489, 417.434937, 100. + 730, 55.5699883, 446.432983, 100. + 731, 80.3450623, 414.615295, 100. + 732, 83.2367935, 442.821136, 100. + 733, 106.907867, 410.521088, 100. + 734, 110.643425, 437.725067, 100. + 735, 133.190231, 405.075562, 100. + 736, 137.674316, 431.118195, 100. + 737, 159.100616, 398.245819, 100. + 738, 164.224823, 423.002411, 100. + 739, 184.549286, 390.027802, 100. + 740, 190.197495, 413.398285, 100. + 741, 209.445496, 380.43866, 100. + 742, 215.499146, 402.339966, 100. + 743, 233.697021, 369.529541, 100. + 744, 240.038361, 389.877319, 100. + 745, 257.088409, 357.470459, 100. + 746, 263.649811, 376.122467, 100. + 747, 279.085999, 344.675598, 100. + 748, 285.982422, 361.318726, 100. + 749, 298.797424, 331.722321, 100. + 750, 306.562988, 345.727051, 100. + 751, 448.509613, 30.5660229, 100. + 752, 418.962402, 28.8342876, 100. + 753, 445.937592, 60.9904785, 100. + 754, 417.099548, 57.4955788, 100. + 755, 441.769073, 90.8010025, 100. + 756, 414.03241, 85.6092377, 100. + 757, 436.391846, 118.202293, 100. + 758, 409.730835, 112.684601, 100. + 759, 429.53244, 144.8508, 100. + 760, 404.071167, 139.068329, 100. + 761, 421.133728, 171.011398, 100. + 762, 396.989044, 164.923126, 100. + 763, 411.21405, 196.618896, 100. + 764, 388.472626, 190.246719, 100. + 765, 399.81778, 221.551819, 100. + 766, 378.544495, 214.959793, 100. + 767, 387.000488, 245.688309, 100. + 768, 367.267426, 238.929092, 100. + 769, 372.842682, 268.875824, 100. + 770, 354.831757, 261.911407, 100. + 771, 357.588928, 290.834473, 100. + 772, 341.694458, 283.422516, 100. + 773, 341.607422, 311.109253, 100. + 774, 328.463165, 302.659302, 100. + 775, 469.57428, 99.4987411, 120. + 776, 475.348907, 66.6589813, 120. + 777, 478.709412, 35.1753426, 120. + 778, 28.7764206, 479.136627, 120. + 779, 57.5026703, 476.543213, 120. + 780, 86.0282516, 472.227844, 120. + 781, 114.222107, 466.21167, 120. + 782, 141.967514, 458.525055, 120. + 783, 169.157501, 449.205688, 120. + 784, 195.692413, 438.297241, 120. + 785, 221.478134, 425.849091, 120. + 786, 246.425385, 411.91568, 120. + 787, 270.449432, 396.556549, 120. + 788, 293.408508, 379.883484, 120. + 789, 315.142822, 362.05661, 120. + 790, 335.586884, 343.192993, 120. + 791, 354.918976, 323.160187, 120. + 792, 373.17923, 301.889496, 120. + 793, 390.288452, 279.418915, 120. + 794, 406.115387, 255.871658, 120. + 795, 420.563751, 231.357162, 120. + 796, 433.562225, 205.970383, 120. + 797, 445.04422, 179.821167, 120. + 798, 454.941071, 153.064117, 120. + 799, 463.192505, 125.907555, 120. + 800, 308.408386, 310.652832, 120. + 801, 26.3173389, 389.829346, 120. + 802, 52.6818466, 388.653961, 120. + 803, 78.9932709, 386.550964, 120. + 804, 105.167198, 383.335907, 120. + 805, 131.121841, 378.891907, 120. + 806, 156.775497, 373.167786, 120. + 807, 182.044037, 366.15274, 120. + 808, 206.836548, 357.865356, 120. + 809, 231.052521, 348.379974, 120. + 810, 254.458435, 337.954651, 120. + 811, 276.378052, 327.271606, 120. + 812, 295.395996, 317.642548, 120. + 813, 315.614349, 297.691711, 120. + 814, 325.730713, 279.077209, 120. + 815, 336.847534, 257.60675, 120. + 816, 347.654297, 234.562912, 120. + 817, 357.458038, 210.555557, 120. + 818, 365.999237, 185.864929, 120. + 819, 373.196716, 160.638474, 120. + 820, 379.032166, 134.960037, 120. + 821, 383.522583, 108.856888, 120. + 822, 386.735046, 82.2811966, 120. + 823, 388.803833, 55.1683884, 120. + 824, 389.893951, 27.6878395, 120. + 825, 324.677673, 329.254822, 120. + 826, 314.936401, 318.413727, 120. + 827, 360.740326, 27.1858025, 120. + 828, 331.282654, 27.3117256, 120. + 829, 301.589142, 27.7472534, 120. + 830, 271.725739, 28.2847176, 120. + 831, 241.741211, 28.8117523, 120. + 832, 211.668472, 29.2747231, 120. + 833, 181.529816, 29.650486, 120. + 834, 151.3414, 29.9301357, 120. + 835, 121.11602, 30.1112804, 120. + 836, 90.8629074, 30.1979103, 120. + 837, 60.5878525, 30.2102776, 120. + 838, 30.2969837, 30.1900826, 120. + 839, 360.044556, 54.3420868, 120. + 840, 330.861023, 54.6320419, 120. + 841, 301.350372, 55.4945908, 120. + 842, 271.607513, 56.5540733, 120. + 843, 241.700684, 57.5973282, 120. + 844, 211.674927, 58.51791, 120. + 845, 181.560577, 59.2677231, 120. + 846, 151.379471, 59.8273239, 120. + 847, 121.14904, 60.1912689, 120. + 848, 90.8833008, 60.369194, 120. + 849, 60.5941391, 60.4040298, 120. + 850, 30.2944794, 60.3757477, 120. + 851, 358.537354, 81.3123169, 120. + 852, 329.752258, 81.8564682, 120. + 853, 300.526794, 83.1806488, 120. + 854, 270.991333, 84.7731705, 120. + 855, 241.23761, 86.3363419, 120. + 856, 211.32576, 87.7171021, 120. + 857, 181.296173, 88.8440933, 120. + 858, 151.178055, 89.6874771, 120. + 859, 120.994743, 90.2388306, 120. + 860, 90.7659607, 90.5144196, 120. + 861, 60.5101814, 90.5805054, 120. + 862, 30.247448, 90.5526047, 120. + 863, 356.004425, 107.98999, 120. + 864, 327.733093, 108.90448, 120. + 865, 298.899048, 110.745125, 120. + 866, 269.668915, 112.897171, 120. + 867, 240.159576, 114.994896, 120. + 868, 210.447906, 116.84613, 120. + 869, 180.585526, 118.359322, 120. + 870, 150.610001, 119.495255, 120. + 871, 120.551254, 120.242989, 120. + 872, 90.4351273, 120.626007, 120. + 873, 60.2864685, 120.733536, 120. + 874, 30.1319237, 120.714272, 120. + 875, 352.309296, 134.344604, 120. + 876, 324.663239, 135.721466, 120. + 877, 296.328308, 138.132324, 120. + 878, 267.507446, 140.876312, 120. + 879, 238.343216, 143.530441, 120. + 880, 208.929749, 145.869308, 120. + 881, 179.330536, 147.783951, 120. + 882, 149.591949, 149.22702, 120. + 883, 119.7509, 150.185593, 120. + 884, 89.8394928, 150.690521, 120. + 885, 59.8885803, 150.853058, 120. + 886, 29.9307117, 150.852463, 120. + 887, 347.375854, 160.354553, 120. + 888, 320.461884, 162.258362, 120. + 889, 292.733185, 165.28627, 120. + 890, 264.428223, 168.65657, 120. + 891, 235.714813, 171.894257, 120. + 892, 206.704163, 174.744217, 120. + 893, 177.471802, 177.082062, 120. + 894, 148.072998, 178.853287, 120. + 895, 118.551918, 180.043442, 120. + 896, 88.9468689, 180.690521, 120. + 897, 59.2944489, 180.926239, 120. + 898, 29.6324768, 180.957291, 120. + 899, 341.169617, 185.971741, 120. + 900, 315.086395, 188.455673, 120. + 901, 288.068298, 192.144928, 120. + 902, 260.386414, 196.178635, 120. + 903, 232.231873, 200.032379, 120. + 904, 203.732239, 203.423309, 120. + 905, 174.974686, 206.212952, 120. + 906, 146.023346, 208.340302, 120. + 907, 116.929596, 209.789658, 120. + 908, 87.7379913, 210.605667, 120. + 909, 58.4906311, 210.938232, 120. + 910, 29.2301559, 211.017517, 120. + 911, 333.686584, 211.116821, 120. + 912, 308.520233, 214.237442, 120. + 913, 282.311462, 218.636719, 120. + 914, 255.358063, 223.376648, 120. + 915, 227.870941, 227.885544, 120. + 916, 199.992142, 231.854523, 120. + 917, 171.819626, 235.131927, 120. + 918, 143.426025, 237.650848, 120. + 919, 114.869835, 239.394516, 120. + 920, 86.201767, 240.413391, 120. + 921, 57.4692993, 240.872528, 120. + 922, 28.7195873, 241.020782, 120. + 923, 324.99762, 235.660843, 120. + 924, 300.833466, 239.486359, 120. + 925, 275.531067, 244.653351, 120. + 926, 249.407806, 250.15097, 120. + 927, 222.69162, 255.363373, 120. + 928, 195.537811, 259.957214, 120. + 929, 168.054626, 263.768616, 120. + 930, 140.322937, 266.725098, 120. + 931, 112.408203, 268.808777, 120. + 932, 84.3672791, 270.074707, 120. + 933, 56.2529373, 270.699554, 120. + 934, 28.1167679, 270.945374, 120. + 935, 315.477386, 259.325012, 120. + 936, 292.469055, 263.937317, 120. + 937, 268.179993, 269.951599, 120. + 938, 242.966888, 276.282837, 120. + 939, 217.090057, 282.270782, 120. + 940, 190.724731, 287.559174, 120. + 941, 163.991791, 291.973114, 120. + 942, 136.981018, 295.435486, 120. + 943, 109.76503, 297.927429, 120. + 944, 82.4069519, 299.507507, 120. + 945, 54.9653358, 300.361877, 120. + 946, 27.4883575, 300.754333, 120. + 947, 306.185425, 281.512512, 120. + 948, 284.595001, 287.069275, 120. + 949, 261.398132, 294.091431, 120. + 950, 237.087585, 301.402466, 120. + 951, 212.013489, 308.294739, 120. + 952, 186.389832, 314.394531, 120. + 953, 160.35466, 319.522034, 120. + 954, 134.007492, 323.599274, 120. + 955, 107.428268, 326.608826, 120. + 956, 80.6871338, 328.616272, 120. + 957, 53.8459587, 329.808197, 120. + 958, 26.9473057, 330.424438, 120. + 959, 298.903351, 301.300476, 120. + 960, 278.958435, 308.21759, 120. + 961, 256.76825, 316.583801, 120. + 962, 233.173798, 325.140198, 120. + 963, 208.694977, 333.149658, 120. + 964, 183.605652, 340.239746, 120. + 965, 158.05777, 346.242981, 120. + 966, 132.157104, 351.091614, 120. + 967, 105.989746, 354.778778, 120. + 968, 79.6325531, 357.374084, 120. + 969, 53.1538162, 359.036652, 120. + 970, 26.6046143, 359.96405, 120. + 971, 26.7750816, 419.111877, 120. + 972, 27.7743015, 448.630798, 120. + 973, 53.5985489, 417.434937, 120. + 974, 55.5699883, 446.432983, 120. + 975, 80.3450623, 414.615295, 120. + 976, 83.2367935, 442.821136, 120. + 977, 106.907867, 410.521088, 120. + 978, 110.643425, 437.725067, 120. + 979, 133.190231, 405.075562, 120. + 980, 137.674316, 431.118195, 120. + 981, 159.100616, 398.245819, 120. + 982, 164.224823, 423.002411, 120. + 983, 184.549286, 390.027802, 120. + 984, 190.197495, 413.398285, 120. + 985, 209.445496, 380.43866, 120. + 986, 215.499146, 402.339966, 120. + 987, 233.697021, 369.529541, 120. + 988, 240.038361, 389.877319, 120. + 989, 257.088409, 357.470459, 120. + 990, 263.649811, 376.122467, 120. + 991, 279.085999, 344.675598, 120. + 992, 285.982422, 361.318726, 120. + 993, 298.797424, 331.722321, 120. + 994, 306.562988, 345.727051, 120. + 995, 448.509613, 30.5660229, 120. + 996, 418.962402, 28.8342876, 120. + 997, 445.937592, 60.9904785, 120. + 998, 417.099548, 57.4955788, 120. + 999, 441.769073, 90.8010025, 120. + 1000, 414.03241, 85.6092377, 120. + 1001, 436.391846, 118.202293, 120. + 1002, 409.730835, 112.684601, 120. + 1003, 429.53244, 144.8508, 120. + 1004, 404.071167, 139.068329, 120. + 1005, 421.133728, 171.011398, 120. + 1006, 396.989044, 164.923126, 120. + 1007, 411.21405, 196.618896, 120. + 1008, 388.472626, 190.246719, 120. + 1009, 399.81778, 221.551819, 120. + 1010, 378.544495, 214.959793, 120. + 1011, 387.000488, 245.688309, 120. + 1012, 367.267426, 238.929092, 120. + 1013, 372.842682, 268.875824, 120. + 1014, 354.831757, 261.911407, 120. + 1015, 357.588928, 290.834473, 120. + 1016, 341.694458, 283.422516, 120. + 1017, 341.607422, 311.109253, 120. + 1018, 328.463165, 302.659302, 120. + 1019, 29.6389694, 499.120758, 53.3333321 + 1020, 59.2194176, 496.480682, 53.3333321 + 1021, 88.5909805, 492.08905, 53.3333321 + 1022, 117.649979, 485.961395, 53.3333321 + 1023, 146.293869, 478.119354, 53.3333321 + 1024, 174.421707, 468.590515, 53.3333321 + 1025, 201.934189, 457.408539, 53.3333321 + 1026, 228.734299, 444.612885, 53.3333321 + 1027, 254.727432, 430.248688, 53.3333321 + 1028, 279.821899, 414.366638, 53.3333321 + 1029, 303.929199, 397.022705, 53.3333321 + 1030, 326.964264, 378.278168, 53.3333321 + 1031, 348.845856, 358.199066, 53.3333321 + 1032, 369.496796, 336.856232, 53.3333321 + 1033, 388.844299, 314.32486, 53.3333321 + 1034, 406.820221, 290.684204, 53.3333321 + 1035, 423.361603, 266.016815, 53.3333321 + 1036, 438.411163, 240.407242, 53.3333321 + 1037, 451.91861, 213.938263, 53.3333321 + 1038, 463.843109, 186.680389, 53.3333321 + 1039, 474.158966, 158.660889, 53.3333321 + 1040, 482.865601, 129.772141, 53.3333321 + 1041, 29.6389694, 499.120758, 26.666666 + 1042, 59.2194176, 496.480682, 26.666666 + 1043, 88.5909805, 492.08905, 26.666666 + 1044, 117.649979, 485.961395, 26.666666 + 1045, 146.293869, 478.119354, 26.666666 + 1046, 174.421707, 468.590515, 26.666666 + 1047, 201.934189, 457.408539, 26.666666 + 1048, 228.734299, 444.612885, 26.666666 + 1049, 254.727432, 430.248688, 26.666666 + 1050, 279.821899, 414.366638, 26.666666 + 1051, 303.929199, 397.022705, 26.666666 + 1052, 326.964264, 378.278168, 26.666666 + 1053, 348.845856, 358.199066, 26.666666 + 1054, 369.496796, 336.856232, 26.666666 + 1055, 388.844299, 314.32486, 26.666666 + 1056, 406.820221, 290.684204, 26.666666 + 1057, 423.361603, 266.016815, 26.666666 + 1058, 438.411163, 240.407242, 26.666666 + 1059, 451.91861, 213.938263, 26.666666 + 1060, 463.843109, 186.680389, 26.666666 + 1061, 474.158966, 158.660889, 26.666666 + 1062, 482.865601, 129.772141, 26.666666 + 1063, 463.192505, 125.907555, 53.3333321 + 1064, 454.941071, 153.064117, 53.3333321 + 1065, 445.04422, 179.821167, 53.3333321 + 1066, 433.562225, 205.970383, 53.3333321 + 1067, 420.563751, 231.357162, 53.3333321 + 1068, 406.115387, 255.871658, 53.3333321 + 1069, 390.288452, 279.418915, 53.3333321 + 1070, 373.17923, 301.889496, 53.3333321 + 1071, 354.918976, 323.160187, 53.3333321 + 1072, 335.586884, 343.192993, 53.3333321 + 1073, 315.142822, 362.05661, 53.3333321 + 1074, 293.408508, 379.883484, 53.3333321 + 1075, 270.449432, 396.556549, 53.3333321 + 1076, 246.425385, 411.91568, 53.3333321 + 1077, 221.478134, 425.849091, 53.3333321 + 1078, 195.692413, 438.297241, 53.3333321 + 1079, 169.157501, 449.205688, 53.3333321 + 1080, 141.967514, 458.525055, 53.3333321 + 1081, 114.222107, 466.21167, 53.3333321 + 1082, 86.0282516, 472.227844, 53.3333321 + 1083, 57.5026703, 476.543213, 53.3333321 + 1084, 28.7764206, 479.136627, 53.3333321 + 1085, 463.192505, 125.907555, 26.666666 + 1086, 454.941071, 153.064117, 26.666666 + 1087, 445.04422, 179.821167, 26.666666 + 1088, 433.562225, 205.970383, 26.666666 + 1089, 420.563751, 231.357162, 26.666666 + 1090, 406.115387, 255.871658, 26.666666 + 1091, 390.288452, 279.418915, 26.666666 + 1092, 373.17923, 301.889496, 26.666666 + 1093, 354.918976, 323.160187, 26.666666 + 1094, 335.586884, 343.192993, 26.666666 + 1095, 315.142822, 362.05661, 26.666666 + 1096, 293.408508, 379.883484, 26.666666 + 1097, 270.449432, 396.556549, 26.666666 + 1098, 246.425385, 411.91568, 26.666666 + 1099, 221.478134, 425.849091, 26.666666 + 1100, 195.692413, 438.297241, 26.666666 + 1101, 169.157501, 449.205688, 26.666666 + 1102, 141.967514, 458.525055, 26.666666 + 1103, 114.222107, 466.21167, 26.666666 + 1104, 86.0282516, 472.227844, 26.666666 + 1105, 57.5026703, 476.543213, 26.666666 + 1106, 28.7764206, 479.136627, 26.666666 +*Element, type=C3D8R +1, 31, 30, 2, 1, 300, 297, 35, 38 +2, 4, 3, 30, 31, 43, 34, 297, 300 +3, 300, 297, 35, 38, 301, 298, 36, 39 +4, 43, 34, 297, 300, 42, 33, 298, 301 +5, 301, 298, 36, 39, 302, 299, 37, 40 +6, 42, 33, 298, 301, 41, 32, 299, 302 +7, 302, 299, 37, 40, 44, 45, 5, 6 +8, 41, 32, 299, 302, 7, 8, 45, 44 + 9, 1, 48, 303, 38, 57, 56, 309, 306 +10, 48, 9, 47, 303, 56, 14, 55, 309 +11, 304, 303, 47, 46, 310, 309, 55, 54 +12, 304, 39, 38, 303, 310, 307, 306, 309 +13, 304, 305, 40, 39, 310, 311, 308, 307 +14, 305, 49, 6, 40, 311, 62, 61, 308 +15, 305, 50, 10, 49, 311, 52, 12, 62 +16, 304, 51, 50, 305, 310, 53, 52, 311 +17, 304, 46, 11, 51, 310, 54, 13, 53 +18, 1, 38, 35, 2, 57, 306, 58, 15 +19, 38, 39, 36, 35, 306, 307, 59, 58 +20, 39, 40, 37, 36, 307, 308, 60, 59 +21, 40, 6, 5, 37, 308, 61, 16, 60 + 22, 9, 47, 55, 14, 23, 176, 533, 177 + 23, 46, 54, 55, 47, 175, 532, 533, 176 + 24, 46, 11, 13, 54, 175, 22, 531, 532 + 25, 336, 339, 103, 102, 580, 583, 181, 180 + 26, 339, 340, 104, 103, 583, 584, 182, 181 + 27, 340, 341, 105, 104, 584, 585, 183, 182 + 28, 341, 342, 106, 105, 585, 586, 184, 183 + 29, 342, 343, 107, 106, 586, 587, 185, 184 + 30, 343, 344, 108, 107, 587, 588, 186, 185 + 31, 344, 345, 109, 108, 588, 589, 187, 186 + 32, 345, 346, 110, 109, 589, 590, 188, 187 + 33, 346, 347, 111, 110, 590, 591, 189, 188 + 34, 347, 348, 112, 111, 591, 592, 190, 189 + 35, 348, 349, 113, 112, 592, 593, 191, 190 + 36, 349, 350, 114, 113, 593, 594, 192, 191 + 37, 350, 63, 17, 114, 594, 137, 20, 192 + 38, 335, 351, 339, 336, 579, 595, 583, 580 + 39, 351, 352, 340, 339, 595, 596, 584, 583 + 40, 352, 353, 341, 340, 596, 597, 585, 584 + 41, 353, 354, 342, 341, 597, 598, 586, 585 + 42, 354, 355, 343, 342, 598, 599, 587, 586 + 43, 355, 356, 344, 343, 599, 600, 588, 587 + 44, 356, 357, 345, 344, 600, 601, 589, 588 + 45, 357, 358, 346, 345, 601, 602, 590, 589 + 46, 358, 359, 347, 346, 602, 603, 591, 590 + 47, 359, 360, 348, 347, 603, 604, 592, 591 + 48, 360, 361, 349, 348, 604, 605, 593, 592 + 49, 361, 362, 350, 349, 605, 606, 594, 593 + 50, 362, 64, 63, 350, 606, 138, 137, 594 + 51, 334, 363, 351, 335, 578, 607, 595, 579 + 52, 363, 364, 352, 351, 607, 608, 596, 595 + 53, 364, 365, 353, 352, 608, 609, 597, 596 + 54, 365, 366, 354, 353, 609, 610, 598, 597 + 55, 366, 367, 355, 354, 610, 611, 599, 598 + 56, 367, 368, 356, 355, 611, 612, 600, 599 + 57, 368, 369, 357, 356, 612, 613, 601, 600 + 58, 369, 370, 358, 357, 613, 614, 602, 601 + 59, 370, 371, 359, 358, 614, 615, 603, 602 + 60, 371, 372, 360, 359, 615, 616, 604, 603 + 61, 372, 373, 361, 360, 616, 617, 605, 604 + 62, 373, 374, 362, 361, 617, 618, 606, 605 + 63, 374, 65, 64, 362, 618, 139, 138, 606 + 64, 333, 375, 363, 334, 577, 619, 607, 578 + 65, 375, 376, 364, 363, 619, 620, 608, 607 + 66, 376, 377, 365, 364, 620, 621, 609, 608 + 67, 377, 378, 366, 365, 621, 622, 610, 609 + 68, 378, 379, 367, 366, 622, 623, 611, 610 + 69, 379, 380, 368, 367, 623, 624, 612, 611 + 70, 380, 381, 369, 368, 624, 625, 613, 612 + 71, 381, 382, 370, 369, 625, 626, 614, 613 + 72, 382, 383, 371, 370, 626, 627, 615, 614 + 73, 383, 384, 372, 371, 627, 628, 616, 615 + 74, 384, 385, 373, 372, 628, 629, 617, 616 + 75, 385, 386, 374, 373, 629, 630, 618, 617 + 76, 386, 66, 65, 374, 630, 140, 139, 618 + 77, 332, 387, 375, 333, 576, 631, 619, 577 + 78, 387, 388, 376, 375, 631, 632, 620, 619 + 79, 388, 389, 377, 376, 632, 633, 621, 620 + 80, 389, 390, 378, 377, 633, 634, 622, 621 + 81, 390, 391, 379, 378, 634, 635, 623, 622 + 82, 391, 392, 380, 379, 635, 636, 624, 623 + 83, 392, 393, 381, 380, 636, 637, 625, 624 + 84, 393, 394, 382, 381, 637, 638, 626, 625 + 85, 394, 395, 383, 382, 638, 639, 627, 626 + 86, 395, 396, 384, 383, 639, 640, 628, 627 + 87, 396, 397, 385, 384, 640, 641, 629, 628 + 88, 397, 398, 386, 385, 641, 642, 630, 629 + 89, 398, 67, 66, 386, 642, 141, 140, 630 + 90, 331, 399, 387, 332, 575, 643, 631, 576 + 91, 399, 400, 388, 387, 643, 644, 632, 631 + 92, 400, 401, 389, 388, 644, 645, 633, 632 + 93, 401, 402, 390, 389, 645, 646, 634, 633 + 94, 402, 403, 391, 390, 646, 647, 635, 634 + 95, 403, 404, 392, 391, 647, 648, 636, 635 + 96, 404, 405, 393, 392, 648, 649, 637, 636 + 97, 405, 406, 394, 393, 649, 650, 638, 637 + 98, 406, 407, 395, 394, 650, 651, 639, 638 + 99, 407, 408, 396, 395, 651, 652, 640, 639 +100, 408, 409, 397, 396, 652, 653, 641, 640 +101, 409, 410, 398, 397, 653, 654, 642, 641 +102, 410, 68, 67, 398, 654, 142, 141, 642 +103, 330, 411, 399, 331, 574, 655, 643, 575 +104, 411, 412, 400, 399, 655, 656, 644, 643 +105, 412, 413, 401, 400, 656, 657, 645, 644 +106, 413, 414, 402, 401, 657, 658, 646, 645 +107, 414, 415, 403, 402, 658, 659, 647, 646 +108, 415, 416, 404, 403, 659, 660, 648, 647 +109, 416, 417, 405, 404, 660, 661, 649, 648 +110, 417, 418, 406, 405, 661, 662, 650, 649 +111, 418, 419, 407, 406, 662, 663, 651, 650 +112, 419, 420, 408, 407, 663, 664, 652, 651 +113, 420, 421, 409, 408, 664, 665, 653, 652 +114, 421, 422, 410, 409, 665, 666, 654, 653 +115, 422, 69, 68, 410, 666, 143, 142, 654 +116, 329, 423, 411, 330, 573, 667, 655, 574 +117, 423, 424, 412, 411, 667, 668, 656, 655 +118, 424, 425, 413, 412, 668, 669, 657, 656 +119, 425, 426, 414, 413, 669, 670, 658, 657 +120, 426, 427, 415, 414, 670, 671, 659, 658 +121, 427, 428, 416, 415, 671, 672, 660, 659 +122, 428, 429, 417, 416, 672, 673, 661, 660 +123, 429, 430, 418, 417, 673, 674, 662, 661 +124, 430, 431, 419, 418, 674, 675, 663, 662 +125, 431, 432, 420, 419, 675, 676, 664, 663 +126, 432, 433, 421, 420, 676, 677, 665, 664 +127, 433, 434, 422, 421, 677, 678, 666, 665 +128, 434, 70, 69, 422, 678, 144, 143, 666 +129, 328, 435, 423, 329, 572, 679, 667, 573 +130, 435, 436, 424, 423, 679, 680, 668, 667 +131, 436, 437, 425, 424, 680, 681, 669, 668 +132, 437, 438, 426, 425, 681, 682, 670, 669 +133, 438, 439, 427, 426, 682, 683, 671, 670 +134, 439, 440, 428, 427, 683, 684, 672, 671 +135, 440, 441, 429, 428, 684, 685, 673, 672 +136, 441, 442, 430, 429, 685, 686, 674, 673 +137, 442, 443, 431, 430, 686, 687, 675, 674 +138, 443, 444, 432, 431, 687, 688, 676, 675 +139, 444, 445, 433, 432, 688, 689, 677, 676 +140, 445, 446, 434, 433, 689, 690, 678, 677 +141, 446, 71, 70, 434, 690, 145, 144, 678 +142, 327, 447, 435, 328, 571, 691, 679, 572 +143, 447, 448, 436, 435, 691, 692, 680, 679 +144, 448, 449, 437, 436, 692, 693, 681, 680 +145, 449, 450, 438, 437, 693, 694, 682, 681 +146, 450, 451, 439, 438, 694, 695, 683, 682 +147, 451, 452, 440, 439, 695, 696, 684, 683 +148, 452, 453, 441, 440, 696, 697, 685, 684 +149, 453, 454, 442, 441, 697, 698, 686, 685 +150, 454, 455, 443, 442, 698, 699, 687, 686 +151, 455, 456, 444, 443, 699, 700, 688, 687 +152, 456, 457, 445, 444, 700, 701, 689, 688 +153, 457, 458, 446, 445, 701, 702, 690, 689 +154, 458, 72, 71, 446, 702, 146, 145, 690 +155, 326, 459, 447, 327, 570, 703, 691, 571 +156, 459, 460, 448, 447, 703, 704, 692, 691 +157, 460, 461, 449, 448, 704, 705, 693, 692 +158, 461, 462, 450, 449, 705, 706, 694, 693 +159, 462, 463, 451, 450, 706, 707, 695, 694 +160, 463, 464, 452, 451, 707, 708, 696, 695 +161, 464, 465, 453, 452, 708, 709, 697, 696 +162, 465, 466, 454, 453, 709, 710, 698, 697 +163, 466, 467, 455, 454, 710, 711, 699, 698 +164, 467, 468, 456, 455, 711, 712, 700, 699 +165, 468, 469, 457, 456, 712, 713, 701, 700 +166, 469, 470, 458, 457, 713, 714, 702, 701 +167, 470, 73, 72, 458, 714, 147, 146, 702 +168, 325, 471, 459, 326, 569, 715, 703, 570 +169, 471, 472, 460, 459, 715, 716, 704, 703 +170, 472, 473, 461, 460, 716, 717, 705, 704 +171, 473, 474, 462, 461, 717, 718, 706, 705 +172, 474, 475, 463, 462, 718, 719, 707, 706 +173, 475, 476, 464, 463, 719, 720, 708, 707 +174, 476, 477, 465, 464, 720, 721, 709, 708 +175, 477, 478, 466, 465, 721, 722, 710, 709 +176, 478, 479, 467, 466, 722, 723, 711, 710 +177, 479, 480, 468, 467, 723, 724, 712, 711 +178, 480, 481, 469, 468, 724, 725, 713, 712 +179, 481, 482, 470, 469, 725, 726, 714, 713 +180, 482, 74, 73, 470, 726, 148, 147, 714 +181, 312, 324, 471, 325, 556, 568, 715, 569 +182, 324, 323, 472, 471, 568, 567, 716, 715 +183, 323, 322, 473, 472, 567, 566, 717, 716 +184, 322, 321, 474, 473, 566, 565, 718, 717 +185, 321, 320, 475, 474, 565, 564, 719, 718 +186, 320, 319, 476, 475, 564, 563, 720, 719 +187, 319, 318, 477, 476, 563, 562, 721, 720 +188, 318, 317, 478, 477, 562, 561, 722, 721 +189, 317, 316, 479, 478, 561, 560, 723, 722 +190, 316, 315, 480, 479, 560, 559, 724, 723 +191, 315, 314, 481, 480, 559, 558, 725, 724 +192, 314, 313, 482, 481, 558, 557, 726, 725 +193, 313, 75, 74, 482, 557, 149, 148, 726 +194, 313, 483, 76, 75, 557, 727, 150, 149 +195, 483, 484, 77, 76, 727, 728, 151, 150 +196, 484, 78, 18, 77, 728, 534, 152, 151 +197, 314, 485, 483, 313, 558, 729, 727, 557 +198, 485, 486, 484, 483, 729, 730, 728, 727 +199, 486, 79, 78, 484, 730, 535, 534, 728 +200, 315, 487, 485, 314, 559, 731, 729, 558 +201, 487, 488, 486, 485, 731, 732, 730, 729 +202, 488, 80, 79, 486, 732, 536, 535, 730 +203, 316, 489, 487, 315, 560, 733, 731, 559 +204, 489, 490, 488, 487, 733, 734, 732, 731 +205, 490, 81, 80, 488, 734, 537, 536, 732 +206, 317, 491, 489, 316, 561, 735, 733, 560 +207, 491, 492, 490, 489, 735, 736, 734, 733 +208, 492, 82, 81, 490, 736, 538, 537, 734 +209, 318, 493, 491, 317, 562, 737, 735, 561 +210, 493, 494, 492, 491, 737, 738, 736, 735 +211, 494, 83, 82, 492, 738, 539, 538, 736 +212, 319, 495, 493, 318, 563, 739, 737, 562 +213, 495, 496, 494, 493, 739, 740, 738, 737 +214, 496, 84, 83, 494, 740, 540, 539, 738 +215, 320, 497, 495, 319, 564, 741, 739, 563 +216, 497, 498, 496, 495, 741, 742, 740, 739 +217, 498, 85, 84, 496, 742, 541, 540, 740 +218, 321, 499, 497, 320, 565, 743, 741, 564 +219, 499, 500, 498, 497, 743, 744, 742, 741 +220, 500, 86, 85, 498, 744, 542, 541, 742 +221, 322, 501, 499, 321, 566, 745, 743, 565 +222, 501, 502, 500, 499, 745, 746, 744, 743 +223, 502, 87, 86, 500, 746, 543, 542, 744 +224, 323, 503, 501, 322, 567, 747, 745, 566 +225, 503, 504, 502, 501, 747, 748, 746, 745 +226, 504, 88, 87, 502, 748, 544, 543, 746 +227, 324, 505, 503, 323, 568, 749, 747, 567 +228, 505, 506, 504, 503, 749, 750, 748, 747 +229, 506, 89, 88, 504, 750, 545, 544, 748 +230, 312, 338, 505, 324, 556, 582, 749, 568 +231, 338, 337, 506, 505, 582, 581, 750, 749 +232, 337, 90, 89, 506, 581, 546, 545, 750 +233, 55, 507, 100, 14, 533, 751, 178, 177 +234, 507, 508, 101, 100, 751, 752, 179, 178 +235, 508, 336, 102, 101, 752, 580, 180, 179 +236, 54, 509, 507, 55, 532, 753, 751, 533 +237, 509, 510, 508, 507, 753, 754, 752, 751 +238, 510, 335, 336, 508, 754, 579, 580, 752 +239, 13, 511, 509, 54, 531, 755, 753, 532 +240, 511, 512, 510, 509, 755, 756, 754, 753 +241, 512, 334, 335, 510, 756, 578, 579, 754 +242, 99, 513, 511, 13, 555, 757, 755, 531 +243, 513, 514, 512, 511, 757, 758, 756, 755 +244, 514, 333, 334, 512, 758, 577, 578, 756 +245, 98, 515, 513, 99, 554, 759, 757, 555 +246, 515, 516, 514, 513, 759, 760, 758, 757 +247, 516, 332, 333, 514, 760, 576, 577, 758 +248, 97, 517, 515, 98, 553, 761, 759, 554 +249, 517, 518, 516, 515, 761, 762, 760, 759 +250, 518, 331, 332, 516, 762, 575, 576, 760 +251, 96, 519, 517, 97, 552, 763, 761, 553 +252, 519, 520, 518, 517, 763, 764, 762, 761 +253, 520, 330, 331, 518, 764, 574, 575, 762 +254, 95, 521, 519, 96, 551, 765, 763, 552 +255, 521, 522, 520, 519, 765, 766, 764, 763 +256, 522, 329, 330, 520, 766, 573, 574, 764 +257, 94, 523, 521, 95, 550, 767, 765, 551 +258, 523, 524, 522, 521, 767, 768, 766, 765 +259, 524, 328, 329, 522, 768, 572, 573, 766 +260, 93, 525, 523, 94, 549, 769, 767, 550 +261, 525, 526, 524, 523, 769, 770, 768, 767 +262, 526, 327, 328, 524, 770, 571, 572, 768 +263, 92, 527, 525, 93, 548, 771, 769, 549 +264, 527, 528, 526, 525, 771, 772, 770, 769 +265, 528, 326, 327, 526, 772, 570, 571, 770 +266, 91, 529, 527, 92, 547, 773, 771, 548 +267, 529, 530, 528, 527, 773, 774, 772, 771 +268, 530, 325, 326, 528, 774, 569, 570, 772 +269, 90, 337, 529, 91, 546, 581, 773, 547 +270, 337, 338, 530, 529, 581, 582, 774, 773 +271, 338, 312, 325, 530, 582, 556, 569, 774 +272, 13, 11, 115, 99, 531, 22, 174, 555 +273, 99, 115, 116, 98, 555, 174, 173, 554 +274, 98, 116, 117, 97, 554, 173, 172, 553 +275, 97, 117, 118, 96, 553, 172, 171, 552 +276, 96, 118, 119, 95, 552, 171, 170, 551 +277, 95, 119, 120, 94, 551, 170, 169, 550 +278, 94, 120, 121, 93, 550, 169, 168, 549 +279, 93, 121, 122, 92, 549, 168, 167, 548 +280, 92, 122, 123, 91, 548, 167, 166, 547 +281, 91, 123, 124, 90, 547, 166, 165, 546 +282, 90, 124, 125, 89, 546, 165, 164, 545 +283, 89, 125, 126, 88, 545, 164, 163, 544 +284, 88, 126, 127, 87, 544, 163, 162, 543 +285, 87, 127, 128, 86, 543, 162, 161, 542 +286, 86, 128, 129, 85, 542, 161, 160, 541 +287, 85, 129, 130, 84, 541, 160, 159, 540 +288, 84, 130, 131, 83, 540, 159, 158, 539 +289, 83, 131, 132, 82, 539, 158, 157, 538 +290, 82, 132, 133, 81, 538, 157, 156, 537 +291, 81, 133, 134, 80, 537, 156, 155, 536 +292, 80, 134, 135, 79, 536, 155, 154, 535 +293, 79, 135, 136, 78, 535, 154, 153, 534 +294, 78, 136, 19, 18, 534, 153, 21, 152 +295, 23, 176, 533, 177, 26, 225, 777, 224 +296, 175, 532, 533, 176, 226, 776, 777, 225 +297, 175, 22, 531, 532, 226, 27, 775, 776 +298, 580, 583, 181, 180, 824, 827, 220, 221 +299, 583, 584, 182, 181, 827, 828, 219, 220 +300, 584, 585, 183, 182, 828, 829, 218, 219 +301, 585, 586, 184, 183, 829, 830, 217, 218 +302, 586, 587, 185, 184, 830, 831, 216, 217 +303, 587, 588, 186, 185, 831, 832, 215, 216 +304, 588, 589, 187, 186, 832, 833, 214, 215 +305, 589, 590, 188, 187, 833, 834, 213, 214 +306, 590, 591, 189, 188, 834, 835, 212, 213 +307, 591, 592, 190, 189, 835, 836, 211, 212 +308, 592, 593, 191, 190, 836, 837, 210, 211 +309, 593, 594, 192, 191, 837, 838, 209, 210 +310, 594, 137, 20, 192, 838, 208, 25, 209 +311, 579, 595, 583, 580, 823, 839, 827, 824 +312, 595, 596, 584, 583, 839, 840, 828, 827 +313, 596, 597, 585, 584, 840, 841, 829, 828 +314, 597, 598, 586, 585, 841, 842, 830, 829 +315, 598, 599, 587, 586, 842, 843, 831, 830 +316, 599, 600, 588, 587, 843, 844, 832, 831 +317, 600, 601, 589, 588, 844, 845, 833, 832 +318, 601, 602, 590, 589, 845, 846, 834, 833 +319, 602, 603, 591, 590, 846, 847, 835, 834 +320, 603, 604, 592, 591, 847, 848, 836, 835 +321, 604, 605, 593, 592, 848, 849, 837, 836 +322, 605, 606, 594, 593, 849, 850, 838, 837 +323, 606, 138, 137, 594, 850, 207, 208, 838 +324, 578, 607, 595, 579, 822, 851, 839, 823 +325, 607, 608, 596, 595, 851, 852, 840, 839 +326, 608, 609, 597, 596, 852, 853, 841, 840 +327, 609, 610, 598, 597, 853, 854, 842, 841 +328, 610, 611, 599, 598, 854, 855, 843, 842 +329, 611, 612, 600, 599, 855, 856, 844, 843 +330, 612, 613, 601, 600, 856, 857, 845, 844 +331, 613, 614, 602, 601, 857, 858, 846, 845 +332, 614, 615, 603, 602, 858, 859, 847, 846 +333, 615, 616, 604, 603, 859, 860, 848, 847 +334, 616, 617, 605, 604, 860, 861, 849, 848 +335, 617, 618, 606, 605, 861, 862, 850, 849 +336, 618, 139, 138, 606, 862, 206, 207, 850 +337, 577, 619, 607, 578, 821, 863, 851, 822 +338, 619, 620, 608, 607, 863, 864, 852, 851 +339, 620, 621, 609, 608, 864, 865, 853, 852 +340, 621, 622, 610, 609, 865, 866, 854, 853 +341, 622, 623, 611, 610, 866, 867, 855, 854 +342, 623, 624, 612, 611, 867, 868, 856, 855 +343, 624, 625, 613, 612, 868, 869, 857, 856 +344, 625, 626, 614, 613, 869, 870, 858, 857 +345, 626, 627, 615, 614, 870, 871, 859, 858 +346, 627, 628, 616, 615, 871, 872, 860, 859 +347, 628, 629, 617, 616, 872, 873, 861, 860 +348, 629, 630, 618, 617, 873, 874, 862, 861 +349, 630, 140, 139, 618, 874, 205, 206, 862 +350, 576, 631, 619, 577, 820, 875, 863, 821 +351, 631, 632, 620, 619, 875, 876, 864, 863 +352, 632, 633, 621, 620, 876, 877, 865, 864 +353, 633, 634, 622, 621, 877, 878, 866, 865 +354, 634, 635, 623, 622, 878, 879, 867, 866 +355, 635, 636, 624, 623, 879, 880, 868, 867 +356, 636, 637, 625, 624, 880, 881, 869, 868 +357, 637, 638, 626, 625, 881, 882, 870, 869 +358, 638, 639, 627, 626, 882, 883, 871, 870 +359, 639, 640, 628, 627, 883, 884, 872, 871 +360, 640, 641, 629, 628, 884, 885, 873, 872 +361, 641, 642, 630, 629, 885, 886, 874, 873 +362, 642, 141, 140, 630, 886, 204, 205, 874 +363, 575, 643, 631, 576, 819, 887, 875, 820 +364, 643, 644, 632, 631, 887, 888, 876, 875 +365, 644, 645, 633, 632, 888, 889, 877, 876 +366, 645, 646, 634, 633, 889, 890, 878, 877 +367, 646, 647, 635, 634, 890, 891, 879, 878 +368, 647, 648, 636, 635, 891, 892, 880, 879 +369, 648, 649, 637, 636, 892, 893, 881, 880 +370, 649, 650, 638, 637, 893, 894, 882, 881 +371, 650, 651, 639, 638, 894, 895, 883, 882 +372, 651, 652, 640, 639, 895, 896, 884, 883 +373, 652, 653, 641, 640, 896, 897, 885, 884 +374, 653, 654, 642, 641, 897, 898, 886, 885 +375, 654, 142, 141, 642, 898, 203, 204, 886 +376, 574, 655, 643, 575, 818, 899, 887, 819 +377, 655, 656, 644, 643, 899, 900, 888, 887 +378, 656, 657, 645, 644, 900, 901, 889, 888 +379, 657, 658, 646, 645, 901, 902, 890, 889 +380, 658, 659, 647, 646, 902, 903, 891, 890 +381, 659, 660, 648, 647, 903, 904, 892, 891 +382, 660, 661, 649, 648, 904, 905, 893, 892 +383, 661, 662, 650, 649, 905, 906, 894, 893 +384, 662, 663, 651, 650, 906, 907, 895, 894 +385, 663, 664, 652, 651, 907, 908, 896, 895 +386, 664, 665, 653, 652, 908, 909, 897, 896 +387, 665, 666, 654, 653, 909, 910, 898, 897 +388, 666, 143, 142, 654, 910, 202, 203, 898 +389, 573, 667, 655, 574, 817, 911, 899, 818 +390, 667, 668, 656, 655, 911, 912, 900, 899 +391, 668, 669, 657, 656, 912, 913, 901, 900 +392, 669, 670, 658, 657, 913, 914, 902, 901 +393, 670, 671, 659, 658, 914, 915, 903, 902 +394, 671, 672, 660, 659, 915, 916, 904, 903 +395, 672, 673, 661, 660, 916, 917, 905, 904 +396, 673, 674, 662, 661, 917, 918, 906, 905 +397, 674, 675, 663, 662, 918, 919, 907, 906 +398, 675, 676, 664, 663, 919, 920, 908, 907 +399, 676, 677, 665, 664, 920, 921, 909, 908 +400, 677, 678, 666, 665, 921, 922, 910, 909 +401, 678, 144, 143, 666, 922, 201, 202, 910 +402, 572, 679, 667, 573, 816, 923, 911, 817 +403, 679, 680, 668, 667, 923, 924, 912, 911 +404, 680, 681, 669, 668, 924, 925, 913, 912 +405, 681, 682, 670, 669, 925, 926, 914, 913 +406, 682, 683, 671, 670, 926, 927, 915, 914 +407, 683, 684, 672, 671, 927, 928, 916, 915 +408, 684, 685, 673, 672, 928, 929, 917, 916 +409, 685, 686, 674, 673, 929, 930, 918, 917 +410, 686, 687, 675, 674, 930, 931, 919, 918 +411, 687, 688, 676, 675, 931, 932, 920, 919 +412, 688, 689, 677, 676, 932, 933, 921, 920 +413, 689, 690, 678, 677, 933, 934, 922, 921 +414, 690, 145, 144, 678, 934, 200, 201, 922 +415, 571, 691, 679, 572, 815, 935, 923, 816 +416, 691, 692, 680, 679, 935, 936, 924, 923 +417, 692, 693, 681, 680, 936, 937, 925, 924 +418, 693, 694, 682, 681, 937, 938, 926, 925 +419, 694, 695, 683, 682, 938, 939, 927, 926 +420, 695, 696, 684, 683, 939, 940, 928, 927 +421, 696, 697, 685, 684, 940, 941, 929, 928 +422, 697, 698, 686, 685, 941, 942, 930, 929 +423, 698, 699, 687, 686, 942, 943, 931, 930 +424, 699, 700, 688, 687, 943, 944, 932, 931 +425, 700, 701, 689, 688, 944, 945, 933, 932 +426, 701, 702, 690, 689, 945, 946, 934, 933 +427, 702, 146, 145, 690, 946, 199, 200, 934 +428, 570, 703, 691, 571, 814, 947, 935, 815 +429, 703, 704, 692, 691, 947, 948, 936, 935 +430, 704, 705, 693, 692, 948, 949, 937, 936 +431, 705, 706, 694, 693, 949, 950, 938, 937 +432, 706, 707, 695, 694, 950, 951, 939, 938 +433, 707, 708, 696, 695, 951, 952, 940, 939 +434, 708, 709, 697, 696, 952, 953, 941, 940 +435, 709, 710, 698, 697, 953, 954, 942, 941 +436, 710, 711, 699, 698, 954, 955, 943, 942 +437, 711, 712, 700, 699, 955, 956, 944, 943 +438, 712, 713, 701, 700, 956, 957, 945, 944 +439, 713, 714, 702, 701, 957, 958, 946, 945 +440, 714, 147, 146, 702, 958, 198, 199, 946 +441, 569, 715, 703, 570, 813, 959, 947, 814 +442, 715, 716, 704, 703, 959, 960, 948, 947 +443, 716, 717, 705, 704, 960, 961, 949, 948 +444, 717, 718, 706, 705, 961, 962, 950, 949 +445, 718, 719, 707, 706, 962, 963, 951, 950 +446, 719, 720, 708, 707, 963, 964, 952, 951 +447, 720, 721, 709, 708, 964, 965, 953, 952 +448, 721, 722, 710, 709, 965, 966, 954, 953 +449, 722, 723, 711, 710, 966, 967, 955, 954 +450, 723, 724, 712, 711, 967, 968, 956, 955 +451, 724, 725, 713, 712, 968, 969, 957, 956 +452, 725, 726, 714, 713, 969, 970, 958, 957 +453, 726, 148, 147, 714, 970, 197, 198, 958 +454, 556, 568, 715, 569, 800, 812, 959, 813 +455, 568, 567, 716, 715, 812, 811, 960, 959 +456, 567, 566, 717, 716, 811, 810, 961, 960 +457, 566, 565, 718, 717, 810, 809, 962, 961 +458, 565, 564, 719, 718, 809, 808, 963, 962 +459, 564, 563, 720, 719, 808, 807, 964, 963 +460, 563, 562, 721, 720, 807, 806, 965, 964 +461, 562, 561, 722, 721, 806, 805, 966, 965 +462, 561, 560, 723, 722, 805, 804, 967, 966 +463, 560, 559, 724, 723, 804, 803, 968, 967 +464, 559, 558, 725, 724, 803, 802, 969, 968 +465, 558, 557, 726, 725, 802, 801, 970, 969 +466, 557, 149, 148, 726, 801, 196, 197, 970 +467, 557, 727, 150, 149, 801, 971, 195, 196 +468, 727, 728, 151, 150, 971, 972, 194, 195 +469, 728, 534, 152, 151, 972, 778, 193, 194 +470, 558, 729, 727, 557, 802, 973, 971, 801 +471, 729, 730, 728, 727, 973, 974, 972, 971 +472, 730, 535, 534, 728, 974, 779, 778, 972 +473, 559, 731, 729, 558, 803, 975, 973, 802 +474, 731, 732, 730, 729, 975, 976, 974, 973 +475, 732, 536, 535, 730, 976, 780, 779, 974 +476, 560, 733, 731, 559, 804, 977, 975, 803 +477, 733, 734, 732, 731, 977, 978, 976, 975 +478, 734, 537, 536, 732, 978, 781, 780, 976 +479, 561, 735, 733, 560, 805, 979, 977, 804 +480, 735, 736, 734, 733, 979, 980, 978, 977 +481, 736, 538, 537, 734, 980, 782, 781, 978 +482, 562, 737, 735, 561, 806, 981, 979, 805 +483, 737, 738, 736, 735, 981, 982, 980, 979 +484, 738, 539, 538, 736, 982, 783, 782, 980 +485, 563, 739, 737, 562, 807, 983, 981, 806 +486, 739, 740, 738, 737, 983, 984, 982, 981 +487, 740, 540, 539, 738, 984, 784, 783, 982 +488, 564, 741, 739, 563, 808, 985, 983, 807 +489, 741, 742, 740, 739, 985, 986, 984, 983 +490, 742, 541, 540, 740, 986, 785, 784, 984 +491, 565, 743, 741, 564, 809, 987, 985, 808 +492, 743, 744, 742, 741, 987, 988, 986, 985 +493, 744, 542, 541, 742, 988, 786, 785, 986 +494, 566, 745, 743, 565, 810, 989, 987, 809 +495, 745, 746, 744, 743, 989, 990, 988, 987 +496, 746, 543, 542, 744, 990, 787, 786, 988 +497, 567, 747, 745, 566, 811, 991, 989, 810 +498, 747, 748, 746, 745, 991, 992, 990, 989 +499, 748, 544, 543, 746, 992, 788, 787, 990 +500, 568, 749, 747, 567, 812, 993, 991, 811 +501, 749, 750, 748, 747, 993, 994, 992, 991 +502, 750, 545, 544, 748, 994, 789, 788, 992 +503, 556, 582, 749, 568, 800, 826, 993, 812 +504, 582, 581, 750, 749, 826, 825, 994, 993 +505, 581, 546, 545, 750, 825, 790, 789, 994 +506, 533, 751, 178, 177, 777, 995, 223, 224 +507, 751, 752, 179, 178, 995, 996, 222, 223 +508, 752, 580, 180, 179, 996, 824, 221, 222 +509, 532, 753, 751, 533, 776, 997, 995, 777 +510, 753, 754, 752, 751, 997, 998, 996, 995 +511, 754, 579, 580, 752, 998, 823, 824, 996 +512, 531, 755, 753, 532, 775, 999, 997, 776 +513, 755, 756, 754, 753, 999, 1000, 998, 997 +514, 756, 578, 579, 754, 1000, 822, 823, 998 +515, 555, 757, 755, 531, 799, 1001, 999, 775 +516, 757, 758, 756, 755, 1001, 1002, 1000, 999 +517, 758, 577, 578, 756, 1002, 821, 822, 1000 +518, 554, 759, 757, 555, 798, 1003, 1001, 799 +519, 759, 760, 758, 757, 1003, 1004, 1002, 1001 +520, 760, 576, 577, 758, 1004, 820, 821, 1002 +521, 553, 761, 759, 554, 797, 1005, 1003, 798 +522, 761, 762, 760, 759, 1005, 1006, 1004, 1003 +523, 762, 575, 576, 760, 1006, 819, 820, 1004 +524, 552, 763, 761, 553, 796, 1007, 1005, 797 +525, 763, 764, 762, 761, 1007, 1008, 1006, 1005 +526, 764, 574, 575, 762, 1008, 818, 819, 1006 +527, 551, 765, 763, 552, 795, 1009, 1007, 796 +528, 765, 766, 764, 763, 1009, 1010, 1008, 1007 +529, 766, 573, 574, 764, 1010, 817, 818, 1008 +530, 550, 767, 765, 551, 794, 1011, 1009, 795 +531, 767, 768, 766, 765, 1011, 1012, 1010, 1009 +532, 768, 572, 573, 766, 1012, 816, 817, 1010 +533, 549, 769, 767, 550, 793, 1013, 1011, 794 +534, 769, 770, 768, 767, 1013, 1014, 1012, 1011 +535, 770, 571, 572, 768, 1014, 815, 816, 1012 +536, 548, 771, 769, 549, 792, 1015, 1013, 793 +537, 771, 772, 770, 769, 1015, 1016, 1014, 1013 +538, 772, 570, 571, 770, 1016, 814, 815, 1014 +539, 547, 773, 771, 548, 791, 1017, 1015, 792 +540, 773, 774, 772, 771, 1017, 1018, 1016, 1015 +541, 774, 569, 570, 772, 1018, 813, 814, 1016 +542, 546, 581, 773, 547, 790, 825, 1017, 791 +543, 581, 582, 774, 773, 825, 826, 1018, 1017 +544, 582, 556, 569, 774, 826, 800, 813, 1018 +545, 531, 22, 174, 555, 775, 27, 227, 799 +546, 555, 174, 173, 554, 799, 227, 228, 798 +547, 554, 173, 172, 553, 798, 228, 229, 797 +548, 553, 172, 171, 552, 797, 229, 230, 796 +549, 552, 171, 170, 551, 796, 230, 231, 795 +550, 551, 170, 169, 550, 795, 231, 232, 794 +551, 550, 169, 168, 549, 794, 232, 233, 793 +552, 549, 168, 167, 548, 793, 233, 234, 792 +553, 548, 167, 166, 547, 792, 234, 235, 791 +554, 547, 166, 165, 546, 791, 235, 236, 790 +555, 546, 165, 164, 545, 790, 236, 237, 789 +556, 545, 164, 163, 544, 789, 237, 238, 788 +557, 544, 163, 162, 543, 788, 238, 239, 787 +558, 543, 162, 161, 542, 787, 239, 240, 786 +559, 542, 161, 160, 541, 786, 240, 241, 785 +560, 541, 160, 159, 540, 785, 241, 242, 784 +561, 540, 159, 158, 539, 784, 242, 243, 783 +562, 539, 158, 157, 538, 783, 243, 244, 782 +563, 538, 157, 156, 537, 782, 244, 245, 781 +564, 537, 156, 155, 536, 781, 245, 246, 780 +565, 536, 155, 154, 535, 780, 246, 247, 779 +566, 535, 154, 153, 534, 779, 247, 248, 778 +567, 534, 153, 21, 152, 778, 248, 24, 193 +568, 13, 99, 115, 11, 53, 1063, 1040, 51 +569, 99, 98, 116, 115, 1063, 1064, 1039, 1040 +570, 98, 97, 117, 116, 1064, 1065, 1038, 1039 +571, 97, 96, 118, 117, 1065, 1066, 1037, 1038 +572, 96, 95, 119, 118, 1066, 1067, 1036, 1037 +573, 95, 94, 120, 119, 1067, 1068, 1035, 1036 +574, 94, 93, 121, 120, 1068, 1069, 1034, 1035 +575, 93, 92, 122, 121, 1069, 1070, 1033, 1034 +576, 92, 91, 123, 122, 1070, 1071, 1032, 1033 +577, 91, 90, 124, 123, 1071, 1072, 1031, 1032 +578, 90, 89, 125, 124, 1072, 1073, 1030, 1031 +579, 89, 88, 126, 125, 1073, 1074, 1029, 1030 +580, 88, 87, 127, 126, 1074, 1075, 1028, 1029 +581, 87, 86, 128, 127, 1075, 1076, 1027, 1028 +582, 86, 85, 129, 128, 1076, 1077, 1026, 1027 +583, 85, 84, 130, 129, 1077, 1078, 1025, 1026 +584, 84, 83, 131, 130, 1078, 1079, 1024, 1025 +585, 83, 82, 132, 131, 1079, 1080, 1023, 1024 +586, 82, 81, 133, 132, 1080, 1081, 1022, 1023 +587, 81, 80, 134, 133, 1081, 1082, 1021, 1022 +588, 80, 79, 135, 134, 1082, 1083, 1020, 1021 +589, 79, 78, 136, 135, 1083, 1084, 1019, 1020 +590, 78, 18, 19, 136, 1084, 252, 250, 1019 +591, 53, 1063, 1040, 51, 52, 1085, 1062, 50 +592, 1063, 1064, 1039, 1040, 1085, 1086, 1061, 1062 +593, 1064, 1065, 1038, 1039, 1086, 1087, 1060, 1061 +594, 1065, 1066, 1037, 1038, 1087, 1088, 1059, 1060 +595, 1066, 1067, 1036, 1037, 1088, 1089, 1058, 1059 +596, 1067, 1068, 1035, 1036, 1089, 1090, 1057, 1058 +597, 1068, 1069, 1034, 1035, 1090, 1091, 1056, 1057 +598, 1069, 1070, 1033, 1034, 1091, 1092, 1055, 1056 +599, 1070, 1071, 1032, 1033, 1092, 1093, 1054, 1055 +600, 1071, 1072, 1031, 1032, 1093, 1094, 1053, 1054 +601, 1072, 1073, 1030, 1031, 1094, 1095, 1052, 1053 +602, 1073, 1074, 1029, 1030, 1095, 1096, 1051, 1052 +603, 1074, 1075, 1028, 1029, 1096, 1097, 1050, 1051 +604, 1075, 1076, 1027, 1028, 1097, 1098, 1049, 1050 +605, 1076, 1077, 1026, 1027, 1098, 1099, 1048, 1049 +606, 1077, 1078, 1025, 1026, 1099, 1100, 1047, 1048 +607, 1078, 1079, 1024, 1025, 1100, 1101, 1046, 1047 +608, 1079, 1080, 1023, 1024, 1101, 1102, 1045, 1046 +609, 1080, 1081, 1022, 1023, 1102, 1103, 1044, 1045 +610, 1081, 1082, 1021, 1022, 1103, 1104, 1043, 1044 +611, 1082, 1083, 1020, 1021, 1104, 1105, 1042, 1043 +612, 1083, 1084, 1019, 1020, 1105, 1106, 1041, 1042 +613, 1084, 252, 250, 1019, 1106, 251, 249, 1041 +614, 52, 1085, 1062, 50, 12, 253, 296, 10 +615, 1085, 1086, 1061, 1062, 253, 254, 295, 296 +616, 1086, 1087, 1060, 1061, 254, 255, 294, 295 +617, 1087, 1088, 1059, 1060, 255, 256, 293, 294 +618, 1088, 1089, 1058, 1059, 256, 257, 292, 293 +619, 1089, 1090, 1057, 1058, 257, 258, 291, 292 +620, 1090, 1091, 1056, 1057, 258, 259, 290, 291 +621, 1091, 1092, 1055, 1056, 259, 260, 289, 290 +622, 1092, 1093, 1054, 1055, 260, 261, 288, 289 +623, 1093, 1094, 1053, 1054, 261, 262, 287, 288 +624, 1094, 1095, 1052, 1053, 262, 263, 286, 287 +625, 1095, 1096, 1051, 1052, 263, 264, 285, 286 +626, 1096, 1097, 1050, 1051, 264, 265, 284, 285 +627, 1097, 1098, 1049, 1050, 265, 266, 283, 284 +628, 1098, 1099, 1048, 1049, 266, 267, 282, 283 +629, 1099, 1100, 1047, 1048, 267, 268, 281, 282 +630, 1100, 1101, 1046, 1047, 268, 269, 280, 281 +631, 1101, 1102, 1045, 1046, 269, 270, 279, 280 +632, 1102, 1103, 1044, 1045, 270, 271, 278, 279 +633, 1103, 1104, 1043, 1044, 271, 272, 277, 278 +634, 1104, 1105, 1042, 1043, 272, 273, 276, 277 +635, 1105, 1106, 1041, 1042, 273, 274, 275, 276 +636, 1106, 251, 249, 1041, 274, 28, 29, 275 +*Nset, nset=Set-8 + 20, 21, 22, 23, 24, 25, 26, 27, 137, 138, 139, 140, 141, 142, 143, 144 + 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160 + 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176 + 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192 + 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208 + 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224 + 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240 + 241, 242, 243, 244, 245, 246, 247, 248, 531, 532, 533, 534, 535, 536, 537, 538 + 539, 540, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554 + 555, 556, 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570 + 571, 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 586 + 587, 588, 589, 590, 591, 592, 593, 594, 595, 596, 597, 598, 599, 600, 601, 602 + 603, 604, 605, 606, 607, 608, 609, 610, 611, 612, 613, 614, 615, 616, 617, 618 + 619, 620, 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631, 632, 633, 634 + 635, 636, 637, 638, 639, 640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 650 + 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666 + 667, 668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 680, 681, 682 + 683, 684, 685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697, 698 + 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714 + 715, 716, 717, 718, 719, 720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 730 + 731, 732, 733, 734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 746 + 747, 748, 749, 750, 751, 752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762 + 763, 764, 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 775, 776, 777, 778 + 779, 780, 781, 782, 783, 784, 785, 786, 787, 788, 789, 790, 791, 792, 793, 794 + 795, 796, 797, 798, 799, 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, 810 + 811, 812, 813, 814, 815, 816, 817, 818, 819, 820, 821, 822, 823, 824, 825, 826 + 827, 828, 829, 830, 831, 832, 833, 834, 835, 836, 837, 838, 839, 840, 841, 842 + 843, 844, 845, 846, 847, 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 858 + 859, 860, 861, 862, 863, 864, 865, 866, 867, 868, 869, 870, 871, 872, 873, 874 + 875, 876, 877, 878, 879, 880, 881, 882, 883, 884, 885, 886, 887, 888, 889, 890 + 891, 892, 893, 894, 895, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905, 906 + 907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919, 920, 921, 922 + 923, 924, 925, 926, 927, 928, 929, 930, 931, 932, 933, 934, 935, 936, 937, 938 + 939, 940, 941, 942, 943, 944, 945, 946, 947, 948, 949, 950, 951, 952, 953, 954 + 955, 956, 957, 958, 959, 960, 961, 962, 963, 964, 965, 966, 967, 968, 969, 970 + 971, 972, 973, 974, 975, 976, 977, 978, 979, 980, 981, 982, 983, 984, 985, 986 + 987, 988, 989, 990, 991, 992, 993, 994, 995, 996, 997, 998, 999, 1000, 1001, 1002 + 1003, 1004, 1005, 1006, 1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 1018 +*Elset, elset=Set-8, generate + 295, 567, 1 +*Nset, nset=Set-9 + 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 + 17, 18, 19, 20, 21, 22, 23, 28, 29, 30, 31, 32, 33, 34, 35, 36 + 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52 + 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68 + 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84 + 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100 + 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116 + 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132 + 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148 + 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164 + 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180 + 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 249, 250, 251, 252 + 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268 + 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284 + 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300 + 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316 + 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332 + 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348 + 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364 + 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380 + 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396 + 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412 + 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428 + 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444 + 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460 + 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476 + 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492 + 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508 + 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524 + 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540 + 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554, 555, 556 + 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570, 571, 572 + 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 586, 587, 588 + 589, 590, 591, 592, 593, 594, 595, 596, 597, 598, 599, 600, 601, 602, 603, 604 + 605, 606, 607, 608, 609, 610, 611, 612, 613, 614, 615, 616, 617, 618, 619, 620 + 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631, 632, 633, 634, 635, 636 + 637, 638, 639, 640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 650, 651, 652 + 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667, 668 + 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 680, 681, 682, 683, 684 + 685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697, 698, 699, 700 + 701, 702, 703, 704, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 715, 716 + 717, 718, 719, 720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 730, 731, 732 + 733, 734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748 + 749, 750, 751, 752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762, 763, 764 + 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 1019, 1020, 1021, 1022, 1023, 1024 + 1025, 1026, 1027, 1028, 1029, 1030, 1031, 1032, 1033, 1034, 1035, 1036, 1037, 1038, 1039, 1040 + 1041, 1042, 1043, 1044, 1045, 1046, 1047, 1048, 1049, 1050, 1051, 1052, 1053, 1054, 1055, 1056 + 1057, 1058, 1059, 1060, 1061, 1062, 1063, 1064, 1065, 1066, 1067, 1068, 1069, 1070, 1071, 1072 + 1073, 1074, 1075, 1076, 1077, 1078, 1079, 1080, 1081, 1082, 1083, 1084, 1085, 1086, 1087, 1088 + 1089, 1090, 1091, 1092, 1093, 1094, 1095, 1096, 1097, 1098, 1099, 1100, 1101, 1102, 1103, 1104 + 1105, 1106 +*Elset, elset=Set-9 + 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 + 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32 + 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48 + 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64 + 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80 + 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96 + 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112 + 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128 + 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144 + 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160 + 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176 + 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192 + 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208 + 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224 + 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240 + 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256 + 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272 + 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288 + 289, 290, 291, 292, 293, 294, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577 + 578, 579, 580, 581, 582, 583, 584, 585, 586, 587, 588, 589, 590, 591, 592, 593 + 594, 595, 596, 597, 598, 599, 600, 601, 602, 603, 604, 605, 606, 607, 608, 609 + 610, 611, 612, 613, 614, 615, 616, 617, 618, 619, 620, 621, 622, 623, 624, 625 + 626, 627, 628, 629, 630, 631, 632, 633, 634, 635, 636 +** Section: Section-2 +*Solid Section, elset=Set-9, material=Material-2 +, +** Section: Section-1 +*Solid Section, elset=Set-8, material=Material-1 +, +*End Part +** +** +** ASSEMBLY +** +*Assembly, name=Assembly +** +*Instance, name=Part-5-1, part=Part-5 +*End Instance +** +*End Assembly +** +** MATERIALS +** +*Material, name=Material-1 +*Material, name=Material-2 diff --git a/Thermo_electro_elasticity/README.md b/Thermo_electro_elasticity/README.md new file mode 100644 index 00000000..a5e659e2 --- /dev/null +++ b/Thermo_electro_elasticity/README.md @@ -0,0 +1,68 @@ +## Overview +Thermo-electro-elastic pump: implementation of a problem with thermo-electro-mechanical coupling. + +The electro-mechanical and thermal problems are solved in a staggered manner. +Two implementations of the electro-mechanical problem are implemented: +- the standard two-field formulation using the electric scalar potential and a + displacement field (valid for compressible materials), and +- a four-field formulation, extending the above with an additional dilatation + and pressure field (valid for both compressible and near-incompressible + materials). + +Concepts related to this implementation are presented in: + +@Article{Mehnert2017a, + author = {M. Mehnert and J.-P. Pelteret and P. Steinmann}, + title = {Numerical modelling of nonlinear thermo-electro-elasticity}, + journal = {Mathematics and Mechanics of Solids}, + year = {2017}, + volume = {22}, + number = {11}, + pages = {2196--2213}, + month = nov, + doi = {10.1177/1081286517729867}, + publisher = {{SAGE} Publications}, +} + +@Article{Pelteret2016a, + author = {Pelteret, J.-P. V. and Davydov, D. and McBride, A. and Vu, D. K. and Steinmann, P.}, + title = {Computational electro-elasticity and magneto-elasticity for quasi-incompressible media immersed in free space}, + journal = {International Journal for Numerical Methods in Engineering}, + year = {2016}, + volume = {108}, + number = {11}, + pages = {1307--1342}, + month = dec, + doi = {10.1002/nme.5254}, +} + +### Related papers + +@Article{Hamkar2012a, + author = {Hamkar, A-W and Hartmann, Stefan}, + title = {Theoretical and numerical aspects in weak-compressible finite strain thermo-elasticity}, + journal = {Journal of Theoretical and Applied Mechanics}, + year = {2012}, + volume = {50}, + pages = {3--22}, + file = {:Articles/Hamkar2012a.pdf:PDF}, +} + +@article{mehnert2018numerical, + title={Numerical modeling of thermo-electro-viscoelasticity with field-dependent material parameters}, + author={Mehnert, Markus and Hossain, Mokarram and Steinmann, Paul}, + journal={International Journal of Non-Linear Mechanics}, + volume={106}, + pages={13--24}, + year={2018}, + publisher={Elsevier} +} + +@article{mehnert2019experimental, + title={Experimental and numerical investigations of the electro-viscoelastic behavior of VHB 4905TM}, + author={Mehnert, Markus and Hossain, Mokarram and Steinmann, Paul}, + journal={European Journal of Mechanics-A/Solids}, + pages={103797}, + year={2019}, + publisher={Elsevier} +} diff --git a/Thermo_electro_elasticity/doc/author b/Thermo_electro_elasticity/doc/author new file mode 100644 index 00000000..7ebe43a5 --- /dev/null +++ b/Thermo_electro_elasticity/doc/author @@ -0,0 +1,2 @@ +Markus Mehnert +Jean-Paul Pelteret diff --git a/Thermo_electro_elasticity/doc/builds-on b/Thermo_electro_elasticity/doc/builds-on new file mode 100644 index 00000000..5c783f36 --- /dev/null +++ b/Thermo_electro_elasticity/doc/builds-on @@ -0,0 +1 @@ +step-44 code_gallery_Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity diff --git a/Thermo_electro_elasticity/doc/dependencies b/Thermo_electro_elasticity/doc/dependencies new file mode 100644 index 00000000..49815215 --- /dev/null +++ b/Thermo_electro_elasticity/doc/dependencies @@ -0,0 +1,3 @@ +DEAL_II_WITH_MPI +DEAL_II_WITH_TRILINOS +DEAL_II_TRILINOS_WITH_SACADO \ No newline at end of file diff --git a/Thermo_electro_elasticity/doc/entry-name b/Thermo_electro_elasticity/doc/entry-name new file mode 100644 index 00000000..8b0ac3d7 --- /dev/null +++ b/Thermo_electro_elasticity/doc/entry-name @@ -0,0 +1 @@ +Nonlinear thermo-electroelasticity diff --git a/Thermo_electro_elasticity/doc/tooltip b/Thermo_electro_elasticity/doc/tooltip new file mode 100644 index 00000000..eb529094 --- /dev/null +++ b/Thermo_electro_elasticity/doc/tooltip @@ -0,0 +1 @@ +A parallelised nonlinear thermo-electroelastic formulation modelling an electro-active pump, using automatic differentiation and MPI from the Trilinos package. diff --git a/Thermo_electro_elasticity/thermo-electro-elastic_pump.cc b/Thermo_electro_elasticity/thermo-electro-elastic_pump.cc new file mode 100644 index 00000000..d7b4490f --- /dev/null +++ b/Thermo_electro_elasticity/thermo-electro-elastic_pump.cc @@ -0,0 +1,2658 @@ +/* --------------------------------------------------------------------- + * Copyright (C) 2020 by the deal.II authors and + * Markus Mehnert and Jean-Paul Pelteret + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + */ + +/* + * Authors: Markus Mehnert, University of Erlangen-Nuremberg, + * Jean-Paul Pelteret, University of Erlangen-Nuremberg, 2020 + */ + +/** + * ========================= + * THERMO-ELECTRO-ELASTICITY + * ========================= + * Problem description: + * Axial deformation (elongation) of a cylinder, with an electric + * field induced in the axial direction and a temperature gradient + * prescribed between the inner and outer radial surfaces. + * + * Summary of features: + * - Staggered one-way coupling of quasi-static electro-elasticity + * and thermal (conductivity) problems + * - Nonlinear, finite deformation quasi-static elasticity + * - Nonlinear quasi-static thermal (conductivity) problem + * - Nonlinear iterative solution scheme (Newton-Raphson) + * that encompasses staggered thermal / coupled EM solution update + * - Parallelisation via Trilinos (and possibly PETSc) + * - Parallel output of solution, residual + * - Choice of direct and indirect solvers + * - Adaptive grid refinement using Kelly error estimator + * - Parameter collection + * - Generic continuum point framework for integrating constitutive models + * - Nonlinear constitutive models: + * - St. Venant Kirchoff + * + Materially linear thermal conductivity + * - Fully decoupled NeoHookean + * + Materially isotropic dielectric material + * + Spatially isotropic thermal conductivity + * - One-way coupled thermo-electro-mechanical model + * based on Markus' paper (Mehnert2015a) + * + Spatially isotropic thermal conductivity + * + * References: + * Wriggers, P. Nonlinear finite element methods. 2008 + * Holzapfel, G. A. Nonlinear solid mechanics. 2007 + * Vu, K., On coupled BEM-FEM simulation of nonlinear electro-elastostatics + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// #include +#include +// #include +#include +#include +// #include + + +#include +#define USE_TRILINOS_LA +namespace LA +{ +#ifdef USE_TRILINOS_LA +using namespace dealii::LinearAlgebraTrilinos; +#else +using namespace dealii::LinearAlgebraPETSc; +#endif +} + +#include +#include +#include + +namespace Coupled_TEE +{ +using namespace dealii; + + +struct Parameters +{ + // Formulation + static const bool use_3_field_formulation = false; + + // Geometry file + static const std::string mesh_file; + + // Boundary ids + static constexpr unsigned int boundary_id_bottom = 0; + static constexpr unsigned int boundary_id_top = 1; + static constexpr unsigned int boundary_id_inner_radius = 2; + static constexpr unsigned int boundary_id_outer_radius = 3; + static constexpr unsigned int boundary_id_cut_bottom = 4; + static constexpr unsigned int boundary_id_cut_left = 5; + static constexpr unsigned int boundary_id_frame = 6; + static constexpr unsigned int boundary_id_outlet_inner_radius = 7; + static constexpr unsigned int boundary_id_outlet_outer_radius = 8; + static constexpr unsigned int boundary_id_cut_outlet = 9; + + // Boundary conditions + static constexpr double temperature_difference = 20.0; + static constexpr double potential_difference = 0.04; // Potential difference in MV + + // Time + static constexpr double dt = 0.1; + static constexpr unsigned int n_timesteps = 10; + static constexpr double time_end = 50.0e-3; + static constexpr double time_delta = time_end/(static_cast(n_timesteps)); + + // Refinement + static constexpr unsigned int n_global_refinements = 0; + static constexpr bool perform_AMR = false; + static constexpr unsigned int n_ts_per_refinement = 3; + static constexpr unsigned int max_grid_level = 4; + static constexpr double frac_refine = 0.3; + static constexpr double frac_coarsen = 0.03; + + // Finite element + static constexpr unsigned int poly_order = 1; + static constexpr unsigned int fe_index_3_field = 0; + static constexpr unsigned int fe_index_1_field = 1; + + // Nonlinear solver + static constexpr unsigned int max_newton_iterations = 10; + static constexpr double max_res_T_norm = 1e-6; + static constexpr double max_res_EM_norm = 1e-9; + static constexpr double max_res_abs = 1e-6; + + // Linear solver: Thermal + static const std::string solver_type_T; + static constexpr double tol_rel_T = 1e-6; + + // Linear solver: Electro-mechanical + static const std::string solver_type_EM; + static constexpr double tol_rel_EM = 1e-6; +}; + +const std::string Parameters::mesh_file = "Pump_um_coarse.inp"; +const std::string Parameters::solver_type_T = "Direct"; +const std::string Parameters::solver_type_EM = "Direct"; + +namespace Material +{ +struct Coefficients +{ + static constexpr double length_scale = 1.0; + + // Parameters in N, mm, V + + // Elastic parameters + static constexpr double g_0 = 0.1e-6;// in N/um^2 + static constexpr double N = 7.84e5; + + // Electro parameters + static constexpr double epsilon_0 = 8.854187817; // F/m = C/(V*m)= (A*s)/(V*m) = N/(EM*EM) + static constexpr double c_1 = epsilon_0; + static constexpr double c_2 = 2000.0*epsilon_0; + + + // Independent of length and voltage units + static constexpr double nu = (Parameters::use_3_field_formulation ? 0.499 : 0.45); // Poisson ratio + static constexpr double mu = g_0; // Small strain shear modulus + + static constexpr double lambda = 2.0*mu*nu/(1.0-2.0*nu); // Lame constant + static constexpr double kappa = 2.0*mu*(1.0+nu)/(3.0*(1.0-2.0*nu)); // mu = 3*10^5 Pa + + // Thermal parameters + static constexpr double c_0 = 460e6; // specific heat capacity in J/(kg*K) + static constexpr double alpha = 20e-6; //thermal expansion coefficient in 1/K + static constexpr double theta_0 = 293.0; // in K + static constexpr double k = 0.50; // Heat conductivity in N/(s*K) +}; + +template +struct Values_AD +{ + using ad_type = Sacado::Fad::DFad; + + Values_AD (const Tensor<2,dim,ad_type> & F, + const Tensor<1,dim,ad_type> & E, + const Tensor<1,dim> & Grad_T, + const double theta, + const ad_type J_tilde, + const ad_type p, + const double alpha, + const double c_2) + : F (F), + E (E), + Grad_T(Grad_T), + theta(theta), + J_tilde(J_tilde), + p(p), + alpha(alpha), + c_2(c_2), + + C (symmetrize(transpose(F)*F)), + C_inv (symmetrize(invert(static_cast >(C)))), + J (determinant(F)), + C_bar(ad_type(std::pow(J ,-2.0/dim))*C), + C_bar_inv (symmetrize(invert(static_cast >(C_bar)))), + + J_theta (ad_type(std::exp(dim*alpha*(theta-Material::Coefficients::theta_0)))), + F_theta (J_theta*unit_symmetric_tensor()), + F_theta_inv (symmetrize(invert(static_cast >(F_theta)))), + + F_EM (F_theta_inv*F), + J_EM (J/J_theta), + C_EM (symmetrize(transpose(F_EM)*F_EM)), + C_EM_inv (symmetrize(invert(static_cast >(C_EM)))), + + C_EM_bar(ad_type(std::pow(J_EM ,-2.0/dim))*C_EM), + C_EM_bar_inv(symmetrize(invert(static_cast >(C_EM_bar)))), + + + I1 (first_invariant(C)), // tr(C) + I1_EM (first_invariant(C_EM)), // tr(C) + I1_bar (first_invariant(C_bar)), + I1_EM_bar (first_invariant(C_EM_bar)), + I4 (E*E), // [ExE].I + I5 (E*(C*E)), // [ExE].C + I5_bar (E*(C_bar*E)), // [ExE].C + I5_EM_bar (E*(C_EM_bar*E)) // [ExE].C + + {} + + // Directly related to solution field + const Tensor<2,dim,ad_type> F; // Deformation gradient + const Tensor<1,dim,ad_type> E; + const Tensor<1,dim> Grad_T; + const double theta; + const ad_type J_tilde; + const ad_type p; + const double alpha; + const double c_2; + + // Commonly used elastic quantities + const SymmetricTensor<2,dim,ad_type> C; // Right Cauchy-Green deformation tensor + const SymmetricTensor<2,dim,ad_type> C_inv; + const ad_type J; + + const SymmetricTensor<2,dim,ad_type> C_bar; + const SymmetricTensor<2,dim,ad_type> C_bar_inv; + + const ad_type J_theta; + const SymmetricTensor<2,dim,ad_type> F_theta; + const SymmetricTensor<2,dim,ad_type> F_theta_inv; + + const Tensor<2,dim,ad_type> F_EM; + const ad_type J_EM; + const SymmetricTensor<2,dim,ad_type> C_EM; + const SymmetricTensor<2,dim,ad_type> C_EM_inv; + + const SymmetricTensor<2,dim,ad_type> C_EM_bar; + const SymmetricTensor<2,dim,ad_type> C_EM_bar_inv; + + // Invariants + const ad_type I1; + const ad_type I1_EM; + const ad_type I1_bar; + const ad_type I1_EM_bar; + const ad_type I4; + const ad_type I5; + const ad_type I5_bar; + const ad_type I5_EM_bar; + + // === First derivatives === + // --- Mechanical --- + SymmetricTensor<2,dim,ad_type> + dI1_dC () const + { + return unit_symmetric_tensor(); + } + + SymmetricTensor<2,dim,ad_type> + dI1_bar_dC_bar () const + { + return unit_symmetric_tensor(); + } + + + SymmetricTensor<2,dim,ad_type> + dI1_EM_bar_dC_EM_bar () const + { + return unit_symmetric_tensor(); + } + + SymmetricTensor<2,dim,ad_type> + dI4_dC () const + { + return SymmetricTensor<2,dim,ad_type>(); + } + + SymmetricTensor<2,dim,ad_type> + dI4_dC_bar () const + { + return SymmetricTensor<2,dim,ad_type>(); + } + + SymmetricTensor<2,dim,ad_type> + dI5_dC () const + { + const SymmetricTensor<2,dim,ad_type> ExE = symmetrize(outer_product(E,E)); + return ExE; + } + + SymmetricTensor<2,dim,ad_type> + dI5_bar_dC_bar () const + { + const SymmetricTensor<2,dim,ad_type> ExE = symmetrize(outer_product(E,E)); + return ExE; + } + + SymmetricTensor<2,dim,ad_type> + dI5_EM_bar_dC_EM_bar () const + { + const SymmetricTensor<2,dim,ad_type> ExE = symmetrize(outer_product(E,E)); + return ExE; + } + + SymmetricTensor<2,dim,ad_type> + dJ_dC () const + { + return ad_type((0.5*J))*C_inv; + } + + SymmetricTensor<2,dim,ad_type> + dJ_EM_dC_EM () const + { + return ad_type((0.5*J_EM))*C_EM_inv; + } + + SymmetricTensor<4,dim,ad_type> + dC_inv_dC () const + { + SymmetricTensor<4,dim,ad_type> dC_inv_dC; + + for (unsigned int A=0; A + dI4_dE () const + { + return 2.0*E; + } + + Tensor<1,dim,ad_type> + dI5_dE () const + { + return 2.0*(C*E); + } + + Tensor<1,dim,ad_type> + dI5_bar_dE () const + { + return 2.0*(C_bar*E); + } + + Tensor<1,dim,ad_type> + dI5_EM_bar_dE () const + { + return 2.0*(C_EM_bar*E); + } + +}; + +template +struct CM_Base_AD +{ + using ad_type = Sacado::Fad::DFad; + + CM_Base_AD (const Tensor<2,dim,ad_type> & F, + const Tensor<1,dim,ad_type> & E, + const Tensor<1,dim> & Grad_T, + const double theta, + const ad_type J_tilde, + const ad_type p, + const double alpha, + const double c_2) + : values_ad (F,E,Grad_T,theta,J_tilde,p,alpha,c_2) + {} + + virtual ~CM_Base_AD () {} + + // --- Kinematic Quantities --- + + const Values_AD values_ad; + + // --- Kinetic Quantities --- + + // Second Piola-Kirchhoff stress tensor + inline SymmetricTensor<2,dim,ad_type> + get_S_3Field () const + { + const double theta_ratio = values_ad.theta/Material::Coefficients::theta_0; + const ad_type &J_theta = this->values_ad.J_theta; + return ad_type(std::pow(J_theta ,-2.0/dim))*(get_S_iso_3Field()+2.0*theta_ratio*get_dPsi_p_dC_EM()); + } + + inline SymmetricTensor<2,dim,ad_type> + get_S_1Field () const + { + return get_S_iso_1Field() + get_S_vol(); + } + + inline SymmetricTensor<2,dim,ad_type> + get_S_iso_3Field () const + { + return 2.0*get_dPsi_iso_dC_EM(); + } + + inline SymmetricTensor<2,dim,ad_type> + get_S_iso_1Field () const + { + return 2.0*get_dPsi_iso_dC(); + } + + inline SymmetricTensor<2,dim,ad_type> + get_S_vol () const + { + return 2.0*get_dPsi_vol_dC(); + } + + // Referential electric displacement vector + inline Tensor<1,dim,ad_type> + get_D () const + { + return -get_dPsi_dE(); + } + + inline ad_type + get_dPsi_dp () const + { + const ad_type &J_EM = this->values_ad.J_EM; + const ad_type &J_tilde = this->values_ad.J_tilde; + + return J_EM - J_tilde; + } + + inline ad_type + get_dPsi_dJ_tilde () const + { + const ad_type &p = this->values_ad.p; + const ad_type &J_tilde = this->values_ad.J_tilde; + + double kappa=Material::Coefficients::kappa; + + return 0.5* kappa * (J_tilde-1.0/J_tilde) - p; + } + + +protected: + // --- Pure mechanical volumetric response --- + + // Derivative of the volumetric free energy with respect to + // $\widetilde{J}$ return $\frac{\partial + // \Psi_{\text{vol}}(\widetilde{J})}{\partial \widetilde{J}}$ + ad_type + get_dW_vol_elastic_dJ (const double kappa, + const double /*g_0*/) const + { + const ad_type &J = values_ad.J; + return (0.5*(kappa)*(J - 1.0/J)); + } + + // Derivative of the volumetric free energy with respect to + // $\widetilde{J}$ return $\frac{\partial + // \Psi_{\text{vol}}(\widetilde{J})}{\partial \widetilde{J}}$ + ad_type + get_dW_J_dJ (const double kappa, + const double /*g_0*/) const + { + const ad_type &J = values_ad.J; + return (0.5*(kappa)*(J - 1.0/J)); + } + + SymmetricTensor<2,dim,ad_type> + get_dW_J_dC (const double kappa, + const double g_0) const + { + // See Wriggers p46 eqs. 3.123, 3.124; Holzapfel p230 + return get_dW_J_dJ(kappa, g_0)*this->values_ad.dJ_dC(); + } + + ad_type + get_dW_vol_elastic_dJ_EM (const double kappa, + const double /*g_0*/) const + { + const ad_type &J_EM = values_ad.J_EM; + return (0.5*(kappa)*(J_EM - 1.0/J_EM)); + } + + SymmetricTensor<2,dim,ad_type> + get_dW_vol_elastic_dC (const double kappa, + const double g_0) const + { + // See Wriggers p46 eqs. 3.123, 3.124; Holzapfel p230 + return get_dW_vol_elastic_dJ(kappa, g_0)*this->values_ad.dJ_dC(); + } + SymmetricTensor<2,dim,ad_type> + get_dW_vol_elastic_dC_EM (const double kappa, + const double g_0) const + { + // See Wriggers p46 eqs. 3.123, 3.124; Holzapfel p230 + return get_dW_vol_elastic_dJ_EM(kappa, g_0)*this->values_ad.dJ_EM_dC_EM(); + } + + SymmetricTensor<2,dim,ad_type> + get_dPsi_p_dC () const + { + const ad_type &p = values_ad.p; + return p*this->values_ad.dJ_dC(); + } + + SymmetricTensor<2,dim,ad_type> + get_dPsi_p_dC_EM () const + { + const ad_type &p = values_ad.p; + return p*this->values_ad.dJ_EM_dC_EM(); + } + + // --- Coupled mechanical response --- + + virtual SymmetricTensor<2,dim,ad_type> + get_dPsi_vol_dC () const = 0; + + virtual SymmetricTensor<2,dim,ad_type> + get_dPsi_vol_dC_EM () const = 0; + + virtual SymmetricTensor<2,dim,ad_type> + get_dPsi_iso_dC () const = 0; + + virtual SymmetricTensor<2,dim,ad_type> + get_dPsi_iso_dC_EM () const = 0; + + // --- Coupled electric response --- + + virtual Tensor<1,dim,ad_type> + get_dPsi_dE () const = 0; + +}; + +// +//template +//struct CM_Incompressible_Uncoupled_8Chain_AD : public CM_Base_AD +//{ +// using ad_type = Sacado::Fad::DFad; +// CM_Incompressible_Uncoupled_8Chain_AD (const Tensor<2,dim,ad_type> & F, +// const Tensor<1,dim,ad_type> & E, +// const Tensor<1,dim> & Grad_T, +// const double theta, +// const ad_type J_tilde, +// const ad_type p, +// const double alpha, +// const double c_2) +// : CM_Base_AD (F,E,Grad_T,theta,J_tilde,p,alpha,c_2) +// {} +// +// +// virtual ~CM_Incompressible_Uncoupled_8Chain_AD () {} +// +//protected: +// +// virtual SymmetricTensor<2,dim,ad_type> +// get_dPsi_iso_dC () const +// { +// +// double mu=Material::Coefficients::g_0; +// double N=Material::Coefficients::N; +// double c_2=this->values_ad.c_2; +// +// return theta_ratio()*(get_dW_FE_iso_elastic_dC(mu,N,c_2)); +// } +// +// virtual SymmetricTensor<2,dim,ad_type> +// get_dPsi_iso_dC_EM () const +// { +// +// +// double mu=Material::Coefficients::g_0; +// double N=Material::Coefficients::N; +// double c_2=this->values_ad.c_2; +// +// return theta_ratio()*(get_dW_FE_iso_elastic_dC_EM(mu,N,c_2)); +// } +// +// +// virtual SymmetricTensor<2,dim,ad_type> +// get_dPsi_vol_dC () const +// { +// return theta_ratio()*this->get_dW_J_dC(Material::Coefficients::kappa, Material::Coefficients::mu) +// - theta_difference()*get_dM_J_dC(Coefficients::kappa,Coefficients::alpha); // Thermal dilatory response M = M(J); +// } +// +// virtual SymmetricTensor<2,dim,ad_type> +// get_dPsi_vol_dC_EM () const +// { +// return unit_symmetric_tensor(); +// } +// +// inline SymmetricTensor<2,dim,ad_type> +// get_dW_FE_iso_elastic_dC (const double g_0, +// const double N, +// const double c_2) const +// { +// +//// const ad_type &lambda=get_lambda(); +//// +//// return ad_type((0.5*g_0/dim)*((dim*N-lambda*lambda)/(N-lambda*lambda)))* +//// (unit_symmetric_tensor()); +// const SymmetricTensor<2,dim,ad_type> &C_inv = this->values_ad.C_inv; +// return ad_type(g_0)*(unit_symmetric_tensor()-C_inv); +// +// +// } +// +// inline SymmetricTensor<2,dim,ad_type> +// get_dW_FE_iso_elastic_dC_EM (const double g_0, +// const double N, +// const double c_2) const +// { +// +// const SymmetricTensor<4,dim,ad_type> P= get_Dev_P(this->values_ad.F_EM);//get_dC_bar_dC(); +// const SymmetricTensor<2,dim,ad_type> dW_FE_dC_EM_bar=get_dW_FE_iso_elastic_dC_EM_bar(g_0,N,c_2); +// +// return (dW_FE_dC_EM_bar*P); +// +// +// } +// +// inline SymmetricTensor<2,dim,ad_type> +// get_dW_FE_iso_elastic_dC_bar (const double g_0, +// const double N, +// const double c_2) const +// { +//// const ad_type &lambda_bar=get_lambda_bar(); +//// +//// return ad_type((0.5*g_0/dim)*((dim*N-lambda_bar*lambda_bar)/(N-lambda_bar*lambda_bar)))* +//// (unit_symmetric_tensor()); +// +// return g_0 /2.0*this->values_ad.dI1_bar_dC_bar(); +// +// } +// +// inline SymmetricTensor<2,dim,ad_type> +// get_dW_FE_iso_elastic_dC_EM_bar (const double g_0, +// const double N, +// const double c_2) const +// { +//// const ad_type &lambda_bar=get_lambda_bar(); +//// +//// return ad_type((0.5*g_0/dim)*((dim*N-lambda_bar*lambda_bar)/(N-lambda_bar*lambda_bar)))* +//// (unit_symmetric_tensor()); +// +// return g_0 /2.0*this->values_ad.dI1_EM_bar_dC_EM_bar(); +// +// } +// +// inline ad_type +// get_dM_J_dJ (const double kappa, const double alpha) const +// { +// const ad_type &J = this->values_ad.J; +// return dim*kappa*alpha/J; +// } +// +// SymmetricTensor<2,dim,ad_type> +// get_dM_J_dC (const double kappa,const double alpha) const +// { +// // See Wriggers p46 eqs. 3.123, 3.124; Holzapfel p230 +// return get_dM_J_dJ(kappa, alpha)*this->values_ad.dJ_dC(); +// } +// +// +// inline SymmetricTensor<4, dim, ad_type> +// get_Dev_P (const Tensor<2, dim, ad_type> &F) const +// { +// const ad_type det_F = determinant(F); +// Assert(det_F > ad_type(0.0), +// ExcMessage("Deformation gradient has a negative determinant.")); +// const Tensor<2,dim,ad_type> C_ns = transpose(F)*F; +// const SymmetricTensor<2,dim,ad_type> C = symmetrize(C_ns); +// const SymmetricTensor<2,dim,ad_type> C_inv = symmetrize(invert(C_ns)); +// +// // See Wriggers p46 equ 3.125 (but transpose indices) +// SymmetricTensor<4,dim,ad_type> Dev_P = outer_product(C,C_inv); // Dev_P = C_x_C_inv +// Dev_P /= -dim; // Dev_P = -[1/dim]C_x_C_inv +// Dev_P += Physics::Elasticity::StandardTensors< dim >::S; // Dev_P = S - [1/dim]C_x_C_inv +// Dev_P *= ad_type(std::pow(det_F, -2.0/dim)); // Dev_P = J^{-2/dim} [S - [1/dim]C_x_C_inv] +// +// return Dev_P; +// } +// +// +// // --- Electric contributions --- +// virtual Tensor<1,dim,ad_type> +// get_dPsi_dE () const +// { +// double c_1=Material::Coefficients::c_1; +// double c_2=Material::Coefficients::c_2; +// +// +// return get_dW_iso_elastic_dE(c_1,c_2) ; +// } +// +// inline Tensor<1,dim,ad_type> +// get_dW_iso_elastic_dE (const double c_1, +// const double c_2) const +// { +// return c_1*this->values_ad.dI4_dE(); +// } +// +// ad_type +// get_lambda() const +// { +// const ad_type &I1=this->values_ad.I1; +// return std::sqrt(I1/dim); +// } +// +// ad_type +// get_lambda_bar() const +// { +// const ad_type &I1_bar=this->values_ad.I1_bar; +// return std::sqrt(I1_bar/dim); +// } +// +// ad_type +// get_lambda_EM_bar() const +// { +// const ad_type &I1_EM_bar=this->values_ad.I1_EM_bar; +// return std::sqrt(I1_EM_bar/dim); +// } +// +// double +// theta_ratio () const +// { +// return this->values_ad.theta/Material::Coefficients::theta_0; +// } +// +// double +// theta_difference () const +// { +// return this->values_ad.theta - Material::Coefficients::theta_0; +// } +// +//}; + + +template +struct CM_Coupled_NeoHooke_AD : public CM_Base_AD +{ + using ad_type = typename CM_Base_AD::ad_type; + + CM_Coupled_NeoHooke_AD (const Tensor<2,dim,ad_type> & F, + const Tensor<1,dim,ad_type> & E, + const Tensor<1,dim> & Grad_T, + const double theta, + const ad_type J_tilde, + const ad_type p, + const double alpha, + const double c_2) + : CM_Base_AD (F,E,Grad_T,theta,J_tilde,p,alpha,c_2) + {} + + virtual ~CM_Coupled_NeoHooke_AD () {} + +protected: + + virtual SymmetricTensor<2,dim,ad_type> + get_dPsi_iso_dC () const + { + double mu=Material::Coefficients::g_0; + double N=Material::Coefficients::N; + double c_2=this->values_ad.c_2; + + return theta_ratio()*(get_dW_FE_iso_elastic_dC(mu,N,c_2)); + } + + virtual SymmetricTensor<2,dim,ad_type> + get_dPsi_iso_dC_EM () const + { + double mu=Material::Coefficients::g_0; + double N=Material::Coefficients::N; + double c_2=this->values_ad.c_2; + + return theta_ratio()*(get_dW_FE_iso_elastic_dC_EM(mu,N,c_2)); + } + + + virtual SymmetricTensor<2,dim,ad_type> + get_dPsi_vol_dC () const + { + return theta_ratio()*this->get_dW_J_dC(Material::Coefficients::kappa, Material::Coefficients::mu) + - theta_difference()*get_dM_J_dC(Coefficients::kappa,Coefficients::alpha); // Thermal dilatory response M = M(J); + } + + virtual SymmetricTensor<2,dim,ad_type> + get_dPsi_vol_dC_EM () const + { + return unit_symmetric_tensor(); + } + + inline SymmetricTensor<2,dim,ad_type> + get_dW_FE_iso_elastic_dC (const double g_0, + const double /*N*/, + const double c_2) const + { + // const SymmetricTensor<2,dim,ad_type> &C_inv = this->values_ad.C_inv; + return (g_0/2.0)*(unit_symmetric_tensor())+ c_2*this->values_ad.dI5_dC(); + } + + inline SymmetricTensor<2,dim,ad_type> + get_dW_FE_iso_elastic_dC_EM (const double g_0, + const double N, + const double c_2) const + { + + const SymmetricTensor<4,dim,ad_type> P= get_Dev_P(this->values_ad.F_EM);//get_dC_bar_dC(); + const SymmetricTensor<2,dim,ad_type> dW_FE_dC_EM_bar=get_dW_FE_iso_elastic_dC_EM_bar(g_0,N,c_2); + + return (dW_FE_dC_EM_bar*P); + } + + inline SymmetricTensor<2,dim,ad_type> + get_dW_FE_iso_elastic_dC_bar (const double g_0, + const double /*N*/, + const double c_2) const + { + return (g_0 /2.0)*this->values_ad.dI1_bar_dC_bar()+ c_2*this->values_ad.dI5_bar_dC_bar(); + } + + inline SymmetricTensor<2,dim,ad_type> + get_dW_FE_iso_elastic_dC_EM_bar (const double g_0, + const double /*N*/, + const double c_2) const + { + return g_0 /2.0*this->values_ad.dI1_EM_bar_dC_EM_bar()+ c_2*this->values_ad.dI5_EM_bar_dC_EM_bar(); + + } + + inline ad_type + get_dM_J_dJ (const double kappa, const double alpha) const + { + const ad_type &J = this->values_ad.J; + return dim*kappa*alpha/J; + } + + SymmetricTensor<2,dim,ad_type> + get_dM_J_dC (const double kappa,const double alpha) const + { + // See Wriggers p46 eqs. 3.123, 3.124; Holzapfel p230 + return get_dM_J_dJ(kappa, alpha)*this->values_ad.dJ_dC(); + } + + + inline SymmetricTensor<4, dim, ad_type> + get_Dev_P (const Tensor<2, dim, ad_type> &F) const + { + const ad_type det_F = determinant(F); + Assert(det_F > ad_type(0.0), + ExcMessage("Deformation gradient has a negative determinant.")); + const Tensor<2,dim,ad_type> C_ns = transpose(F)*F; + const SymmetricTensor<2,dim,ad_type> C = symmetrize(C_ns); + const SymmetricTensor<2,dim,ad_type> C_inv = symmetrize(invert(C_ns)); + + // See Wriggers p46 equ 3.125 (but transpose indices) + SymmetricTensor<4,dim,ad_type> Dev_P = outer_product(C,C_inv); // Dev_P = C_x_C_inv + Dev_P /= -dim; // Dev_P = -[1/dim]C_x_C_inv + Dev_P += Physics::Elasticity::StandardTensors< dim >::S; // Dev_P = S - [1/dim]C_x_C_inv + Dev_P *= ad_type(std::pow(det_F, -2.0/dim)); // Dev_P = J^{-2/dim} [S - [1/dim]C_x_C_inv] + + return Dev_P; + } + + + // --- Electric contributions --- + virtual Tensor<1,dim,ad_type> + get_dPsi_dE () const + { + double c_1=Material::Coefficients::c_1; + double c_2=this->values_ad.c_2; + + return get_dW_iso_elastic_dE(c_1,c_2) ; + } + + inline Tensor<1,dim,ad_type> + get_dW_iso_elastic_dE (const double c_1, + const double c_2) const + { + return c_1*this->values_ad.dI4_dE()+c_2*this->values_ad.dI5_dE(); + } + + ad_type + get_lambda() const + { + const ad_type &I1=this->values_ad.I1; + return std::sqrt(I1/dim); + } + + ad_type + get_lambda_bar() const + { + const ad_type &I1_bar=this->values_ad.I1_bar; + return std::sqrt(I1_bar/dim); + } + + ad_type + get_lambda_EM_bar() const + { + const ad_type &I1_EM_bar=this->values_ad.I1_EM_bar; + return std::sqrt(I1_EM_bar/dim); + } + + double + theta_ratio () const + { + return this->values_ad.theta/Material::Coefficients::theta_0; + } + + double + theta_difference () const + { + return this->values_ad.theta - Material::Coefficients::theta_0; + } + +}; + + +} + +template +class CoupledProblem +{ + using Continuum_Point_Coupled_NeoHooke_AD = Material::CM_Coupled_NeoHooke_AD; + // using Continuum_Point_8_Chain_uncoupled_AD = Material::CM_Incompressible_Uncoupled_8Chain_AD; + +public: + CoupledProblem (); + ~CoupledProblem (); + void + run (); + +private: + void + make_grid (); + void + set_active_fe_indices(); + void + setup_system (); + void + make_constraints (const unsigned int newton_iteration, const unsigned int timestep); + void + assemble_system_thermo (); + void + solve_thermo (LA::MPI::BlockVector & solution_update); + void + assemble_system_mech (); + void + solve_mech (LA::MPI::BlockVector & solution_update); + void + solve_nonlinear_timestep (const int ts); + void + output_results (const unsigned int timestep) const; + + const unsigned int n_blocks; + const unsigned int first_u_component; // Displacement + const unsigned int V_component; // Voltage / Potential difference + const unsigned int J_component; // Temperature + const unsigned int p_component; // Temperature + const unsigned int T_component; // Temperature + const unsigned int n_components; + + enum + { + EM_block = 0, + T_block = 1 + }; + + enum + { + u_dof = 0, + V_dof = 1, + J_dof = 2, + p_dof = 3, + T_dof = 4 + }; + + enum + { + coupled_material_id = 1, + uncoupled_material_id = 2 + }; + + const FEValuesExtractors::Vector displacement; + const FEValuesExtractors::Scalar x_displacement; + const FEValuesExtractors::Scalar y_displacement; + const FEValuesExtractors::Scalar z_displacement; + const FEValuesExtractors::Scalar voltage; + const FEValuesExtractors::Scalar dilatation; + const FEValuesExtractors::Scalar pressure; + const FEValuesExtractors::Scalar temperature; + + MPI_Comm mpi_communicator; + const unsigned int n_mpi_processes; + const unsigned int this_mpi_process; + mutable ConditionalOStream pcout; + mutable TimerOutput computing_timer; + + parallel::shared::Triangulation triangulation; + hp::FECollection fe_collection; + hp::DoFHandler dof_handler; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + std::vector locally_owned_partitioning; + std::vector locally_relevant_partitioning; + + const unsigned int poly_order; + + hp::QCollection q_collection; + + FESystem fe_cell_3Field; + FESystem fe_cell_1Field; + + QGauss qf_cell_3Field; + QGauss qf_cell_1Field; + + AffineConstraints hanging_node_constraints; + AffineConstraints dirichlet_constraints; + AffineConstraints periodicity_constraints; + AffineConstraints all_constraints; + + LA::MPI::BlockSparseMatrix system_matrix; + LA::MPI::BlockVector system_rhs; + LA::MPI::BlockVector locally_relevant_solution; + LA::MPI::BlockVector locally_relevant_solution_update; + // LA::MPI::BlockVector completely_distributed_solution_update; +}; + +template +CoupledProblem::CoupledProblem () +: +n_blocks (2), +first_u_component (0), // Displacement +V_component (first_u_component + dim), // Voltage / Potential difference +J_component (V_component+1), +p_component (J_component+1), +T_component (p_component+1), // Temperature +n_components (T_component+1), + +displacement(first_u_component), +x_displacement(first_u_component), +y_displacement(first_u_component+1), +z_displacement(dim==3 ? first_u_component+2 : first_u_component+1), +voltage(V_component), +dilatation(J_component), +pressure(p_component), +temperature(T_component), + +mpi_communicator (MPI_COMM_WORLD), +n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)), +this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)), +pcout(std::cout, this_mpi_process == 0), +computing_timer(mpi_communicator, + pcout, + TimerOutput::summary, + TimerOutput::wall_times), +triangulation(mpi_communicator, Triangulation::maximum_smoothing), + +dof_handler(triangulation), + +poly_order (Parameters::poly_order), +fe_cell_3Field(FE_Q (poly_order), dim, + FE_Q (poly_order), 1, // Voltage + FE_DGPMonomial(poly_order - 1), 1, // Dilatation + FE_DGPMonomial(poly_order - 1), 1, // Pressure + FE_Q (poly_order), 1), // Temperature +fe_cell_1Field(FE_Q (poly_order), dim, + FE_Q (poly_order), 1, // Voltage + FE_Nothing(), 1, // Dilatation + FE_Nothing(), 1, // Pressure + FE_Q (poly_order), 1), // Temperature + +qf_cell_3Field(poly_order+1), +qf_cell_1Field(poly_order+1) +{ + fe_collection.push_back(fe_cell_3Field); + fe_collection.push_back(fe_cell_1Field); + q_collection.push_back(qf_cell_3Field); + q_collection.push_back(qf_cell_1Field); +} + + +template +CoupledProblem::~CoupledProblem () +{ + dof_handler.clear(); +} + + +template +void +CoupledProblem::make_grid () //Generate thick walled cylinder +{ + TimerOutput::Scope timer_scope (computing_timer, "Make grid"); + + GridIn grid_in; + grid_in.attach_triangulation (triangulation); + std::ifstream input_file(Parameters::mesh_file.c_str()); + + grid_in.read_abaqus (input_file); + + static CylindricalManifold manifold_cylinder_X (0); // Manifold id 1 + static CylindricalManifold manifold_cylinder_Z (2); // Manifold id 2 + + // Set boundary and manifold ID's for this tricky geometry. + // Note: X-aligned cylinder manifold > Z-aligned cylinder manifold > Straight/Planar manifold + + // Set straight/planar boundary and manifold IDs + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + for (; cell != endc; ++cell) + { + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + { + if (cell->face(f)->at_boundary()) + { + const Point face_center = cell->face(f)->center(); + if (face_center[2] == 120) + { // Faces at cylinder bottom + cell->face(f)->set_boundary_id(Parameters::boundary_id_bottom); + } + else if (face_center[2] == 80) + { // Faces at cylinder top + cell->face(f)->set_boundary_id(Parameters::boundary_id_top); + } + else if (face_center[2] == 0) + { // Faces at cylinder top + cell->face(f)->set_boundary_id(Parameters::boundary_id_frame); + } + else if (face_center[0] == 0) + { // Faces at cylinder top + cell->face(f)->set_boundary_id(Parameters::boundary_id_cut_left); + } + else if (face_center[1] == 0) + { // Faces at cylinder top + cell->face(f)->set_boundary_id(Parameters::boundary_id_cut_bottom); + } + else if (face_center[0] >= 539) + { // Faces at cylinder top + cell->face(f)->set_boundary_id(Parameters::boundary_id_cut_outlet); + } + else + { + // Catch all, mainly for faces at external surface of the cylinder + // Some of the erroneously set boundary ID's will be corrected later. + cell->face(f)->set_boundary_id(Parameters::boundary_id_outer_radius); + } + } + } + } + + // Set Z-aligned cylinder boundary and manifold IDs + cell = triangulation.begin_active(); + for (; cell != endc; ++cell) + { + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + { + if (cell->face(f)->at_boundary()) + { + const Point face_center = cell->face(f)->center(); + if (face_center[0] > 0 && face_center[1] > 0 && face_center[2] > 0) + for (unsigned int v=0; v::vertices_per_face; ++v) + { + const Point &pt = cell->face(f)->vertex(v); + if (abs(sqrt(pt[0]*pt[0] + pt[1]*pt[1]) - 460) < 1e-3 && + face_center[2] < 80) + { + // Faces at internal surface of the cylinder + cell->face(f)->set_boundary_id(Parameters::boundary_id_inner_radius); + cell->face(f)->set_all_manifold_ids(2); + break; + } + else if (abs(sqrt(pt[0]*pt[0] + pt[1]*pt[1]) - 490) < 1e-6) + { + // Faces at external surface of the cylinder + cell->face(f)->set_boundary_id(Parameters::boundary_id_outer_radius); + cell->face(f)->set_all_manifold_ids(2); + break; + } + } + + } + } + } + + // Set X-aligned cylinder boundary and manifold IDs + cell = triangulation.begin_active(); + for (; cell != endc; ++cell) + { + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + { + if (cell->face(f)->at_boundary()) + { + const Point face_center = cell->face(f)->center(); + if (face_center[0] > 0 && face_center[1] > 0 && face_center[2] > 0) // Face not on cartesian planes + if (face_center[0] > 460 && face_center[0] < 540 && + sqrt(face_center[1]*face_center[1] + face_center[2]*face_center[2]) <= 70) + { + for (unsigned int v=0; v::vertices_per_face; ++v) + { + const Point &pt = cell->face(f)->vertex(v); + if (abs(sqrt(pt[1]*pt[1] + pt[2]*pt[2]) - 50) < 1e-6) + { + // Faces at internal surface of the cylinder + cell->face(f)->set_boundary_id(Parameters::boundary_id_outlet_inner_radius); + cell->face(f)->set_all_manifold_ids(1); + break; + } + else if (abs(sqrt(pt[1]*pt[1] + pt[2]*pt[2]) - 70) < 1e-6 && face_center[0] > 490) + { + // Faces at external surface of the cylinder + cell->face(f)->set_boundary_id(Parameters::boundary_id_outlet_outer_radius); + cell->face(f)->set_all_manifold_ids(1); + break; + } + } + } + } + } + } + + triangulation.set_manifold (1, manifold_cylinder_X); + triangulation.set_manifold (2, manifold_cylinder_Z); + + triangulation.refine_global (Parameters::n_global_refinements); +} + + +template +void +CoupledProblem::set_active_fe_indices() +{ + for (typename hp::DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(); + cell != dof_handler.end(); + ++cell) + { + if (!cell->is_locally_owned()) + continue; + + if (Parameters::use_3_field_formulation) + cell->set_active_fe_index(Parameters::fe_index_3_field); + else if (Parameters::use_3_field_formulation==false) + cell->set_active_fe_index(Parameters::fe_index_1_field); + else + Assert(false, ExcNotImplemented()); + } +} + + +template +void +CoupledProblem::setup_system () +{ + TimerOutput::Scope timer_scope (computing_timer, "System setup"); + + pcout << "Setting up the thermo-electro-mechanical system..." << std::endl; + set_active_fe_indices(); + dof_handler.distribute_dofs(fe_collection); + + std::vector block_component(n_components, EM_block); // Displacement + block_component[V_component] = EM_block; // Voltage + block_component[J_component] = EM_block; // Dilatation + block_component[p_component] = EM_block; // Pressure + block_component[T_component] = T_block; // Temperature + + DoFRenumbering::Cuthill_McKee(dof_handler); + DoFRenumbering::component_wise(dof_handler, block_component); + + std::vector dofs_per_block(n_blocks); + DoFTools::count_dofs_per_block(dof_handler, dofs_per_block, block_component); + const types::global_dof_index &n_u_V = dofs_per_block[0]; + const types::global_dof_index &n_th = dofs_per_block[1]; + + pcout + << "Number of active cells: " + << triangulation.n_active_cells() + << std::endl + << "Total number of cells: " + << triangulation.n_cells() + << std::endl + << "Number of degrees of freedom: " + << dof_handler.n_dofs() + << " (" << n_u_V << '+' << n_th << ')' + << std::endl; + + locally_owned_dofs.clear(); + locally_relevant_dofs.clear(); + locally_owned_dofs = dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + locally_owned_partitioning.clear(); + locally_relevant_partitioning.clear(); + locally_owned_partitioning.reserve(n_blocks); + locally_relevant_partitioning.reserve(n_blocks); + for (unsigned int b=0; b +void +CoupledProblem::make_constraints (const unsigned int newton_iteration, const unsigned int timestep) +{ + TimerOutput::Scope timer_scope (computing_timer, "Make constraints"); + + if (newton_iteration >= 2) + { + pcout << std::string(14, ' ') << std::flush; + return; + } + if (newton_iteration == 0) + { + dirichlet_constraints.clear(); + dirichlet_constraints.reinit (locally_relevant_dofs); + + pcout << " CST T" << std::flush; + + // const double temperature_difference_per_ts = Parameters::temperature_difference/static_cast(Parameters::n_timesteps); + if (timestep==1) + { + // Prescribed temperature at inner radius + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_inner_radius, + ConstantFunction(Material::Coefficients::theta_0+Parameters::temperature_difference,n_components), + dirichlet_constraints, + fe_collection.component_mask(temperature)); + + // Prescribed temperature at inner radius + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_outlet_inner_radius, + ConstantFunction(Material::Coefficients::theta_0+Parameters::temperature_difference,n_components), + dirichlet_constraints, + fe_collection.component_mask(temperature)); + + // Prescribed temperature at top + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_top, + ConstantFunction(Material::Coefficients::theta_0+Parameters::temperature_difference,n_components), + dirichlet_constraints, + fe_collection.component_mask(temperature)); + + // Prescribed temperature at bottom + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_bottom, + ConstantFunction(Material::Coefficients::theta_0,n_components), + dirichlet_constraints, + fe_collection.component_mask(temperature)); + + // Prescribed temperature at outer radius + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_outer_radius, + ConstantFunction(Material::Coefficients::theta_0,n_components), + dirichlet_constraints, + fe_collection.component_mask(temperature)); + + // Prescribed temperature at outer radius + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_outlet_outer_radius, + ConstantFunction(Material::Coefficients::theta_0,n_components), + dirichlet_constraints, + fe_collection.component_mask(temperature)); + } + else + { + // Prescribed temperature at top + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_top, + ZeroFunction(n_components), + dirichlet_constraints, + fe_collection.component_mask(temperature)); + + // Prescribed temperature at bottom + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_bottom, + ZeroFunction(n_components), + dirichlet_constraints, + fe_collection.component_mask(temperature)); + + // Prescribed temperature at inner radius + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_inner_radius, + ZeroFunction(n_components), + dirichlet_constraints, + fe_collection.component_mask(temperature)); + + // Prescribed temperature at outer radius + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_outer_radius, + ZeroFunction(n_components), + dirichlet_constraints, + fe_collection.component_mask(temperature)); + + // Prescribed temperature at outer radius + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_outlet_outer_radius, + ZeroFunction(n_components), + dirichlet_constraints, + fe_collection.component_mask(temperature)); + + // Prescribed temperature at inner radius + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_outlet_inner_radius, + ZeroFunction(n_components), + dirichlet_constraints, + fe_collection.component_mask(temperature)); + } + + + pcout << " CST M" << std::flush; + { + const double potential_difference_per_ts = Parameters::potential_difference/(static_cast(Parameters::n_timesteps)); + + // Y-Cut Surface + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_cut_bottom, + ZeroFunction(n_components), + dirichlet_constraints, + fe_collection.component_mask(y_displacement)); + + // X-Cut Surface + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_cut_left, + ZeroFunction(n_components), + dirichlet_constraints, + fe_collection.component_mask(x_displacement)); + + // Frame Cut Surface + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_frame, + ZeroFunction(n_components), + dirichlet_constraints, + fe_collection.component_mask(z_displacement)); + + // Frame Cut Surface + /* VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_outer_radius, + ZeroFunction(n_components), + dirichlet_constraints, + fe_cell.component_mask(x_displacement)|fe_cell.component_mask(y_displacement)); + */ + + // Fixed outlet + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_cut_outlet, + ZeroFunction(n_components), + dirichlet_constraints, + fe_collection.component_mask(x_displacement) | + fe_collection.component_mask(y_displacement) | + fe_collection.component_mask(z_displacement)); + + // Prescribed voltage at lower surface + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_top, + //ZeroFunction(n_components), + ConstantFunction(+potential_difference_per_ts/2,n_components), + dirichlet_constraints, + fe_collection.component_mask(voltage)); + + // Prescribed voltage at upper surface + VectorTools::interpolate_boundary_values(dof_handler, + Parameters::boundary_id_bottom, + //ZeroFunction(n_components), + ConstantFunction(-potential_difference_per_ts/2,n_components), + dirichlet_constraints, + fe_collection.component_mask(voltage)); + } + + dirichlet_constraints.close(); + } + else + { + pcout << " CST ZERO " << std::flush; + // Remove inhomogenaities + for (types::global_dof_index d=0; d::left_object_wins); + all_constraints.close(); +} + +template +void +CoupledProblem::assemble_system_thermo() +{ + using ad_type = Sacado::Fad::DFad; + + TimerOutput::Scope timer_scope (computing_timer, "Assembly: Thermal"); + pcout << " ASM T" << std::flush; + + hp::FEValues hp_fe_values(fe_collection, + q_collection, + update_values | update_quadrature_points | + update_JxW_values | update_gradients); + + FullMatrix cell_matrix; + Vector cell_rhs; + + std::vector local_dof_indices; + + typename hp::DoFHandler::active_cell_iterator cell = dof_handler + .begin_active(), + endc = dof_handler.end(); + for (; cell != endc; ++cell) + { + if (cell->is_locally_owned() == false) continue; + + hp_fe_values.reinit(cell); + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + const unsigned int n_q_points = fe_values.n_quadrature_points; + + cell_matrix.reinit(cell->get_fe().dofs_per_cell, + cell->get_fe().dofs_per_cell); + cell_rhs.reinit(cell->get_fe().dofs_per_cell); + cell_matrix = 0; + cell_rhs = 0; + + local_dof_indices.resize(cell->get_fe().dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + + // Values at integration points + std::vector< Tensor<2,dim> > Grad_u(n_q_points); // Material gradient of displacement + std::vector< Tensor<1,dim> > Grad_V(n_q_points); // Material gradient of voltage + std::vector< Tensor<1,dim> > Grad_T(n_q_points); // Material gradient of temperature + std::vector theta(n_q_points); // Temperature + + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double &JxW = fe_values.JxW(q_point); + + fe_values[displacement].get_function_gradients(locally_relevant_solution, Grad_u); + fe_values[voltage].get_function_gradients(locally_relevant_solution, Grad_V); + fe_values[temperature].get_function_gradients(locally_relevant_solution, Grad_T); + fe_values[temperature].get_function_values(locally_relevant_solution, theta); + + //Deformation gradient at quadrature point + const Tensor<2,dim> F_q_point = (static_cast< Tensor<2,dim> >(unit_symmetric_tensor()) + Grad_u[q_point]); + + const Tensor<2,dim> F_inv = invert(F_q_point); + Tensor<2,dim> K; + + + + K=Material::Coefficients::k*symmetrize(F_inv*transpose(F_inv)); + + + const Tensor<1,dim> Q = K*(-Grad_T[q_point]); + + for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i) + { + const unsigned int i_group = fe_cell_3Field.system_to_base_index(i).first.first; + + const Tensor<1,dim> &Grad_Nx_i_T = fe_values[temperature].gradient(i, q_point); + + + for (unsigned int j = 0; j < cell->get_fe().dofs_per_cell; ++j) + { + + const unsigned int j_group = fe_cell_3Field.system_to_base_index(j).first.first; + + const Tensor<1,dim> &Grad_Nx_j_T = fe_values[temperature].gradient(j, q_point); + + + if ((i_group == T_dof) && (j_group == T_dof)) + { + // T-T terms + cell_matrix(i, j) -= (Grad_Nx_i_T*K*Grad_Nx_j_T) * JxW; + } + } + + // RHS = -Residual + if (i_group == T_dof) + { + // T terms + cell_rhs(i) -= (Grad_Nx_i_T*Q) * JxW; + + } + } + } + + + all_constraints.distribute_local_to_global(cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } + + system_matrix.compress (VectorOperation::add); + system_rhs.compress (VectorOperation::add); +} + +template +void +CoupledProblem::assemble_system_mech () + +{ + using ad_type = Sacado::Fad::DFad; + + TimerOutput::Scope timer_scope (computing_timer, "Assembly: Mechanical"); + pcout << " ASM M" << std::flush; + + hp::FEValues hp_fe_values(fe_collection, + q_collection, + update_values | update_quadrature_points | + update_JxW_values | update_gradients); + + + FullMatrix cell_matrix; + Vector cell_rhs; + + std::vector local_dof_indices; + + typename hp::DoFHandler::active_cell_iterator cell = dof_handler + .begin_active(), + endc = dof_handler.end(); + for (; cell != endc; ++cell) + { + if (cell->is_locally_owned() == false) continue; + + hp_fe_values.reinit(cell); + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + const unsigned int n_q_points = fe_values.n_quadrature_points; + cell_matrix.reinit(cell->get_fe().dofs_per_cell, + cell->get_fe().dofs_per_cell); + cell_rhs.reinit(cell->get_fe().dofs_per_cell); + + cell_matrix = 0; + cell_rhs = 0; + + local_dof_indices.resize(cell->get_fe().dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + + + const types::material_id mat_id = cell->material_id(); + + { + const unsigned int n_independent_variables = local_dof_indices.size(); + std::vector local_dof_values(n_independent_variables); + cell->get_dof_values(locally_relevant_solution, local_dof_values.begin(), local_dof_values.end()); + + // We now retreive a set of degree-of-freedom values that + // have the operations that are performed with them tracked. + std::vector local_dof_values_ad (n_independent_variables); + for (unsigned int i=0; i > Grad_u_ad (n_q_points, Tensor<2,dim,ad_type>()); + fe_values[displacement].get_function_gradients_from_local_dof_values(local_dof_values_ad, Grad_u_ad); + std::vector< Tensor<1,dim,ad_type> > Grad_V_ad (n_q_points, Tensor<1,dim,ad_type>()); + fe_values[voltage].get_function_gradients_from_local_dof_values(local_dof_values_ad, Grad_V_ad); + std::vector< ad_type> J_tilde (n_q_points, 1.0); + fe_values[dilatation].get_function_values_from_local_dof_values(local_dof_values_ad, J_tilde); + std::vector< ad_type> p (n_q_points, 1.0); + fe_values[pressure].get_function_values_from_local_dof_values(local_dof_values_ad, p); + std::vector theta(n_q_points); // Temperature + fe_values[temperature].get_function_values(locally_relevant_solution, theta); + std::vector< Tensor<1,dim> > Grad_T(n_q_points); + fe_values[temperature].get_function_gradients(locally_relevant_solution, Grad_T); + + std::vector cell_residual_ad(cell->get_fe().dofs_per_cell, ad_type(0.0)); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const Tensor<2,dim,ad_type> F_ad = Physics::Elasticity::Kinematics::F(Grad_u_ad[q_point]); + + SymmetricTensor<2,dim,ad_type> S_ad; + Tensor<1,dim,ad_type> D_ad; + ad_type dPsi_dJ_tilde; + ad_type dPsi_dp; + + const double JxW = fe_values.JxW(q_point); + + const double alpha = Material::Coefficients::alpha; + const double c_2 = Material::Coefficients::c_2; + + if (Parameters::use_3_field_formulation) + { + if(mat_id==coupled_material_id) + { + const Continuum_Point_Coupled_NeoHooke_AD cp_ad (F_ad, -Grad_V_ad[q_point],Grad_T[q_point], theta[q_point],J_tilde[q_point],p[q_point],alpha,c_2); + S_ad=cp_ad.get_S_3Field(); + D_ad=cp_ad.get_D(); + dPsi_dJ_tilde= cp_ad.get_dPsi_dJ_tilde(); + dPsi_dp= cp_ad.get_dPsi_dp(); + } + else + { + const Continuum_Point_Coupled_NeoHooke_AD cp_ad (F_ad, -Grad_V_ad[q_point],Grad_T[q_point], theta[q_point],J_tilde[q_point],p[q_point],alpha,0.0); + S_ad=cp_ad.get_S_3Field(); + D_ad=cp_ad.get_D(); + dPsi_dJ_tilde= cp_ad.get_dPsi_dJ_tilde(); + dPsi_dp= cp_ad.get_dPsi_dp(); + } + + for(unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i) + { + const unsigned int i_group = fe_cell_3Field.system_to_base_index(i).first.first; + + if (i_group == u_dof) + { + const SymmetricTensor<2,dim,ad_type> dE_ad_I = symmetrize(transpose(F_ad)*fe_values[displacement].gradient(i, q_point)); + cell_residual_ad[i] += (dE_ad_I*S_ad) * JxW; // residual + + } + else if (i_group == V_dof) + { + const Tensor<1,dim> &Grad_Nx_i_V = fe_values[voltage].gradient(i, q_point); + cell_residual_ad[i] -= (Grad_Nx_i_V*D_ad) * JxW; + } + else if (i_group == J_dof) + { + const double &NJ_i_value = fe_values[dilatation].value(i,q_point); + cell_residual_ad[i] -= NJ_i_value * dPsi_dJ_tilde * JxW; + } + else if (i_group == p_dof) + { + const double &Np_i_value = fe_values[pressure].value(i,q_point); + cell_residual_ad[i] -= Np_i_value * dPsi_dp * JxW; + } + } + } + else + { + const ad_type det_F_qp = determinant(F_ad); + const ad_type p_qp = p[q_point]; // TODO: FixME! 1/dim trace(sigma) + std::cout << "det_F_qp: " << det_F_qp << " p_qp: " << p_qp << std::endl; + + if(mat_id==coupled_material_id) + { + const Continuum_Point_Coupled_NeoHooke_AD cp_ad (F_ad, -Grad_V_ad[q_point],Grad_T[q_point], theta[q_point],det_F_qp,p_qp,alpha,c_2); + S_ad=cp_ad.get_S_1Field(); + D_ad=cp_ad.get_D(); + } + else + { + const Continuum_Point_Coupled_NeoHooke_AD cp_ad (F_ad, -Grad_V_ad[q_point],Grad_T[q_point], theta[q_point],det_F_qp,p_qp,alpha,0.0); + S_ad=cp_ad.get_S_1Field(); + D_ad=cp_ad.get_D(); + } + + for(unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i) + { + const unsigned int i_group = fe_cell_1Field.system_to_base_index(i).first.first; + + if (i_group == u_dof) + { + const SymmetricTensor<2,dim,ad_type> dE_ad_I = symmetrize(transpose(F_ad)*fe_values[displacement].gradient(i, q_point)); + cell_residual_ad[i] += (dE_ad_I*S_ad) * JxW; // residual + + } + else if (i_group == V_dof) + { + const Tensor<1,dim> &Grad_Nx_i_V = fe_values[voltage].gradient(i, q_point); + cell_residual_ad[i] -= (Grad_Nx_i_V*D_ad) * JxW; + } + } + } + + } + + if (Parameters::use_3_field_formulation) + { + for (unsigned int I=0; I +void +CoupledProblem::solve_thermo (LA::MPI::BlockVector &locally_relevant_solution_update) +{ + TimerOutput::Scope timer_scope (computing_timer, "Solve: Thermal"); + pcout << " SLV T" << std::flush; + + // const std::string solver_type = "Iterative"; + const std::string solver_type = "Direct"; + + LA::MPI::BlockVector + completely_distributed_solution_update (locally_owned_partitioning, + mpi_communicator); + + { // Direct solver +#ifdef USE_TRILINOS_LA + SolverControl solver_control(1, 1e-12); + TrilinosWrappers::SolverDirect solver (solver_control); + + + solver.solve(system_matrix.block(T_block, T_block), + completely_distributed_solution_update.block(T_block), + system_rhs.block(T_block)); +#else + AssertThrow(false, ExcNotImplemented()); +#endif + } + + all_constraints.distribute(completely_distributed_solution_update); + locally_relevant_solution_update.block(T_block) = completely_distributed_solution_update.block(T_block); +} + + +template +void +CoupledProblem::solve_mech (LA::MPI::BlockVector &locally_relevant_solution_update) +{ + TimerOutput::Scope timer_scope (computing_timer, "Solve: Mechanical"); + pcout << " SLV M" << std::flush; + + LA::MPI::BlockVector + completely_distributed_solution_update (locally_owned_partitioning, + mpi_communicator); + + { // Direct solver +#ifdef USE_TRILINOS_LA + SolverControl solver_control(1, 1e-12); + TrilinosWrappers::SolverDirect solver (solver_control); + + solver.solve(system_matrix.block(EM_block, EM_block), + completely_distributed_solution_update.block(EM_block), + system_rhs.block(EM_block)); +#else + AssertThrow(false, ExcNotImplemented()); +#endif + } + + all_constraints.distribute(completely_distributed_solution_update); + locally_relevant_solution_update.block(EM_block) = completely_distributed_solution_update.block(EM_block); + +} + +template +void +CoupledProblem::output_results (const unsigned int timestep) const +{ + using ad_type = Sacado::Fad::DFad; + TimerOutput::Scope timer_scope (computing_timer, "Post-processing"); + + unsigned int scalar_components; + unsigned int vector_components; + unsigned int tensor_components; + + switch(dim) + { + case 2: + { + scalar_components = 1; + vector_components = 2; + tensor_components = 4; + } + break; + + case 3: + { + scalar_components = 1; + vector_components = 3; + tensor_components = 9; + } + break; + default: + { + AssertThrow(false, ExcMessage("Magnetoelasticity::Output_Results -> vector_components -> dimension missmatch!")); + } + break; + } + + std::vector > deformation_gradient (tensor_components, Vector (triangulation.n_active_cells())); + std::vector > Piola_Kirchhoff (tensor_components, Vector (triangulation.n_active_cells())); + + std::vector > electric_field (vector_components, Vector (triangulation.n_active_cells())); + std::vector > electric_displacement (vector_components, Vector (triangulation.n_active_cells())); + + // FEValues fe_values(fe_cell, + // qf_cell, + // update_values | + // update_gradients | + // update_quadrature_points | + // update_JxW_values); + + hp::FEValues hp_fe_values(fe_collection, + q_collection, + update_values | update_quadrature_points | + update_JxW_values | update_gradients); + + std::vector local_dof_indices; + + typename hp::DoFHandler::active_cell_iterator cell = dof_handler + .begin_active(), + endc = dof_handler.end(); + for (unsigned int i = 0; cell != endc; ++cell, ++i) + { + if (cell->is_locally_owned() == false) continue; + + hp_fe_values.reinit(cell); + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + const unsigned int n_q_points = fe_values.n_quadrature_points; + local_dof_indices.resize(cell->get_fe().dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + + const types::material_id mat_id = cell->material_id(); + + { + const unsigned int n_independent_variables = local_dof_indices.size(); + std::vector local_dof_values(n_independent_variables); + cell->get_dof_values(locally_relevant_solution, local_dof_values.begin(), local_dof_values.end()); + + // We now retreive a set of degree-of-freedom values that + // have the operations that are performed with them tracked. + std::vector local_dof_values_ad (n_independent_variables); + for (unsigned int i=0; i > Grad_u_ad (n_q_points, Tensor<2,dim,ad_type>()); + fe_values[displacement].get_function_gradients_from_local_dof_values(local_dof_values_ad, Grad_u_ad); + std::vector< Tensor<1,dim,ad_type> > Grad_V_ad (n_q_points, Tensor<1,dim,ad_type>()); + fe_values[voltage].get_function_gradients_from_local_dof_values(local_dof_values_ad, Grad_V_ad); + std::vector< ad_type> J_tilde (n_q_points, 1.0); + fe_values[dilatation].get_function_values_from_local_dof_values(local_dof_values_ad, J_tilde); + std::vector< ad_type> p (n_q_points, 1.0); + fe_values[pressure].get_function_values_from_local_dof_values(local_dof_values_ad, p); + std::vector theta(n_q_points); // Temperature + fe_values[temperature].get_function_values(locally_relevant_solution, theta); + std::vector< Tensor<1,dim> > Grad_T(n_q_points); + fe_values[temperature].get_function_gradients(locally_relevant_solution, Grad_T); + + + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const Tensor<2,dim,ad_type> F_ad = Physics::Elasticity::Kinematics::F(Grad_u_ad[q_point]); + SymmetricTensor<2,dim,ad_type> S_ad; + Tensor<1,dim,ad_type> D_ad; + ad_type dPsi_dJ_tilde; + ad_type dPsi_dp; + + + const double alpha = Material::Coefficients::alpha; + const double c_2 = Material::Coefficients::c_2; + + if(mat_id==coupled_material_id) + { + + const Continuum_Point_Coupled_NeoHooke_AD cp_ad (F_ad, -Grad_V_ad[q_point],Grad_T[q_point], theta[q_point],J_tilde[q_point],p[q_point],alpha,c_2); + S_ad=cp_ad.get_S_3Field(); + D_ad=cp_ad.get_D(); + dPsi_dJ_tilde= cp_ad.get_dPsi_dJ_tilde(); + dPsi_dp= cp_ad.get_dPsi_dp(); + } + else + { + const Continuum_Point_Coupled_NeoHooke_AD cp_ad (F_ad, -Grad_V_ad[q_point],Grad_T[q_point], theta[q_point],J_tilde[q_point],p[q_point],alpha,0.0); + S_ad=cp_ad.get_S_3Field(); + D_ad=cp_ad.get_D(); + dPsi_dJ_tilde= cp_ad.get_dPsi_dJ_tilde(); + dPsi_dp= cp_ad.get_dPsi_dp(); + } + + + deformation_gradient[0][i] += (F_ad[0][0].val())/n_q_points; //F_xx + deformation_gradient[1][i] += (F_ad[0][1].val())/n_q_points; //F_xy + deformation_gradient[2][i] += (F_ad[1][0].val())/n_q_points; //F_yx + deformation_gradient[3][i] += (F_ad[1][1].val())/n_q_points; //F_yy + + Piola_Kirchhoff[0][i] += (S_ad[0][0].val())/n_q_points; //S_xx + Piola_Kirchhoff[1][i] += (S_ad[0][1].val())/n_q_points; //S_xy + Piola_Kirchhoff[2][i] += (S_ad[1][0].val())/n_q_points; //S_yx + Piola_Kirchhoff[3][i] += (S_ad[1][1].val())/n_q_points; //S_yy + + electric_field[0][i] += (-Grad_V_ad[q_point][0].val())/n_q_points; //E_x + electric_field[1][i] += (-Grad_V_ad[q_point][1].val())/n_q_points; //E_y + + electric_displacement[0][i] += (D_ad[0].val())/n_q_points; //D_x + electric_displacement[1][i] += (D_ad[1].val())/n_q_points; //D_x + + + } + } + + } + + // Write out main data file + struct Filename + { + static std::string get_filename_vtu (unsigned int process, + unsigned int cycle, + const unsigned int n_digits = 4) + { + std::ostringstream filename_vtu; + filename_vtu + << "solution-" + << (std::to_string(dim) + "d") + << "." + << Utilities::int_to_string (process, n_digits) + << "." + << Utilities::int_to_string(cycle, n_digits) + << ".vtu"; + return filename_vtu.str(); + } + + static std::string get_filename_pvtu (unsigned int timestep, + const unsigned int n_digits = 4) + { + std::ostringstream filename_vtu; + filename_vtu + << "solution-" + << (std::to_string(dim) + "d") + << "." + << Utilities::int_to_string(timestep, n_digits) + << ".pvtu"; + return filename_vtu.str(); + } + + static std::string get_filename_pvd (void) + { + std::ostringstream filename_vtu; + filename_vtu + << "solution-" + << (std::to_string(dim) + "d") + << ".pvd"; + return filename_vtu.str(); + } + }; + + // DataOut data_out; + DataOut> data_out; + data_out.attach_dof_handler (dof_handler); + + std::vector solution_names (n_components, "displacement"); + solution_names[V_component] = "voltage"; + solution_names[J_component] = "dilatation"; + solution_names[p_component] = "pressure"; + solution_names[T_component] = "temperature"; + + std::vector residual_names (solution_names); + + std::vector reaction_forces_names (solution_names); + + for (unsigned int i=0; i < solution_names.size(); ++i) + { + solution_names[i].insert(0, "soln_"); + residual_names[i].insert(0, "res_"); + } + + std::vector + data_component_interpretation(dim, + DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar); + data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar); + data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar); + data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar); + + data_out.add_data_vector(locally_relevant_solution, solution_names, + DataOut>::type_dof_data, + data_component_interpretation); + + data_out.add_data_vector (deformation_gradient[0], "F_xx"); + data_out.add_data_vector (deformation_gradient[1], "F_xy"); + data_out.add_data_vector (deformation_gradient[2], "F_yx"); + data_out.add_data_vector (deformation_gradient[3], "F_yy"); + + data_out.add_data_vector (Piola_Kirchhoff[0], "S_xx"); + data_out.add_data_vector (Piola_Kirchhoff[1], "S_xy"); + data_out.add_data_vector (Piola_Kirchhoff[2], "S_yx"); + data_out.add_data_vector (Piola_Kirchhoff[3], "S_yy"); + + data_out.add_data_vector (electric_field[0], "E_x"); + data_out.add_data_vector (electric_field[1], "E_y"); + + data_out.add_data_vector (electric_displacement[0], "D_x"); + data_out.add_data_vector (electric_displacement[1], "D_y"); + + + data_out.build_patches (poly_order); + + const std::string filename_vtu = Filename::get_filename_vtu(this_mpi_process, + timestep); + std::ofstream output(filename_vtu.c_str()); + data_out.write_vtu(output); + + // Collection of files written in parallel + // This next set of steps should only be performed + // by master process + if (this_mpi_process == 0) + { + // List of all files written out at this timestep by all processors + std::vector parallel_filenames_vtu; + for (unsigned int p=0; p < n_mpi_processes; ++p) + { + parallel_filenames_vtu.push_back(Filename::get_filename_vtu(p, + timestep)); + } + + const std::string filename_pvtu (Filename::get_filename_pvtu(timestep)); + std::ofstream pvtu_master(filename_pvtu.c_str()); + data_out.write_pvtu_record(pvtu_master, + parallel_filenames_vtu); + + // Time dependent data master file + static std::vector > time_and_name_history; + time_and_name_history.push_back (std::make_pair (timestep, + filename_pvtu)); + const std::string filename_pvd (Filename::get_filename_pvd()); + std::ofstream pvd_output (filename_pvd.c_str()); + DataOutBase::write_pvd_record (pvd_output, time_and_name_history); + } +} + + + +struct L2_norms +{ + L2_norms (const unsigned int block, + const std::vector &locally_owned_partitioning, + const std::vector &locally_relevant_partitioning, + const MPI_Comm &mpi_communicator) + : block (block), + locally_owned_partitioning (locally_owned_partitioning), + locally_relevant_partitioning (locally_relevant_partitioning), + mpi_communicator (mpi_communicator) + {} + + const unsigned int block; + const std::vector &locally_owned_partitioning; + const std::vector &locally_relevant_partitioning; + const MPI_Comm &mpi_communicator; + + double value = 1.0; + double value_norm = 1.0; + + void + set (const LA::MPI::BlockVector & vector, + const AffineConstraints & all_constraints) + { + LA::MPI::BlockVector vector_zeroed; + vector_zeroed.reinit (locally_owned_partitioning, + locally_relevant_partitioning, + mpi_communicator, + true); + vector_zeroed = vector; + all_constraints.set_zero(vector_zeroed); + + value = vector_zeroed.block(block).l2_norm(); + + // Reset if unsensible values + if (value == 0.0) value = 1.0; + value_norm = value; + } + + void + normalise (const L2_norms & norm_0) + { + value_norm/=norm_0.value; + } +}; + +template +void +CoupledProblem::solve_nonlinear_timestep (const int ts) +{ + L2_norms ex_T (T_block, + locally_owned_partitioning, + locally_relevant_partitioning, + mpi_communicator); + L2_norms ex_EM (EM_block, + locally_owned_partitioning, + locally_relevant_partitioning, + mpi_communicator); + + // locally_relevant_solution_t1 = locally_relevant_solution; + + + L2_norms res_T_0(ex_T), update_T_0(ex_T); + L2_norms res_T(ex_T), update_T(ex_T); + L2_norms res_EM_0(ex_EM), update_EM_0(ex_EM); + L2_norms res_EM(ex_EM), update_EM(ex_EM); + + pcout + << std::string(52,' ') + << "|" + << " RES_T " << std::string(2,' ') + << " NUP_T " << std::string(2,' ') + << " RES_EM " << std::string(2,' ') + << " NUP_EM " + << std::endl; + + for (unsigned int n=0; n < Parameters::max_newton_iterations; ++n) + { + pcout << "IT " << n << std::flush; + + LA::MPI::BlockVector locally_relevant_solution_update; + locally_relevant_solution_update.reinit (locally_relevant_partitioning, + mpi_communicator); + + make_constraints(n, ts); + + // === THERMAL PROBLEM === + + system_matrix = 0; + system_rhs = 0; + locally_relevant_solution_update = 0; + + assemble_system_thermo(); + solve_thermo(locally_relevant_solution_update); + locally_relevant_solution.block(T_block) += locally_relevant_solution_update.block(T_block); + // locally_relevant_solution.compress (VectorOperation::add); + + // Compute temperature residual + { + res_T.set(system_rhs, all_constraints); + update_T.set(locally_relevant_solution_update, + all_constraints); + + if (n == 0 || n == 1) + { + res_T_0.set(system_rhs, all_constraints); + update_T_0.set(locally_relevant_solution_update, + all_constraints); + } + + res_T.normalise(res_T_0); + update_T.normalise(update_T_0); + } + + // === ELECTRO-MECHANICAL PROBLEM === + + system_matrix = 0; + system_rhs = 0; + locally_relevant_solution_update = 0; + + assemble_system_mech(); + solve_mech(locally_relevant_solution_update); + locally_relevant_solution.block(EM_block) += locally_relevant_solution_update.block(EM_block); + // locally_relevant_solution.compress (VectorOperation::add); + + // To analyse the residual, we must reassemble both + // systems since they depend on one another + // assemble_system_thermo(); + // assemble_system_mech(); + + // Compute electro-mechanical residual + { + res_EM.set(system_rhs, all_constraints); + update_EM.set(locally_relevant_solution_update, + all_constraints); + + if (n == 0 || n == 1) + { + res_EM_0.set(system_rhs, all_constraints); + update_EM_0.set(locally_relevant_solution_update, + all_constraints); + } + + res_EM.normalise(res_EM_0); + update_EM.normalise(update_EM_0); + } + + pcout + << std::fixed + << std::setprecision(3) + << std::setw(7) + << std::scientific + << "|" + << " " << res_T.value_norm + << " " << update_T.value_norm + << " " << res_EM.value_norm + << " " << update_EM.value_norm + << std::endl; + + bool converged_abs=false; + bool converged_rel=false; + + { + if((res_T.value < Parameters::max_res_abs) && + (res_EM.value < Parameters::max_res_abs)) + { + converged_abs = true; + } + + if((res_T.value_norm < Parameters::max_res_T_norm) && + (res_EM.value_norm < Parameters::max_res_EM_norm)) + { + converged_rel = true; + } + } + + if (converged_abs || converged_rel) + { + pcout + << "Converged." + << std::endl; + break; + } + + if (n == (Parameters::max_newton_iterations-1)) + { + pcout + << "No convergence... :-/" + << std::endl; + } + } + + pcout + << "Absolute values of residuals and Newton update:" + << std::endl + << "res_T: " << res_T.value + << "\t update_T: " << update_T.value + << std::endl + << "res_EM: " << res_EM.value + << "\t update_EM: " << update_EM.value + << std::endl; + +} + + +template +void +CoupledProblem::run () +{ + make_grid(); + setup_system(); + + if (Parameters::use_3_field_formulation) + { + // Currently, this function is not implemented for any Triangulation in the + // parallel namespace in deal.II 9.1 and 9.2. + // So for this reason we do the projection for the dilatation field manually. + const bool use_deal_II_project = false; + if (use_deal_II_project) + { + AffineConstraints constraints; + constraints.close(); + + const ComponentSelectFunction J_mask (J_component, n_components); + + LA::MPI::BlockVector tmp (locally_owned_partitioning, + mpi_communicator); + + VectorTools::project (dof_handler, + constraints, + q_collection, + J_mask, + tmp); + + locally_relevant_solution = tmp; + } + else + { + FE_DGPMonomial fe_dilatation(poly_order - 1); + QGauss qf_cell_dilatation(poly_order+1); + + Vector initial_value_dof(fe_dilatation.dofs_per_cell); + Vector initial_value_qp(qf_cell_dilatation.size()); + for (unsigned int q=0; q qpoint_to_dof_matrix (fe_dilatation.dofs_per_cell, + qf_cell_dilatation.size()); + FETools::compute_projection_from_quadrature_points_matrix + (fe_dilatation, + qf_cell_dilatation, + qf_cell_3Field, + qpoint_to_dof_matrix); + + Vector cell_rhs (fe_cell_3Field.dofs_per_cell); + std::vector local_dof_indices (fe_cell_3Field.dofs_per_cell); + + LA::MPI::BlockVector soln_projected_dilatation (locally_owned_partitioning, + mpi_communicator); + + for (typename hp::DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(); + cell != dof_handler.end(); + ++cell) + { + if (!cell->is_locally_owned()) + continue; + Assert(cell->active_fe_index() == Parameters::fe_index_3_field, ExcInternalError()); + + cell_rhs = 0; + cell->get_dof_indices(local_dof_indices); + + qpoint_to_dof_matrix.vmult (initial_value_dof, + initial_value_qp); + + for (unsigned int i = 0, i_proj = 0; i < cell->get_fe().dofs_per_cell; ++i) + { + const unsigned int i_group = fe_cell_3Field.system_to_base_index(i).first.first; + if (i_group == J_dof) + { + Assert(i_proj < initial_value_dof.size(), ExcInternalError()); + cell_rhs(i) = initial_value_dof(i_proj++); + } + } + + all_constraints.distribute_local_to_global(cell_rhs, + local_dof_indices, + soln_projected_dilatation); + } + + soln_projected_dilatation.compress(VectorOperation::add); + locally_relevant_solution = soln_projected_dilatation; + } + } + + output_results(0); + + double time = Parameters::dt; + for (unsigned int ts = 1; + ts<=Parameters::n_timesteps; + ++ts, time += Parameters::dt) + { + pcout + << std::endl + << std::string(100,'=') + << std::endl + << "Timestep: " << ts + << "\t Time: " << time + << std::endl + << std::string(100,'=') + << std::endl; + solve_nonlinear_timestep(ts); + output_results(ts); + + } + + +} +} + + +int +main (int argc, char *argv[]) +{ + try + { + using namespace dealii; + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + deallog.depth_console (0); + + Coupled_TEE::CoupledProblem<3> coupled_thermo_electro_elastic_problem_3d; + coupled_thermo_electro_elastic_problem_3d.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; + std::cerr << "Exception on processing: " + << std::endl + << exc.what() + << std::endl + << "Aborting!" + << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; + std::cerr << "Unknown exception!" + << std::endl + << "Aborting!" + << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +}