1696 {
1697 constexpr MInt nDim = 3;
1698
1699 MInt hasConnectionInfo = 0;
1702 MBool attributeExists = pio.hasAttribute(
"hasConnectionInfo",
"Connectivity");
1703 if(attributeExists) {
1704 m_log <<
"Grid file has connection info!" << endl;
1705 pio.getAttribute(&hasConnectionInfo, "hasConnectionInfo", "Connectivity");
1706 }
1708 } else {
1710 }
1711
1712 MBool readInConnectionInfo =
true;
1713
1714 if(!hasConnectionInfo) {
1715 readInConnectionInfo = false;
1716 } else {
1718 readInConnectionInfo =
1719 Context::getSolverProperty<MBool>(
"readInConnectionInfo",
solverId(), AT_, &readInConnectionInfo);
1720 }
1721 }
1722
1723 if(readInConnectionInfo) {
1724 m_log <<
"Connection info available, reading connection information from grid file!" << endl;
1726 } else {
1728 cout << "Starting grid connection info search..." << endl;
1729 }
1730 m_log <<
" Starting grid connection info search..." << endl;
1731 MInt offset[3], size[3], count = 0;
1732
1733
1734
1735 MInt totalGridCells = 0;
1736 MInt** totalGridBlockDim =
nullptr;
1737 MInt* totalGridBlockCells =
nullptr;
1738
1741
1744 for(
MInt dim = 0; dim < 3; ++dim) {
1745 temp *=
m_grid->getBlockNoCells(i, dim) + 1;
1746 totalGridBlockDim[i][dim] =
m_grid->getBlockNoCells(i, dim) + 1;
1747 }
1748 totalGridBlockCells[i] = temp;
1749 if(temp != 1) totalGridCells += temp;
1750 }
1751
1754 countNode = countNode + 2 * totalGridBlockDim[i][1] * totalGridBlockDim[i][2]
1755 + 2 * totalGridBlockDim[i][2] * (totalGridBlockDim[i][0] - 2)
1756 + 2 * (totalGridBlockDim[i][0] - 2) * (totalGridBlockDim[i][1] - 2);
1757 }
1758
1760 coordinates.fill(-1.01010101);
1761 std::vector<pointType<nDim>> nodeMap;
1762 nodeMap.resize(countNode);
1763
1764
1765
1766 const MInt converter[] = {2, 1, 0};
1767 auto getCurrentWindow = [
this, converter](
const MInt blockId,
const MInt*
const offset_,
const MInt*
const size_) {
1768 std::vector<MInt> currentWindows;
1770 if(
inputWindows[windowId]->blockId != blockId)
continue;
1772
1773 MBool takeIt =
true;
1774 for(
MInt dim = 0; dim < nDim; ++dim) {
1775 if(
inputWindows[windowId]->startindex[dim] >= offset_[converter[dim]] + size_[converter[dim]]
1776 ||
inputWindows[windowId]->endindex[dim] < offset_[converter[dim]]) {
1777 takeIt = false;
1778 break;
1779 }
1780 }
1781 if(takeIt) currentWindows.push_back(windowId);
1782 }
1783 return currentWindows;
1784 };
1786 auto getBCFromWindow = [
this, &pointsFoundPerWindow](
const MInt idx[nDim],
1787 const std::vector<MInt>& currentWindows) {
1788 std::vector<MInt> BCset;
1789
1790
1791
1792
1793 for(auto windowId : currentWindows) {
1794 MBool takeIt =
true;
1795 for(
MInt dim = 0; dim < nDim; ++dim) {
1797 takeIt = false;
1798 break;
1799 }
1800 }
1801 if(takeIt) {
1802 ++pointsFoundPerWindow[windowId];
1804
1805
1806 }
1807 }
1808
1809
1810
1811
1812 if(BCset.size() > 3)
mTerm(1,
"");
1813 if(BCset.size() == 0) BCset.push_back(6000);
1814 MInt BC[3] = {0, 0, 0};
1816 for(auto it = BCset.begin(); it != BCset.end(); ++it)
1817 BC[cnt++] = *it;
1818 return std::tuple<MInt, MInt, MInt>(BC[0], BC[1], BC[2]);
1819 };
1820
1822
1824 cout << "Reading in all coordinates of all block faces" << endl;
1825 }
1826
1828
1831 stringstream number;
1832 number << i << "/";
1833 bName += number.str();
1834
1835 ParallelIo::size_type ioOffset[3] = {0, 0, 0};
1836 ParallelIo::size_type ioSize[3] = {0, 0, 0};
1837
1841
1842
1844 ioOffset[2] = 0;
1845 ioSize[2] = 1;
1846 ioOffset[1] = 0;
1847 ioSize[1] =
m_grid->getBlockNoCells(i, 1) + 1;
1848 ioOffset[0] = 0;
1849 ioSize[0] =
m_grid->getBlockNoCells(i, 0) + 1;
1850
1851 pio.readArray(&coordinates(0, memsize), bName, "x", 3, ioOffset, ioSize);
1852 pio.readArray(&coordinates(1, memsize), bName, "y", 3, ioOffset, ioSize);
1853 pio.readArray(&coordinates(2, memsize), bName, "z", 3, ioOffset, ioSize);
1854 } else {
1855 ioOffset[2] = 0;
1856 ioSize[2] = 0;
1857 ioOffset[1] = 0;
1858 ioSize[1] = 0;
1859 ioOffset[0] = 0;
1860 ioSize[0] = 0;
1861
1863 pio.readArray(&empty, bName, "x", 3, ioOffset, ioSize);
1864 pio.readArray(&empty, bName, "y", 3, ioOffset, ioSize);
1865 pio.readArray(&empty, bName, "z", 3, ioOffset, ioSize);
1866 }
1867 for(int dim = 0; dim < 3; ++dim) {
1868 offset[dim] = ioOffset[dim];
1869 size[dim] = ioSize[dim];
1870 }
1871
1872 std::vector<MInt> currentWindows = getCurrentWindow(i, offset, size);
1873
1874
1875
1876 for(
MInt a = offset[2];
a < offset[2] + size[2]; ++
a) {
1877 for(
MInt c = offset[0]; c < offset[0] + size[0]; ++c) {
1878 for(
MInt b = offset[1];
b < offset[1] + size[1]; ++
b) {
1879 nodeMap[count].blockId = i;
1880 nodeMap[count].pos[0] =
a;
1881 nodeMap[count].pos[1] =
b;
1882 nodeMap[count].pos[2] = c;
1883 std::tie(nodeMap[count].BC[0], nodeMap[count].BC[1], nodeMap[count].BC[2]) =
1884 getBCFromWindow(nodeMap[count].pos, currentWindows);
1885 count++;
1886 }
1887 }
1888 }
1889 memsize += size[0] * size[1] * size[2];
1890
1894
1896 ioOffset[2] =
m_grid->getBlockNoCells(i, 2);
1897 ioSize[2] = 1;
1898 ioOffset[1] = 0;
1899 ioSize[1] =
m_grid->getBlockNoCells(i, 1) + 1;
1900 ioOffset[0] = 0;
1901 ioSize[0] =
m_grid->getBlockNoCells(i, 0) + 1;
1902 pio.readArray(&coordinates(0, memsize), bName, "x", 3, ioOffset, ioSize);
1903 pio.readArray(&coordinates(1, memsize), bName, "y", 3, ioOffset, ioSize);
1904 pio.readArray(&coordinates(2, memsize), bName, "z", 3, ioOffset, ioSize);
1905 } else {
1906 ioOffset[2] = 0;
1907 ioSize[2] = 0;
1908 ioOffset[1] = 0;
1909 ioSize[1] = 0;
1910 ioOffset[0] = 0;
1911 ioSize[0] = 0;
1913 pio.readArray(&empty, bName, "x", 3, ioOffset, ioSize);
1914 pio.readArray(&empty, bName, "y", 3, ioOffset, ioSize);
1915 pio.readArray(&empty, bName, "z", 3, ioOffset, ioSize);
1916 }
1917 for(int dim = 0; dim < 3; ++dim) {
1918 offset[dim] = ioOffset[dim];
1919 size[dim] = ioSize[dim];
1920 }
1921
1922 currentWindows = getCurrentWindow(i, offset, size);
1923
1924 for(
MInt a = offset[2];
a < offset[2] + size[2]; ++
a) {
1925 for(
MInt c = offset[0]; c < offset[0] + size[0]; ++c) {
1926 for(
MInt b = offset[1];
b < offset[1] + size[1]; ++
b) {
1927 nodeMap[count].blockId = i;
1928 nodeMap[count].pos[0] =
a;
1929 nodeMap[count].pos[1] =
b;
1930 nodeMap[count].pos[2] = c;
1931 std::tie(nodeMap[count].BC[0], nodeMap[count].BC[1], nodeMap[count].BC[2]) =
1932 getBCFromWindow(nodeMap[count].pos, currentWindows);
1933 count++;
1934 }
1935 }
1936 }
1937 memsize += size[0] * size[1] * size[2];
1938
1942
1944 ioOffset[2] = 1;
1945 ioSize[2] =
m_grid->getBlockNoCells(i, 2) - 1;
1946 ioOffset[1] = 0;
1947 ioSize[1] = 1;
1948 ioOffset[0] = 0;
1949 ioSize[0] =
m_grid->getBlockNoCells(i, 0) + 1;
1950 pio.readArray(&coordinates(0, memsize), bName, "x", 3, ioOffset, ioSize);
1951 pio.readArray(&coordinates(1, memsize), bName, "y", 3, ioOffset, ioSize);
1952 pio.readArray(&coordinates(2, memsize), bName, "z", 3, ioOffset, ioSize);
1953 } else {
1954 ioOffset[2] = 0;
1955 ioSize[2] = 0;
1956 ioOffset[1] = 0;
1957 ioSize[1] = 0;
1958 ioOffset[0] = 0;
1959 ioSize[0] = 0;
1961 pio.readArray(&empty, bName, "x", 3, ioOffset, ioSize);
1962 pio.readArray(&empty, bName, "y", 3, ioOffset, ioSize);
1963 pio.readArray(&empty, bName, "z", 3, ioOffset, ioSize);
1964 }
1965 for(int dim = 0; dim < 3; ++dim) {
1966 offset[dim] = ioOffset[dim];
1967 size[dim] = ioSize[dim];
1968 }
1969
1970 currentWindows = getCurrentWindow(i, offset, size);
1971
1972 for(
MInt b = offset[1];
b < offset[1] + size[1]; ++
b) {
1973 for(
MInt c = offset[0]; c < offset[0] + size[0]; ++c) {
1974 for(
MInt a = offset[2];
a < offset[2] + size[2]; ++
a) {
1975 nodeMap[count].blockId = i;
1976 nodeMap[count].pos[0] =
a;
1977 nodeMap[count].pos[1] =
b;
1978 nodeMap[count].pos[2] = c;
1979 std::tie(nodeMap[count].BC[0], nodeMap[count].BC[1], nodeMap[count].BC[2]) =
1980 getBCFromWindow(nodeMap[count].pos, currentWindows);
1981 count++;
1982 }
1983 }
1984 }
1985 memsize += size[0] * size[1] * size[2];
1986
1990
1992 ioOffset[2] = 1;
1993 ioSize[2] =
m_grid->getBlockNoCells(i, 2) - 1;
1994 ioOffset[1] =
m_grid->getBlockNoCells(i, 1);
1995 ioSize[1] = 1;
1996 ioOffset[0] = 0;
1997 ioSize[0] =
m_grid->getBlockNoCells(i, 0) + 1;
1998 pio.readArray(&coordinates(0, memsize), bName, "x", 3, ioOffset, ioSize);
1999 pio.readArray(&coordinates(1, memsize), bName, "y", 3, ioOffset, ioSize);
2000 pio.readArray(&coordinates(2, memsize), bName, "z", 3, ioOffset, ioSize);
2001 } else {
2002 ioOffset[2] = 0;
2003 ioSize[2] = 0;
2004 ioOffset[1] = 0;
2005 ioSize[1] = 0;
2006 ioOffset[0] = 0;
2007 ioSize[0] = 0;
2009 pio.readArray(&empty, bName, "x", 3, ioOffset, ioSize);
2010 pio.readArray(&empty, bName, "y", 3, ioOffset, ioSize);
2011 pio.readArray(&empty, bName, "z", 3, ioOffset, ioSize);
2012 }
2013 for(int dim = 0; dim < 3; ++dim) {
2014 offset[dim] = ioOffset[dim];
2015 size[dim] = ioSize[dim];
2016 }
2017
2018 currentWindows = getCurrentWindow(i, offset, size);
2019
2020 for(
MInt b = offset[1];
b < offset[1] + size[1]; ++
b) {
2021 for(
MInt c = offset[0]; c < offset[0] + size[0]; ++c) {
2022 for(
MInt a = offset[2];
a < offset[2] + size[2]; ++
a) {
2023 nodeMap[count].blockId = i;
2024 nodeMap[count].pos[0] =
a;
2025 nodeMap[count].pos[1] =
b;
2026 nodeMap[count].pos[2] = c;
2027 std::tie(nodeMap[count].BC[0], nodeMap[count].BC[1], nodeMap[count].BC[2]) =
2028 getBCFromWindow(nodeMap[count].pos, currentWindows);
2029 count++;
2030 }
2031 }
2032 }
2033 memsize += size[0] * size[1] * size[2];
2034
2038
2040 ioOffset[2] = 1;
2041 ioSize[2] =
m_grid->getBlockNoCells(i, 2) - 1;
2042 ioOffset[1] = 1;
2043 ioSize[1] =
m_grid->getBlockNoCells(i, 1) - 1;
2044 ioOffset[0] = 0;
2045 ioSize[0] = 1;
2046 pio.readArray(&coordinates(0, memsize), bName, "x", 3, ioOffset, ioSize);
2047 pio.readArray(&coordinates(1, memsize), bName, "y", 3, ioOffset, ioSize);
2048 pio.readArray(&coordinates(2, memsize), bName, "z", 3, ioOffset, ioSize);
2049 } else {
2050 ioOffset[2] = 0;
2051 ioSize[2] = 0;
2052 ioOffset[1] = 0;
2053 ioSize[1] = 0;
2054 ioOffset[0] = 0;
2055 ioSize[0] = 0;
2057 pio.readArray(&empty, bName, "x", 3, ioOffset, ioSize);
2058 pio.readArray(&empty, bName, "y", 3, ioOffset, ioSize);
2059 pio.readArray(&empty, bName, "z", 3, ioOffset, ioSize);
2060 }
2061 for(int dim = 0; dim < 3; ++dim) {
2062 offset[dim] = ioOffset[dim];
2063 size[dim] = ioSize[dim];
2064 }
2065
2066 currentWindows = getCurrentWindow(i, offset, size);
2067
2068 for(
MInt c = offset[0]; c < offset[0] + size[0]; ++c) {
2069 for(
MInt b = offset[1];
b < offset[1] + size[1]; ++
b) {
2070 for(
MInt a = offset[2];
a < offset[2] + size[2]; ++
a) {
2071 nodeMap[count].blockId = i;
2072 nodeMap[count].pos[0] =
a;
2073 nodeMap[count].pos[1] =
b;
2074 nodeMap[count].pos[2] = c;
2075 std::tie(nodeMap[count].BC[0], nodeMap[count].BC[1], nodeMap[count].BC[2]) =
2076 getBCFromWindow(nodeMap[count].pos, currentWindows);
2077 count++;
2078 }
2079 }
2080 }
2081 memsize += size[0] * size[1] * size[2];
2082
2086
2088 ioOffset[2] = 1;
2089 ioSize[2] =
m_grid->getBlockNoCells(i, 2) - 1;
2090 ioOffset[1] = 1;
2091 ioSize[1] =
m_grid->getBlockNoCells(i, 1) - 1;
2092 ioOffset[0] =
m_grid->getBlockNoCells(i, 0);
2093 ioSize[0] = 1;
2094 pio.readArray(&coordinates(0, memsize), bName, "x", 3, ioOffset, ioSize);
2095 pio.readArray(&coordinates(1, memsize), bName, "y", 3, ioOffset, ioSize);
2096 pio.readArray(&coordinates(2, memsize), bName, "z", 3, ioOffset, ioSize);
2097 } else {
2098 ioOffset[2] = 0;
2099 ioSize[2] = 0;
2100 ioOffset[1] = 0;
2101 ioSize[1] = 0;
2102 ioOffset[0] = 0;
2103 ioSize[0] = 0;
2105 pio.readArray(&empty, bName, "x", 3, ioOffset, ioSize);
2106 pio.readArray(&empty, bName, "y", 3, ioOffset, ioSize);
2107 pio.readArray(&empty, bName, "z", 3, ioOffset, ioSize);
2108 }
2109 for(int dim = 0; dim < 3; ++dim) {
2110 offset[dim] = ioOffset[dim];
2111 size[dim] = ioSize[dim];
2112 }
2113
2114 currentWindows = getCurrentWindow(i, offset, size);
2115
2116 for(
MInt c = offset[0]; c < offset[0] + size[0]; ++c) {
2117 for(
MInt b = offset[1];
b < offset[1] + size[1]; ++
b) {
2118 for(
MInt a = offset[2];
a < offset[2] + size[2]; ++
a) {
2119 nodeMap[count].blockId = i;
2120 nodeMap[count].pos[0] =
a;
2121 nodeMap[count].pos[1] =
b;
2122 nodeMap[count].pos[2] = c;
2123 std::tie(nodeMap[count].BC[0], nodeMap[count].BC[1], nodeMap[count].BC[2]) =
2124 getBCFromWindow(nodeMap[count].pos, currentWindows);
2125 count++;
2126 }
2127 }
2128 }
2129 memsize += size[0] * size[1] * size[2];
2130 }
2131
2132
2135 MInt noWindowPoints = 1;
2136 for(
MInt dim = 0; dim < nDim; ++dim) {
2138 }
2139
2140 if(noWindowPoints > pointsFoundPerWindow[windowId])
mTerm(1,
"noWindowPOints > pointsFoundPerWindow");
2141 }
2142
2143
2149
2150
2151 MPI_Datatype matrix;
2155
2160 cout << "Building up connections for multiblock connection search" << endl;
2161 }
2162 vector<Point<3>> pts;
2163 for(
MInt j = 0; j < count; ++j) {
2164 Point<3> a(coordinates(0, j), coordinates(1, j), coordinates(2, j));
2166 nodeMap[j].found = false;
2167 }
2168
2169 MFloat m_gridEps = 0.0000001;
2170
2173 MInt numConnections = 0, numSingularConnections = 0, nfound, tempnum = 0, tempcount;
2175 for(
MInt i = 0; i < count; ++i) {
2176 if(!nodeMap[i].found) {
2177 Point<3> a(coordinates(0, i), coordinates(1, i), coordinates(2, i));
2178 nfound = tree.locatenear(
a, m_gridEps, results, 10,
false);
2179 for(
MInt countNode2 = 0; countNode2 < nfound; ++countNode2) {
2180 for(
MInt countNode3 = countNode2 + 1; countNode3 < nfound; ++countNode3) {
2181 MBool check =
false;
2182 for(
MInt bc2 = 0; bc2 < nDim; ++bc2) {
2183 if(nodeMap[results[countNode2]].BC[bc2] == 0) break;
2184 for(
MInt bc3 = 0; bc3 < nDim; ++bc3) {
2185 if(nodeMap[results[countNode2]].BC[bc2] == nodeMap[results[countNode3]].BC[bc3]) {
2186 check = true;
2187
2189 nodeMap[results[countNode2]].BC[bc2],
2190 nodeMap[results[countNode2]].blockId,
2191 nodeMap[results[countNode2]].pos,
2192 nodeMap[results[countNode3]].blockId,
2193 nodeMap[results[countNode3]].pos)) {
2194 numConnections = numConnections + 1;
2195 }
2196 }
2197 }
2198 }
2199 if(!check) {
2200 cout << "DANGER: x|y|z=" << coordinates(0, i) << "|" << coordinates(1, i) << "|" << coordinates(2, i)
2201 << " nodeMap2: inpuBlockId=" << nodeMap[results[countNode2]].blockId << " BC =";
2202 for(
MInt d = 0; d < nDim; ++d)
2203 cout << " " << nodeMap[results[countNode2]].BC[d];
2204 cout << " nodeMap3: blockId=" << nodeMap[results[countNode3]].blockId << " BC =";
2205 for(
MInt d = 0; d < nDim; ++d)
2206 cout << " " << nodeMap[results[countNode3]].BC[d];
2207 cout << endl;
2208
2209 }
2210
2211
2212
2213
2214
2215
2216
2217
2218
2219
2220
2221
2222
2223
2224 }
2225 nodeMap[results[countNode2]].found = true;
2226 }
2227
2228 std::set<MInt> BCs;
2229 for(
MInt countNode2 = 0; countNode2 < nfound; ++countNode2) {
2230 for(
MInt d = 0; d < nDim; ++d) {
2231 if(nodeMap[results[countNode2]].BC[d] != 0) {
2232 BCs.insert(nodeMap[results[countNode2]].BC[d]);
2233 }
2234 }
2235 }
2236
2237
2238
2239 if(nfound == 3) {
2240 add = false;
2241 tempcount = 0;
2242
2243
2244 for(
MInt countNode2 = 0; countNode2 < nfound; ++countNode2) {
2246 if(nodeMap[results[countNode2]].blockId == j) {
2247 tempnum = j;
2248 break;
2249 }
2250 }
2251 for(
MInt j = 0; j < 3; ++j) {
2252
2253 if(nodeMap[results[countNode2]].pos[j] == 0
2254 || nodeMap[results[countNode2]].pos[j] ==
m_grid->getBlockNoCells(tempnum, 2 - j)) {
2255 ++tempcount;
2256 }
2257 }
2258 }
2259
2260 if(tempcount == 9 || tempcount == 6) {
2261 add = true;
2262 }
2263
2264 if(add) {
2265 for(
MInt countNode2 = 0; countNode2 < nfound; ++countNode2) {
2266 if(BCs.size() > 1) cout << "CAUTION: This singular point has more than one BC!" << endl;
2267 for(const auto& BC : BCs) {
2268
2270 nodeMap[results[countNode2]].blockId,
2271 nodeMap[results[countNode2]].pos,
2272 nodeMap[results[countNode2]].blockId,
2273 nodeMap[results[countNode2]].pos,
2274 3)) {
2275 numSingularConnections = numSingularConnections + 1;
2276 }
2277 }
2278 }
2279 }
2280 }
2281
2282
2283
2284 if(nfound == 5) {
2285 for(
MInt countNode2 = 0; countNode2 < nfound; ++countNode2) {
2286 if(BCs.size() > 1) cout << "CAUTION: This singular point has more than one BC!" << endl;
2287 for(const auto& BC : BCs) {
2288
2290 nodeMap[results[countNode2]].blockId,
2291 nodeMap[results[countNode2]].pos,
2292 nodeMap[results[countNode2]].blockId,
2293 nodeMap[results[countNode2]].pos,
2294 5)) {
2295 numSingularConnections = numSingularConnections + 1;
2296 }
2297 }
2298 }
2299 }
2300
2301 if(nfound == 6) {
2302 for(
MInt countNode2 = 0; countNode2 < nfound; ++countNode2) {
2303 if(BCs.size() > 1) cout << "CAUTION: This singular point has more than one BC!" << endl;
2304 for(const auto& BC : BCs) {
2305
2307 nodeMap[results[countNode2]].blockId,
2308 nodeMap[results[countNode2]].pos,
2309 nodeMap[results[countNode2]].blockId,
2310 nodeMap[results[countNode2]].pos,
2311 6)) {
2312 numSingularConnections = numSingularConnections + 1;
2313 }
2314 }
2315 }
2316 }
2317 }
2318 }
2319
2320
2322 MInt pcount = 0, plocation = 0, numofpoints[3];
2323
2326 for(
MInt j = 0; j < 3; j++) {
2328 }
2329 pcount += numofpoints[0] * numofpoints[1] * numofpoints[2];
2330 }
2331 }
2332
2333 if(pcount > 0) {
2335 periodicCoordinates.fill(-1.01010101);
2336 std::vector<pointType<nDim>> nodeMapP;
2337 nodeMapP.resize(pcount);
2338
2340 cout << "Loading periodic face coordinates" << endl;
2341 }
2346 stringstream number;
2348 bName += number.str();
2349
2350 ParallelIo::size_type ioOffset[3] = {0, 0, 0};
2351 ParallelIo::size_type ioSize[3] = {0, 0, 0};
2353 for(
MInt j = 0; j < 3; j++) {
2356 }
2357
2358 pio.readArray(&periodicCoordinates(0, plocation), bName, "x", 3, ioOffset, ioSize);
2359 pio.readArray(&periodicCoordinates(1, plocation), bName, "y", 3, ioOffset, ioSize);
2360 pio.readArray(&periodicCoordinates(2, plocation), bName, "z", 3, ioOffset, ioSize);
2361 } else {
2362 ioOffset[2] = 0;
2363 ioSize[2] = 0;
2364 ioOffset[1] = 0;
2365 ioSize[1] = 0;
2366 ioOffset[0] = 0;
2367 ioSize[0] = 0;
2369 pio.readArray(&empty, bName, "x", 3, ioOffset, ioSize);
2370 pio.readArray(&empty, bName, "y", 3, ioOffset, ioSize);
2371 pio.readArray(&empty, bName, "z", 3, ioOffset, ioSize);
2372 }
2373 for(int dim = 0; dim < 3; ++dim) {
2374 offset[dim] = ioOffset[dim];
2375 size[dim] = ioSize[dim];
2376 }
2377
2378
2379 for(
MInt c = offset[0]; c < offset[0] + size[0]; ++c) {
2380 for(
MInt b = offset[1];
b < offset[1] + size[1]; ++
b) {
2381 for(
MInt a = offset[2];
a < offset[2] + size[2]; ++
a) {
2383 nodeMapP[plocation].blockId =
inputWindows[i]->blockId;
2384 nodeMapP[plocation].pos[0] =
a;
2385 nodeMapP[plocation].pos[1] =
b;
2386 nodeMapP[plocation].pos[2] = c;
2387 plocation++;
2388 }
2389 }
2390 }
2391 pmemsize += size[0] * size[1] * size[2];
2392 }
2393 }
2394
2397
2401
2403
2405 cout << "Computing periodic window displacements" << endl;
2406 }
2407 m_log <<
"Computing periodic window displacements" << endl;
2413 cout.precision(10);
2414 periodicWindowCenter.fill(F0);
2415 periodicWindowNoNodes.fill(0);
2416
2417
2418 const MInt periodicOffset = 4401;
2419 for(
MInt i = 0; i < pcount; ++i) {
2420 if(nodeMapP[i].BC[0] < 4401 || nodeMapP[i].BC[0] > 4406) {
2421 continue;
2422 }
2423
2424
2425 for(
MInt dim = 0; dim < nDim; dim++) {
2426 periodicWindowCenter(nodeMapP[i].blockId, dim, nodeMapP[i].BC[0] - periodicOffset) +=
2427 periodicCoordinates(dim, i);
2428 }
2429
2430 periodicWindowNoNodes(nodeMapP[i].blockId, nodeMapP[i].BC[0] - periodicOffset) += 1;
2431 }
2432
2433
2434
2435
2436
2437 multimap<MInt, pair<vector<MFloat>,
MInt>> oddPeriodicSides;
2438 multimap<MInt, pair<vector<MFloat>,
MInt>> evenPeriodicSides;
2439
2440
2441
2442
2443
2444
2445
2446
2448 for(
MInt periodicWindowId = 0; periodicWindowId < 2 * nDim; periodicWindowId++) {
2449 if(periodicWindowNoNodes(blockId, periodicWindowId) <= 0) {
2450 continue;
2451 }
2452 vector<MFloat> centerDimensions(nDim);
2453 for(
MInt dim = 0; dim < nDim; dim++) {
2454 periodicWindowCenter(blockId, dim, periodicWindowId) /=
2455 (
MFloat)periodicWindowNoNodes(blockId, periodicWindowId);
2456 centerDimensions[dim] =
2457 periodicWindowCenter(blockId, dim,
2458 periodicWindowId);
2459 cout << "periodic Window center " << periodicWindowCenter(blockId, dim, periodicWindowId) << endl;
2460 }
2461 if(periodicWindowId % 2 == 0) {
2462
2463 evenPeriodicSides.insert(pair<
MInt, pair<vector<MFloat>,
MInt>>(
2464 blockId, pair<vector<MFloat>,
MInt>(centerDimensions, periodicWindowId)));
2465 cout << "dimensions of centerDimensions (even)" << centerDimensions[0] << " " << centerDimensions[1] << " "
2466 << centerDimensions[2] << endl;
2467 } else {
2468
2469 oddPeriodicSides.insert(pair<
MInt, pair<vector<MFloat>,
MInt>>(
2470 blockId, pair<vector<MFloat>,
MInt>(centerDimensions, periodicWindowId)));
2471 cout << "dimensions of centerDimensions (odd)" << centerDimensions[0] << " " << centerDimensions[1] << " "
2472 << centerDimensions[2] << endl;
2473 }
2474 }
2475 }
2476
2477
2478
2479
2480
2481 pair<multimap<MInt, pair<vector<MFloat>,
MInt>>::iterator, multimap<MInt, pair<vector<MFloat>,
MInt>>::iterator>
2482 evenIt;
2483 pair<multimap<MInt, pair<vector<MFloat>,
MInt>>::iterator, multimap<MInt, pair<vector<MFloat>,
MInt>>::iterator>
2484 oddIt;
2485
2486
2488 "periodic Displacements per block");
2489
2491 MInt noEvenPerodicWindowsBlockIds = evenPeriodicSides.count(blockId);
2492 MInt noOddPeriodicWindowsBlockIds = oddPeriodicSides.count(blockId);
2493 cout << "we have evenSides " << noEvenPerodicWindowsBlockIds << " and we have oddSides "
2494 << noOddPeriodicWindowsBlockIds << " in BlockId " << blockId << endl;
2495 if(noEvenPerodicWindowsBlockIds == 0) continue;
2496 if(noOddPeriodicWindowsBlockIds == 0)
2497 continue;
2498
2499 evenIt = evenPeriodicSides.equal_range(blockId);
2500 multimap<MInt, pair<vector<MFloat>,
MInt>>::iterator eit = evenIt.first;
2501 while(eit != evenIt.second) {
2502 MBool tobeErased =
false;
2504 ceil(((
MFloat)((*eit).second.second + 1) / 2.0) - 1);
2505 MInt evenSide = (*eit).second.second;
2506 MInt oddSide = evenSide + 1;
2507 oddIt = oddPeriodicSides.equal_range(blockId);
2508 multimap<MInt, pair<vector<MFloat>,
MInt>>::iterator oit = oddIt.first;
2509
2510 while(oit != oddIt.second) {
2511 if((*oit).second.second == oddSide) {
2512
2513 for(
MInt dim = 0; dim < nDim; dim++) {
2514 cout << "first element is " << (*oit).second.first[dim] << " and we subtract "
2515 << (*eit).second.first[dim] << endl;
2516 periodicDisplacementsBlockIds(blockId, evenSideDim + nDim * dim) =
2517 abs((*oit).second.first[dim] - (*eit).second.first[dim]);
2518 }
2519 oit = oddPeriodicSides.erase(oit);
2520 tobeErased = true;
2521 break;
2522 } else {
2523 oit++;
2524 }
2525 }
2526 if(tobeErased) {
2527 eit = evenPeriodicSides.erase(eit);
2528 } else {
2529 eit++;
2530 }
2531 }
2532 }
2533
2534 MFloat periodicEpsilon = 0.000001;
2535
2536
2537 if(evenPeriodicSides.size() != 0
2538 && oddPeriodicSides.size() != 0) {
2540 MInt noEvenPerodicWindowsBlockId = evenPeriodicSides.count(blockId);
2541 if(noEvenPerodicWindowsBlockId == 0) continue;
2542 evenIt = evenPeriodicSides.equal_range(blockId);
2543 multimap<MInt, pair<vector<MFloat>,
MInt>>::iterator eit = evenIt.first;
2544 while(eit != evenIt.second) {
2545 MBool tobeErased =
false;
2547 ceil(((
MFloat)((*eit).second.second + 1) / 2.0) - 1);
2548
2549 multimap<MInt, pair<vector<MFloat>,
MInt>>::iterator oit = oddPeriodicSides.begin();
2550 while(oit != oddPeriodicSides.end()) {
2551
2552 MInt equalConnectionDims = 0;
2553 for(
MInt d = 0; d < nDim; d++) {
2554 if(abs((*oit).second.first[d] - (*eit).second.first[d]) < periodicEpsilon) equalConnectionDims++;
2555 }
2556 if(equalConnectionDims == 2) {
2557 for(
MInt dim = 0; dim < nDim; dim++) {
2558 periodicDisplacementsBlockIds(blockId, evenSideDim + nDim * dim) =
2559 abs((*oit).second.first[dim] - (*eit).second.first[dim]);
2560
2561 periodicDisplacementsBlockIds((*oit).first, evenSideDim + nDim * dim) =
2562 periodicDisplacementsBlockIds(blockId, evenSideDim + nDim * dim);
2563 }
2564 tobeErased = true;
2565 oit = oddPeriodicSides.erase(oit);
2566 } else {
2567 oit++;
2568 }
2569 }
2570 if(tobeErased) {
2571 eit = evenPeriodicSides.erase(eit);
2572 } else {
2573 eit++;
2574 }
2575 }
2576 }
2577 }
2578
2579
2580 for(
MInt periodicWindowId = 0; periodicWindowId < 2 * nDim; periodicWindowId++) {
2581 if(periodicWindowId % 2 == 1) {
2582 for(
MInt dim = 0; dim < nDim; dim++) {
2583 const MInt displacementId = (
MFloat)(periodicWindowId + 1) / 2.0 - 1;
2584 periodicDisplacements[dim * nDim + displacementId] =
2585 periodicWindowCenter(
m_blockId, dim, periodicWindowId)
2586 - periodicWindowCenter(
m_blockId, dim, periodicWindowId - 1);
2587 m_log <<
m_blockId <<
" periodicWindowCenter for window: " << periodicWindowId + periodicOffset
2588 << " dim: " << dim << " displacementId: " << displacementId
2589 << " displacement: " << periodicDisplacements[dim * nDim + displacementId] << endl;
2590 }
2591 }
2592 }
2593
2594 cout << "old approach" << endl;
2595 for(
MInt periodicWindowId = 0; periodicWindowId < 2 * nDim; periodicWindowId++) {
2596 if(periodicWindowId % 2 == 1) {
2597 for(
MInt dim = 0; dim < nDim; dim++) {
2598 const MInt displacementId = (
MFloat)(periodicWindowId + 1) / 2.0 - 1;
2599 periodicDisplacements[dim * nDim + displacementId] =
2600 periodicWindowCenter(
m_blockId, dim, periodicWindowId)
2601 - periodicWindowCenter(
m_blockId, dim, periodicWindowId - 1);
2603 << "periodicWindowCenter for window: " << periodicWindowId + periodicOffset << " dim: " << dim
2604 << " displacementId: " << displacementId
2605 << " displacement: " << periodicDisplacements[dim * nDim + displacementId] << endl;
2606 }
2607 }
2608 }
2609 cout << "new approach " << endl;
2610 for(
MInt periodicWindowId = 0; periodicWindowId < 2 * nDim; periodicWindowId++) {
2611 if(periodicWindowId % 2 == 1) {
2612 const MInt displacementId = (
MFloat)(periodicWindowId + 1) / 2.0 - 1;
2613 for(
MInt dim = 0; dim < nDim; dim++) {
2615 << "periodicWindowCenter for window: " << periodicWindowId + periodicOffset << " dim: " << dim
2616 << " displacementId: " << displacementId
2617 <<
" displacement: " << periodicDisplacementsBlockIds(
m_blockId, dim * nDim + displacementId) << endl;
2618 }
2619 }
2620 }
2624 cout <<
"BlockId " <<
b <<
": " << endl;
2625 cout << "\t";
2626 for(
MInt i = 0; i < nDim * nDim; i++) {
2627 cout << periodicDisplacementsBlockIds(
b, i) <<
" ";
2628 }
2629 cout << endl;
2630 }
2631 }
2635 for(
MInt i = 0; i < nDim * nDim; i++) {
2636 periodicDisplacements[i] = periodicDisplacementsBlockIds(
b, i);
2637 }
2638 }
2639 }
2640
2641
2644
2646
2647 for(
MInt i = 0; i < pcount; ++i) {
2648 tmppoint[0] = periodicCoordinates(0, i);
2649 tmppoint[1] = periodicCoordinates(1, i);
2650 tmppoint[2] = periodicCoordinates(2, i);
2652 Point<3> aaa(tmppoint[0], tmppoint[1], tmppoint[2]);
2653 nfound = tree.locatenear(aaa, m_gridEps, results, 10, false);
2654 for(
MInt countNode2 = 0; countNode2 < nfound; ++countNode2) {
2655
2657 nodeMapP[i].blockId,
2658 nodeMapP[i].pos,
2659 nodeMap[results[countNode2]].blockId,
2660 nodeMap[results[countNode2]].pos)) {
2661 numConnections = numConnections + 1;
2662 }
2663 }
2664
2665 if(nfound == 5) {
2666
2667 if(
addConnection(nodeMapP[i].BC[0], nodeMapP[i].blockId, nodeMapP[i].pos, nodeMapP[i].blockId,
2668 nodeMapP[i].pos, 5)) {
2669 numSingularConnections = numSingularConnections + 1;
2670 }
2671 }
2672
2673 if(nfound == 6) {
2674
2675 if(
addConnection(nodeMapP[i].BC[0], nodeMapP[i].blockId, nodeMapP[i].pos, nodeMapP[i].blockId,
2676 nodeMapP[i].pos, 6)) {
2677 numSingularConnections = numSingularConnections + 1;
2678 }
2679 }
2680 }
2681 }
2682
2684 cout << "Assemble multiblock connections" << endl;
2685 }
2691
2692
2699
2700
2701
2704 break;
2705 }
2706
2707 if((*it)->BC >= 4400 && (*it)->BC < 4410) continue;
2708 MInt noSameWindows = 0;
2709 (*it)->BCsingular[noSameWindows++] = (*it)->BC;
2710 (*it)->BC = 0;
2713 break;
2714 }
2715 const unique_ptr<StructuredWindowMap<nDim>>& map1 = (*it);
2716 const unique_ptr<StructuredWindowMap<nDim>>& map2 = (*it2);
2718 if(test) {
2719 if((*it)->BCsingular[0] == (*it2)->BC || (*it)->Nstar != (*it2)->Nstar) {
2720 mTerm(1,
"There is no hope anymore!");
2721 }
2722
2723 if(noSameWindows >= (*it)->Nstar) {
2725 }
2726
2727 (*it)->BCsingular[noSameWindows++] = (*it2)->BC;
2729 } else {
2730 ++it2;
2731 }
2732 }
2733 ++it;
2734 }
2735
2736
2737
2738
2739
2740
2741
2742
2743
2744
2745
2746
2747
2748
2749
2750
2751
2752
2753
2754
2755
2756
2757
2758
2759
2760
2761
2762
2763
2764
2765
2766
2767
2768
2769
2770
2771
2772
2774
2776
2777 unique_ptr<StructuredWindowMap<nDim>> windowMap1 = make_unique<StructuredWindowMap<nDim>>();
2780 windowMap1->Nstar = -1;
2781 windowMap1->SingularId = -1;
2782 windowMap1->dc1 =
window2d[i]->dc1;
2783 windowMap1->dc2 =
window2d[i]->dc2;
2784
2785
2786 if(windowMap1->dc1 * windowMap1->dc2 > 0) {
2787 MInt tempdim = abs(windowMap1->dc2) - 1;
2788 windowMap1->step2[tempdim] = -1;
2789 }
2790
2791
2793
2794
2796 unique_ptr<StructuredWindowMap<nDim>> windowMap = make_unique<StructuredWindowMap<nDim>>();
2799 windowMap->Nstar = -1;
2800 windowMap->SingularId = -1;
2803
2804 if(windowMap->dc1 * windowMap->dc2 > 0) {
2805 MInt tempdim = abs(windowMap->dc2) - 1;
2806 windowMap->step2[tempdim] = -1;
2807 }
2808
2812 }
2813 }
2814
2815
2817
2819
2820 unique_ptr<StructuredWindowMap<nDim>> windowMap1 = make_unique<StructuredWindowMap<nDim>>();
2823 windowMap1->Nstar = -1;
2824 windowMap1->SingularId = -1;
2825 windowMap1->dc1 =
window1d[i]->dc1;
2826 windowMap1->dc2 =
window1d[i]->dc2;
2827
2828
2829 for(
MInt j = 0; j < 3; ++j) {
2830 if(windowMap1->start1[j] == windowMap1->end1[j]) {
2831 if(windowMap1->start1[j] == 0 && windowMap1->start2[windowMap1->order[j]] == 0) {
2832 windowMap1->step2[windowMap1->order[j]] = -1;
2833 }
2834 if(windowMap1->start1[j] > 0 && windowMap1->start2[windowMap1->order[j]] > 0) {
2835 windowMap1->step2[windowMap1->order[j]] = -1;
2836 }
2837 }
2838 }
2839
2840
2841
2844 if(test) {
2847 }
2848 }
2849
2850
2852
2854 unique_ptr<StructuredWindowMap<nDim>> windowMap = make_unique<StructuredWindowMap<nDim>>();
2857 windowMap->Nstar = -1;
2858 windowMap->SingularId = -1;
2861
2862
2863 for(
MInt j = 0; j < 3; ++j) {
2864 if(windowMap->start1[j] == windowMap->end1[j]) {
2865 if(windowMap->start1[j] == 0 && windowMap->start2[windowMap->order[j]] == 0) {
2866 windowMap->step2[windowMap->order[j]] = -1;
2867 }
2868 if(windowMap->start1[j] > 0 && windowMap->start2[windowMap->order[j]] > 0) {
2869 windowMap->step2[windowMap->order[j]] = -1;
2870 }
2871 }
2872 }
2873
2876
2877
2878
2881 if(test) {
2884 }
2885 }
2886
2888 }
2889 }
2890
2891
2893
2895
2896
2897 unique_ptr<StructuredWindowMap<nDim>> windowMap1 = make_unique<StructuredWindowMap<nDim>>();
2900 windowMap1->Nstar = -1;
2901 windowMap1->SingularId = -1;
2902 windowMap1->dc1 =
window0d[i]->dc1;
2903 windowMap1->dc2 =
window0d[i]->dc2;
2904
2905 for(
MInt j = 0; j < 3; ++j) {
2906 if(windowMap1->start1[j] == windowMap1->end1[j]) {
2907 if(windowMap1->start1[j] == 0 && windowMap1->start2[windowMap1->order[j]] == 0) {
2908 windowMap1->step2[windowMap1->order[j]] = -1;
2909 }
2910 if(windowMap1->start1[j] > 0 && windowMap1->start2[windowMap1->order[j]] > 0) {
2911 windowMap1->step2[windowMap1->order[j]] = -1;
2912 }
2913 }
2914 }
2915
2916
2917
2920 if(test) {
2923 }
2924 }
2925
2926
2928
2930 unique_ptr<StructuredWindowMap<nDim>> windowMap = make_unique<StructuredWindowMap<nDim>>();
2933 windowMap->Nstar = -1;
2934 windowMap->SingularId = -1;
2937
2938
2939 for(
MInt j = 0; j < 3; ++j) {
2940 if(windowMap->start1[j] == windowMap->end1[j]) {
2941 if(windowMap->start1[j] == 0 && windowMap->start2[windowMap->order[j]] == 0) {
2942 windowMap->step2[windowMap->order[j]] = -1;
2943 }
2944 if(windowMap->start1[j] > 0 && windowMap->start2[windowMap->order[j]] > 0) {
2945 windowMap->step2[windowMap->order[j]] = -1;
2946 }
2947 }
2948 }
2949
2952
2953
2954
2957 if(test) {
2960 }
2961 }
2963 }
2964 }
2965
2966
2968
2969#ifndef NDEBUG
2971 cout << "======= GLOBALSTRUCTUREDBNDRYCNDMAPS ======" << endl;
2974 }
2975
2976 cout << " ======= SINGULARWINDOW =============" << endl;
2979 }
2980
2981 cout << " ======== window0d ==============" << endl;
2984 }
2985
2986 cout << " ======== window1d ==============" << endl;
2989 }
2990
2991 cout << " ======== window2d ==============" << endl;
2994 }
2995 }
2996#endif
2997
2998 if(!hasConnectionInfo) {
2999 cout << "Writing connection info" << endl;
3001 } else {
3003 cout << "##########################################################################" << endl;
3004 cout << "WARNING: OLD CONNECTION INFO EXISTS, DELETE BEFORE NEW ONE CAN BE WRITTEN!" << endl;
3005 cout << "##########################################################################" << endl;
3006 }
3007 }
3008
3010 cout << "Connection identification and window creation finished! NUM_SINGULAR=" << numSingularConnections << endl;
3011 }
3012 }
3013}
void mAlloc(T *&a, const MLong N, const MString &objectName, MString function)
allocates memory for one-dimensional array 'a' of size N
static MBool propertyExists(const MString &name, MInt solver=m_noSolvers)
This function checks if a property exists in general.
void writeConnectionWindowInformation3D(MFloat *periodicDisplacements)
void deleteDuplicateWindows(std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > &, std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > &)
void readConnectionWindowInformation3D(MFloat *periodicDisplacments)
void multiBlockAssembling()
void singularityAssembling()
void deleteDuplicateCommMaps()
MBool addConnection(MInt connectiontype, MInt b1, MInt *p1, MInt b2, MInt *p2)
void periodicPointsChange(MFloat *pt, MInt type, MFloat *periodicDisplacements)
int MPI_Type_commit(MPI_Datatype *datatype, const MString &name)
same as MPI_Type_commit
int MPI_Type_contiguous(int count, MPI_Datatype old_type, MPI_Datatype *new_type_p, const MString &name)
same as MPI_Type_contiguous
1.9.5