/usr/include/lorene/C++/Include/scalar.h is in liblorene-dev 0.0.0~cvs20161116+dfsg-1ubuntu4.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 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 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 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 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618 1619 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681 1682 1683 1684 1685 1686 1687 1688 1689 1690 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1731 1732 1733 1734 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1750 1751 1752 1753 1754 1755 1756 1757 1758 1759 1760 1761 1762 1763 1764 1765 1766 1767 1768 1769 1770 1771 1772 1773 1774 1775 1776 1777 1778 1779 1780 1781 1782 1783 1784 1785 1786 1787 1788 1789 1790 1791 1792 1793 1794 1795 1796 1797 1798 1799 1800 1801 1802 1803 1804 1805 1806 1807 1808 1809 1810 1811 1812 1813 1814 1815 1816 1817 1818 1819 1820 1821 1822 1823 1824 1825 1826 1827 1828 1829 1830 1831 1832 1833 1834 1835 1836 1837 1838 1839 1840 1841 1842 1843 1844 1845 1846 1847 1848 1849 1850 1851 1852 1853 1854 1855 1856 1857 1858 1859 1860 1861 1862 1863 1864 1865 1866 1867 1868 1869 1870 1871 1872 1873 1874 1875 1876 1877 1878 1879 1880 1881 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 1895 1896 1897 1898 1899 1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 1914 1915 1916 1917 1918 1919 1920 1921 1922 1923 1924 1925 1926 1927 1928 1929 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 | /*
* Definition of Lorene class Scalar
*
*/
/*
* Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
*
* Copyright (c) 1999-2000 Jean-Alain Marck (for previous class Cmp)
* Copyright (c) 1999-2002 Eric Gourgoulhon (for previous class Cmp)
* Copyright (c) 1999-2001 Philippe Grandclement (for previous class Cmp)
* Copyright (c) 2000-2002 Jerome Novak (for previous class Cmp)
* Copyright (c) 2000-2001 Keisuke Taniguchi (for previous class Cmp)
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
#ifndef __SCALAR_H_
#define __SCALAR_H_
/*
* $Id: scalar.h,v 1.94 2015/12/18 15:52:51 j_novak Exp $
* $Log: scalar.h,v $
* Revision 1.94 2015/12/18 15:52:51 j_novak
* New method is_nan() for class Scalar.
*
* Revision 1.93 2015/09/10 13:28:05 j_novak
* New methods for the class Hot_Eos
*
* Revision 1.92 2014/10/13 08:52:36 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.91 2013/06/05 15:06:10 j_novak
* Legendre bases are treated as standard bases, when the multi-grid
* (Mg3d) is built with BASE_LEG.
*
* Revision 1.90 2013/01/11 15:44:54 j_novak
* Addition of Legendre bases (part 2).
*
* Revision 1.89 2012/12/19 13:59:56 j_penner
* added a few lines to the documentation for scalar_match_tau function
*
* Revision 1.88 2012/01/17 15:05:46 j_penner
* *** empty log message ***
*
* Revision 1.87 2012/01/17 10:16:27 j_penner
* functions added: sarra_filter_r, sarra_filter_r_all_domains, Heaviside
*
* Revision 1.86 2011/04/08 13:13:09 e_gourgoulhon
* Changed the comment of function val_point to indicate specifically the
* division by r^dzpuis in the compactified external domain.
*
* Revision 1.85 2008/09/29 13:23:51 j_novak
* Implementation of the angular mapping associated with an affine
* mapping. Things must be improved to take into account the domain index.
*
* Revision 1.84 2008/09/22 19:08:01 j_novak
* New methods to deal with boundary conditions
*
* Revision 1.83 2008/05/24 15:05:22 j_novak
* New method Scalar::match_tau to match the output of an explicit time-marching scheme with the tau method.
*
* Revision 1.82 2007/12/21 16:06:16 j_novak
* Methods to filter Tensor, Vector and Sym_tensor objects.
*
* Revision 1.81 2007/10/31 10:33:11 j_novak
* Added exponential filters to smooth Gibbs-type phenomena.
*
* Revision 1.80 2007/06/21 19:56:36 k_taniguchi
* Introduction of another filtre_r.
*
* Revision 1.79 2007/05/06 10:48:08 p_grandclement
* Modification of a few operators for the vorton project
*
* Revision 1.78 2007/01/16 15:05:59 n_vasset
* New constructor (taking a Scalar in mono-domain angular grid for
* boundary) for function sol_elliptic_boundary
*
* Revision 1.77 2006/05/26 09:00:09 j_novak
* New members for multiplication or division by cos(theta).
*
* Revision 1.76 2005/11/30 13:48:06 e_gourgoulhon
* Replaced M_PI/2 by 1.57... in argument list of sol_elliptic_sin_zec
* (in order not to require the definition of M_PI).
*
* Revision 1.75 2005/11/30 11:09:03 p_grandclement
* Changes for the Bin_ns_bh project
*
* Revision 1.74 2005/11/17 15:29:46 e_gourgoulhon
* Added arithmetics with Mtbl.
*
* Revision 1.73 2005/10/25 08:56:34 p_grandclement
* addition of std_spectral_base in the case of odd functions near the origin
*
* Revision 1.72 2005/09/07 13:10:47 j_novak
* Added a filter setting to zero all mulitpoles in a given range.
*
* Revision 1.71 2005/08/26 14:02:38 p_grandclement
* Modification of the elliptic solver that matches with an oscillatory exterior solution
* small correction in Poisson tau also...
*
* Revision 1.70 2005/08/25 12:14:07 p_grandclement
* Addition of a new method to solve the scalar Poisson equation, based on a multi-domain Tau-method
*
* Revision 1.69 2005/06/09 07:56:25 f_limousin
* Implement a new function sol_elliptic_boundary() and
* Vector::poisson_boundary(...) which solve the vectorial poisson
* equation (method 6) with an inner boundary condition.
*
* Revision 1.68 2005/06/08 12:35:18 j_novak
* New method for solving divergence-like ODEs.
*
* Revision 1.67 2005/05/20 14:42:30 j_novak
* Added the method Scalar::get_spectral_base().
*
* Revision 1.66 2005/04/04 21:28:57 e_gourgoulhon
* Added argument lambda to method poisson_angu
* to deal with the generalized angular Poisson equation:
* Lap_ang u + lambda u = source.
*
* Revision 1.65 2004/12/14 09:09:39 f_limousin
* Modif. comments.
*
* Revision 1.64 2004/11/23 12:41:53 f_limousin
* Intoduce function poisson_dir_neu(...) to solve a scalar poisson
* equation with a mixed boundary condition (Dirichlet + Neumann).
*
* Revision 1.63 2004/10/11 15:09:00 j_novak
* The radial manipulation functions take Scalar as arguments, instead of Cmp.
* Added a conversion operator from Scalar to Cmp.
* The Cmp radial manipulation function make conversion to Scalar, call to the
* Map_radial version with a Scalar argument and back.
*
* Revision 1.62 2004/08/24 09:14:40 p_grandclement
* Addition of some new operators, like Poisson in 2d... It now requieres the
* GSL library to work.
*
* Also, the way a variable change is stored by a Param_elliptic is changed and
* no longer uses Change_var but rather 2 Scalars. The codes using that feature
* will requiere some modification. (It should concern only the ones about monopoles)
*
* Revision 1.61 2004/07/27 08:24:26 j_novak
* Modif. comments
*
* Revision 1.60 2004/07/26 16:02:21 j_novak
* Added a flag to specify whether the primitive should be zero either at r=0
* or at r going to infinity.
*
* Revision 1.59 2004/07/06 13:36:27 j_novak
* Added methods for desaliased product (operator |) only in r direction.
*
* Revision 1.58 2004/06/22 08:49:57 p_grandclement
* Addition of everything needed for using the logarithmic mapping
*
* Revision 1.57 2004/06/14 15:24:23 e_gourgoulhon
* Added method primr (radial primitive).
*
* Revision 1.56 2004/06/11 14:29:56 j_novak
* Scalar::multipole_spectrum() is now a const method.
*
* Revision 1.55 2004/05/24 14:07:31 e_gourgoulhon
* Method set_domain now includes a call to del_deriv() for safety.
*
* Revision 1.54 2004/05/07 11:26:10 f_limousin
* New method filtre_r(int*)
*
* Revision 1.53 2004/03/22 13:12:43 j_novak
* Modification of comments to use doxygen instead of doc++
*
* Revision 1.52 2004/03/17 15:58:47 p_grandclement
* Slight modification of sol_elliptic_no_zec
*
* Revision 1.51 2004/03/11 12:07:30 e_gourgoulhon
* Added method visu_section_anim.
*
* Revision 1.50 2004/03/08 15:45:38 j_novak
* Modif. comment
*
* Revision 1.49 2004/03/05 15:09:40 e_gourgoulhon
* Added method smooth_decay.
*
* Revision 1.48 2004/03/01 09:57:02 j_novak
* the wave equation is solved with Scalars. It now accepts a grid with a
* compactified external domain, which the solver ignores and where it copies
* the values of the field from one time-step to the next.
*
* Revision 1.47 2004/02/27 09:43:58 f_limousin
* New methods filtre_phi(int) and filtre_theta(int).
*
* Revision 1.46 2004/02/26 22:46:26 e_gourgoulhon
* Added methods derive_cov, derive_con and derive_lie.
*
* Revision 1.45 2004/02/21 17:03:49 e_gourgoulhon
* -- Method "point" renamed "val_grid_point".
* -- Method "set_point" renamed "set_grid_point".
*
* Revision 1.44 2004/02/19 22:07:35 e_gourgoulhon
* Added argument "comment" in method spectral_display.
*
* Revision 1.43 2004/02/11 09:47:44 p_grandclement
* Addition of a new elliptic solver, matching with the homogeneous solution
* at the outer shell and not solving in the external domain (more details
* coming soon ; check your local Lorene dealer...)
*
* Revision 1.42 2004/01/28 16:46:22 p_grandclement
* Addition of the sol_elliptic_fixe_der_zero stuff
*
* Revision 1.41 2004/01/28 13:25:40 j_novak
* The ced_mult_r arguments have been suppressed from the Scalar::*dsd* methods.
* In the div/mult _r_dzpuis, there is no more default value.
*
* Revision 1.40 2004/01/28 10:39:17 j_novak
* Comments modified.
*
* Revision 1.39 2004/01/27 15:10:01 j_novak
* New methods Scalar::div_r_dzpuis(int) and Scalar_mult_r_dzpuis(int)
* which replace div_r_inc*. Tried to clean the dzpuis handling.
* WARNING: no testing at this point!!
*
* Revision 1.38 2004/01/23 13:25:44 e_gourgoulhon
* Added methods set_inner_boundary and set_outer_boundary.
* Methods set_val_inf and set_val_hor, which are particular cases of
* the above, have been suppressed.
*
* Revision 1.37 2004/01/22 16:10:09 e_gourgoulhon
* Added (provisory) method div_r_inc().
*
* Revision 1.36 2003/12/16 06:32:20 e_gourgoulhon
* Added method visu_box.
*
* Revision 1.35 2003/12/14 21:46:35 e_gourgoulhon
* Added argument start_dx in visu_section.
*
* Revision 1.34 2003/12/11 16:19:38 e_gourgoulhon
* Added method visu_section for visualization with OpenDX.
*
* Revision 1.33 2003/12/11 14:48:47 p_grandclement
* Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
*
* Revision 1.32 2003/11/13 13:43:53 p_grandclement
* Addition of things needed for Bhole::update_metric (const Etoile_bin&, double, double)
*
* Revision 1.31 2003/11/06 14:43:37 e_gourgoulhon
* Gave a name to const arguments in certain method prototypes (e.g.
* constructors) to correct a bug of DOC++.
*
* Revision 1.30 2003/11/04 22:55:50 e_gourgoulhon
* Added new methods mult_cost(), mult_sint() and div_sint().
*
* Revision 1.29 2003/10/29 13:09:11 e_gourgoulhon
* -- Added integer argument to derivative functions dsdr, etc...
* so that one can choose the dzpuis of the result (default=2).
* -- Change of method name: laplacien --> laplacian.
*
* Revision 1.28 2003/10/29 11:00:42 e_gourgoulhon
* Virtual functions dec_dzpuis and inc_dzpuis have now an integer argument to
* specify by which amount dzpuis is to be increased.
* Accordingly virtual methods dec2_dzpuis and inc2_dzpuis have been suppressed.
*
* Revision 1.27 2003/10/20 14:26:02 j_novak
* New assignement operators.
*
* Revision 1.26 2003/10/19 19:46:33 e_gourgoulhon
* -- Method spectral_display now virtual (from Tensor), list of argument
* changed.
*
* Revision 1.25 2003/10/17 13:46:14 j_novak
* The argument is now between 1 and 3 (instead of 0->2)
*
* Revision 1.24 2003/10/16 15:23:41 e_gourgoulhon
* Name of method div_r_ced() changed to div_r_inc2().
* Name of method div_rsint_ced() changed to div_rsint_inc2().
*
* Revision 1.23 2003/10/15 21:10:11 e_gourgoulhon
* Added method poisson_angu().
*
* Revision 1.22 2003/10/15 16:03:35 j_novak
* Added the angular Laplace operator for Scalar.
*
* Revision 1.21 2003/10/15 10:29:05 e_gourgoulhon
* Added new members p_dsdt and p_stdsdp.
* Added new methods dsdt(), stdsdp() and div_tant().
*
* Revision 1.20 2003/10/13 13:52:39 j_novak
* Better managment of derived quantities.
*
* Revision 1.19 2003/10/10 15:57:27 j_novak
* Added the state one (ETATUN) to the class Scalar
*
* Revision 1.18 2003/10/08 14:24:08 j_novak
* replaced mult_r_zec with mult_r_ced
*
* Revision 1.17 2003/10/06 16:16:02 j_novak
* New constructor from a Tensor.
*
* Revision 1.16 2003/10/06 13:58:45 j_novak
* The memory management has been improved.
* Implementation of the covariant derivative with respect to the exact Tensor
* type.
*
* Revision 1.15 2003/10/05 21:06:31 e_gourgoulhon
* - Added new methods div_r_ced() and div_rsint_ced().
* - Added new virtual method std_spectral_base()
* - Removed method std_spectral_base_scal()
*
* Revision 1.14 2003/10/01 13:02:58 e_gourgoulhon
* Suppressed the constructor from Map* .
*
* Revision 1.13 2003/09/29 12:52:56 j_novak
* Methods for changing the triad are implemented.
*
* Revision 1.12 2003/09/25 09:33:36 j_novak
* Added methods for integral calculation and various manipulations
*
* Revision 1.11 2003/09/25 09:11:21 e_gourgoulhon
* Added functions for radial operations (divr, etc...)
*
* Revision 1.10 2003/09/25 08:55:23 e_gourgoulhon
* Added members raccord*.
*
* Revision 1.9 2003/09/25 08:50:11 j_novak
* Added the members import
*
* Revision 1.8 2003/09/25 08:13:51 j_novak
* Added method for calculating derivatives
*
* Revision 1.7 2003/09/25 07:59:26 e_gourgoulhon
* Added prototypes for PDE resolutions.
*
* Revision 1.6 2003/09/25 07:17:58 j_novak
* Method asymptot implemented.
*
* Revision 1.5 2003/09/24 20:53:38 e_gourgoulhon
* Added -- constructor by conversion from a Cmp
* -- assignment from Cmp
*
* Revision 1.4 2003/09/24 15:10:54 j_novak
* Suppression of the etat flag in class Tensor (still present in Scalar)
*
* Revision 1.3 2003/09/24 12:01:44 j_novak
* Added friend functions for math.
*
* Revision 1.2 2003/09/24 10:22:01 e_gourgoulhon
* still in progress...
*
* Revision 1.1 2003/09/22 12:50:47 e_gourgoulhon
* First version: not ready yet!
*
*
* $Header: /cvsroot/Lorene/C++/Include/scalar.h,v 1.94 2015/12/18 15:52:51 j_novak Exp $
*
*/
// Headers Lorene
#include "valeur.h"
#include "tensor.h"
namespace Lorene {
class Param ;
class Cmp ;
class Param_elliptic ;
/**
* Tensor field of valence 0 (or component of a tensorial field).
* \ingroup (tensor)
*
*
*/
class Scalar : public Tensor {
// Data :
// -----
protected:
/** The logical state \c ETATNONDEF (undefined), \c ETATZERO (null),
* \c ETATUN (one), or \c ETATQCQ (ordinary).
*/
int etat ;
/**
* Power of \e r by which the quantity represented by \c this
* must be divided in the compactified external domain (CED) in order
* to get the correct physical values
*/
int dzpuis ;
Valeur va ; ///< The numerical value of the \c Scalar
// Derived data :
// ------------
protected:
/// Pointer on \f$\partial/\partial r\f$ of \c *this (0x0 if not up to date)
mutable Scalar* p_dsdr ;
/** Pointer on \f$1/r \partial/\partial \theta\f$ of \c *this
* (0x0 if not up to date)
*/
mutable Scalar* p_srdsdt ;
/** Pointer on \f$1/(r\sin\theta) \partial/\partial \phi\f$ of \c *this
* (0x0 if not up to date)
*/
mutable Scalar* p_srstdsdp ;
/// Pointer on \f$\partial/\partial \theta\f$ of \c *this (0x0 if not up to date)
mutable Scalar* p_dsdt ;
/** Pointer on \f$1/\sin\theta \partial/\partial \phi\f$ of \c *this
* (0x0 if not up to date)
*/
mutable Scalar* p_stdsdp ;
/** Pointer on \f$\partial/\partial x\f$ of \c *this ,
* where \f$x=r\sin\theta \cos\phi\f$ (0x0 if not up to date)
*/
mutable Scalar* p_dsdx ;
/** Pointer on \f$\partial/\partial y\f$ of \c *this ,
* where \f$y=r\sin\theta \sin\phi\f$(0x0 if not up to date)
*/
mutable Scalar* p_dsdy ;
/** Pointer on \f$\partial/\partial z\f$ of \c *this ,
* where \f$z=r\cos\theta\f$ (0x0 if not up to date)
*/
mutable Scalar* p_dsdz ;
/** Pointer on the Laplacian of \c *this (0x0 if not up to date)
*/
mutable Scalar* p_lap ;
/** Pointer on the Laplacian of \c *this (0x0 if not up to date)
*/
mutable Scalar* p_lapang ;
/// Pointer on \f$\partial/\partial radial \f$ of \c *this
mutable Scalar* p_dsdradial ;
/// Pointer on \f$\partial/\partial \rho \f$ of \c *this
mutable Scalar* p_dsdrho ;
/** Power of \e r by which the last computed Laplacian has been
* multiplied in the compactified external domain.
*/
mutable int ind_lap ;
/** Pointer on the space integral of \c *this (values in each
* domain) (0x0 if not up to date)
*/
mutable Tbl* p_integ ;
// Constructors - Destructor
// -------------------------
public:
explicit Scalar(const Map& mpi) ; ///< Constructor from mapping
/// Constructor from a Tensor (must be of valence 0)
Scalar(const Tensor& a) ;
Scalar(const Scalar& a) ; ///< Copy constructor
Scalar(const Cmp& a) ; ///< Constructor by conversion of a Cmp
/// Constructor from a file (see \c sauve(FILE*) )
Scalar(const Map&, const Mg3d&, FILE* ) ;
virtual ~Scalar() ; ///< Destructor
// Memory management
// -----------------
protected:
void del_t() ; ///< Logical destructor
virtual void del_deriv() const; ///< Logical destructor of the derivatives
void set_der_0x0() const; ///< Sets the pointers for derivatives to 0x0
public:
/**
* Sets the logical state to \c ETATNONDEF (undefined).
* Calls the logical destructor of the \c Valeur \c va
* deallocates the memory occupied by all the derivatives.
*/
virtual void set_etat_nondef() ;
/**
* Sets the logical state to \c ETATZERO (zero).
* Calls the logical destructor of the \c Valeur \c va and
* deallocates the memory occupied by all the derivatives.
*/
virtual void set_etat_zero() ;
/**
* Sets the logical state to \c ETATQCQ (ordinary state).
* If the state is already \c ETATQCQ , this function does nothing.
* Otherwise, it calls the logical destructor of the \c Valeur \c va and
* deallocates the memory occupied by all the derivatives.
*/
virtual void set_etat_qcq() ;
/**
* Sets the logical state to \c ETATUN (one).
* Fills the \c Valeur \c va with ones and
* deallocates the memory occupied by all the derivatives.
*/
void set_etat_one() ;
/**
* Sets the logical state to \c ETATQCQ (ordinary state)
* and performs the memory allocation of all the
* elements, down to the \c double arrays of the \c Tbl s.
* This function performs in fact recursive calls to \c set_etat_qcq()
* on each element of the chain \c Scalar ->
* \c Valeur -> \c Mtbl -> \c Tbl .
*/
virtual void allocate_all() ;
/**
* Sets the \c Scalar to zero in a hard way.
* 1/ Sets the logical state to \c ETATQCQ , i.e. to an ordinary state.
* 2/ Fills the \c Valeur \c va with zeros.
* NB: this function must be used for debugging purposes only.
* For other operations, the functions \c set_etat_zero()
* or \c annule(int,int) must be perferred.
*/
void annule_hard() ;
// Extraction of information
// -------------------------
public:
/** Returns the logical state \c ETATNONDEF (undefined),
* \c ETATZERO (null) or \c ETATQCQ (ordinary).
*/
int get_etat() const {return etat;} ;
/// Returns \c dzpuis
int get_dzpuis() const {return dzpuis;} ;
/** Returns \c true if the last domain is compactified and
* \c *this is not zero in this domain
*/
bool dz_nonzero() const ;
/** Returns \c false if the last domain is compactified
* and \c *this is not zero in this domain and \c dzpuis
* is not equal to \c dzi , otherwise return true.
*/
bool check_dzpuis(int dzi) const ;
/** Looks for NaNs (not a number) in the scalar field.
* If at least one NaN is found, it returns \c true. If the flag \c verbose
* is set to \c true, it outputs to the standard output the indices where NaNs
* have been found.
*/
bool is_nan(bool verbose=false) const ;
// Assignment
// -----------
public:
/// Assignment to another \c Scalar defined on the same mapping
void operator=(const Scalar& a) ;
/// Assignment to a \c Tensor (of valence 0)
virtual void operator=(const Tensor& a) ;
void operator=(const Cmp& a) ; ///< Assignment to a \c Cmp
void operator=(const Valeur& a) ; ///< Assignment to a \c Valeur
void operator=(const Mtbl& a) ; ///< Assignment to a \c Mtbl
void operator=(double ) ; ///< Assignment to a \c double
void operator=(int ) ; ///< Assignment to an \c int
// Conversion oprators
//---------------------
operator Cmp() const ; ///< Conversion to a \c Cmp
// Access to individual elements
// -----------------------------
public:
/// Returns \c va (read only version)
const Valeur& get_spectral_va() const {return va;} ;
/// Returns \c va (read/write version)
Valeur& set_spectral_va() {return va;} ;
/** Read/write of the value in a given domain.
* This method should be used only to set the value in a given
* domain (it performs a call to \c del_deriv); for reading the
* value in a domain without changing it, the method \c domain(int )
* is preferable.
*
* @param l [input] domain index
* @return writable \c Tbl containing the value of the field in domain \c l .
*/
Tbl& set_domain(int l) {
assert(etat == ETATQCQ) ;
del_deriv() ;
return va.set(l) ;
};
/** Read-only of the value in a given domain.
* @param l [input] domain index
* @return \c Tbl containing the value of the field in domain \c l .
*/
const Tbl& domain(int l) const {
assert( (etat == ETATQCQ) || (etat == ETATUN) ) ;
return va(l) ;
};
/** Returns the value of the field at a specified grid point.
* @param l [input] domain index
* @param k [input] \f$\phi\f$ index
* @param j [input] \f$\theta\f$ index
* @param i [input] \e r (\f$\xi\f$) index
*/
double val_grid_point(int l, int k, int j, int i) const {
assert(etat != ETATNONDEF) ;
if (etat == ETATZERO) {
double zero = 0. ;
return zero ;
}
else {
if (etat == ETATUN) {
double one = 1. ;
return one ;
}
else{
return va(l, k, j, i) ;
}
}
};
/** Computes the value of the field at an
* arbitrary point \f$(r, \theta, \phi)\f$, by means of the spectral
* expansion.
* NB: if \f$(r, \theta, \phi)\f$ is a point of the spectral grid,
* the method \c val_grid_point is to be preferred,
* being much more efficient.
* @param r [input] value of the coordinate \e r
* @param theta [input] value of the coordinate \f$\theta\f$
* @param phi [input] value of the coordinate \f$\phi\f$
* @return value at the point \f$(r, \theta, \phi)\f$
* of the field represented by \c *this . NB: in the compactified
* external domain, the returned value is the actual value of the field,
* i.e. the stored value divided by \f$ r^{\rm dzpuis} \f$.
*/
double val_point(double r, double theta, double phi) const ;
/** Setting the value of the field at a given grid point.
* CAUTION: to gain in efficiency (especially when this method is
* invoqued inside a loop), the method \c del_deriv() (to delete
* the derived members) is not called by \c set_grid_point .
* It must thus be invoqued by the user, after all the calls
* to \c set_grid_point have been performed.
*
* @param l [input] domain index
* @param k [input] \f$\phi\f$ index
* @param j [input] \f$\theta\f$ index
* @param i [input] \e r (\f$\xi\f$) index
* @return writable value of the field at the specified grid point
*/
double& set_grid_point(int l, int k, int j, int i) {
assert(etat == ETATQCQ) ;
return va.set(l, k, j, i) ;
};
/**
* Sets the \c Scalar to zero in several domains.
* @param l_min [input] The \c Scalar will be set (logically) to zero
* in the domains whose indices are in the range
* \c [l_min,l_max] .
* @param l_max [input] see the comments for \c l_min .
*
* Note that \c annule(0,va.mg->get_nzone()-1) is equivalent to
* \c set_etat_zero() .
*/
virtual void annule(int l_min, int l_max) ;
/** Sets the value of the \c Scalar at the inner boundary of a given
* domain.
* @param l [input] domain index
* @param x [input] (constant) value at the inner boundary of domain no. \c l
*/
void set_inner_boundary(int l, double x) ;
/** Sets the value of the \c Scalar at the outer boundary of a given
* domain.
* @param l [input] domain index
* @param x [input] (constant) value at the outer boundary of domain no. \c l
*/
void set_outer_boundary(int l, double x) ;
/**
* Gives the spectrum in terms of multipolar modes \e l .
* @return a \c Tbl of size (nzone, lmax), where lmax is the
* maximal multipolar momentum over all domains. The \e l -th
* element contains the L1 norm of the \e l -th multipole
* (\e i.e. a sum over all \e m of the norms (coefficient space)
* of the component of a given \f$Y_l^m\f$.
*/
Tbl multipole_spectrum () const ;
/** Returns the \c Tbl containing the values of angular coefficients
* at the outer boundary.
* @param l_dom [input] domain index
* @param leave_ylm [input] flag to decide whether the coefficients
* are expressed in spherical harmonics or Fourier base
*/
Tbl tbl_out_bound(int l_dom, bool leave_ylm = false) ;
/** Returns the \c Tbl containing the values of angular coefficients
* at the inner boundary.
* @param l_dom [input] domain index
* @param leave_ylm [input] flag to decide whether the coefficients
* are expressed in spherical harmonics or Fourier base
*/
Tbl tbl_in_bound(int n, bool leave_ylm = false) ;
/** Returns the \c Scalar containing the values of angular coefficients
* at the outer boundary.
* @param l_dom [input] domain index
* @param leave_ylm [input] flag to decide whether the coefficients
* are expressed in spherical harmonics or Fourier base
*/
Scalar scalar_out_bound(int n, bool leave_ylm = false) ;
// Differential operators and others
// ---------------------------------
public:
/** Returns \f$\partial / \partial r\f$ of \c *this .
* If \c dzpuis is zero, then the returned \c Scalar has
* \c dzpuis = 2. It is increased by 1 otherwise.
*/
const Scalar& dsdr() const ;
/** Returns \f$1/r \partial / \partial \theta\f$ of \c *this .
* If \c dzpuis is zero, then the returned \c Scalar has
* \c dzpuis = 2. It is increased by 1 otherwise.
*/
const Scalar& srdsdt() const ;
/** Returns \f$1/(r\sin\theta) \partial / \partial \phi\f$ of \c *this .
* If \c dzpuis is zero, then the returned \c Scalar has
* \c dzpuis = 2. It is increased by 1 otherwise.
*/
const Scalar& srstdsdp() const ;
/** Returns \f$\partial / \partial \theta\f$ of \c *this .
*/
const Scalar& dsdt() const ;
/** Returns \f$\partial / \partial r\f$ of \c *this if the mapping is
* affine (class \c Map_af) and \f$\partial / \partial \ln r\f$
* of \c *this if the mapping
* is logarithmic (class \c Map_log).
* If \c dzpuis is zero, then the returned \c Scalar has
* \c dzpuis = 2. It is increased by 1 otherwise.
*/
const Scalar& dsdradial() const ;
/** Returns \f$\partial / \partial \rho \f$ of \c *this .
* If \c dzpuis is zero, then the returned \c Scalar has
* \c dzpuis = 2. It is increased by 1 otherwise.
*/
const Scalar& dsdrho() const ;
/** Returns \f$1/\sin\theta \partial / \partial \phi\f$ of \c *this .
*/
const Scalar& stdsdp() const ;
/** Returns \f$\partial/\partial x\f$ of \c *this ,
* where \f$x=r\sin\theta \cos\phi\f$.
* If \c dzpuis is zero, then the returned \c Scalar has
* \c dzpuis = 2. It is increased by 1 otherwise.
*/
const Scalar& dsdx() const ;
/** Returns \f$\partial/\partial y\f$ of \c *this ,
* where \f$y=r\sin\theta \sin\phi\f$.
* If \c dzpuis is zero, then the returned \c Scalar has
* \c dzpuis = 2. It is increased by 1 otherwise.
*/
const Scalar& dsdy() const ;
/** Returns \f$\partial/\partial z\f$ of \c *this ,
* where \f$z=r\cos\theta\f$.
* If \c dzpuis is zero, then the returned \c Scalar has
* \c dzpuis = 2. It is increased by 1 otherwise.
*/
const Scalar& dsdz() const ;
/** Returns \f$\partial/\partial x_i\f$ of \c *this ,
* where \f$x_i = (x, y, z)\f$.
* If \c dzpuis is zero, then the returned \c Scalar has
* \c dzpuis = 2. It is increased by 1 otherwise.
* @param i [input] i=1 for \e x , i=2 for \e y , i=3 for \e z .
*/
const Scalar& deriv(int i) const ;
/** Returns the gradient (1-form = covariant vector) of \c *this
* @param gam metric components only used to get the triad with
* respect to which the components of the result are defined
*/
const Vector& derive_cov(const Metric& gam) const ;
/** Returns the "contravariant" derivative of \c *this with respect
* to some metric \f$\gamma\f$, by raising the index of the
* gradient (cf. method \c derive_cov() ) with
* \f$\gamma\f$.
*/
const Vector& derive_con(const Metric& gam) const ;
/// Computes the derivative of \c this along a vector field \c v
Scalar derive_lie(const Vector& v) const ;
/** Returns the Laplacian of \c *this
* @param ced_mult_r [input] Determines the quantity computed in
* the compactified external domain (CED)
* (\e u in the field represented by \c *this ) :
* \li ced_mult_r = 0 : \f$\Delta u\f$
* \li ced_mult_r = 2 : \f$r^2 \, \Delta u\f$
* \li ced_mult_r = 4 (default) : \f$r^4 \, \Delta u\f$
*/
const Scalar& laplacian(int ced_mult_r = 4) const ;
/** Returns the angular Laplacian \f$\Delta_{\theta\varphi}\f$ of \c *this ,
* where \f$\Delta_{\theta\varphi} f = \frac{\partial^2 f}
* {\partial \theta^2} + \frac{1}{\tan \theta} \frac{\partial f}
* {\partial \theta} +\frac{1}{\sin^2 \theta}\frac{\partial^2 f}
* {\partial \varphi^2}\f$
*
*/
const Scalar& lapang() const ;
/// Division by \e r everywhere; \c dzpuis is not changed.
void div_r() ;
/** Division by \e r everywhere but with the output flag \c dzpuis
* set to \c ced_mult_r .
* @param ced_mult_r [input] value of \c dzpuis of the result.
*/
void div_r_dzpuis(int ced_mult_r) ;
/** Division by \e r in the compactified external domain (CED), the
* \c dzpuis flag is not changed.
*/
void div_r_ced() ;
/// Multiplication by \e r everywhere; \c dzpuis is not changed.
void mult_r() ;
/** Multiplication by \e r everywhere but with the output flag \c dzpuis
* set to \c ced_mult_r .
* @param ced_mult_r [input] value of \c dzpuis of the result.
*/
void mult_r_dzpuis(int ced_mult_r) ;
/** Multiplication by \e r in the compactified external domain (CED), the
* \c dzpuis flag is not changed.
*/
void mult_r_ced() ;
/// Multiplication by \f$r\sin\theta\f$ everywhere; \c dzpuis is not changed.
void mult_rsint() ;
/** Multiplication by \f$r\sin\theta\f$ but with the output flag \c dzpuis
* set to \c ced_mult_r .
* @param ced_mult_r [input] value of \c dzpuis of the result.
*/
void mult_rsint_dzpuis(int ced_mult_r) ;
/// Division by \f$r\sin\theta\f$ everywhere; \c dzpuis is not changed.
void div_rsint() ;
/** Division by \f$r\sin\theta\f$ but with the output flag \c dzpuis
* set to \c ced_mult_r .
* @param ced_mult_r [input] value of \c dzpuis of the result.
*/
void div_rsint_dzpuis(int ced_mult_r) ;
void mult_cost() ; ///< Multiplication by \f$\cos\theta\f$
void div_cost() ; ///< Division by \f$\cos\theta\f$
void mult_sint() ; ///< Multiplication by \f$\sin\theta\f$
void div_sint() ; ///< Division by \f$\sin\theta\f$
void div_tant() ; ///< Division by \f$\tan\theta\f$
/** Computes the radial primitive which vanishes for \f$r\to \infty\f$.
* i.e. the function
* \f$ F(r,\theta,\varphi) = \int_r^\infty f(r',\theta,\varphi) \, dr' \f$
* where \e f is the function represented by \c *this
* (and must have a \c dzpuis = 2).
* @param null_infty if true (default), the primitive is null
* at infinity (or on the grid boundary). \e F is null at the
* center otherwise
* @return function \e F
*/
Scalar primr(bool null_infty = true) const ;
/** Computes the integral over all space of \c *this .
* The computed quantity is (\e u being the field represented by
* \c *this )
* \f$\int u \, r^2 \sin\theta \, dr\, d\theta \, d\phi\f$.
* Note that in the compactified external domain (CED), \c dzpuis
* must be 4 for the computation to take place.
*/
double integrale() const ;
/** Computes the integral in each domain of \c *this .
* The computed quantity is (\e u being the field represented by
* \c *this )
* \f$\int u \, r^2 \sin\theta \, dr\, d\theta \, d\phi\f$
* in each domain. The result is returned a \c Tbl on the
* various domains.
* Note that in the compactified external domain (CED), \c dzpuis
* must be 4 for the computation to take place.
*/
const Tbl& integrale_domains() const ;
/** Decreases by \c dec units the value of \c dzpuis and
* changes accordingly the values of the \c Scalar in the
* compactified external domain (CED).
*/
virtual void dec_dzpuis(int dec = 1) ;
/** Increases by \c inc units the value of \c dzpuis and
* changes accordingly the values of the \c Scalar in the
* compactified external domain (CED).
*/
virtual void inc_dzpuis(int inc = 1) ;
/** Sets a new vectorial basis (triad) of decomposition and modifies
* the components accordingly.
*/
virtual void change_triad(const Base_vect& new_triad) ;
/**
* Sets the \c n lasts coefficients in \e r to 0 in the
* external domain.
*/
void filtre (int n) ;
/**
* Sets the \c n lasts coefficients in \e r to 0 in all
* domains.
*/
void filtre_r (int* nn) ;
/**
* Sets the \c n last coefficients in \e r to 0 in the domain \c nzone .
*/
void filtre_r (int n, int nzone) ;
/**
* Applies an exponential filter to the spectral coefficients in the radial direction.
* The filter is of the type: \f$ \forall n\leq N,\, b_n = \sigma(n/N ) a_n\f$, with
* \f$ \sigma(x) = \exp\left( -\ln (10^\alpha ) x^{2p} \right) \f$ and \e N the number
* of radial coefficients.
* @param lzmin, lzmax [input] the indices of the domain where the filter is applied
* (in [\c lzmin , \c lzmax ])
* @param p [input] the order of the filter
* @param alpha [input] \f$\alpha\f$ appearing in the above formula.
*/
// virtual void exponential_filter_r(int lzmin, int lzmax, int p,
// double alpha= -16.) ;
virtual void exponential_filter_r(int lzmin, int lzmax, int p,
double alpha= -16.) ;
/**
* Applies an exponential filter to the spectral coefficients in the radial direction.
* The filter is of the type: \f$ \forall n\leq N,\, b_n = \sigma(n/N ) a_n\f$, with
* \f$ \sigma(x) = \exp\left( \alpha x^{p} \right) \f$ and \e N the number
* of radial coefficients.
* @param lzmin, lzmax [input] the indices of the domain where the filter is applied
* (in [\c lzmin , \c lzmax ])
* @param p [input] the order of the filter
* @param alpha [input] \f$\alpha\f$ appearing in the above formula.
*/
void sarra_filter_r(int lzmin, int lzmax, double p,
double alpha= -1E-16) ;
/**
* Applies an exponential filter in radial direction in all domains.
* (see \c Scalar:exponential_filter_r ). Note that this may cause
* regularity problems at the origin if applied in a nucleus.
*/
void exp_filter_r_all_domains(Scalar &ss, int p, double alpha=-16.) ;
/**
* Applies an exponential filter in radial direction in all domains
* for the case where p is a double
* (see \c Scalar:sarra_filter_r ). Note that this may cause
* regularity problems at the origin if applied in a nucleus.
*/
void sarra_filter_r_all_domains(double p, double alpha=1E-16) ;
/**
* Applies an exponential filter to the spectral coefficients in the angular directions.
* The filter is of the type:
* \f$ \forall \ell \leq \ell_{\rm max},\, \forall m,\, b_{\ell m} = \sigma(\ell/\ell_{\rm max} ) a_{\ell m}\f$, with
* \f$ \sigma(x) \f$ defined for \c Scalar::exponential_filter_r and
* \f$\ell_{\rm max}\f$ the number of spherical harmonics used.
*/
virtual void exponential_filter_ylm(int lzmin, int lzmax, int p,
double alpha= -16.) ;
/**
* Sets all the multipolar components between \c l_min and \c l_max
* to zero. This is done for [ \c l_min , \c l_max ] and all relevant
* \c m in the spherical harmonics expansion basis. If \c ylm_output
* is set to \c true , then the spectral expansion basis of \c this is
* left to be that of spherical harmonics.
*/
void annule_l (int l_min, int l_max, bool ylm_output= false ) ;
/**
* Sets the \c n lasts coefficients in \f$\Phi\f$ to 0 in the
* domain \c zone .
*/
void filtre_phi (int n, int zone) ;
/**
* Sets the \c n lasts coefficients in \f$\theta\f$ to 0 in the
* domain \c nz1 to \c nz2 when expressed in spherical harmonics.
*/
void filtre_tp(int nn, int nz1, int nz2) ;
/**
* Substracts all the components behaving like \f$r^{-n}\f$ in the external
* domain, with \e n strictly lower than \c puis , so that \c *this
* decreases at least like \f$r^{\tt puis} \f$ at infinity.
*/
void fixe_decroissance (int puis) ;
/** Performs a \f$C^k\f$ matching of the last non-compactified shell with
* a decaying function \f$\sum_{j=0}^k {\alpha_j \over r^{\ell+n+j}}\f$ where
* \f$\ell\f$ is the spherical harmonic index and \e n is some
* specifiable parameter.
*/
void smooth_decay(int k, int n) ;
/**
* Performs the \f$C^n\f$ matching of the nucleus with respect to the
* first shell.
*/
void raccord(int n) ;
/**
* Performs the \f$C^1\f$ matching of the external domain with respect to
* the last shell using function like \f$\frac{1}{r^i}\f$ with
* \f${\tt puis} \leq i \leq {\tt puis+nbre}\f$ for each spherical harmonics
* with \f$l \leq {\tt lmax}\f$.
*/
void raccord_c1_zec(int puis, int nbre, int lmax) ;
/**
* Matching of the external domain with the outermost shell
*/
void raccord_externe(int puis, int nbre, int lmax) ;
/**
* Method for matching accross domains and imposing boundary condition.
* Matching of the field represented by \c this accross domains and imposition of the
* boundary condition using the tau method.
* @param par_bc [input] \c Param to control the boundary conditions
* par_bc must contain (at a minimum)
* a modifiable Tbl which specifies a physical boundary
* two integers, one specifying the domain that has the boundary
* the other specifying the number of conditions
* 1 -> Dirichlet
* 2 -> Robin (which may reduce to von Neumann, see below)
* two doubles, specifying the Robin BC parameters. If the first is zero, we see that
* Robin will reduce to von Neumann
* @param par_mat [input/output] optional \c Param in which the matching matrices are
* stored (together with their LU decomposition).
*/
void match_tau(Param& par_bc, Param* par_mat=0x0) ;
/**
* Method for matching accross domains and imposing boundary condition.
* Matching of the field represented by \c this accross domains and imposition of the
* boundary condition using the tau method.
* @param par_bc [input] \c Param to control the boundary conditions
* @param par_mat [input/output] optional \c Param in which the matching matrices are
* stored (together with their LU decomposition).
*/
void match_tau_shell(Param& par_bc, Param* par_mat=0x0) ;
/**
* Method for matching accross domains and imposing boundary condition.
* Matching of the field represented by \c this accross domains and imposition of the
* boundary condition using the collocation method.
* @param par_bc [input] \c Param to control the boundary conditions
* @param par_mat [input/output] optional \c Param in which the matching matrices are
* stored (together with their LU decomposition).
*/
void match_collocation(Param& par_bc, Param* par_mat=0x0) ;
// Outputs
// -------
public:
virtual void sauve(FILE *) const ; ///< Save in a file
/** Displays the spectral coefficients and the associated
* basis functions. This function shows only the values greater than a
* given threshold.
* @param comment comment to be printed at top of the display
* (default: 0x0 = nothing printed)
* @param threshold [input] Value above which a coefficient is printed
* (default: 1.e-7)
* @param precision [input] Number of printed digits (default: 4)
* @param ostr [input] Output stream used for the printing (default: cout)
*/
virtual void spectral_display(const char* comment = 0x0,
double threshold = 1.e-7, int precision = 4,
ostream& ostr = cout) const ;
/// Display
friend ostream& operator<<(ostream& , const Scalar & ) ;
/** 3D visualization via a plane section.
* Prepares files for visualization by OpenDX of the values of the field in
* a plane x=const, y=const or z=const
*
* @param section_type [input] defines the type of section :
* \li 'x' for a plane x = a with a = const (parameter \c aa )
* \li 'y' for a plane y = a with a = const (parameter \c aa )
* \li 'z' for a plane z = a with a = const (parameter \c aa )
* @param aa [input] constant a defining the section plane
* @param umin [input] defines with \c umax the range of the plane coordinate u
* @param umax [input] defines with \c umin the range of the plane coordinate u
* @param vmin [input] defines with \c vmax the range of the plane coordinate v
* @param vmax [input] defines with \c vmin the range of the plane coordinate v
* @param title [input] title for the graph (for OpenDX legend)
* @param filename [input] name for the file which will be the input for
* OpenDX; the default 0x0 is transformed into "scalar_section"
* @param start_dx [input] determines whether OpenDX must be launched (as a
* subprocess) to view the field; if set to \c false , only input files
* for future usage of OpenDX are created
* @param nu [input] number of points in the u direction (uniform sampling)
* @param nv [input] number of points in the v direction (uniform sampling)
*
*/
void visu_section(const char section_type, double aa, double umin, double umax, double vmin,
double vmax, const char* title = 0x0, const char* filename = 0x0,
bool start_dx = true, int nu = 200, int nv = 200) const ;
/** 3D visualization via a plane section.
* Prepares files for visualization by OpenDX of the values of the field in
* any given plane.
*
* @param plane [input] : 2D \c Tbl defining the section plane: \c plane
* must of dimension 3x3 with the following content:
* \li \c plane(0,i) : absolute Cartesian coordinates (xa0,ya0,za0) of some
* point in the plane considered as the origin for the plane coordinates
* (u,v): \c plane(0,0) = xa0 , \c plane(0,1) = ya0 ,
* \li \c plane(0,2) = za0
* \li \c plane(1,i) : components w.r.t. absolute Cartesian coordinates
* of the u-coordinate unit vector in the section plane
* \li \c plane(2,i) : components w.r.t. absolute Cartesian coordinates
* of the v-coordinate unit vector in the section plane
* @param umin [input] defines with \c umax the range of the plane coordinate u
* @param umax [input] defines with \c umin the range of the plane coordinate u
* @param vmin [input] defines with \c vmax the range of the plane coordinate v
* @param vmax [input] defines with \c vmin the range of the plane coordinate v
* @param title [input] title for the graph (for OpenDX legend)
* @param filename [input] name for the file which will be the input for
* OpenDX; the default 0x0 is transformed into "scalar_section"
* @param start_dx [input] determines whether OpenDX must be launched (as a
* subprocess) to view the field; if set to \c false , only input files
* for future usage of OpenDX are created
* @param nu [input] number of points in the u direction (uniform sampling)
* @param nv [input] number of points in the v direction (uniform sampling)
*
*/
void visu_section(const Tbl& plane, double umin, double umax, double vmin,
double vmax, const char* title = 0x0, const char* filename = 0x0,
bool start_dx = true, int nu = 200, int nv = 200) const ;
/** 3D visualization via time evolving plane section (animation).
* Prepares files for visualization by OpenDX of the values of the field in
* a plane x=const, y=const or z=const at successive time steps
*
* @param section_type [input] defines the type of section :
* \li 'x' for a plane x = a with a = const (parameter \c aa )
* \li 'y' for a plane y = a with a = const (parameter \c aa )
* \li 'z' for a plane z = a with a = const (parameter \c aa )
* @param aa [input] constant a defining the section plane
* @param umin [input] defines with \c umax the range of the plane coordinate u
* @param umax [input] defines with \c umin the range of the plane coordinate u
* @param vmin [input] defines with \c vmax the range of the plane coordinate v
* @param vmax [input] defines with \c vmin the range of the plane coordinate v
* @param jtime [input] time step label
* @param ttime [input] time t corresponding to \c jtime
* @param jgraph [input] number of time steps between two graphs: the graph
* will be generated only if \c jtime is a multiple of \c jgraph
* @param title [input] title for the graph (for OpenDX legend)
* @param filename_root [input] beginning of the names for the files which will
* be the input for OpenDX (the end of names will be automatically generated
* from the time steps); the default 0x0 is transformed into "anim"
* @param start_dx [input] determines whether OpenDX must be launched (as a
* subprocess) to view the field; if set to \c false , only input files
* for future usage of OpenDX are created
* @param nu [input] number of points in the u direction (uniform sampling)
* @param nv [input] number of points in the v direction (uniform sampling)
*
*/
void visu_section_anim(const char section_type, double aa, double umin,
double umax, double vmin, double vmax, int jtime, double ttime,
int jgraph = 1, const char* title = 0x0, const char* filename_root = 0x0,
bool start_dx = false, int nu = 200, int nv = 200) const ;
/** 3D visualization (volume rendering) via OpenDX.
* Prepares files for visualization by OpenDX of the values of the field in
* some rectangular box.
*
* @param xmin [input] defines with \c xmax the x range of the visualization box
* @param xmax [input] defines with \c xmin the x range of the visualization box
* @param ymin [input] defines with \c ymax the y range of the visualization box
* @param ymax [input] defines with \c ymin the y range of the visualization box
* @param zmin [input] defines with \c zmax the z range of the visualization box
* @param zmax [input] defines with \c zmin the z range of the visualization box
* @param title [input] title for the graph (for OpenDX legend)
* @param filename [input] name for the file which will be the input for
* OpenDX; the default 0x0 is transformed into "scalar_box"
* @param start_dx [input] determines whether OpenDX must be launched (as a
* subprocess) to view the field; if set to \c false , only input files
* for future usage of OpenDX are created
* @param nx [input] number of points in the x direction (uniform sampling)
* @param ny [input] number of points in the y direction (uniform sampling)
* @param nz [input] number of points in the z direction (uniform sampling)
*
*/
void visu_box(double xmin, double xmax, double ymin, double ymax,
double zmin, double zmax, const char* title0 = 0x0,
const char* filename0 = 0x0, bool start_dx = true, int nx = 40, int ny = 40,
int nz = 40) const ;
// Member arithmetics
// ------------------
public:
void operator+=(const Scalar &) ; ///< += Scalar
void operator-=(const Scalar &) ; ///< -= Scalar
void operator*=(const Scalar &) ; ///< *= Scalar
// Manipulation of spectral bases
// ------------------------------
/** Sets the spectral bases of the \c Valeur \c va to the standard ones
* for a scalar field
*/
virtual void std_spectral_base() ;
/** Sets the spectral bases of the \c Valeur \c va to the standard odd ones
* for a scalar field
*/
virtual void std_spectral_base_odd() ;
/** Sets the spectral bases of the \c Valeur \c va
*/
void set_spectral_base(const Base_val& ) ;
/// Returns the spectral bases of the \c Valeur \c va
const Base_val& get_spectral_base( ) const {return va.base ;} ;
/** Modifies the \c dzpuis flag.
* NB: this method does not change the field values stored in
* the compactified external domain (use methods \c dec_dzpuis() ,
* etc... for this purpose).
*/
void set_dzpuis(int ) ;
/** Asymptotic expansion at r = infinity.
*
* Determines the coefficients \f$a_k(\theta, \phi)\f$ of the expansion
* \f[
* \sum_{k=0}^n {a_k(\theta, \phi) \over r^k}
* \f]
* of \c *this when \f$r \rightarrow \infty\f$.
*
* @param n order of the expansion
* @param flag : output
* @return Array of \c n +1 \c Valeur s on \c mg->angu
* describing the coefficients \f$a_k(\theta, \phi)\f$.
* This array is allocated by the routine.
*
*/
Valeur** asymptot(int n, const int flag = 0) const ;
// PDE resolution
// --------------
public:
/** Solves the scalar Poisson equation with \c *this as a source.
* The source \f$\sigma\f$ of the equation \f$\Delta u = \sigma\f$ is
* represented by the \c Scalar \c *this .
* Note that \c dzpuis must be equal to 2 or 4, i.e. that the
* quantity stored in \c *this is in fact \f$r^2 \sigma\f$ or
* \f$r^4 \sigma\f$ in the compactified external domain.
* The solution \e u with the boundary condition \e u =0 at spatial
* infinity is the returned \c Scalar.
*/
Scalar poisson() const ;
/** Solves the scalar Poisson equation with \c *this as a source
* (version with parameters to control the resolution).
* The source \f$\sigma\f$ of the equation \f$\Delta u = \sigma\f$ is
* represented by the \c Scalar \c *this .
* Note that \c dzpuis must be equal to 2 or 4, i.e. that the
* quantity stored in \c *this is in fact \f$r^2 \sigma\f$ or
* \f$r^4 \sigma\f$ in the compactified external domain.
* @param par [input/output] possible parameters
* @param uu [input/output] solution \e u with the boundary condition
* \e u =0 at spatial infinity.
*/
void poisson(Param& par, Scalar& uu) const ;
/** Solves the scalar Poisson equation with \c *this as a source using a real Tau method
* The source \f$\sigma\f$ of the equation \f$\Delta u = \sigma\f$ is
* represented by the \c Scalar \c *this .
* Note that \c dzpuis must be equal to 2, 3 or 4, i.e. that the
* quantity stored in \c *this is in fact \f$r^2 \sigma\f$, \f$r^3 \sigma\f$ or
* \f$r^4 \sigma\f$ in the compactified external domain.
* The solution \e u with the boundary condition \e u =0 at spatial
* infinity is the returned \c Scalar.
*/
Scalar poisson_tau() const ;
/** Solves the scalar Poisson equation with \c *this as a source using a real Tau method
* (version with parameters to control the resolution)
* The source \f$\sigma\f$ of the equation \f$\Delta u = \sigma\f$ is
* represented by the \c Scalar \c *this .
* Note that \c dzpuis must be equal to 2, 3 or 4, i.e. that the
* quantity stored in \c *this is in fact \f$r^2 \sigma\f$, \f$r^3 \sigma\f$ or
* \f$r^4 \sigma\f$ in the compactified external domain.
* The solution \e u with the boundary condition \e u =0 at spatial
* infinity is the returned \c Scalar.
*/
void poisson_tau(Param& par, Scalar& uu) const ;
/**
* Is identicall to \c Scalar::poisson() . The regularity condition at the
* origin is replace by a boundary condition of the Dirichlet type.
*
* @param limite [input] : angular function. The boundary condition is
* given by \c limite[num] .
* @param num [input] : index of the boudary at which the condition is to
* be fullfilled.
*
* More precisely we impose the solution is equal to \c limite[num] at the
* boundary between the domains \c num and \c num+1 (the latter one being
* a shell).
*
*/
Scalar poisson_dirichlet (const Valeur& limite, int num) const ;
/**
* Idem as \c Scalar::poisson_dirichlet , the boundary condition being on
* the radial derivative of the solution.
*/
Scalar poisson_neumann (const Valeur&, int) const ;
/**
* Is identicall to \c Scalar::poisson() . The regularity condition at the
* origin is replace by a mixed boundary condition (Dirichlet + Neumann).
*
* @param limite [input] : angular function. The boundary condition is
* given by \c limite[num] .
* @param num [input] : index of the boudary at which the condition is to
* be fullfilled.
* @param fact_dir [input] : double in front of \f$\Psi\f$ (if \f$\Psi\f$
* is the variable solved).
* @param fact_neu [input] : double in front of the radial derivative
* of \f$\Psi\f$.
*
* More precisely we impose \f$ fact\_dir.\Psi + fact\_neu.\frac{\partial
* \Phi}{\partial r}\f$ is equal to the source at the
* boundary between the domains \c num and \c num+1 (the latter one being
* a shell).
*/
Scalar poisson_dir_neu (const Valeur& limite , int num,
double fact_dir, double fact_neu) const ;
/**
* Idem as \c Scalar::poisson_dirichlet , the boundary condition being on
* both the function and its radial derivative. The boundary condition
* at infinity is relaxed.
*/
Scalar poisson_frontiere_double (const Valeur&, const Valeur&, int) const ;
/** Solves the scalar Poisson equation with \c *this as a source
* (version with parameters to control the resolution).
* The source \f$\sigma\f$ of the equation \f$\Delta u = \sigma\f$ is
* represented by the \c Scalar \c *this .
* The regularized source
* \f$\sigma_{\rm regu} = \sigma - \sigma_{\rm div}\f$
* is constructed and solved.
* Note that \c dzpuis must be equal to 2 or 4, i.e. that the
* quantity stored in \c *this is in fact \f$r^2 \sigma\f$ or
* \f$r^4 \sigma\f$ in the compactified external domain.
* @param k_div [input] regularization degree of the procedure
* @param nzet [input] number of domains covering the star
* @param unsgam1 [input] parameter \f$1/(\gamma-1)\f$ where \f$\gamma\f$
* denotes the adiabatic index
* @param par [input/output] possible parameters
@param uu [input/output] solution
* @param uu_regu [output] solution of the regular part of
* the source.
* @param uu_div [output] solution of the diverging part of
* the source.
* @param duu_div [output] derivative of the diverging potential.
* @param source_regu [output] regularized source
* @param source_div [output] diverging part of the source
*/
void poisson_regular(int k_div, int nzet, double unsgam1, Param& par,
Scalar& uu, Scalar& uu_regu, Scalar& uu_div,
Tensor& duu_div,
Scalar& source_regu, Scalar& source_div) const ;
/** Checks if a Poisson equation with \c *this as a source
* has been correctly solved.
*
* @param uu [input] Solution \e u of the Poisson equation
* \f$\Delta u = \sigma\f$, \f$\sigma\f$ being
* represented by the \c Scalar \c *this .
*
* @param ostr [input/output] Output stream used for displaying
* \c err .
*
* @param detail [input] \li if \c true displays \c err(0,*) ,
* \c err(1,*) and \c err(2,*)
* \li if \c false (default), displays only
* the relative error \c err(0,*) .
*
* @return 2-D \c Tbl \c err decribing the errors in each
* domain:
* \li \c err(0,l) : Relative error in domain no. \c l ,
* defined as the maximum value of
* \f$|\Delta u - \sigma|\f$ in that domain divided by \e M ,
* where \e M is the maximum value of \f$|\sigma|\f$
* over all domains if \c dzpuis = 0 or \f$\sigma\f$ is
* zero in the compactified external domain (CED). If
* \c dzpuis != 0 and \f$\sigma\f$ does not vanish in the
* CED, the value of \e M used in the
* non-compactified domains is the maximum value over
* these domains, whereas the value of \e M used in the
* compactified external domain is the maximum value
* on that particular domain.
* \li \c err(1,l) : Maximum value of the absolute error
* \f$|\Delta u - \sigma|\f$ in domain no. \c l
* \li \c err(2,l) : Maximum value of \f$|\sigma|\f$ in domain
* no. \c l
*/
Tbl test_poisson(const Scalar& uu, ostream& ostr,
bool detail = false) const ;
/** Solves the (generalized) angular Poisson equation with \c *this
* as source.
* The generalized angular Poisson equation is
* \f$\Delta_{\theta\varphi} u + \lambda u = \sigma\f$,
* where \f$\Delta_{\theta\varphi} u := \frac{\partial^2 u}
* {\partial \theta^2} + \frac{1}{\tan \theta} \frac{\partial u}
* {\partial \theta} +\frac{1}{\sin^2 \theta}\frac{\partial^2 u}
* {\partial \varphi^2}\f$.
*
* @param lambda [input] coefficient \f$\lambda\f$ in the above equation
* (default value = 0)
* @return solution \e u .
*/
Scalar poisson_angu(double lambda =0) const ;
/** Performs one time-step integration (from \f$t=J \to J+1\f$) of the
* scalar d'Alembert equation with \c *this being the value of
* the function \e f at time \e J .
*
* Works only with an affine mapping (class \c Map_af ) and,
* if the last domain is a compactified one, it simply copies
* the value of the field in this last domain at the time-step \e J
* to the last domain of the returned solution.
* @param par [input/output] possible parameters to control the
* resolution of the d'Alembert equation:
* \li \c par.get_double(0) : [input] the time step \e dt ,
* \li \c par.get_int(0) : [input] the type of boundary conditions
* set at the outer boundary (0 : reflexion, 1 : Sommerfeld
* outgoing wave, valid only for \e l=0 components, 2 : Bayliss
* \& Turkel outgoing wave, valid for \e l=0, 1, 2 components)
* \li \c par.get_int_mod(0) : [input/output] set to 0 at first
* call, is used as a working flag after (must not be modified after
* first call)
* \li \c par.get_tensor_mod(0) : [input] (optional) if the wave
* equation is on a curved space-time, this is the potential in front
* of the Laplace operator. It has to be a Scalar and updated at
* every time-step (for a potential depending on time).
* Note: there are many other working objects attached to this
* \c Param , so one should not modify it.
* There should be exactly one \c Param for each wave equation to be
* solved.
* @param fJm1 [input] solution \f$f^{J-1}\f$ at time \e J-1
* @param source [input] source \f$\sigma\f$ of the d'Alembert equation
* \f$\diamond u = \sigma\f$.
* @return solution \f$f^{J+1}\f$ at time \e J+1
* with boundary conditions defined by \c par.get_int(0) .
*/
Scalar avance_dalembert(Param& par, const Scalar& fJm1, const Scalar& source)
const ;
/**
* Resolution of a general elliptic equation, putting zero at infinity.
* @param params [input] the operators and variables to be used.
**/
Scalar sol_elliptic(Param_elliptic& params) const ;
/**
* Resolution of a general elliptic equation, putting zero at infinity
* and with inner boundary conditions.
* @param params [input] the operators and variables to be used.
* @param bound [input] : the boundary condition
* @param fact_dir : 1 Dirchlet condition, 0 Neumann condition
* @param fact_neu : 0 Dirchlet condition, 1 Neumann condition
**/
Scalar sol_elliptic_boundary(Param_elliptic& params, const Mtbl_cf& bound,
double fact_dir, double fact_neu) const ;
/** Resolution of general elliptic equation, with inner boundary conditions as Scalars
* on mono-domain angulare grids
**/
Scalar sol_elliptic_boundary(Param_elliptic& params, const Scalar& bound,
double fact_dir, double fact_neu) const ;
/** Solves the scalar 2-dimensional elliptic equation
* with \c *this as a source.
* Note that \c dzpuis must be equal to 2, 3 or 4, i.e.
* The solution \e u with the boundary condition \e u =0 at spatial
* infinity is the returned \c Scalar.
*/
Scalar sol_elliptic_2d(Param_elliptic&) const ;
/** Solves a pseudo-1d elliptic equation
* with \c *this as a source.
*
*/
Scalar sol_elliptic_pseudo_1d(Param_elliptic&) const ;
/**
* Resolution of a general elliptic equation, putting a given value
* at the outermost
* shell and not solving in the compactified domain.
* @param params [input] the operators and variables to be used.
* @param val [input] value at the last shell.
**/
Scalar sol_elliptic_no_zec(Param_elliptic& params, double val = 0) const ;
/**
* Resolution of a general elliptic equation solving in the
* compactified domain and putting a given value at the inner boundary.
* @param params [input] the operators and variables to be used.
* @param val [input] value at the inner boundary of the external domain.
**/
Scalar sol_elliptic_only_zec(Param_elliptic& params, double val) const ;
/**
* General elliptic solver.
* The equation is not solved in the compactified domain and the
* matching is done with an homogeneous solution.
* @param params [input] the operators and variables to be used.
* @param coef [output] : coefficient of the oscillatory solution
* in the external domain.
* @param phases [output] : phases (i.e. choice of the homogeneous solution to match with).
**/
Scalar sol_elliptic_sin_zec(Param_elliptic& params, double* coefs, double* phases) const ;
/**
* Resolution of a general elliptic equation fixing the dericative at
* the origin and relaxing one continuity condition.
*
* @param val [input] value of the derivative.
* @param params [input] the operators and variables to be used.
**/
Scalar sol_elliptic_fixe_der_zero(double val,
Param_elliptic& params) const ;
/**
* Resolution of a divergence-like equation.
* The equation solved reads: \f$ \frac{\partial \phi}{\partial r} + \frac{n}{r}
* \phi = \sigma\f$, with \f$\phi\f$ the unknown and \f$\sigma\f$ the source
* represented by \c this.
*
* @param n [input] the coefficient in front of the 1/r term.
*
* @return the solution to the equation.
*/
Scalar sol_divergence(int n) const ;
// Import from other mapping
// -------------------------
/** Assignment to another \c Scalar defined on a different mapping.
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param ci [input] \c Scalar to be imported.
*/
void import(const Scalar& ci) ;
/** Assignment to another \c Scalar defined on a different mapping.
* Case where the \c Scalar is symmetric with respect to the plane y=0.
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param ci [input] \c Scalar to be imported.
*/
void import_symy(const Scalar& ci) ;
/** Assignment to another \c Scalar defined on a different mapping.
* Case where the \c Scalar is antisymmetric with respect to the
* plane y=0.
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param ci [input] \c Scalar to be imported.
*/
void import_asymy(const Scalar& ci) ;
/** Assignment to another \c Scalar defined on a different mapping.
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param nzet [input] Number of domains of the destination
* mapping (i.e. \c this->mp ) where the
* importation is performed: the assignment
* is done for the domains whose indices are
* between 0 and \c nzet-1 . In the other
* domains, \c *this is set to zero.
* @param ci [input] \c Scalar to be imported.
*/
void import(int nzet, const Scalar& ci) ;
/** Assignment to another \c Scalar defined on a different mapping.
* Case where the \c Scalar is symmetric with respect to the plane y=0.
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param nzet [input] Number of domains of the destination
* mapping (i.e. \c this->mp ) where the
* importation is performed: the assignment
* is done for the domains whose indices are
* between 0 and \c nzet-1 . In the other
* domains, \c *this is set to zero.
* @param ci [input] \c Scalar to be imported.
*/
void import_symy(int nzet, const Scalar& ci) ;
/** Assignment to another \c Scalar defined on a different mapping.
* Case where the \c Scalar is antisymmetric with respect to the
* plane y=0.
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param nzet [input] Number of domains of the destination
* mapping (i.e. \c this->mp ) where the
* importation is performed: the assignment
* is done for the domains whose indices are
* between 0 and \c nzet-1 . In the other
* domains, \c *this is set to zero.
* @param ci [input] \c Scalar to be imported.
*/
void import_asymy(int nzet, const Scalar& ci) ;
protected:
/** Assignment to another \c Scalar defined on a different mapping,
* when the two mappings do not have a particular relative orientation.
*
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param nzet [input] Number of domains of the destination
* mapping (i.e. \c this->mp ) where the
* importation is performed: the assignment
* is done for the domains whose indices are
* between 0 and \c nzet-1 . In the other
* domains, \c *this is set to zero.
* @param ci [input] \c Scalar to be imported.
*/
void import_gal(int nzet, const Scalar& ci) ;
/** Assignment to another \c Scalar defined on a different mapping,
* when the two mappings have aligned Cartesian axis.
*
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param nzet [input] Number of domains of the destination
* mapping (i.e. \c this->mp ) where the
* importation is performed: the assignment
* is done for the domains whose indices are
* between 0 and \c nzet-1 . In the other
* domains, \c *this is set to zero.
* @param ci [input] \c Scalar to be imported.
*/
void import_align(int nzet, const Scalar& ci) ;
/** Assignment to another \c Scalar defined on a different mapping,
* when the two mappings have anti-aligned Cartesian axis (i.e.
* \f$x_1 = - x_2\f$, \f$y_1 = - y_2\f$, \f$z_1 = z_2\f$).
*
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param nzet [input] Number of domains of the destination
* mapping (i.e. \c this->mp ) where the
* importation is performed: the assignment
* is done for the domains whose indices are
* between 0 and \c nzet-1 . In the other
* domains, \c *this is set to zero.
* @param ci [input] \c Scalar to be imported.
*/
void import_anti(int nzet, const Scalar& ci) ;
/** Assignment to another \c Scalar defined on a different mapping,
* when the two mappings have aligned Cartesian axis.
* Case where the \c Scalar is symmetric with respect to the plane y=0.
*
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param nzet [input] Number of domains of the destination
* mapping (i.e. \c this->mp ) where the
* importation is performed: the assignment
* is done for the domains whose indices are
* between 0 and \c nzet-1 . In the other
* domains, \c *this is set to zero.
* @param ci [input] \c Scalar to be imported.
*/
void import_align_symy(int nzet, const Scalar& ci) ;
/** Assignment to another \c Scalar defined on a different mapping,
* when the two mappings have anti-aligned Cartesian axis (i.e.
* \f$x_1 = - x_2\f$, \f$y_1 = - y_2\f$, \f$z_1 = z_2\f$).
* Case where the \c Scalar is symmetric with respect to the plane y=0.
*
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param nzet [input] Number of domains of the destination
* mapping (i.e. \c this->mp ) where the
* importation is performed: the assignment
* is done for the domains whose indices are
* between 0 and \c nzet-1 . In the other
* domains, \c *this is set to zero.
* @param ci [input] \c Scalar to be imported.
*/
void import_anti_symy(int nzet, const Scalar& ci) ;
/** Assignment to another \c Scalar defined on a different mapping,
* when the two mappings have aligned Cartesian axis.
* Case where the \c Scalar is antisymmetric with respect to the
* plane y=0.
*
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param nzet [input] Number of domains of the destination
* mapping (i.e. \c this->mp ) where the
* importation is performed: the assignment
* is done for the domains whose indices are
* between 0 and \c nzet-1 . In the other
* domains, \c *this is set to zero.
* @param ci [input] \c Scalar to be imported.
*/
void import_align_asymy(int nzet, const Scalar& ci) ;
/** Assignment to another \c Scalar defined on a different mapping,
* when the two mappings have anti-aligned Cartesian axis (i.e.
* \f$x_1 = - x_2\f$, \f$y_1 = - y_2\f$, \f$z_1 = z_2\f$).
* Case where the \c Scalar is antisymmetric with respect to the
* plane y=0.
*
* This assignment is performed point to point by means of the
* spectral expansion of the original \c Scalar.
* @param nzet [input] Number of domains of the destination
* mapping (i.e. \c this->mp ) where the
* importation is performed: the assignment
* is done for the domains whose indices are
* between 0 and \c nzet-1 . In the other
* domains, \c *this is set to zero.
* @param ci [input] \c Scalar to be imported.
*/
void import_anti_asymy(int nzet, const Scalar& ci) ;
friend Scalar operator-(const Scalar& ) ;
friend Scalar operator+(const Scalar&, const Scalar &) ;
friend Scalar operator+(const Scalar&, const Mtbl&) ;
friend Scalar operator+(const Scalar&, double ) ;
friend Scalar operator-(const Scalar &, const Scalar &) ;
friend Scalar operator-(const Scalar&, const Mtbl&) ;
friend Scalar operator-(const Scalar&, double ) ;
friend Scalar operator*(const Scalar &, const Scalar &) ;
friend Scalar operator%(const Scalar &, const Scalar &) ;
friend Scalar operator|(const Scalar &, const Scalar &) ;
friend Scalar operator*(const Mtbl&, const Scalar &) ;
friend Scalar operator*(double, const Scalar &) ;
friend Scalar operator/(const Scalar &, const Scalar &) ;
friend Scalar operator/(const Scalar &, const Mtbl &) ;
friend Scalar operator/(const Mtbl &, const Scalar &) ;
friend Scalar operator/(const Scalar&, double ) ;
friend Scalar operator/(double, const Scalar &) ;
friend Scalar sin(const Scalar& ) ;
friend Scalar cos(const Scalar& ) ;
friend Scalar tan(const Scalar& ) ;
friend Scalar asin(const Scalar& ) ;
friend Scalar acos(const Scalar& ) ;
friend Scalar atan(const Scalar& ) ;
friend Scalar exp(const Scalar& ) ;
friend Scalar Heaviside(const Scalar& ) ;
friend Scalar log(const Scalar& ) ;
friend Scalar log10(const Scalar& ) ;
friend Scalar sqrt(const Scalar& ) ;
friend Scalar racine_cubique (const Scalar& ) ;
friend Scalar pow(const Scalar& , int ) ;
friend Scalar pow(const Scalar& , double ) ;
friend Scalar abs(const Scalar& ) ;
friend double totalmax(const Scalar& ) ;
friend double totalmin(const Scalar& ) ;
friend Tbl max(const Scalar& ) ;
friend Tbl min(const Scalar& ) ;
friend Tbl norme(const Scalar& ) ;
friend Tbl diffrel(const Scalar& a, const Scalar& b) ;
friend Tbl diffrelmax(const Scalar& a, const Scalar& b) ;
};
ostream& operator<<(ostream& , const Scalar & ) ;
// Prototypage de l'arithmetique
/**
* \defgroup sal_mat Scalar mathematics
* \ingroup (tensor)
* @{
*/
Scalar operator+(const Scalar& ) ; ///< + Scalar
Scalar operator-(const Scalar& ) ; ///< \c - Scalar
Scalar operator+(const Scalar&, const Scalar &) ; ///< Scalar + Scalar
Scalar operator+(const Scalar&, const Mtbl&) ; ///< Scalar + Mbtl
Scalar operator+(const Mtbl&, const Scalar&) ; ///< Mtbl + Scalar
Scalar operator+(const Scalar&, double ) ; ///< Scalar + double
Scalar operator+(double, const Scalar& ) ; ///< double + Scalar
Scalar operator+(const Scalar&, int ) ; ///< Scalar + int
Scalar operator+(int, const Scalar& ) ; ///< int + Scalar
Scalar operator-(const Scalar &, const Scalar &) ; ///< Scalar - Scalar
Scalar operator-(const Scalar&, const Mtbl&) ; ///< Scalar - Mbtl
Scalar operator-(const Mtbl&, const Scalar&) ; ///< Mtbl - Scalar
Scalar operator-(const Scalar&, double ) ; ///< Scalar - double
Scalar operator-(double, const Scalar& ) ; ///< double - Scalar
Scalar operator-(const Scalar&, int ) ; ///< Scalar - int
Scalar operator-(int, const Scalar& ) ; ///< int - Scalar
Scalar operator*(const Scalar &, const Scalar &) ; ///< Scalar * Scalar
/// Scalar * Scalar with desaliasing
Scalar operator%(const Scalar &, const Scalar &) ;
/// Scalar * Scalar with desaliasing only in \e r
Scalar operator|(const Scalar &, const Scalar &) ;
Scalar operator*(const Mtbl&, const Scalar&) ; ///< Mtbl * Scalar
Scalar operator*(const Scalar&, const Mtbl&) ; ///< Scalar * Mtbl
Scalar operator*(const Scalar&, double ) ; ///< Scalar * double
Scalar operator*(double, const Scalar &) ; ///< double * Scalar
Scalar operator*(const Scalar&, int ) ; ///< Scalar * int
Scalar operator*(int, const Scalar& ) ; ///< int * Scalar
Scalar operator/(const Scalar &, const Scalar &) ; ///< Scalar / Scalar
Scalar operator/(const Scalar&, double ) ; ///< Scalar / double
Scalar operator/(double, const Scalar &) ; ///< double / Scalar
Scalar operator/(const Scalar&, int ) ; ///< Scalar / int
Scalar operator/(int, const Scalar &) ; ///< int / Scalar
Scalar operator/(const Scalar &, const Mtbl&) ; ///< Scalar / Mtbl
Scalar operator/(const Mtbl&, const Scalar &) ; ///< Mtbl / Scalar
Scalar sin(const Scalar& ) ; ///< Sine
Scalar cos(const Scalar& ) ; ///< Cosine
Scalar tan(const Scalar& ) ; ///< Tangent
Scalar asin(const Scalar& ) ; ///< Arcsine
Scalar acos(const Scalar& ) ; ///< Arccosine
Scalar atan(const Scalar& ) ; ///< Arctangent
Scalar exp(const Scalar& ) ; ///< Exponential
Scalar Heaviside(const Scalar& ) ; ///< Heaviside function
Scalar log(const Scalar& ) ; ///< Neperian logarithm
Scalar log10(const Scalar& ) ; ///< Basis 10 logarithm
Scalar sqrt(const Scalar& ) ; ///< Square root
Scalar racine_cubique (const Scalar& ) ; ///< Cube root
Scalar pow(const Scalar& , int ) ; ///< Power \f${\tt Scalar}^{\tt int}\f$
Scalar pow(const Scalar& , double ) ; ///< Power \f${\tt Scalar}^{\tt double}\f$
Scalar abs(const Scalar& ) ; ///< Absolute value
/**
* Maximum values of a \c Scalar in each domain.
* @return 1-D \c Tbl of size the number of domains, the elements of which
* are the set of the maximum values in each domain.
*/
double totalmax(const Scalar& ) ;
/**
* Minimum values of a \c Scalar in each domain.
* @return 1-D \c Tbl of size the number of domains, the elements of which
* are the set of the minimum values in each domain.
*/
double totalmin(const Scalar& ) ;
/**
* Maximum values of a \c Scalar in each domain.
* @return 1-D \c Tbl of size the number of domains, the elements of which
* are the set of the maximum values in each domain.
*/
Tbl max(const Scalar& ) ;
/**
* Minimum values of a \c Scalar in each domain.
* @return 1-D \c Tbl of size the number of domains, the elements of which
* are the set of the minimum values in each domain.
*/
Tbl min(const Scalar& ) ;
/**
* Sums of the absolute values of all the values of the \c Scalar
* in each domain.
* @return 1-D \c Tbl of size the number of domains, the elements of which
* are the set of the sums of the absolute values in each domain.
*/
Tbl norme(const Scalar& ) ;
/**
* Relative difference between two \c Scalar (norme version).
* @return 1-D \c Tbl of size the number of domains, the elements of which
* are \c norme[a(l)-b(l)]/norme[b(l)] if \c b(l)!=0 and
* \c norme[a(l)-b(l)] if \c b(l)=0 , where \c a(l) and
* \c b(l) denote symbolically the values of \c a and \c b
* in domain no. \c l .
*/
Tbl diffrel(const Scalar& a, const Scalar& b) ;
/**
* Relative difference between two \c Scalar (max version).
* @return 1-D \c Tbl of size the number of domains, the elements of which
* are \c max[abs(a(l)-b(l))]/max[abs(b(l))] if \c b(l)!=0 and
* \c max[abs(a(l)-b(l))] if \c b(l)=0 , where \c a(l) and
* \c b(l) denote symbolically the values of \c a and \c b
* in domain no. \c l .
*/
Tbl diffrelmax(const Scalar& a, const Scalar& b) ;
/**
* Applies an exponential filter in angular directions in all domains.
* (see \c Scalar:exponential_filter_ylm ).
*/
void exp_filter_ylm_all_domains(Scalar& ss, int p, double alpha=-16.) ;
/** @} */
}
#endif
|