1293 {
1295 const MUint noCutCells = cutCellData.size();
1296
1297
1298
1299 static constexpr MInt cornersSOLVERtoMC[8] = {1, 2, 0, 3, 5, 6, 4, 7};
1300 static constexpr MInt edgesMCtoSOLVER[12] = {0, 2, 1, 3, 4, 6, 5, 7, 10, 8, 9, 11};
1301 static constexpr MInt facesMCtoSOLVER[6] = {0, 2, 1, 3, 4, 5};
1302 static constexpr MFloat signStencil[8][3] = {{-F1, -F1, -F1}, {F1, -F1, -F1}, {-F1, F1, -F1}, {F1, F1, -F1},
1303 {-F1, -F1, F1}, {F1, -F1, F1}, {-F1, F1, F1}, {F1, F1, F1}};
1304 static constexpr MInt edgeCornerCode[12][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}, {4, 6}, {5, 7},
1305 {4, 5}, {6, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
1306 static constexpr MInt faceEdgeCode[6][4] = {{0, 10, 4, 8}, {1, 9, 5, 11}, {2, 8, 6, 9},
1307 {3, 11, 7, 10}, {2, 1, 3, 0}, {6, 4, 7, 5}};
1308 static constexpr MInt edgeFaceCode[12][2] = {{0, 4}, {1, 4}, {2, 4}, {3, 4}, {0, 5}, {1, 5},
1309 {2, 5}, {3, 5}, {0, 2}, {1, 2}, {0, 3}, {1, 3}};
1310 static constexpr MInt faceCornerCode[6][4] = {{0, 2, 6, 4}, {3, 1, 5, 7}, {1, 0, 4, 5},
1311 {2, 3, 7, 6}, {0, 1, 3, 2}, {5, 4, 6, 7}};
1312 const MInt maxNoSets = 6;
1314
1317 for(
MInt i = 0; i < 15; i++) {
1318 presentCases(s, i) = 0;
1319 }
1320 }
1321
1322 const MInt maxNoFaces = 100;
1323 const MInt maxNoVertices = 300;
1324 const MInt maxNoEdges = 300;
1325
1326 std::vector<polyVertex> vertices[maxNoSets];
1327 std::vector<polyFace> faces[maxNoSets];
1328 std::vector<polyEdge> edges[maxNoSets];
1329 std::vector<polyCutCell> cutCells;
1330
1331 stack<MInt> faceStack;
1332
1333 std::vector<polyMultiFace> multiFaces;
1334 MIntScratchSpace edgeCutCellPointer(maxNoEdges, AT_,
"edgeCutCellPointer");
1335 std::list<MInt> pureBodyEdges;
1337 MIntScratchSpace multiFaceConnection(maxNoFaces, AT_,
"multiFaceConnection");
1338
1342
1343 std::vector<CsgPolygon> csg_polygons;
1344 std::vector<CsgVertex> csg_vertices;
1345 std::vector<Csg> csg;
1346 std::vector<CsgPolygon> result;
1347 std::vector<polyVertex> vertices_result;
1349 std::list<std::pair<MInt, MInt>> openEdges;
1350
1352 std::vector<MInt> faceVertices;
1355 "normalDotProduct");
1358 MBoolScratchSpace isPartOfBodySurface(maxNoVertices, AT_,
"isPartOfBodySurface");
1360
1361
1362
1363
1367
1368 NEW_SUB_TIMER_STATIC(tCutFace_0, "cutFaceNew_0", tCutGroup);
1369 NEW_SUB_TIMER_STATIC(tCutFace_1a, "cutFaceNew_1a", tCutGroup);
1370 NEW_SUB_TIMER_STATIC(tCutFace_1, "cutFaceNew_1", tCutGroup);
1371 NEW_SUB_TIMER_STATIC(tCutFace_2, "cutFaceNew_2", tCutGroup);
1372 NEW_SUB_TIMER_STATIC(tCutFace_3, "cutFaceNew_3", tCutGroup);
1373 NEW_SUB_TIMER_STATIC(tCutFace_3b, "cutFaceNew_3b", tCutGroup);
1374 NEW_SUB_TIMER_STATIC(tCutFace_4, "cutFaceNew_4", tCutGroup);
1375 NEW_SUB_TIMER_STATIC(tCutFace_5a, "cutFaceNew_5a", tCutGroup);
1376 NEW_SUB_TIMER_STATIC(tCutFace_5b, "cutFaceNew_5b", tCutGroup);
1377 NEW_SUB_TIMER_STATIC(tCutFace_6, "cutFaceNew_6", tCutGroup);
1378 NEW_SUB_TIMER_STATIC(tCutFace_7, "cutFaceNew_7", tCutGroup);
1379
1380#ifdef CutCell_DEBUG
1381 const MInt debugDomainId = 5;
1382 const MInt debugTimeStep = 0;
1383 const MInt debugCellId = 21866;
1384#endif
1385
1386
1387
1388
1389
1390
1391
1394
1395
1396 for(
MUint cutc = 0; cutc < noCutCells; cutc++) {
1397
1398 RECORD_TIMER_START(tCutFace_0);
1399
1400
1401
1402
1403 const MInt cellId = cutCellData[cutc].cellId;
1404
1405
1407 const MFloat cellHalfLength = F1B2 * cellLength0;
1408
1409 for(
MBool& externalFace : cutCellData[cutc].externalFaces) {
1410 externalFace = false;
1411 }
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1423 vertices[s].clear();
1424 faces[s].clear();
1425 edges[s].clear();
1426 cutSetPointers[s] = -1;
1427 noCutPointsFromSet[s] = 0;
1428 }
1429 cutCells.clear();
1430 const MBool isGapCell = cutCellData[cutc].isGapCell;
1431
1435 startSet = 1;
1437 }
1438
1439
1440
1442 MBool isCompletelyOutside =
false;
1443 MBool hasAmbiguousFaces =
false;
1444 for(
MInt s = startSet; s < endSet; s++) {
1445
1446
1447
1448 unsigned char outcode_MC = 0;
1450
1451 MBool currentOutcode = (!cutCellData[cutc].cornerIsInsideGeometry[s][c]);
1452 if(currentOutcode) {
1453 outcode_MC = outcode_MC | (1 << cornersSOLVERtoMC[c]);
1454 }
1455 }
1456 outcode_MC_set[s] = outcode_MC;
1457 if(outcode_MC == 0) {
1458 isCompletelyOutside = true;
1459 }
1460 }
1461
1462 for(
MInt cutPoint = 0; cutPoint < cutCellData[cutc].noCutPoints; cutPoint++) {
1465
1467 }
1468
1469 noCutPointsFromSet[set]++;
1470 }
1471 for(
MInt set = startSet; set < endSet; set++) {
1472 if(noCutPointsFromSet[set] && !isCompletelyOutside) {
1473 cutSetPointers[set] = noCutSets++;
1474 }
1475 }
1476 MInt maxCutsPerFace = 0;
1478 for(
MInt cutPoint = 0; cutPoint < cutCellData[cutc].noCutPoints; cutPoint++) {
1479 MInt edge = cutCellData[cutc].cutEdges[cutPoint];
1480 cutsPerFace[edgeFaceCode[edge][0]]++;
1481 cutsPerFace[edgeFaceCode[edge][1]]++;
1482 }
1483 for(
MInt k = 0; k < CC::noFaces; k++) {
1484 maxCutsPerFace =
mMax(maxCutsPerFace, cutsPerFace[k]);
1485 }
1486
1487 RECORD_TIMER_STOP(tCutFace_0);
1488
1489
1490 for(
MInt set = startSet; set < endSet; set++) {
1491 MInt setIndex = cutSetPointers[set];
1492 if(setIndex < 0) {
1493 continue;
1494 }
1495 MInt outcode_MC = outcode_MC_set[set];
1496
1497 MInt currentCase = cases[outcode_MC][0];
1498 if(!caseStatesLs[currentCase][0]) {
1499 cerr << "Error: Case not implemented: " << currentCase << endl;
1500 mTerm(1, AT_,
"Case not implemented, see error file for deatails.");
1501 }
1502
1503 if(!caseStatesLs[currentCase][1]) {
1504 hasAmbiguousFaces = true;
1505 }
1506 }
1507
1508
1509
1510 MBool simpleMarchingCubesSucceeds =
false;
1511
1512 if(!hasAmbiguousFaces && noCutSets == 1 && maxCutsPerFace == 2 && !isCompletelyOutside
1514
1515 RECORD_TIMER_START(tCutFace_1a);
1516
1518 RECORD_TIMER_STOP(tCutFace_1a);
1519 }
1520
1521 if(simpleMarchingCubesSucceeds) {
1522 continue;
1523 }
1524
1527 noTriangles[i] = 0;
1528 }
1529
1530
1531
1532 for(
MInt set = startSet; set < endSet; set++) {
1533 MInt setIndex = cutSetPointers[set];
1534 if(setIndex < 0) {
1535 continue;
1536 }
1537
1538#ifdef CutCell_DEBUG
1539 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
1540 cerr << "Cutting-Cell with " << set << " Index " << setIndex << " bodyId "
1541 << cutCellData[cutc].associatedBodyIds[set] << endl;
1542 }
1543#endif
1544
1545 RECORD_TIMER_START(tCutFace_1);
1546
1547
1548 MInt bodyId = cutCellData[cutc].associatedBodyIds[set];
1549 if(bodyId < 0) continue;
1550
1551
1552
1553
1554 MInt cutPoints[12] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1};
1555 MInt cutPointToVertexMap[12] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1};
1556 MInt noCutPoints = 0;
1557
1558 for(
MInt cutPoint = 0; cutPoint < cutCellData[cutc].noCutPoints; cutPoint++) {
1560
1562 continue;
1563 }
1564 }
1565 noCutPoints++;
1566
1567 cutPoints[cutCellData[cutc].cutEdges[cutPoint]] = cutPoint;
1568
1569 vertices[setIndex].emplace_back(cutCellData[cutc].cutPoints[cutPoint], cutPoint, 1);
1570 cutPointToVertexMap[cutPoint] = vertices[setIndex].size() - 1;
1571 }
1572 MInt cornerToVertexMap[8] = {-1, -1, -1, -1, -1, -1, -1, -1};
1574
1575 if(!cutCellData[cutc].cornerIsInsideGeometry[set][c]) {
1576 std::array<MFloat, nDim> tmp_coords{};
1578 tmp_coords[i] =
a_coordinate(cellId, i) + signStencil[c][i] * cellHalfLength;
1579 }
1580 vertices[setIndex].emplace_back(tmp_coords, c, 0);
1581 cornerToVertexMap[c] = vertices[setIndex].size() - 1;
1582 }
1583 }
1584
1585
1586
1587
1588
1589
1590
1591 MInt outcode_MC = outcode_MC_set[set];
1592
1593
1594 MInt currentCase = cases[outcode_MC][0];
1595 presentCases(set, currentCase)++;
1596 MInt currentSubCase = cases[outcode_MC][1];
1598 if(!caseStatesLs[currentCase][0]) {
1599 cerr << "Error:Case not implemented: " << currentCase << endl;
1600 mTerm(1, AT_,
"Case not implemented, see error file for deatails.");
1601 }
1602
1603#ifdef CutCell_DEBUG
1605#endif
1606
1607#ifdef CutCell_DEBUG
1608 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
1609 cerr << "cases: " << outcode_MC << " case " << currentCase << " subcase " << currentSubCase << " ambiguous "
1610 << caseStatesLs[currentCase][1] << endl;
1611 }
1612#endif
1613
1614
1615 if(!caseStatesLs[currentCase][1]) {
1616
1617 MInt faceMax = noAmbiguousFaces[currentCase];
1618 const MInt* testPointer;
1619 switch(currentCase) {
1620 case 3:
1621 testPointer = test3[currentSubCase];
1622 break;
1623 case 6:
1624 testPointer = test6[currentSubCase];
1625 break;
1626 case 7:
1627 testPointer = test7[currentSubCase];
1628 break;
1629 case 10:
1630 testPointer = test10[currentSubCase];
1631 break;
1632 case 12:
1633 testPointer = test12[currentSubCase];
1634 break;
1635 case 13:
1636 testPointer = test13[currentSubCase];
1637 break;
1638 default:
1639 mTerm(1, AT_,
"this should not happen!");
1640 break;
1641 }
1642 for(
MInt i = 0; i < faceMax; i++) {
1643 MInt face = testPointer[i];
1644 MBool revert =
false;
1645 if(face < F0) {
1646 if(face == -6) {
1647 face = 0;
1648 } else {
1649 face = -face;
1650 }
1651 revert = true;
1652 }
1653 face = facesMCtoSOLVER[face];
1654
1655 MBool testResult = cutCellData[cutc].faceCentroidIsInsideGeometry[set][face];
1656
1657 if(revert) testResult = !testResult;
1658
1659 if(testResult) {
1660 subConfig +=
IPOW2(i);
1661 }
1662 }
1663 }
1664
1665
1666
1667 const MInt* tilingPointer =
nullptr;
1668 MBool addVertex =
false;
1669 noTriangles[set] = noTriangles_simpleCases[currentCase];
1670 switch(currentCase) {
1671 case 0:
1672 break;
1673
1674 case 1:
1675 tilingPointer = tiling1Ls[currentSubCase];
1676 break;
1677 case 2:
1678 tilingPointer = tiling2Ls[currentSubCase];
1679 break;
1680 case 4:
1681 tilingPointer = tiling4Ls[currentSubCase];
1682 break;
1683 case 5:
1684 tilingPointer = tiling5Ls[currentSubCase];
1685 break;
1686 case 8:
1687 tilingPointer = tiling8Ls[currentSubCase];
1688 break;
1689 case 9:
1690 tilingPointer = tiling9Ls[currentSubCase];
1691 break;
1692 case 11:
1693 tilingPointer = tiling11Ls[currentSubCase];
1694 break;
1695 case 14:
1696 tilingPointer = tiling14Ls[currentSubCase];
1697 break;
1698
1699 case 3:
1700 switch(subConfig) {
1701 case 0:
1702 tilingPointer = tiling3_1Ls[currentSubCase];
1703 noTriangles[set] = 2;
1704 break;
1705 case 1:
1706 tilingPointer = tiling3_2Ls[currentSubCase];
1707 noTriangles[set] = 4;
1708 break;
1709 default:
1710 mTerm(1, AT_,
"This should not happen - we have case 3 with subConfig > 1... exiting");
1711 break;
1712 }
1713 break;
1714
1715 case 6:
1716 switch(subConfig) {
1717 case 0:
1718 tilingPointer = tiling6_1_1Ls[currentSubCase];
1719 noTriangles[set] = 3;
1720 break;
1721 case 1:
1722 tilingPointer = tiling6_2Ls[currentSubCase];
1723 noTriangles[set] = 5;
1724 break;
1725 default:
1726 mTerm(1, AT_,
"This should not happen - we have case 6 with subConfig > 1... exiting");
1727 break;
1728 }
1729 break;
1730
1731 case 7:
1732 switch(subConfig) {
1733 case 0:
1734 tilingPointer = tiling7_1Ls[currentSubCase];
1735 noTriangles[set] = 3;
1736 break;
1737 case 1:
1738 tilingPointer = tiling7_2Ls[currentSubCase][0];
1739 noTriangles[set] = 5;
1740 break;
1741 case 2:
1742 tilingPointer = tiling7_2Ls[currentSubCase][1];
1743 noTriangles[set] = 5;
1744 break;
1745 case 3:
1746 tilingPointer = tiling7_3Ls[currentSubCase][0];
1747 noTriangles[set] = 9;
1748 addVertex = true;
1749 break;
1750 case 4:
1751 tilingPointer = tiling7_2Ls[currentSubCase][2];
1752 noTriangles[set] = 5;
1753 break;
1754 case 5:
1755 tilingPointer = tiling7_3Ls[currentSubCase][1];
1756 noTriangles[set] = 9;
1757 addVertex = true;
1758 break;
1759 case 6:
1760 tilingPointer = tiling7_3Ls[currentSubCase][2];
1761 noTriangles[set] = 9;
1762 addVertex = true;
1763 break;
1764 case 7:
1765 tilingPointer = tiling7_4_1Ls[currentSubCase];
1766 noTriangles[set] = 5;
1767 break;
1768 default:
1769 mTerm(1, AT_,
"This should not happen - we have case 7 with subConfig > 7... exiting");
1770 break;
1771 }
1772 break;
1773
1774 case 10:
1775 switch(subConfig) {
1776 case 0:
1777 tilingPointer = tiling10_1_1Ls[currentSubCase];
1778 noTriangles[set] = 4;
1779 break;
1780 case 1:
1781 tilingPointer = tiling10_2Ls[currentSubCase];
1782 noTriangles[set] = 8;
1783 addVertex = true;
1784 break;
1785 case 2:
1786 tilingPointer = tiling10_2_Ls[currentSubCase];
1787 noTriangles[set] = 8;
1788 addVertex = true;
1789 break;
1790 case 3:
1791 tilingPointer = tiling10_1_1_Ls[currentSubCase];
1792 noTriangles[set] = 4;
1793 break;
1794 default:
1795 mTerm(1, AT_,
"This should not happen - we have case 10 with subConfig > 3... exiting");
1796 break;
1797 }
1798 break;
1799
1800 case 12:
1801 switch(subConfig) {
1802 case 0:
1803 tilingPointer = tiling12_1_1Ls[currentSubCase];
1804 noTriangles[set] = 4;
1805 break;
1806 case 1:
1807 tilingPointer = tiling12_2Ls[currentSubCase];
1808 noTriangles[set] = 8;
1809 addVertex = true;
1810 break;
1811 case 2:
1812 tilingPointer = tiling12_2_Ls[currentSubCase];
1813 noTriangles[set] = 8;
1814 addVertex = true;
1815 break;
1816 case 3:
1817 tilingPointer = tiling12_1_1_Ls[currentSubCase];
1818 noTriangles[set] = 4;
1819 break;
1820 default:
1821 mTerm(1, AT_,
"This should not happen - we have case 12 with subConfig > 3... exiting");
1822 break;
1823 }
1824 break;
1825
1826 case 13: {
1827 MInt config_identificator = 0;
1828 subConfig = subconfig13[subConfig];
1829 config_identificator = 0;
1830 switch(subConfig) {
1831 case 0:
1832 tilingPointer = tiling13_1Ls[currentSubCase];
1833 noTriangles[set] = 4;
1834 break;
1835
1836 case 1:
1837 case 2:
1838 case 3:
1839 case 4:
1840 case 5:
1841 case 6:
1842 config_identificator = subConfig - 1;
1843 tilingPointer = tiling13_2Ls[currentSubCase][config_identificator];
1844 noTriangles[set] = 6;
1845 break;
1846
1847 case 7:
1848 case 8:
1849 case 9:
1850 case 10:
1851 case 11:
1852 case 12:
1853 case 13:
1854 case 14:
1855 case 15:
1856 case 16:
1857 case 17:
1858 case 18:
1859 config_identificator = subConfig - 7;
1860 tilingPointer = tiling13_3Ls[currentSubCase][config_identificator];
1861 noTriangles[set] = 10;
1862 addVertex = true;
1863 break;
1864
1865 case 19:
1866 case 20:
1867 case 21:
1868 case 22:
1869 config_identificator = subConfig - 19;
1870 tilingPointer = tiling13_4Ls[currentSubCase][config_identificator];
1871 noTriangles[set] = 12;
1872 addVertex = true;
1873 break;
1874
1875 case 23:
1876 case 24:
1877 case 25:
1878 case 26:
1879 config_identificator = subConfig - 23;
1880 tilingPointer = tiling13_5_1Ls[currentSubCase][config_identificator];
1881 noTriangles[set] = 6;
1882 break;
1883
1884 case 27:
1885 case 28:
1886 case 29:
1887 case 30:
1888 case 31:
1889 case 32:
1890 case 33:
1891 case 34:
1892 case 35:
1893 case 36:
1894 case 37:
1895 case 38:
1896 config_identificator = subConfig - 27;
1897 tilingPointer = tiling13_3_Ls[currentSubCase][config_identificator];
1898 noTriangles[set] = 10;
1899 addVertex = true;
1900 break;
1901
1902 case 39:
1903 case 40:
1904 case 41:
1905 case 42:
1906 case 43:
1907 case 44:
1908 config_identificator = subConfig - 39;
1909 tilingPointer = tiling13_2_Ls[currentSubCase][config_identificator];
1910 noTriangles[set] = 6;
1911 break;
1912
1913 case 45:
1914 tilingPointer = tiling13_1_Ls[currentSubCase];
1915 noTriangles[set] = 4;
1916 break;
1917
1918 default:
1919 mTerm(1, AT_,
"impossible MC case 13 subConfig, how could this happen? exiting...");
1920 break;
1921 }
1922 } break;
1923
1924 default:
1925 mTerm(1, AT_,
"invalid MC case, how could this happen? exiting...");
1926 break;
1927 }
1928 RECORD_TIMER_STOP(tCutFace_1);
1929
1930
1931 RECORD_TIMER_START(tCutFace_2);
1932
1933 MInt additionalVertexId = -1;
1934 if(addVertex) {
1935 additionalVertexId = vertices[setIndex].size();
1936 vertices[setIndex].emplace_back(-1, 2);
1937 addPoint(&vertices[setIndex], cutPointToVertexMap, noCutPoints,
1938 vertices[setIndex][additionalVertexId].coordinates);
1939 }
1940
1941 ASSERT(noCutPoints >=
nDim,
"");
1942 ASSERT(noCutPoints == caseCutPoints[currentCase], to_string(
grid().tree().globalId(cellId)));
1943
1944 for(
MInt t = 0; t < noTriangles[set]; t++) {
1945
1946
1947 MInt currentFace = faces[setIndex].size();
1948 faces[setIndex].emplace_back(-1, 1, bodyId);
1949
1951 for(
MInt pt = 0; pt < 3; pt++) {
1952 MInt cutEdge = tilingPointer[t * 3 + pt];
1953 if(cutEdge == 12) {
1954 p[pt] = additionalVertexId;
1955 } else {
1956 p[pt] = cutPointToVertexMap[cutPoints[edgesMCtoSOLVER[cutEdge]]];
1957 }
1958 ASSERT(p[pt] >= 0, to_string(cutEdge) + " " + to_string(isGapCell) + " " + to_string(cellId) + " "
1959 + to_string(
grid().tree().globalId(cellId)) +
" " + to_string(edgesMCtoSOLVER[cutEdge])
1960 + " " + to_string(cutPoints[edgesMCtoSOLVER[cutEdge]]) + " "
1961 + to_string(cutPointToVertexMap[cutPoints[edgesMCtoSOLVER[cutEdge]]]) + " "
1962 + to_string(cutCellData[cutc].noCutPoints));
1963 }
1964 for(
MInt te = 0; te < 3; te++) {
1966 ASSERT(point > -1, "");
1967 vertices[setIndex][point].surfaceIdentificators.insert(bodyId +
m_noDirs);
1968 MInt nextPoint =
p[(te + 1) % 3];
1969 MBool edgeExists =
false;
1970 for(
MInt e = 0; (unsigned)e < vertices[setIndex][p[te]].edges.size(); e++) {
1971 MInt edge = vertices[setIndex][point].edges[e];
1972 if(edges[setIndex][edge].vertices[1] == nextPoint) {
1973
1974 faces[setIndex][currentFace].edges.emplace_back(edge, 1);
1975 const MInt vertexId = edges[setIndex][edge].vertices[0];
1976 faces[setIndex][currentFace].vertices.push_back(vertexId);
1977 edges[setIndex][edge].face[0] = currentFace;
1978 edgeExists = true;
1979 } else if(edges[setIndex][edge].vertices[0] == nextPoint) {
1980
1981 faces[setIndex][currentFace].edges.emplace_back(edge, -1);
1982 const MInt vertexId = edges[setIndex][edge].vertices[1];
1983 faces[setIndex][currentFace].vertices.push_back(vertexId);
1984
1985 edges[setIndex][edge].face[1] = currentFace;
1986 edgeExists = true;
1987 }
1988 }
1989 if(!edgeExists) {
1990
1991 MInt newEdge = edges[setIndex].size();
1993 if(vertices[setIndex][point].
pointType == 2 || vertices[setIndex][nextPoint].
pointType == 2) {
1994 edgeType = 3;
1995 }
1997 if(edgeType == 2) {
1998 MInt p0 = vertices[setIndex][point].pointId;
1999 MInt p1 = vertices[setIndex][nextPoint].pointId;
2000
2001 MInt edge0 = cutCellData[cutc].cutEdges[p0];
2002
2003 MInt edge1 = cutCellData[cutc].cutEdges[p1];
2004 if(edgeFaceCode[edge0][0] == edgeFaceCode[edge1][0] || edgeFaceCode[edge0][0] == edgeFaceCode[edge1][1])
2005 edgeId = edgeFaceCode[edge0][0];
2006 else if(edgeFaceCode[edge0][1] == edgeFaceCode[edge1][0]
2007 || edgeFaceCode[edge0][1] == edgeFaceCode[edge1][1])
2008 edgeId = edgeFaceCode[edge0][1];
2009 }
2010 edges[setIndex].emplace_back(point, nextPoint, edgeId, edgeType);
2011 vertices[setIndex][point].edges.push_back(newEdge);
2012 vertices[setIndex][nextPoint].edges.push_back(newEdge);
2013 faces[setIndex][currentFace].edges.emplace_back(newEdge, 1);
2014 const MInt vertexId = edges[setIndex][newEdge].vertices[0];
2015 faces[setIndex][currentFace].vertices.push_back(vertexId);
2016
2017 edges[setIndex][newEdge].face[0] = currentFace;
2018 }
2019 }
2020
2021 computeNormal(vertices[setIndex][p[0]].coordinates, vertices[setIndex][p[1]].coordinates,
2022 vertices[setIndex][p[2]].coordinates, faces[setIndex][currentFace].normal,
2023 faces[setIndex][currentFace].w);
2024 }
2025
2026 for(auto& f : faces[setIndex]) {
2027 ASSERT(f.vertices.size() == 3, "");
2028 }
2029
2030#ifdef CutCell_DEBUG
2031 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2032 cerr << "Body surfaces: " << endl;
2033 cerr << "Set: " << set << " noTriangles " << noTriangles[set] << " noFaces " << faces[setIndex].size()
2034 << " noEdges " << edges[setIndex].size() << " noVertices " << vertices[setIndex].size() << endl;
2035 }
2036#endif
2037
2038
2039 RECORD_TIMER_STOP(tCutFace_2);
2040
2041 RECORD_TIMER_START(tCutFace_3);
2042
2043
2044
2045
2046
2047
2049
2050 MBool p0Fluid = !cutCellData[cutc].cornerIsInsideGeometry[set][edgeCornerCode[e][0]];
2051
2052 MBool p1Fluid = !cutCellData[cutc].cornerIsInsideGeometry[set][edgeCornerCode[e][1]];
2053 if(p0Fluid && p1Fluid) {
2054 MInt v0 = cornerToVertexMap[edgeCornerCode[e][0]];
2055 MInt v1 = cornerToVertexMap[edgeCornerCode[e][1]];
2056 edges[setIndex].emplace_back(v0, v1, e, 0);
2057 vertices[setIndex][v0].edges.push_back(edges[setIndex].size() - 1);
2058 vertices[setIndex][v1].edges.push_back(edges[setIndex].size() - 1);
2059 } else if(p0Fluid) {
2060 MInt v0 = cornerToVertexMap[edgeCornerCode[e][0]];
2061 MInt v1 = cutPointToVertexMap[cutPoints[e]];
2062 edges[setIndex].emplace_back(v0, v1, e, 1);
2063 vertices[setIndex][v0].edges.push_back(edges[setIndex].size() - 1);
2064 vertices[setIndex][v1].edges.push_back(edges[setIndex].size() - 1);
2065 } else if(p1Fluid) {
2066 MInt v0 = cornerToVertexMap[edgeCornerCode[e][1]];
2067 MInt v1 = cutPointToVertexMap[cutPoints[e]];
2068 edges[setIndex].emplace_back(v0, v1, e, 1);
2069 vertices[setIndex][v0].edges.push_back(edges[setIndex].size() - 1);
2070 vertices[setIndex][v1].edges.push_back(edges[setIndex].size() - 1);
2071 }
2072 }
2073
2074
2075
2076 for(
MInt cartFace = 0; cartFace <
m_noDirs; cartFace++) {
2077
2078 MBool edgeTouched[4] = {
false,
false,
false,
false};
2079 for(
MInt e = 0; e < 4; e++) {
2080 if(edgeTouched[e]) continue;
2081 MInt cartEdge = faceEdgeCode[cartFace][e];
2082 MInt startVertex = cornerToVertexMap[faceCornerCode[cartFace][e]];
2083
2084 if(startVertex == -1) continue;
2085
2086
2087 MInt currentFace = faces[setIndex].size();
2088 faces[setIndex].emplace_back(cartFace, 0, -1);
2090 faces[setIndex][currentFace].normal[i] = F0;
2091 }
2092 MInt normalDir = cartFace / 2;
2094 if(cartFace % 2 == 0) {
2095 normalSign = F1;
2096 }
2097 faces[setIndex][currentFace].w = -vertices[setIndex][startVertex].coordinates[normalDir] * normalSign;
2098 faces[setIndex][currentFace].normal[normalDir] = normalSign;
2099
2100
2102 MInt currentVertex = startVertex;
2103 do {
2104 vertices[setIndex][currentVertex].surfaceIdentificators.insert(cartFace);
2105 edgeTouched[currentE] = true;
2106
2107
2108 MBool edgeFound =
false;
2109 for(
MInt ve = 0; (unsigned)ve < vertices[setIndex][currentVertex].edges.size(); ve++) {
2110 MInt edge = vertices[setIndex][currentVertex].edges[ve];
2111 if(edges[setIndex][edge].edgeType > 1)
2112 continue;
2113 if(edges[setIndex][edge].edgeId != cartEdge) continue;
2114
2115 edgeFound = true;
2116 MInt edgeDirection = 1;
2117 if(edges[setIndex][edge].vertices[1] == currentVertex) {
2118 edgeDirection = -1;
2119 }
2120 faces[setIndex][currentFace].edges.emplace_back(edge, edgeDirection);
2121 MInt vertexId = edges[setIndex][edge].vertices[0];
2122 if(edgeDirection == -1) vertexId = edges[setIndex][edge].vertices[1];
2123 faces[setIndex][currentFace].vertices.push_back(vertexId);
2124
2125 if(edges[setIndex][edge].face[0] > -1)
2126 edges[setIndex][edge].face[1] = currentFace;
2127 else
2128 edges[setIndex][edge].face[0] = currentFace;
2129 if(edgeDirection == 1)
2130 currentVertex = edges[setIndex][edge].vertices[1];
2131 else
2132 currentVertex = edges[setIndex][edge].vertices[0];
2133 break;
2134 }
2135 ASSERT(edgeFound, "");
2136
2137
2138 if(vertices[setIndex][currentVertex].
pointType == 0) {
2139 currentE = (currentE + 1) % 4;
2140 cartEdge = faceEdgeCode[cartFace][currentE];
2141 }
else if(vertices[setIndex][currentVertex].
pointType == 1) {
2142 edgeFound = false;
2143 for(
MInt ve = 0; (unsigned)ve < vertices[setIndex][currentVertex].edges.size(); ve++) {
2144 MInt edge = vertices[setIndex][currentVertex].edges[ve];
2145 if(edges[setIndex][edge].edgeType != 2)
2146 continue;
2147 if(edges[setIndex][edge].edgeId != cartFace) continue;
2148
2149 edgeFound = true;
2150 MInt edgeDirection = 1;
2151 if(edges[setIndex][edge].vertices[1] == currentVertex) {
2152 edgeDirection = -1;
2153 }
2154 faces[setIndex][currentFace].edges.emplace_back(edge, edgeDirection);
2155 MInt vertexId = edges[setIndex][edge].vertices[0];
2156 if(edgeDirection == -1) vertexId = edges[setIndex][edge].vertices[1];
2157 faces[setIndex][currentFace].vertices.push_back(vertexId);
2158
2159 edges[setIndex][edge].face[1] = currentFace;
2160 if(edgeDirection == 1)
2161 currentVertex = edges[setIndex][edge].vertices[1];
2162 else
2163 currentVertex = edges[setIndex][edge].vertices[0];
2164
2165 cartEdge = cutCellData[cutc].cutEdges[vertices[setIndex][currentVertex].pointId];
2166 while(faceEdgeCode[cartFace][currentE] != cartEdge) {
2167 currentE = (currentE + 1) % 4;
2168 }
2169 break;
2170 }
2171 ASSERT(edgeFound, "");
2172 }
2173 } while(currentVertex != startVertex);
2174
2175 }
2176 }
2177
2178#ifdef CutCell_DEBUG
2179 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2180 cerr << "Including cartesian faces: " << endl;
2181 cerr << "Set: " << set << " SetIndex:" << setIndex << " faces : " << faces[setIndex].size()
2182 << " edges : " << edges[setIndex].size() << " certices : " << vertices[setIndex].size() << endl;
2183 for(
MUint i = 0; i < vertices[setIndex].size(); i++) {
2184 cerr << " vertices " << setprecision(19) << vertices[setIndex][i].coordinates[0] << " "
2185 << vertices[setIndex][i].coordinates[1] << " " << vertices[setIndex][i].coordinates[2] << endl;
2186 }
2187 }
2188#endif
2189
2190
2191 RECORD_TIMER_STOP(tCutFace_3);
2192 }
2193
2194 RECORD_TIMER_START(tCutFace_3b);
2195
2196 MInt noIndividualCuts = 0;
2197 for(
MInt set = startSet; set < endSet; set++) {
2198 if(cutSetPointers[set] > -1) {
2199 noIndividualCuts++;
2200 }
2201 }
2202
2203#ifdef CutCell_DEBUG
2204 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2205 cerr << "noIndividualCuts: " << noIndividualCuts << endl;
2206 }
2207#endif
2208
2209 if(noIndividualCuts == 0) {
2210 cutCellData[cutc].volume = F0;
2211 RECORD_TIMER_STOP(tCutFace_3b);
2212 continue;
2213 }
2214 ASSERT(noIndividualCuts > 0, "");
2215 ASSERT(noCutSets > 0, " ");
2216 ASSERT(noCutSets == noIndividualCuts, "");
2217
2218
2219
2220 MBool error =
false;
2221 const MInt startSetIndex = 0;
2222 MInt referenceSet = startSet;
2223 for(
MInt set = startSet; set < endSet; set++) {
2224 if(cutSetPointers[set] > -1) {
2225 referenceSet = set;
2226 break;
2227 }
2228 }
2229
2230
2231
2232
2233
2234
2235
2236
2237 for(
MInt set = startSet; set < endSet; set++) {
2238 MInt setIndex = cutSetPointers[set];
2239 if(setIndex < 0) continue;
2240
2241 for(
MInt faceId = 0; faceId < (signed)faces[setIndex].size(); faceId++) {
2242 ASSERT(faces[setIndex][faceId].edges.size() == faces[setIndex][faceId].vertices.size(), "");
2243 for(
MInt edgeId = 0; edgeId < (signed)faces[setIndex][faceId].edges.size(); edgeId++) {
2244 const MInt edge = faces[setIndex][faceId].edges[edgeId].first;
2245 const MInt newVertexId = faces[setIndex][faceId].vertices[edgeId];
2246 const MInt direction = faces[setIndex][faceId].edges[edgeId].second;
2247 MInt vertexId = edges[setIndex][edge].vertices[0];
2248 if(direction == -1) vertexId = edges[setIndex][edge].vertices[1];
2249 ASSERT(newVertexId == vertexId, "");
2250 }
2251 }
2252 }
2253
2254
2255
2256
2257
2258 MInt noInitialFaces = 0;
2259
2260
2261
2262
2263
2264 if(noIndividualCuts > 1) {
2265
2268
2269
2270
2271 for(
MInt set = referenceSet + 1; set < endSet; set++) {
2272 MInt setIndex = cutSetPointers[set];
2273 if(setIndex < 0) continue;
2274
2275 csg.clear();
2276 csg_polygons.clear();
2277
2278
2279
2280 for(
MInt face = faces[startSetIndex].size() - 1; face > -1; face--) {
2281 csg_vertices.clear();
2282 polyFace* faceP = &faces[startSetIndex][face];
2283
2284 for(
MInt e = 0; (unsigned)e < faceP->vertices.size(); e++) {
2285 const MInt vertexId = faceP->vertices[e];
2286 csg_vertices.emplace_back(CsgVector(vertices[startSetIndex][vertexId].coordinates), vertexId,
2287 startSetIndex);
2288 }
2289 if(faceP->faceType == 0) {
2290 CsgVector normal(faceP->normal);
2291 CsgPlane plane(normal, -faceP->w);
2292 csg_polygons.emplace_back(csg_vertices, startSetIndex, faceP->faceId, faceP->faceType, faceP->bodyId,
2293 plane);
2294 } else {
2295 csg_polygons.emplace_back(csg_vertices, startSetIndex, faceP->faceId, faceP->faceType, faceP->bodyId);
2296 }
2297 }
2298
2299 csg.emplace_back(csg_polygons);
2300 csg_polygons.clear();
2301
2302
2303 for(
MInt face = faces[setIndex].size() - 1; face > -1; face--) {
2304 csg_vertices.clear();
2305 polyFace* faceP = &faces[setIndex][face];
2306
2307 for(
MInt e = 0; (unsigned)e < faceP->vertices.size(); e++) {
2308 const MInt vertexId = faceP->vertices[e];
2309 csg_vertices.emplace_back(CsgVector(vertices[setIndex][vertexId].coordinates), vertexId, setIndex);
2310 }
2311 if(faceP->faceType == 0) {
2312 CsgVector normal(faceP->normal);
2313 CsgPlane plane(normal, -faceP->w);
2314 csg_polygons.emplace_back(csg_vertices, setIndex, faceP->faceId, faceP->faceType, faceP->bodyId, plane);
2315 } else {
2316 csg_polygons.emplace_back(csg_vertices, setIndex, faceP->faceId, faceP->faceType, faceP->bodyId);
2317 }
2318 }
2319 csg.emplace_back(csg_polygons);
2320
2321#ifdef CutCell_DEBUG
2322 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2323 cerr << "Combine: reference " << referenceSet << " with " << set << "Index " << startSetIndex << "face-size "
2324 << faces[startSetIndex].size() << "Index " << setIndex << "face-size " << faces[setIndex].size() << endl;
2325 }
2326#endif
2327
2328
2329
2330
2331 result.clear();
2332 result = csg[0].intersect(csg[1]);
2333 vertices_result.clear();
2334
2335
2336#ifdef CutCell_DEBUG
2337 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2338 cerr << "result-size: " << result.size() << "vertices_result " << vertices_result.size() << endl;
2339
2340 for(
MInt poligon = 0; (unsigned)poligon < result.size(); poligon++) {
2341 cerr << "resultId " << poligon << " body " << result[poligon].bodyId << " type " << result[poligon].faceType
2342 << " setIndex " << setIndex << endl;
2343 for(
MInt vertex = 0; (unsigned)vertex < result[poligon].vertices.size(); vertex++) {
2344 cerr << "vertex " << vertex << "coordinate: " << setprecision(21)
2345 << result[poligon].vertices[vertex].pos.xx[0] << " " << result[poligon].vertices[vertex].pos.xx[1]
2346 << " " << result[poligon].vertices[vertex].pos.xx[2] << " "
2347 << result[poligon].vertices[vertex].vertexId << " " << result[poligon].vertices[vertex].setIndex
2348 << endl;
2349 }
2350 cerr << "faceId " << result[poligon].faceId << endl;
2351 }
2352 }
2353#endif
2354
2355
2356
2357
2358
2359
2360
2361
2363 for(
MInt j = 0; j < maxNoVertices; j++) {
2364 vertices_renamed(i, j) = -1;
2365 }
2366 }
2367
2368
2369 MInt noVertices = 0;
2370
2371
2372
2373
2374
2375 for(
MInt p = 0; (unsigned)p < result.size();
p++) {
2376 ASSERT(result[p].vertices.size() <= (unsigned)maxNoVertices, "");
2377 ASSERT((signed)result[p].vertices.size() >= 3, "");
2378
2379
2380 for(
MInt v = 0; (unsigned)v < result[p].vertices.size(); v++) {
2381 MInt vertexId = result[
p].vertices[v].vertexId;
2382 MInt vertexSetIndex = result[
p].vertices[v].setIndex;
2383 if(vertexSetIndex < 0) continue;
2384
2385 if(vertexId < 0) continue;
2386
2387
2388
2389 if(vertices_renamed(vertexSetIndex, vertexId) == -1) {
2390
2391 MBool vertexFound =
false;
2392 MInt vertexIndex = -1;
2393
2394
2395 for(
MInt i = 0; (unsigned)i < vertices_result.size(); i++) {
2397 coord_diff += pow(result[p].vertices[v].pos.xx[0] - vertices_result[i].coordinates[0], 2);
2398 coord_diff += pow(result[p].vertices[v].pos.xx[1] - vertices_result[i].coordinates[1], 2);
2399 coord_diff += pow(result[p].vertices[v].pos.xx[2] - vertices_result[i].coordinates[2], 2);
2400
2402
2403
2404 if(coord_diff < difBound1) {
2405#ifdef CutCell_DEBUG
2406 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2407 cerr <<
"Small diff " <<
p <<
" " << vertexId <<
" " << vertices_result[i].coordinates[0] <<
" "
2408 << vertices_result[i].coordinates[1] << " " << vertices_result[i].coordinates[2] << " eps "
2409 << difBound1 <<
" dif " << coord_diff <<
" " << result[
p].vertices[v].pos.xx[0] <<
" "
2410 << result[
p].vertices[v].pos.xx[1] <<
" " << result[
p].vertices[v].pos.xx[2] <<
" " << endl;
2411 }
2412#endif
2413
2414 vertexFound = true;
2415 vertexIndex = i;
2416 break;
2417 }
2418 }
2419 if(vertexFound) {
2420 vertices_renamed(vertexSetIndex, vertexId) = vertexIndex;
2421 } else {
2422
2423 vertices_renamed(vertexSetIndex, vertexId) = noVertices;
2424 vertices_result.emplace_back(vertices[vertexSetIndex][vertexId].pointId,
2425 vertices[vertexSetIndex][vertexId].
pointType);
2426
2427 ASSERT(vertices[vertexSetIndex][vertexId].
pointType == 0
2428 || vertices[vertexSetIndex][vertexId].
pointType == 1
2429 || vertices[vertexSetIndex][vertexId].
pointType == 2,
2430 "");
2431 vertices_result[noVertices].coordinates[0] = vertices[vertexSetIndex][vertexId].coordinates[0];
2432 vertices_result[noVertices].coordinates[1] = vertices[vertexSetIndex][vertexId].coordinates[1];
2433 vertices_result[noVertices].coordinates[2] = vertices[vertexSetIndex][vertexId].coordinates[2];
2434 noVertices++;
2435 }
2436 }
2437
2438 result[
p].vertices[v].vertexId = vertices_renamed(vertexSetIndex, vertexId);
2439 result[
p].vertices[v].setIndex = -1;
2440 }
2441 }
2442
2443
2444
2445
2446 for(
MInt p = 0; (unsigned)p < result.size();
p++) {
2447 for(
MInt v = 0; (unsigned)v < result[p].vertices.size(); v++) {
2448 MInt vertexId = result[
p].vertices[v].vertexId;
2449 if(vertexId == -1) {
2450
2451 MBool vertexFound =
false;
2452 MInt vertexIndex = -1;
2453 for(
MInt i = 0; (unsigned)i < vertices_result.size(); i++) {
2455 coord_diff += pow(result[p].vertices[v].pos.xx[0] - vertices_result[i].coordinates[0], 2);
2456 coord_diff += pow(result[p].vertices[v].pos.xx[1] - vertices_result[i].coordinates[1], 2);
2457 coord_diff += pow(result[p].vertices[v].pos.xx[2] - vertices_result[i].coordinates[2], 2);
2458
2460
2461
2462 MFloat difBound = difBound1;
2464 difBound = difBound2;
2465 }
2466
2467 if(coord_diff < difBound) {
2468#ifdef CutCell_DEBUG
2469 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2470 cerr <<
"2 Small diff " <<
p <<
" " << vertexId <<
" " << vertices_result[i].coordinates[0] <<
" "
2471 << vertices_result[i].coordinates[1] << " " << vertices_result[i].coordinates[2] << " eps "
2472 << coord_diff <<
" " << result[
p].vertices[v].pos.xx[0] <<
" "
2473 << result[
p].vertices[v].pos.xx[1] <<
" " << result[
p].vertices[v].pos.xx[2] <<
" " << endl;
2474 }
2475#endif
2476
2477 vertexFound = true;
2478 vertexIndex = i;
2479 break;
2480 }
2481 }
2482 if(vertexFound) {
2483 result[
p].vertices[v].vertexId = vertexIndex;
2484 result[
p].vertices[v].setIndex = -1;
2485 } else {
2486
2487
2488 vertices_result.emplace_back(-1, 3);
2489 vertices_result[noVertices].coordinates[0] = result[
p].vertices[v].pos.xx[0];
2490 vertices_result[noVertices].coordinates[1] = result[
p].vertices[v].pos.xx[1];
2491 vertices_result[noVertices].coordinates[2] = result[
p].vertices[v].pos.xx[2];
2492 result[
p].vertices[v].vertexId = noVertices;
2493 result[
p].vertices[v].setIndex = -1;
2494 noVertices++;
2495 }
2496 }
2497 }
2498 }
2499
2500
2501 vertices[startSetIndex].swap(vertices_result);
2502 vertices_result.clear();
2503 edges[startSetIndex].clear();
2504 faces[startSetIndex].clear();
2505
2506#ifdef CutCell_DEBUG
2507 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2508 cerr << "VertexId-check: " << endl;
2509
2510 for(
MInt poligon = 0; (unsigned)poligon < result.size(); poligon++) {
2511 cerr << "resultId " << poligon << " body " << result[poligon].bodyId << " type " << result[poligon].faceType
2512 << " setIndex " << setIndex << endl;
2513 for(
MInt vertex = 0; (unsigned)vertex < result[poligon].vertices.size(); vertex++) {
2514 cerr << "vertex " << vertex << "coordinate: " << setprecision(21)
2515 << result[poligon].vertices[vertex].pos.xx[0] << " " << result[poligon].vertices[vertex].pos.xx[1]
2516 << " " << result[poligon].vertices[vertex].pos.xx[2] << " "
2517 << result[poligon].vertices[vertex].vertexId << " " << result[poligon].vertices[vertex].setIndex
2518 << endl;
2519 }
2520 cerr << "faceId " << result[poligon].faceId << endl;
2521 }
2522
2523 for(
MInt poligon = 0; (unsigned)poligon < result.size(); poligon++) {
2524 for(
MInt vertex = 0; (unsigned)vertex < result[poligon].vertices.size(); vertex++) {
2525 MInt vertexId = result[poligon].vertices[vertex].vertexId;
2527 if(fabs(result[poligon].vertices[vertex].pos.xx[i] - vertices[startSetIndex][vertexId].coordinates[i])
2529 cerr <<
" VertexId-missmatch " <<
cellId <<
" " << vertexId <<
" " << vertex <<
" " << i
2530 << setprecision(21) << result[poligon].vertices[vertex].pos.xx[i] << " "
2531 << vertices[startSetIndex][vertexId].coordinates[i] << endl;
2532 }
2533 }
2534 }
2535 }
2536 }
2537
2538#endif
2539
2540
2541
2542
2543
2545 noInitialFaces = 0;
2547
2548
2549 for(
MInt p = 0; (unsigned)p < result.size();
p++) {
2551 MBool polygonValid =
true;
2552 MInt validVertices = 0;
2553 faceMapping[
p] = -1;
2554
2555 for(
MInt v = 0; (unsigned)v < result[p].vertices.size(); v++) {
2556 const MInt j = (v + 1) % result[p].vertices.size();
2557 const MInt vertexId = result[
p].vertices[v].vertexId;
2558 const MInt vertexIdNext = result[
p].vertices[j].vertexId;
2559 if(vertexId == vertexIdNext) {
2560 vertexValid[v] = false;
2561 } else {
2562 vertexValid[v] = true;
2563 validVertices++;
2564 }
2565
2566 if(vertexValid[v]) {
2567 MInt surfaceIndicator = -1;
2568 if(result[p].faceType == 0) {
2569 surfaceIndicator = result[
p].faceId;
2570 } else {
2571 surfaceIndicator = result[
p].bodyId +
m_noDirs;
2572 }
2573 vertices[startSetIndex][vertexId].surfaceIdentificators.insert(surfaceIndicator);
2574 }
2575 }
2576 if(validVertices < 3) {
2577#ifdef CutCell_DEBUG
2578
2579
2580
2581
2582
2583
2584
2585
2586
2587
2588
2589#endif
2590 polygonValid = false;
2591 }
2592
2593 if(!polygonValid) continue;
2594
2595 faces[startSetIndex].emplace_back(result[p].faceId, result[p].faceType, result[p].bodyId);
2596 faces[startSetIndex][noFaces].normal[0] = result[
p].plane.normal.xx[0];
2597 faces[startSetIndex][noFaces].normal[1] = result[
p].plane.normal.xx[1];
2598 faces[startSetIndex][noFaces].normal[2] = result[
p].plane.normal.xx[2];
2599 ASSERT(
2600 !(std::isnan(result[p].plane.normal.xx[0] + result[p].plane.normal.xx[1] + result[p].plane.normal.xx[2])),
2601 "");
2602 faces[startSetIndex][noFaces].tmpSetIndex = result[
p].setIndex;
2603 faces[startSetIndex][noFaces].w = -result[
p].plane.w;
2604 faceMapping[
p] = noFaces;
2605 faces[startSetIndex][noFaces].edges.clear();
2606 faces[startSetIndex][noFaces].vertices.clear();
2607
2608
2609 for(
MInt v = 0; (unsigned)v < result[p].vertices.size(); v++) {
2610 const MInt vertexId = result[
p].vertices[v].vertexId;
2611 MBool vertexUsed =
false;
2612 for(
MInt vertice : faces[startSetIndex][noFaces].vertices) {
2613 if(vertexId == vertice) {
2614 vertexUsed = true;
2615 break;
2616 }
2617 }
2618
2619 if(!vertexUsed) {
2620 faces[startSetIndex][noFaces].vertices.push_back(vertexId);
2621 vertices[startSetIndex][vertexId].faceIds.push_back(noFaces);
2622 }
2623 }
2624 ASSERT((signed)faces[startSetIndex][noFaces].vertices.size() >= 3, "");
2625 noFaces++;
2626 noInitialFaces++;
2627 }
2628
2629#ifdef CutCell_DEBUG
2630 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2631 cerr << "Current-Face-Count " << noFaces << endl;
2632 }
2633#endif
2634
2635
2636
2638#if defined CutCell_DEBUG || !defined NDEBUG
2639 MBool onlySingleEdges =
true;
2640#endif
2641 for(
MInt p = 0; (unsigned)p < result.size();
p++) {
2642 MInt faceId = faceMapping[
p];
2643 if(faceId < 0) continue;
2644
2645
2646 for(
MInt v = 0; (unsigned)v < result[p].vertices.size(); v++) {
2647 MInt j = (v + 1) % result[p].vertices.size();
2648 MInt vertexId = result[
p].vertices[v].vertexId;
2649 MInt vertexIdNext = result[
p].vertices[j].vertexId;
2650 if(vertexId == vertexIdNext) {
2651 continue;
2652 }
2653
2655 MBool edgeFound =
false;
2657 for(
MInt e = 0; (unsigned)e < edges[startSetIndex].size(); e++) {
2658 MInt v0 = edges[startSetIndex][e].vertices[0];
2659 MInt v1 = edges[startSetIndex][e].vertices[1];
2660 if(vertexId == v0 && vertexIdNext == v1) {
2661 edgeFound = true;
2662 edgeId = e;
2663 direction = 1;
2664 break;
2665 }
2666 if(vertexId == v1 && vertexIdNext == v0) {
2667 edgeFound = true;
2668 edgeId = e;
2669 direction = -1;
2670 break;
2671 }
2672 }
2673
2674 if(!edgeFound) {
2675 edges[startSetIndex].emplace_back(vertexId, vertexIdNext, -1, -1);
2676 edges[startSetIndex][noEdges].face[0] = faceId;
2677 edgeId = noEdges;
2678 direction = 1;
2679 vertices[startSetIndex][vertexId].edges.push_back(noEdges);
2680 vertices[startSetIndex][vertexIdNext].edges.push_back(noEdges);
2681 noEdges++;
2682
2683#ifdef CutCell_DEBUG
2685 cerr << "Adding Edge " << noEdges << " " << vertexId << " " << vertexIdNext << " " << faceId << endl;
2686 }
2687#endif
2688
2689
2690 } else if(edges[startSetIndex][edgeId].face[1] < 0) {
2691 edges[startSetIndex][edgeId].face[1] = faceId;
2692
2693#ifdef CutCell_DEBUG
2695 cerr << "Adding 2nd face to egde " << edgeId << " " << vertexId << " " << vertexIdNext << " " << faceId
2696 << endl;
2697 cerr << "edge vertices are: " << edges[startSetIndex][edgeId].vertices[0] << " "
2698 << edges[startSetIndex][edgeId].vertices[1] << endl;
2699 }
2700#endif
2701 } else {
2702
2703
2704#if defined CutCell_DEBUG || !defined NDEBUG
2705 onlySingleEdges = false;
2706#endif
2707 edges[startSetIndex].emplace_back(vertexId, vertexIdNext, -1, -1);
2708 edges[startSetIndex][noEdges].face[0] = faceId;
2709 edgeId = noEdges;
2710 direction = 1;
2711 vertices[startSetIndex][vertexId].edges.push_back(noEdges);
2712 vertices[startSetIndex][vertexIdNext].edges.push_back(noEdges);
2713 noEdges++;
2714#ifdef CutCell_DEBUG
2715 cerr <<
cellId <<
"Has a duplicate edge!" << endl;
2717 cerr << "Adding Duplicate Edge " << noEdges << " " << vertexId << " " << vertexIdNext << " " << faceId
2718 << endl;
2719 }
2720#endif
2721 }
2722 faces[startSetIndex][faceId].edges.emplace_back(edgeId, direction);
2723 }
2724 }
2725
2726
2727 for(
MInt faceId = 0; faceId < static_cast<MInt>(faces[startSetIndex].size()); faceId++) {
2728 ASSERT(faces[startSetIndex][faceId].edges.size() == faces[startSetIndex][faceId].vertices.size(), "");
2729 for(
MInt edgeId = 0; edgeId < static_cast<MInt>(faces[startSetIndex][faceId].edges.size()); edgeId++) {
2730 const MInt edge = faces[startSetIndex][faceId].edges[edgeId].first;
2731 const MInt newVertexId = faces[startSetIndex][faceId].vertices[edgeId];
2732 const MInt direction = faces[startSetIndex][faceId].edges[edgeId].second;
2733 MInt vertexId = edges[startSetIndex][edge].vertices[0];
2734 if(direction == -1) {
2735 vertexId = edges[startSetIndex][edge].vertices[1];
2736 }
2737 ASSERT(newVertexId == vertexId, "");
2738 }
2739 }
2740
2741
2742
2743
2744
2745
2746
2747 openEdges.clear();
2748
2749 for(
MInt e = 0; e < static_cast<MInt>(edges[startSetIndex].size()); e++) {
2750 ASSERT(edges[startSetIndex][e].face[0] > -1, "");
2751 if(edges[startSetIndex][e].face[1] == -1) {
2752 MInt face = edges[startSetIndex][e].face[0];
2754 for(
MInt fe = 0; fe < (signed)faces[startSetIndex][face].edges.size(); fe++) {
2755 if(faces[startSetIndex][face].edges[fe].first == e) {
2756 dir = faces[startSetIndex][face].edges[fe].second;
2757 }
2758 }
2759 ASSERT(dir, "");
2760 if(dir == 1) {
2761 openEdges.emplace_back(e, -1);
2762 } else {
2763 openEdges.emplace_back(e, 1);
2764 }
2765 }
2766 }
2767
2768 vector<std::pair<MInt, MInt>> openEdgeList;
2770 openEdgeId.fill(-1);
2771 for(auto& openEdge : openEdges) {
2772 openEdgeId(openEdge.first) =
static_cast<MInt>(openEdgeList.size());
2773 openEdgeList.push_back(openEdge);
2774 }
2775
2776#ifdef CutCell_DEBUG
2778 cerr << openEdges.size() << " initially open edges! " << endl;
2779 cerr << "Count: " << faces[startSetIndex].size() << " faces " << edges[startSetIndex].size() << " edges "
2780 << vertices[startSetIndex].size() << " vertices " << endl;
2781
2782 for(auto it2 = openEdges.begin(); it2 != openEdges.end(); it2++) {
2783 const MInt edgdeId = (*it2).first;
2784 cerr << " openEdge with Id " << (*it2).first << " v1 " << edges[startSetIndex][edgdeId].vertices[0] << " "
2785 << edges[startSetIndex][edgdeId].vertices[1] << endl;
2786 }
2787 }
2788#endif
2789
2790
2791
2792
2793
2794
2795
2796
2797
2798
2799#if defined CutCell_DEBUG || !defined NDEBUG
2800 for(
MInt i = 0; i < (signed)openEdgeList.size(); i++) {
2801 const MInt e = openEdgeList[i].first;
2802 const MInt dir = openEdgeList[i].second;
2803 const MInt id0 = (dir == 1) ? 0 : 1;
2804 const MInt id1 = (dir == 1) ? 1 : 0;
2805 const MInt v0 = edges[startSetIndex][e].vertices[id0];
2806 const MInt v1 = edges[startSetIndex][e].vertices[id1];
2807
2808 ASSERT(edges[startSetIndex][e].face[0] > -1, "");
2809
2810
2811 for(
MInt j = 0; j < (signed)vertices[startSetIndex][v0].edges.size(); j++) {
2812 const MInt e2 = vertices[startSetIndex][v0].edges[j];
2813
2814 if(e2 == e) continue;
2815
2816 const MInt face = edges[startSetIndex][e2].face[0];
2818 for(
MInt fe = 0; fe < (signed)faces[startSetIndex][face].edges.size(); fe++) {
2819 if(faces[startSetIndex][face].edges[fe].first == e2) {
2820 dir2 = faces[startSetIndex][face].edges[fe].second;
2821 }
2822 }
2823 const MInt id00 = (dir2 == 1) ? 0 : 1;
2824 const MInt id11 = (dir2 == 1) ? 1 : 0;
2825 const MInt v00 = edges[startSetIndex][e2].vertices[id00];
2826 const MInt v11 = edges[startSetIndex][e2].vertices[id11];
2827
2828 ASSERT(edges[startSetIndex][e2].face[0] > -1, "");
2829 if(onlySingleEdges) {
2830 ASSERT(v0 != v00 || v1 != v11, "");
2831 ASSERT(v0 != v11 || v1 != v00, "");
2832 }
2833 }
2834 }
2835#endif
2836
2837
2838
2839
2840 MBool somethingChanged = !openEdgeList.empty();
2841 while(somethingChanged) {
2842 somethingChanged = false;
2843
2844 MBool alreadyResolved =
false;
2845
2846 for(
MUint k1 = 0; k1 < openEdgeList.size(); k1++) {
2847
2848 const MInt e1 = openEdgeList[k1].first;
2849 if(openEdgeId(e1) < 0) continue;
2850 const MInt dir1 = openEdgeList[k1].second;
2851 const MInt rdir1 = (dir1 == 1) ? -1 : 1;
2852 const MInt id0 = (dir1 == 1) ? 0 : 1;
2853 const MInt id1 = (dir1 == 1) ? 1 : 0;
2854 const MInt v0 = edges[startSetIndex][e1].vertices[id0];
2855 const MInt v1 = edges[startSetIndex][e1].vertices[id1];
2856 const MInt face = edges[startSetIndex][e1].face[0];
2857
2858
2860 for(fe = 0; (unsigned)fe < faces[startSetIndex][face].edges.size(); fe++) {
2861 if(faces[startSetIndex][face].edges[fe].first == e1) {
2862 break;
2863 }
2864 }
2865 ASSERT(fe < (signed)faces[startSetIndex][face].edges.size(), "");
2866
2867
2868
2870 for(fv = 0; (unsigned)fv < faces[startSetIndex][face].vertices.size(); fv++) {
2871 if(faces[startSetIndex][face].vertices[fv] == v0 || faces[startSetIndex][face].vertices[fv] == v1) {
2872 if(fv + 1 < (signed)faces[startSetIndex][face].vertices.size()
2873 && (faces[startSetIndex][face].vertices[fv + 1] == v0
2874 || faces[startSetIndex][face].vertices[fv + 1] == v1)) {
2875 break;
2876 }
2877 }
2878 }
2879 fv = fv + 1;
2880
2884 a[i] = vertices[startSetIndex][v1].coordinates[i] - vertices[startSetIndex][v0].coordinates[i];
2885 a_abs +=
a[i] *
a[i];
2886 }
2887 a_abs = sqrt(a_abs);
2888
2889
2890 for(
MInt k2 = 0; (unsigned)k2 < vertices[startSetIndex][v0].edges.size(); k2++) {
2891
2892 const MInt e2 = vertices[startSetIndex][v0].edges[k2];
2893 if(e2 == e1) continue;
2894 if(openEdgeId(e2) < 0) continue;
2895 ASSERT(openEdgeList[openEdgeId(e2)].first == e2, "");
2896
2897 const MInt dir2 = openEdgeList[openEdgeId(e2)].second;
2898 const MInt id00 = (dir2 == 1) ? 0 : 1;
2899 const MInt id11 = (dir2 == 1) ? 1 : 0;
2900 const MInt v00 = edges[startSetIndex][e2].vertices[id00];
2901 const MInt v11 = edges[startSetIndex][e2].vertices[id11];
2902
2903 if(v11 != v0) continue;
2904
2909 b[i] = vertices[startSetIndex][v11].coordinates[i] - vertices[startSetIndex][v00].coordinates[i];
2910 b_abs +=
b[i] *
b[i];
2911 dotProduct +=
a[i] *
b[i];
2912 }
2913 b_abs = sqrt(b_abs);
2914
2915 dotProduct = fabs(F1 + (dotProduct / (a_abs * b_abs)));
2916
2917
2918 if(dotProduct < 1e-10 && b_abs < a_abs) {
2919
2920
2921
2922
2923
2924
2925
2932 MBool foundMatchingEdge =
false;
2933 for(
MUint k3 = 0; k3 < openEdgeList.size(); k3++) {
2934 e3 = openEdgeList[k3].first;
2935 if(openEdgeId(e3) < 0) continue;
2936 dir3 = openEdgeList[k3].second;
2937 id000 = (dir3 == 1) ? 0 : 1;
2938 id111 = (dir3 == 1) ? 1 : 0;
2939 v000 = edges[startSetIndex][e3].vertices[id000];
2940 v111 = edges[startSetIndex][e3].vertices[id111];
2941 if(v111 == v00 && v000 == v1) {
2942 foundMatchingEdge = true;
2943 break;
2944 }
2945 if(v000 == v00 && v111 == v1) {
2946 foundMatchingEdge = true;
2947 break;
2948 }
2949 }
2950
2951
2952
2953
2954
2955
2956
2957
2958
2959
2960
2961
2962
2963
2964
2965
2966
2967
2968
2969
2970
2971 edges[startSetIndex][e1].vertices[id0] = v00;
2972
2973
2974 edges[startSetIndex][e2].face[1] = face;
2975
2976
2977 MInt otherDir = (edges[startSetIndex][e2].vertices[id0] == v0) ? rdir1 : dir1;
2978 auto pos = faces[startSetIndex][face].edges.begin() + fe;
2979 pos += (dir1 == 1) ? 0 : 1;
2980 faces[startSetIndex][face].edges.insert(pos, make_pair(e2, otherDir));
2981
2982
2984
2985 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v0].edges.size(); edg++) {
2986 if(vertices[startSetIndex][v0].edges[edg] == e1) {
2987 vertices[startSetIndex][v0].edges[edg] = e2;
2988 }
2989 }
2990 } else {
2991
2992 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v0].edges.size(); edg++) {
2993 if(vertices[startSetIndex][v0].edges[edg] == e1) {
2994 auto epos = vertices[startSetIndex][v0].edges.begin() + edg;
2995 vertices[startSetIndex][v0].edges.erase(epos);
2996 }
2997 }
2998 }
2999
3000 if(
m_multiCutCell) vertices[startSetIndex][v00].edges.push_back(e1);
3001
3002
3003 if(foundMatchingEdge) {
3004 MInt otherDir2 = (edges[startSetIndex][e2].vertices[id0] == v0) ? rdir1 : dir1;
3005
3006
3007 for(fe = 0; fe < (signed)faces[startSetIndex][face].edges.size(); fe++) {
3008 if(faces[startSetIndex][face].edges[fe].first == e1) {
3009 faces[startSetIndex][face].edges[fe].first = e3;
3010 faces[startSetIndex][face].edges[fe].second = otherDir2;
3011 }
3012 }
3013
3014
3015 edges[startSetIndex][e3].face[1] = face;
3016
3017
3018 edges[startSetIndex][e1].face[0] = -1;
3019 edges[startSetIndex][e1].face[1] = -1;
3020
3021#ifdef CutCell_DEBUG
3023 cerr << "Replacing edge " << e1 << " with " << e3 << " and adding vertex " << v00 << " to face "
3024 << face << endl;
3025 cerr << "Intermediate edge is " << e2 << " with v " << v1 << " v2 " << v0 << endl;
3026 }
3027#endif
3028
3030
3031 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v00].edges.size(); edg++) {
3032 if(vertices[startSetIndex][v00].edges[edg] == e1) {
3033 vertices[startSetIndex][v00].edges[edg] = e3;
3034 }
3035 }
3036 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v1].edges.size(); edg++) {
3037 if(vertices[startSetIndex][v1].edges[edg] == e1) {
3038 vertices[startSetIndex][v1].edges[edg] = e3;
3039 }
3040 }
3041 } else {
3042
3043 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v00].edges.size(); edg++) {
3044 if(vertices[startSetIndex][v00].edges[edg] == e1) {
3045 auto epos = vertices[startSetIndex][v00].edges.begin() + edg;
3046 vertices[startSetIndex][v00].edges.erase(epos);
3047 }
3048 }
3049 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v1].edges.size(); edg++) {
3050 if(vertices[startSetIndex][v1].edges[edg] == e1) {
3051 auto epos = vertices[startSetIndex][v1].edges.begin() + edg;
3052 vertices[startSetIndex][v1].edges.erase(epos);
3053 }
3054 }
3055 }
3056
3057
3058 openEdgeId(e1) = -1;
3059 openEdgeId(e3) = -1;
3060 }
3061
3062
3063 auto vpos = faces[startSetIndex][face].vertices.begin() + fv;
3064 if(vpos > faces[startSetIndex][face].vertices.end()) {
3065 faces[startSetIndex][face].vertices.push_back(v00);
3066 } else {
3067 faces[startSetIndex][face].vertices.insert(vpos, v00);
3068 }
3069
3070
3071 openEdgeId(e2) = -1;
3072
3073 somethingChanged = true;
3074 alreadyResolved = true;
3075 break;
3076 }
3077 }
3078
3079 if(alreadyResolved) continue;
3080
3081
3082 for(
MInt k2 = 0; (unsigned)k2 < vertices[startSetIndex][v1].edges.size(); k2++) {
3083 const MInt e2 = vertices[startSetIndex][v1].edges[k2];
3084 if(e2 == e1) continue;
3085 if(openEdgeId(e2) < 0) continue;
3086 ASSERT(openEdgeList[openEdgeId(e2)].first == e2, "");
3087 const MInt dir2 = openEdgeList[openEdgeId(e2)].second;
3088 const MInt id00 = (dir2 == 1) ? 0 : 1;
3089 const MInt id11 = (dir2 == 1) ? 1 : 0;
3090 const MInt v00 = edges[startSetIndex][e2].vertices[id00];
3091 const MInt v11 = edges[startSetIndex][e2].vertices[id11];
3092
3093 if(v00 != v1) continue;
3094
3099 b[i] = vertices[startSetIndex][v11].coordinates[i] - vertices[startSetIndex][v00].coordinates[i];
3100 b_abs +=
b[i] *
b[i];
3101 dotProduct +=
a[i] *
b[i];
3102 }
3103 b_abs = sqrt(b_abs);
3104
3105 dotProduct = fabs(F1 + (dotProduct / (a_abs * b_abs)));
3106
3107
3108 if(dotProduct < 1e-10 && b_abs < a_abs) {
3109
3110
3111
3112
3113
3114
3115
3122 MBool foundMatchingEdge =
false;
3123 for(auto& k3 : openEdgeList) {
3124 e3 = k3.first;
3125 if(openEdgeId(e3) < 0) continue;
3126 dir3 = k3.second;
3127 id000 = (dir3 == 1) ? 0 : 1;
3128 id111 = (dir3 == 1) ? 1 : 0;
3129 v000 = edges[startSetIndex][e3].vertices[id000];
3130 v111 = edges[startSetIndex][e3].vertices[id111];
3131 if(v000 == v0 && v111 == v11) {
3132 foundMatchingEdge = true;
3133 break;
3134 } else if(v000 == v11 && v111 == v0) {
3135 foundMatchingEdge = true;
3136 break;
3137 }
3138 }
3139
3140
3141
3142
3143
3144
3145
3146
3147
3148
3149
3150
3151
3152
3153
3154
3155
3156
3157
3158
3159 MInt otherDir = (edges[startSetIndex][e2].vertices[id1] == v1) ? rdir1 : dir1;
3160 auto pos = faces[startSetIndex][face].edges.begin() + fe;
3161 pos += (dir1 == 1) ? 1 : 0;
3162 faces[startSetIndex][face].edges.insert(pos, make_pair(e2, otherDir));
3163
3164
3165 auto vpos = faces[startSetIndex][face].vertices.begin() + fv;
3166 if(vpos > faces[startSetIndex][face].vertices.end()) {
3167 faces[startSetIndex][face].vertices.push_back(v11);
3168 } else {
3169 faces[startSetIndex][face].vertices.insert(vpos, v11);
3170 }
3171
3172
3173 edges[startSetIndex][e2].face[1] = face;
3174
3175
3177 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v1].edges.size(); edg++) {
3178 if(vertices[startSetIndex][v1].edges[edg] == e1) {
3179 vertices[startSetIndex][v1].edges[edg] = e2;
3180 }
3181 }
3182 } else {
3183
3184 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v1].edges.size(); edg++) {
3185 if(vertices[startSetIndex][v1].edges[edg] == e1) {
3186 auto epos = vertices[startSetIndex][v1].edges.begin() + edg;
3187 vertices[startSetIndex][v1].edges.erase(epos);
3188 }
3189 }
3190 }
3191
3192 if(
m_multiCutCell) vertices[startSetIndex][v11].edges.push_back(e1);
3193
3194
3195 edges[startSetIndex][e1].vertices[id1] = v11;
3196
3197
3198 openEdgeId(e2) = -1;
3199
3200
3201 if(foundMatchingEdge) {
3202
3203
3204
3205 const MInt otherDir2 = (edges[startSetIndex][e3].vertices[id1] == v1) ? rdir1 : dir1;
3206
3207
3208 for(fe = 0; fe < (signed)faces[startSetIndex][face].edges.size(); fe++) {
3209 if(faces[startSetIndex][face].edges[fe].first == e1) {
3210 faces[startSetIndex][face].edges[fe].first = e3;
3211 faces[startSetIndex][face].edges[fe].second = otherDir2;
3212 break;
3213 }
3214 }
3215
3216#ifdef CutCell_DEBUG
3218 cerr << "Replacing edge " << e1 << " with " << e3 << " and adding vertex " << v00 << " to face "
3219 << face << endl;
3220 cerr << "Intermediate edge is " << e2 << " with v " << v11 << " v2 " << v1 << endl;
3221 }
3222#endif
3223
3224
3225 edges[startSetIndex][e1].face[0] = -1;
3226 edges[startSetIndex][e1].face[1] = -1;
3227
3228
3230 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v0].edges.size(); edg++) {
3231 if(vertices[startSetIndex][v0].edges[edg] == e1) {
3232 vertices[startSetIndex][v0].edges[edg] = e3;
3233 }
3234 }
3235 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v11].edges.size(); edg++) {
3236 if(vertices[startSetIndex][v11].edges[edg] == e1) {
3237 vertices[startSetIndex][v11].edges[edg] = e3;
3238 }
3239 }
3240 } else {
3241
3242 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v0].edges.size(); edg++) {
3243 if(vertices[startSetIndex][v0].edges[edg] == e1) {
3244 auto epos = vertices[startSetIndex][v0].edges.begin() + edg;
3245 vertices[startSetIndex][v0].edges.erase(epos);
3246 }
3247 }
3248 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v11].edges.size(); edg++) {
3249 if(vertices[startSetIndex][v11].edges[edg] == e1) {
3250 auto epos = vertices[startSetIndex][v11].edges.begin() + edg;
3251 vertices[startSetIndex][v11].edges.erase(epos);
3252 }
3253 }
3254 }
3255
3256 openEdgeId(e1) = -1;
3257 openEdgeId(e3) = -1;
3258
3259
3260 edges[startSetIndex][e3].face[1] = face;
3261 }
3262
3263 somethingChanged = true;
3264 break;
3265 }
3266 }
3267 }
3268 }
3269
3270
3271
3272
3273
3274#ifdef CutCell_DEBUG
3275 for(
MUint i = 0; i < openEdgeList.size(); i++) {
3276 const MInt e = openEdgeList[i].first;
3277 const MInt dir = openEdgeList[i].second;
3278 const MInt id0 = (dir == 1) ? 0 : 1;
3279 const MInt id1 = (dir == 1) ? 1 : 0;
3280 const MInt v0 = edges[startSetIndex][e].vertices[id0];
3281 const MInt v1 = edges[startSetIndex][e].vertices[id1];
3282
3283
3284 if(edges[startSetIndex][e].face[0] == -1 && edges[startSetIndex][e].face[1] == -1) continue;
3285
3286
3287 for(
MInt j = 0; (unsigned)j < vertices[startSetIndex][v0].edges.size(); j++) {
3288 const MInt e2 = vertices[startSetIndex][v0].edges[j];
3289 if(e2 == e) continue;
3290 const MInt dir2 = openEdgeList[openEdgeId(e2)].second;
3291 const MInt id00 = (dir2 == 1) ? 0 : 1;
3292 const MInt id11 = (dir2 == 1) ? 1 : 0;
3293 const MInt v00 = edges[startSetIndex][e2].vertices[id00];
3294 const MInt v11 = edges[startSetIndex][e2].vertices[id11];
3295
3296
3297 if(edges[startSetIndex][e2].face[0] == -1 && edges[startSetIndex][e2].face[1] == -1) continue;
3298 if(onlySingleEdges) {
3299 ASSERT(v0 != v00 || v1 != v11, "");
3300 ASSERT(v0 != v11 || v1 != v00, "");
3301 }
3302 }
3303 }
3304#endif
3305
3306
3307 for(
MInt faceId = 0; faceId < (signed)faces[startSetIndex].size(); faceId++) {
3308 for(
MInt edgeId = 0; edgeId < (signed)faces[startSetIndex][faceId].edges.size(); edgeId++) {
3309 const MInt edge = faces[startSetIndex][faceId].edges[edgeId].first;
3310 ASSERT(edges[startSetIndex][edge].face[0] > -1, "");
3311 }
3312 }
3313
3314
3315 for(
MInt v = 0; v < (signed)vertices[startSetIndex].size(); v++) {
3316 for(
MInt edgeId = 0; (unsigned)edgeId < vertices[startSetIndex][v].edges.size(); edgeId++) {
3317 const MInt edge = vertices[startSetIndex][v].edges[edgeId];
3318 ASSERT(edges[startSetIndex][edge].face[0] > -1, "");
3319 }
3320 }
3321
3322
3323 for(
MInt faceId = 0; faceId < (signed)faces[startSetIndex].size(); faceId++) {
3324 ASSERT(faces[startSetIndex][faceId].edges.size() == faces[startSetIndex][faceId].vertices.size(), "");
3325 for(
MInt edgeId = 0; edgeId < (signed)faces[startSetIndex][faceId].edges.size(); edgeId++) {
3326 const MInt edge = faces[startSetIndex][faceId].edges[edgeId].first;
3327 const MInt newVertexId = faces[startSetIndex][faceId].vertices[edgeId];
3328 const MInt direction = faces[startSetIndex][faceId].edges[edgeId].second;
3329 MInt vertexId = edges[startSetIndex][edge].vertices[0];
3330 if(direction == -1) vertexId = edges[startSetIndex][edge].vertices[1];
3331 ASSERT(newVertexId == vertexId, "");
3332 }
3333 }
3334
3335
3336 std::list<std::pair<MInt, MInt>> openEdgesOld(openEdges);
3337 openEdges.clear();
3338 for(auto& it1 : openEdgesOld) {
3339 if(openEdgeId(it1.first) < 0) continue;
3340 openEdges.push_back(it1);
3341 }
3342
3343
3344#ifdef CutCell_DEBUG
3345 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
3346 cerr << openEdges.size() << " open edge remaining: " << endl;
3347 cerr << "Count: " << faces[startSetIndex].size() << " faces " << edges[startSetIndex].size() << " edges "
3348 << vertices[startSetIndex].size() << " vertices " << endl;
3349
3350 for(auto it2 = openEdges.begin(); it2 != openEdges.end(); it2++) {
3351 const MInt edgeId = (*it2).first;
3352 cerr << " openEdge with Id " << edgeId << " v1 " << edges[startSetIndex][edgeId].vertices[0] << " "
3353 << edges[startSetIndex][edgeId].vertices[1] << endl;
3354 }
3355 }
3356
3357
3358
3359
3360
3361
3362
3363
3364
3365 if(openEdges.size() == 1) {
3366 cerr << " Only 1 open edge remaining: This means trouble! " << endl;
3367 cerr << "Re-checking coordinate-distances! " << endl;
3368 for(
MInt i = 0; (unsigned)i < vertices[startSetIndex].size(); i++) {
3369 for(
MInt j = 0; (unsigned)j < vertices[startSetIndex].size(); j++) {
3370 if(i == j) continue;
3372 coord_diff +=
3373 pow(vertices[startSetIndex][i].coordinates[0] - vertices[startSetIndex][j].coordinates[0], 2);
3374 coord_diff +=
3375 pow(vertices[startSetIndex][i].coordinates[1] - vertices[startSetIndex][j].coordinates[1], 2);
3376 coord_diff +=
3377 pow(vertices[startSetIndex][i].coordinates[2] - vertices[startSetIndex][j].coordinates[2], 2);
3378
3379 MFloat dif_old = coord_diff;
3380 coord_diff = pow(coord_diff, 0.5);
3381
3382 if(coord_diff < cellLength0) {
3383 cerr << " Coord-dif " << coord_diff << " i " << i << " " << j << endl;
3384 cerr <<
" dif-old " << dif_old <<
" " <<
m_eps * 10 << endl;
3385 }
3386 }
3387 }
3388 }
3389
3390#endif
3391
3392
3393
3394
3395
3396
3397
3398 for(auto it1 = openEdges.begin(); it1 != openEdges.end(); it1++) {
3399
3400 const MInt e1 = (*it1).first;
3401 const MInt dir1 = (*it1).second;
3402 MInt v11 = edges[startSetIndex][e1].vertices[0];
3403 MInt v12 = edges[startSetIndex][e1].vertices[1];
3404 if(dir1 == -1) {
3406 v11 = v12;
3407 v12 = tmp;
3408 }
3412 a[i] = vertices[startSetIndex][v12].coordinates[i] - vertices[startSetIndex][v11].coordinates[i];
3413 a_abs +=
a[i] *
a[i];
3414 }
3415 a_abs = sqrt(a_abs);
3416
3417 for(auto& openEdge : openEdges) {
3418 const MInt e2 = openEdge.first;
3419 const MInt dir2 = openEdge.second;
3420
3421 if(e2 == e1) {
3422 dotProdMatrix(e1, e2) = 1000.0;
3423 dotProdMatrix(e2, e1) = 1000.0;
3424 continue;
3425 }
3426 MInt v21 = edges[startSetIndex][e2].vertices[0];
3427 MInt v22 = edges[startSetIndex][e2].vertices[1];
3428 if(dir2 == -1) {
3430 v21 = v22;
3431 v22 = tmp;
3432 }
3433
3434 if(v12 == v21) {
3435
3440 b[i] = vertices[startSetIndex][v22].coordinates[i] - vertices[startSetIndex][v21].coordinates[i];
3441 b_abs +=
b[i] *
b[i];
3442 dotProduct +=
a[i] *
b[i];
3443 }
3444 b_abs = sqrt(b_abs);
3445 dotProduct = abs(F1 - abs(dotProduct / (a_abs * b_abs)));
3446 dotProdMatrix(e1, e2) = dotProduct;
3447 } else {
3448
3449 dotProdMatrix(e1, e2) = 1000.0;
3450 }
3451 }
3452 }
3453
3454
3455
3456 while(!openEdges.empty()) {
3457 faceVertices.clear();
3458
3460 MInt eLast = openEdges.front().first;
3461 MInt dir = openEdges.front().second;
3466 surfaceIdentificatorCounters[i] = 0;
3467 }
3468
3469
3470 faces[startSetIndex].emplace_back(faceId, faceType, bodyId);
3471
3472 faces[startSetIndex][noFaces].edges.emplace_back(eLast, dir);
3473 faces[startSetIndex][noFaces].isLine = 1;
3474
3475
3476 openEdges.pop_front();
3477
3478 MInt startVertex = edges[startSetIndex][eLast].vertices[0];
3479 MInt vertex = edges[startSetIndex][eLast].vertices[1];
3480
3481 if(dir == -1) {
3482 MInt tmp = startVertex;
3483 startVertex = vertex;
3484 vertex = tmp;
3485 }
3486
3487 faces[startSetIndex][noFaces].vertices.push_back(startVertex);
3488
3489#ifdef CutCell_DEBUG
3490 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
3491 cerr << "Adding face line starting with edge " << eLast << " with vertex " << startVertex << " and "
3492 << vertex << endl;
3493 }
3494#endif
3495
3496
3497 faceVertices.push_back(vertex);
3498
3499 for(
MInt surfaceIdentificator : vertices[startSetIndex][vertex].surfaceIdentificators) {
3500 surfaceIdentificatorCounters[surfaceIdentificator]++;
3501 }
3502
3503 edges[startSetIndex][eLast].face[1] = noFaces;
3504
3505
3506 while(vertex != startVertex) {
3507
3508 auto bestEdgeIt = openEdges.end();
3509 MFloat minDotProduct = 1001.0;
3510
3511
3512
3513
3514 for(auto it = openEdges.begin(); it != openEdges.end(); it++) {
3515
3516 e = (*it).first;
3517 if(dotProdMatrix(eLast, e) < minDotProduct) {
3518
3519 dir = (*bestEdgeIt).second;
3520 MInt vStart = edges[startSetIndex][e].vertices[0];
3521 MInt vEnd = edges[startSetIndex][e].vertices[1];
3522 if(dir == -1) {
3524 vStart = vEnd;
3525 vEnd = tmp;
3526 }
3527 if(vStart == vertex) {
3528 minDotProduct = dotProdMatrix(eLast, e);
3529 bestEdgeIt = it;
3530 } else if(vEnd == vertex) {
3531 minDotProduct = dotProdMatrix(eLast, e);
3532 bestEdgeIt = it;
3533 }
3534 }
3535 }
3536
3537 if(bestEdgeIt == openEdges.end()) {
3539 cerr <<
grid().domainId() <<
"open edges for cell " <<
cellId <<
": ";
3540 cerr << "Open edges size: " << openEdges.size() << endl;
3541 for(auto& openEdge : openEdges) {
3542 cerr << openEdge.first;
3543 }
3544 mTerm(1, AT_,
"Open edge for cutCell!");
3545 }
3546
3547
3548 e = (*bestEdgeIt).first;
3549 dir = (*bestEdgeIt).second;
3550 MInt vStart = edges[startSetIndex][e].vertices[0];
3551 MInt vEnd = edges[startSetIndex][e].vertices[1];
3552 if(dir == -1) {
3554 vStart = vEnd;
3555 vEnd = tmp;
3556 }
3557 if(vStart != vertex) {
3558 if(vEnd == vertex) {
3559
3560 if(dir == -1) {
3561 dir = 1;
3562 } else {
3563 dir = -1;
3564 }
3565 MInt vEndOld = vEnd;
3566 vEnd = vStart;
3567 vStart = vEndOld;
3568 } else {
3569 cerr <<
grid().domainId() <<
" missmatching edges for " <<
cellId <<
" "
3570 <<
grid().tree().globalId(cellId) <<
" vertices: " << vStart <<
" " << vEnd <<
" " << vertex
3571 << endl;
3572 }
3573 }
3574
3575 eLast = e;
3576
3577
3578 faces[startSetIndex][noFaces].edges.emplace_back(e, dir);
3579 faces[startSetIndex][noFaces].vertices.push_back(vStart);
3580
3581#ifdef CutCell_DEBUG
3582 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
3583 cerr << "Adding edge " << e << " with vertex " << vStart << " and " << vEnd << endl;
3584 cerr << "The new face now has " << faces[startSetIndex][noFaces].vertices.size() << " vertices!" << endl;
3585 }
3586#endif
3587
3588 openEdges.erase(bestEdgeIt);
3589 vertex = vEnd;
3590 faceVertices.push_back(vertex);
3591 for(
MInt surfaceIdentificator : vertices[startSetIndex][vertex].surfaceIdentificators)
3592 surfaceIdentificatorCounters[surfaceIdentificator]++;
3593 edges[startSetIndex][e].face[1] = noFaces;
3594 }
3595
3596 if(faceVertices.size() <= 2) {
3597 cerr <<
grid().domainId() <<
" face with only two vertices for " <<
cellId <<
" "
3598 <<
grid().tree().globalId(cellId) <<
" face : " << noFaces << endl;
3599 }
3600
3601
3602 MInt noDoubleVertices = 0;
3603 ASSERT(vertices[startSetIndex].size() <= (unsigned)maxNoVertices, "");
3604 for(
MInt v = 0; (unsigned)v < vertices[startSetIndex].size(); v++) {
3605 vertexTouches[v] = 0;
3606 }
3607 for(
MInt fv = 0; (unsigned)fv < faceVertices.size(); fv++) {
3608 vertexTouches[faceVertices[fv]]++;
3609 }
3610 for(
MInt v = 0; (unsigned)v < vertices[startSetIndex].size(); v++) {
3611 if(vertexTouches[v] > 1) {
3612 noDoubleVertices++;
3613 }
3614 }
3615
3616 if(noDoubleVertices != 0) {
3617 cerr <<
"double vertices for " <<
cellId <<
" on " <<
grid().domainId() << endl;
3618 error = true;
3619 }
3620
3621
3622
3623 MInt bestIdentificator = -1;
3626 if(surfaceIdentificatorCounters[i] > maxHits) {
3627 maxHits = surfaceIdentificatorCounters[i];
3628 bestIdentificator = i;
3629 }
3630 }
3631 ASSERT(bestIdentificator > -1, "");
3633 faceType = (bestIdentificator <
m_noDirs ? 0 : 1);
3634 if(faceType == 0) {
3635 faceId = bestIdentificator;
3636 } else {
3637 bodyId = bestIdentificator -
m_noDirs;
3638 }
3639 faces[startSetIndex][noFaces].faceType = faceType;
3640 faces[startSetIndex][noFaces].faceId = faceId;
3641 faces[startSetIndex][noFaces].bodyId = bodyId;
3642
3643
3644 MBool partnerFound =
false;
3645 for(
MInt e2 = 0; (unsigned)e2 < faces[startSetIndex][noFaces].edges.size(); e2++) {
3646 MInt edge = faces[startSetIndex][noFaces].edges[e2].first;
3647 MInt otherFace = edges[startSetIndex][edge].face[0];
3648 MInt otherFaceType = faces[startSetIndex][otherFace].faceType;
3649 MInt otherFaceId = faces[startSetIndex][otherFace].faceId;
3650 MInt otherBodyId = faces[startSetIndex][otherFace].bodyId;
3651 if(otherFaceType == faceType && otherFaceId == faceId && otherBodyId == bodyId) {
3652 partnerFound = true;
3653 faces[startSetIndex][noFaces].normal[0] = faces[startSetIndex][otherFace].normal[0];
3654 faces[startSetIndex][noFaces].normal[1] = faces[startSetIndex][otherFace].normal[1];
3655 faces[startSetIndex][noFaces].normal[2] = faces[startSetIndex][otherFace].normal[2];
3656 faces[startSetIndex][noFaces].w = faces[startSetIndex][otherFace].w;
3657
3658 break;
3659 }
3660 }
3661
3662 ASSERT(partnerFound, "");
3663
3664 noFaces++;
3665 }
3666
3667 }
3668
3669 }
3670
3671#ifdef CutCell_DEBUG
3672 if(
grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
3673 cerr << "Final count: " << faces[startSetIndex].size() << " faces " << edges[startSetIndex].size() << " edges "
3674 << vertices[startSetIndex].size() << " vertices " << endl;
3675 }
3676#endif
3677
3678
3679
3680 for(
MInt faceId = 0; faceId < (signed)faces[startSetIndex].size(); faceId++) {
3681 ASSERT(faces[startSetIndex][faceId].edges.size() == faces[startSetIndex][faceId].vertices.size(), "");
3682 for(
MInt edgeId = 0; edgeId < (signed)faces[startSetIndex][faceId].edges.size(); edgeId++) {
3683 const MInt edge = faces[startSetIndex][faceId].edges[edgeId].first;
3684 const MInt newVertexId = faces[startSetIndex][faceId].vertices[edgeId];
3685 const MInt direction = faces[startSetIndex][faceId].edges[edgeId].second;
3686 MInt vertexId = edges[startSetIndex][edge].vertices[0];
3687 if(direction == -1) vertexId = edges[startSetIndex][edge].vertices[1];
3688 ASSERT(newVertexId == vertexId, "");
3689 }
3690 }
3691
3692 RECORD_TIMER_STOP(tCutFace_3b);
3693 RECORD_TIMER_START(tCutFace_4);
3694
3695
3696
3697 ASSERT(faceStack.empty(), "");
3698
3699
3700 if(faces[startSetIndex].size() == 0) {
3701 cutCells.emplace_back(cellId, &
a_coordinate(cellId, 0));
3702 }
3703
3704
3705 for(
MInt faceCounter = 0; (unsigned)faceCounter < faces[startSetIndex].size(); faceCounter++) {
3706 if(faces[startSetIndex][faceCounter].cutCell > -1) continue;
3707
3708
3709 if(faces[startSetIndex][faceCounter].isLine) continue;
3710
3711#ifdef CutCell_DEBUG
3713 cerr << "faceCounter: " << faceCounter << endl;
3714 }
3715#endif
3716
3717 faceStack.push(faceCounter);
3718 const MInt currentCutCell = cutCells.size();
3719 cutCells.emplace_back(cellId, &
a_coordinate(cellId, 0));
3720 faces[startSetIndex][faceCounter].cutCell = currentCutCell;
3721 cutCells[currentCutCell].faces.push_back(faceCounter);
3722
3723 while(!faceStack.empty()) {
3724 MInt currentFace = faceStack.top();
3725 faceStack.pop();
3726
3727 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][currentFace].edges.size(); e++) {
3728 MInt edge = faces[startSetIndex][currentFace].edges[e].first;
3729 MInt otherFace = edges[startSetIndex][edge].face[0];
3730 if(otherFace == currentFace) {
3731 otherFace = edges[startSetIndex][edge].face[1];
3732 }
3733
3734
3735 if(faces[startSetIndex][otherFace].cutCell == -1) {
3736#ifdef CutCell_DEBUG
3738 cerr << "faceCounter: " << faceCounter << " cutCell " << currentCutCell
3739 << " is adding neighbor: " << otherFace << endl;
3740 }
3741#endif
3742 cutCells[currentCutCell].faces.push_back(otherFace);
3743 faces[startSetIndex][otherFace].cutCell = currentCutCell;
3744 faceStack.push(otherFace);
3745 }
3746 }
3747 }
3748 }
3749 ASSERT(faceStack.empty(), "");
3750
3751 RECORD_TIMER_STOP(tCutFace_4);
3752
3753
3754 RECORD_TIMER_START(tCutFace_5a);
3756 RECORD_TIMER_STOP(tCutFace_5a);
3757
3758#ifdef CutCell_DEBUG
3759 if(
globalTimeStep == debugTimeStep && cellId == debugCellId &&
grid().domainId() == debugDomainId) {
3761 }
3762#endif
3763
3764 if(error) {
3765 cerr << "[" << -1 << "]"
3766 << ": Warning: removal/suppression of doubly touched vertices seems to have failed!" << endl;
3768 }
3769
3770
3771 RECORD_TIMER_START(tCutFace_5b);
3772
3773
3774
3775
3776
3777
3778
3779
3780
3781
3783
3784
3785 for(
MInt e = 0; (unsigned)e < edges[startSetIndex].size(); e++) {
3786 MInt face = edges[startSetIndex][e].face[0];
3787 if(face < 0)
3788 edgeCutCellPointer[e] = -1;
3789 else
3790 edgeCutCellPointer[e] = faces[startSetIndex][face].cutCell;
3791 }
3792 for(
MInt f = 0; (unsigned)f < faces[startSetIndex].size(); f++) {
3793 multiFaceConnection[f] = -1;
3794 }
3795 for(
MInt c = 0; (unsigned)c < cutCells.size(); c++) {
3797 case 1:
3798 case 3:
3799 case 4:
3800 {
3801 list<MInt> bodyVertices;
3802 set<MInt> locBodySrfcs;
3804 isBodyEdge.fill(-2);
3805
3806
3807 for(auto& vert : vertices[startSetIndex]) {
3808 sort(vert.edges.begin(), vert.edges.end());
3809 auto last = unique(vert.edges.begin(), vert.edges.end());
3810 vert.edges.erase(last, vert.edges.end());
3811 }
3812
3813 for(
MInt e = 0; (unsigned)e < edges[startSetIndex].size(); e++) {
3814 if(edgeCutCellPointer[e] != c) continue;
3815 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].faceType == 0
3816 && faces[startSetIndex][edges[startSetIndex][e].face[1]].faceType == 0) {
3817 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].faceId
3818 == faces[startSetIndex][edges[startSetIndex][e].face[1]].faceId) {
3819 pureBodyEdges.push_back(e);
3820 }
3822 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].faceType == 1
3823 && faces[startSetIndex][edges[startSetIndex][e].face[1]].faceType == 1) {
3824 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].bodyId
3825 == faces[startSetIndex][edges[startSetIndex][e].face[1]].bodyId) {
3826 pureBodyEdges.push_back(e);
3827 isBodyEdge(e) = faces[startSetIndex][edges[startSetIndex][e].face[0]].bodyId;
3828 bodyVertices.push_back(edges[startSetIndex][e].vertices[0]);
3829 bodyVertices.push_back(edges[startSetIndex][e].vertices[1]);
3830 locBodySrfcs.insert(faces[startSetIndex][edges[startSetIndex][e].face[0]].bodyId);
3831 } else {
3832 isBodyEdge(e) = -1;
3833 }
3834 }
3835 }
3836 }
3837
3839 bodyVertices.sort();
3840 bodyVertices.unique();
3841 const MInt noLocBodySrfcs = (signed)locBodySrfcs.size();
3842 MInt concaveCnt = 0;
3843 for(
MInt vert : bodyVertices) {
3844 if(vertices[startSetIndex][vert].
pointType > 1) {
3845 MInt noOuterEdges = 0;
3846 MInt noBodyEdges = 0;
3847 MInt outerEdges[6] = {-1};
3848 MInt bodyEdges[6] = {-1};
3849 for(
MInt edge : vertices[startSetIndex][vert].edges) {
3850 if(isBodyEdge(edge) == -1) {
3851 outerEdges[noOuterEdges] = edge;
3852 noOuterEdges++;
3853 }
3854 if(isBodyEdge(edge) > -1) {
3855 bodyEdges[noBodyEdges] = edge;
3856 noBodyEdges++;
3857 }
3858 }
3859 if(noOuterEdges == 2 && noBodyEdges > 0) {
3860 for(
MInt k = 0; k < noBodyEdges; k++) {
3861 MInt edge = bodyEdges[k];
3862 MInt otherVert = (edges[startSetIndex][edge].vertices[0] == vert) ? 1 : 0;
3863 otherVert = edges[startSetIndex][edge].vertices[otherVert];
3864 MInt face0 = edges[startSetIndex][edge].face[0];
3865 MInt face1 = edges[startSetIndex][edge].face[1];
3866
3871 MInt noEdges0 = (signed)faces[startSetIndex][face0].edges.size();
3872 MInt noEdges1 = (signed)faces[startSetIndex][face1].edges.size();
3876 for(
MInt e = 0; e < noEdges0; e++) {
3877 if(faces[startSetIndex][face0].edges[e].first == edge) {
3878 edge0 = e;
3879 dir0 = faces[startSetIndex][face0].edges[e].second;
3880 if(prevEdge < 0 || nextEdge < 0) {
3881 MInt vA = edges[startSetIndex][edge].vertices[dir0 == 1 ? 1 : 0];
3882 MInt otherEdge = (vA == vert) ? (edge0 + 1) % noEdges0 : (edge0 + noEdges0 - 1) % noEdges0;
3883 otherEdge = faces[startSetIndex][face0].edges[otherEdge].first;
3884 if(otherEdge == outerEdges[0] || otherEdge == outerEdges[1]) {
3885 MInt otherEdge2 = (otherEdge == outerEdges[0]) ? outerEdges[1] : outerEdges[0];
3886 nextEdge = (vA == vert) ? otherEdge : otherEdge2;
3887 prevEdge = (vA == vert) ? otherEdge2 : otherEdge;
3888 side = 0;
3889 }
3890 }
3891 }
3892 }
3893 for(
MInt e = 0; e < noEdges1; e++) {
3894 if(faces[startSetIndex][face1].edges[e].first == edge) {
3895 edge1 = e;
3896 dir1 = faces[startSetIndex][face1].edges[e].second;
3897 if(prevEdge < 0 || nextEdge < 0) {
3898 MInt vA = edges[startSetIndex][edge].vertices[dir1 == 1 ? 1 : 0];
3899 MInt otherEdge = (vA == vert) ? (edge1 + 1) % noEdges1 : (edge1 + noEdges1 - 1) % noEdges1;
3900 otherEdge = faces[startSetIndex][face1].edges[otherEdge].first;
3901 if(otherEdge == outerEdges[0] || otherEdge == outerEdges[1]) {
3902 MInt otherEdge2 = (otherEdge == outerEdges[0]) ? outerEdges[1] : outerEdges[0];
3903 nextEdge = (vA == vert) ? otherEdge : otherEdge2;
3904 prevEdge = (vA == vert) ? otherEdge2 : otherEdge;
3905 side = 1;
3906 }
3907 }
3908 }
3909 }
3910 if(edge0 < 0 || edge1 < 0)
mTerm(1, AT_,
"ERROR: pure body edge not found.");
3911 if(dir0 == dir1)
mTerm(1, AT_,
"ERROR: pure body edge not found (2).");
3912 if(faces[startSetIndex][face0].edges[edge0].first != edge
3913 || faces[startSetIndex][face1].edges[edge1].first != edge)
3914 mTerm(1, AT_,
"ERROR: pure body edge not found (3).");
3915 ASSERT(outerEdges[0] > -1 && outerEdges[1] > -1, "");
3916
3917 if(prevEdge > -1 && nextEdge > -1) {
3918 ASSERT(side > -1, "");
3919 MInt id_p0 = edges[startSetIndex][prevEdge].vertices[0] == vert ? 0 : 1;
3920 MInt id_p1 = edges[startSetIndex][nextEdge].vertices[0] == vert ? 0 : 1;
3921 if(edges[startSetIndex][prevEdge].vertices[id_p0] != vert)
3922 mTerm(1, AT_,
"ERROR: pure body edge not found (4).");
3923 if(edges[startSetIndex][nextEdge].vertices[id_p1] != vert)
3924 mTerm(1, AT_,
"ERROR: pure body edge not found (5).");
3925 MInt otherId[2] = {1, 0};
3926 MInt vA2 = edges[startSetIndex][prevEdge].vertices[otherId[id_p0]];
3927 MInt vB2 = edges[startSetIndex][prevEdge].vertices[id_p0];
3928 MInt vA = edges[startSetIndex][nextEdge].vertices[id_p1];
3929 MInt vB = edges[startSetIndex][nextEdge].vertices[otherId[id_p1]];
3930
3932 MInt vO = otherVert;
3933
3945 MFloat normal[3] = {F0, F0, F0};
3946 MFloat area0 = faces[startSetIndex][face0].area;
3947 MFloat area1 = faces[startSetIndex][face1].area;
3948 ASSERT(area0 + area1 > F0, "");
3950 normal[i] = (area0 * faces[startSetIndex][face0].normal[i]
3951 + area1 * faces[startSetIndex][face1].normal[i])
3952 / (area0 + area1);
3953 }
3956 vertices[startSetIndex][vB2].coordinates[i] - vertices[startSetIndex][vA2].coordinates[i];
3958 vertices[startSetIndex][vB].coordinates[i] - vertices[startSetIndex][vA].coordinates[i];
3960 vertices[startSetIndex][vO].coordinates[i] - vertices[startSetIndex][vV].coordinates[i];
3961 vec1[i] = d1;
3962 vec2[i] = d2;
3963 vec3[i] = d3;
3964 dotp += d1 * d2;
3968 }
3969 abs1 = sqrt(abs1);
3970 abs2 = sqrt(abs2);
3971 abs3 = sqrt(abs3);
3972 dotp /= (abs1 * abs2);
3974 vec1[i] /= abs1;
3975 vec2[i] /= abs2;
3976 vec3[i] /= abs3;
3977 }
3980 dotp2 += ori[i] * normal[i];
3981 }
3984 dotp3 += ori[i] * normal[i];
3985 }
3988 dotp4 += ori[i] * normal[i];
3989 }
3990
3991 if(dotp2 < -0.1 && fabs(F1 - dotp) > 1e-2
3992 &&
mMin(area0, area1) > 1e-5 * pow(cellLength0, (
MFloat)(
nDim - 1))) {
3993 if(noLocBodySrfcs + concaveCnt > maxNoSurfaces) {
3994 cerr << "Warning: removal of concave polygons was skipped, since "
3995 "FvBndryCell<nDim>::m_maxNoSurfaces is not large enough."
3996 << endl;
3997 continue;
3998 }
3999
4000 if(dotp3 > -1e-3 && dotp4 > -1e-3) {
4001 auto it2 = pureBodyEdges.begin();
4002 while(it2 != pureBodyEdges.end()) {
4003
4004
4005 MInt bedge = (*it2);
4006 if(edgeCutCellPointer[bedge] != c) {
4007 it2++;
4008 continue;
4009 }
4010 MInt faceA = edges[startSetIndex][bedge].face[0];
4011 MInt faceB = edges[startSetIndex][bedge].face[1];
4012 if((face0 == faceA && face1 == faceB) || (face1 == faceA && face0 == faceB)) {
4013 it2 = pureBodyEdges.erase(it2);
4014 concaveCnt++;
4015 } else {
4016 it2++;
4017 }
4018 }
4019 }
4020 }
4021 }
4022 }
4023 }
4024 }
4025 }
4026 }
4027
4028 multiFaces.clear();
4029
4030
4031 MBool somethingChanged =
true;
4032 MInt currentMultiFace = 0;
4033 while(!pureBodyEdges.empty()) {
4034 MInt pbedge = pureBodyEdges.front();
4035 pureBodyEdges.pop_front();
4036 multiFaces.emplace_back();
4037 MInt pbface0 = edges[startSetIndex][pbedge].face[0];
4038 multiFaceConnection[pbface0] = currentMultiFace;
4039 multiFaces[currentMultiFace].faces.push_back(pbface0);
4040 multiFaces[currentMultiFace].bodyId = faces[startSetIndex][pbface0].bodyId;
4041 MInt pbface1 = edges[startSetIndex][pbedge].face[1];
4042 multiFaceConnection[pbface1] = currentMultiFace;
4043 multiFaces[currentMultiFace].faces.push_back(pbface1);
4044 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][pbface0].edges.size(); e++) {
4045 if(faces[startSetIndex][pbface0].edges[e].first != pbedge) {
4046 multiFaces[currentMultiFace].edges.emplace_back(faces[startSetIndex][pbface0].edges[e].first,
4047 faces[startSetIndex][pbface0].edges[e].second);
4048 } else {
4049
4051 for(
MInt ee = 0; (unsigned)ee < faces[startSetIndex][pbface1].edges.size(); ee++) {
4052 if(faces[startSetIndex][pbface1].edges[ee].first == pbedge) {
4053 edgeIndex = ee;
4054 break;
4055 }
4056 }
4057 for(
MInt ee = 1; (unsigned)ee < faces[startSetIndex][pbface1].edges.size(); ee++) {
4058 MInt eee = (edgeIndex + ee) % faces[startSetIndex][pbface1].edges.size();
4059 multiFaces[currentMultiFace].edges.emplace_back(faces[startSetIndex][pbface1].edges[eee].first,
4060 faces[startSetIndex][pbface1].edges[eee].second);
4061 }
4062 }
4063 }
4064 somethingChanged = true;
4065
4066 while(somethingChanged) {
4067 somethingChanged = false;
4068 auto it = pureBodyEdges.begin();
4069 while(it != pureBodyEdges.end()) {
4070
4072 if(edgeCutCellPointer[edge] != c) {
4073 it++;
4074 continue;
4075 }
4076 MInt face0 = edges[startSetIndex][edge].face[0];
4077 MInt face1 = edges[startSetIndex][edge].face[1];
4078 if(multiFaceConnection[face0] == currentMultiFace && multiFaceConnection[face1] == currentMultiFace) {
4079 it = pureBodyEdges.erase(it);
4080 } else {
4081 if(multiFaceConnection[face0] == currentMultiFace) {
4082 multiFaceConnection[face1] = currentMultiFace;
4083 it = pureBodyEdges.erase(it);
4084 multiFaces[currentMultiFace].faces.push_back(face1);
4085 auto it2 = multiFaces[currentMultiFace].edges.begin();
4086 for(it2 = multiFaces[currentMultiFace].edges.begin();
4087 it2 != multiFaces[currentMultiFace].edges.end();
4088 it2++) {
4089 MInt edge2 = (*it2).first;
4090 if(edge2 == edge) {
4091 break;
4092 }
4093 }
4094
4096 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face1].edges.size(); e++) {
4097 if(faces[startSetIndex][face1].edges[e].first == edge) {
4098 edgeIndex = e;
4099 break;
4100 }
4101 }
4102 for(
MInt e = 1; (unsigned)e < faces[startSetIndex][face1].edges.size(); e++) {
4103 MInt ee = (edgeIndex + e) % faces[startSetIndex][face1].edges.size();
4104 multiFaces[currentMultiFace].edges.insert(
4105 it2, make_pair(faces[startSetIndex][face1].edges[ee].first,
4106 faces[startSetIndex][face1].edges[ee].second));
4107 }
4108 multiFaces[currentMultiFace].edges.erase(it2);
4109 somethingChanged = true;
4110 } else if(multiFaceConnection[face1] == currentMultiFace) {
4111 multiFaceConnection[face0] = currentMultiFace;
4112 it = pureBodyEdges.erase(it);
4113 multiFaces[currentMultiFace].faces.push_back(face0);
4114 auto it2 = multiFaces[currentMultiFace].edges.begin();
4115 for(it2 = multiFaces[currentMultiFace].edges.begin();
4116 it2 != multiFaces[currentMultiFace].edges.end();
4117 it2++) {
4118 MInt edge2 = (*it2).first;
4119 if(edge2 == edge) {
4120 break;
4121 }
4122 }
4123
4125 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face0].edges.size(); e++) {
4126 if(faces[startSetIndex][face0].edges[e].first == edge) {
4127 edgeIndex = e;
4128 break;
4129 }
4130 }
4131 for(
MInt e = 1; (unsigned)e < faces[startSetIndex][face0].edges.size(); e++) {
4132 MInt ee = (edgeIndex + e) % faces[startSetIndex][face0].edges.size();
4133 multiFaces[currentMultiFace].edges.insert(
4134 it2, make_pair(faces[startSetIndex][face0].edges[ee].first,
4135 faces[startSetIndex][face0].edges[ee].second));
4136 }
4137 multiFaces[currentMultiFace].edges.erase(it2);
4138 somethingChanged = true;
4139 } else {
4140 it++;
4141 }
4142 }
4143 }
4144 }
4145 currentMultiFace++;
4146 }
4147 } break;
4148 case 2:
4149 {
4150 for(
MInt e = 0; (unsigned)e < edges[startSetIndex].size(); e++) {
4151 if(edgeCutCellPointer[e] != c) continue;
4152 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].faceType == 1
4153 && faces[startSetIndex][edges[startSetIndex][e].face[1]].faceType == 1) {
4154 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].bodyId
4155 == faces[startSetIndex][edges[startSetIndex][e].face[1]].bodyId) {
4156 pureBodyEdges.push_back(e);
4157 }
4158 } else if(faces[startSetIndex][edges[startSetIndex][e].face[0]].faceType == 0
4159 && faces[startSetIndex][edges[startSetIndex][e].face[1]].faceType == 0) {
4160 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].faceId
4161 == faces[startSetIndex][edges[startSetIndex][e].face[1]].faceId) {
4162 pureBodyEdges.push_back(e);
4163 }
4164 }
4165 }
4166
4167 const MInt faceNum = faces[startSetIndex].size();
4168 ASSERT(faceNum <= maxNoFaces, "");
4169 for(
MInt i = 0; i < faceNum; i++)
4170 for(
MInt j = 0; j < faceNum; j++)
4171 normalDotProduct(i, j) = F0;
4172 MInt noBodyFaces = 0;
4173 for(
MInt f = 0; (unsigned)f < cutCells[c].faces.size(); f++) {
4174 MInt face = cutCells[c].faces[f];
4175 if(faces[startSetIndex][face].faceType == 1) bodyFaces[noBodyFaces++] = face;
4176 }
4177 for(
MInt f = 0; f < noBodyFaces; f++) {
4178 MInt face1 = bodyFaces[f];
4179 for(
MInt ff = 0; ff <= f; ff++) {
4180 MInt face2 = bodyFaces[ff];
4182 normalDotProduct(face1, face2) +=
4183 faces[startSetIndex][face1].normal[i] * faces[startSetIndex][face2].normal[i];
4184 normalDotProduct(face1, face2) = F1 - normalDotProduct(face1, face2);
4185 normalDotProduct(face2, face1) = normalDotProduct(face1, face2);
4186 }
4187 }
4188
4189 multiFaces.clear();
4190
4191
4192 MBool somethingChanged =
true;
4193 MInt currentMultiFace = 0;
4194 while(!pureBodyEdges.empty()) {
4195 MInt pbedge = pureBodyEdges.front();
4196 pureBodyEdges.pop_front();
4197 MInt pbface0 = edges[startSetIndex][pbedge].face[0];
4198 MInt pbface1 = edges[startSetIndex][pbedge].face[1];
4200 multiFaces.emplace_back();
4201 multiFaceConnection[pbface0] = currentMultiFace;
4202 multiFaces[currentMultiFace].faces.push_back(pbface0);
4203 multiFaces[currentMultiFace].bodyId = faces[startSetIndex][pbface0].bodyId;
4204 multiFaceConnection[pbface1] = currentMultiFace;
4205 multiFaces[currentMultiFace].faces.push_back(pbface1);
4206 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][pbface0].edges.size(); e++) {
4207 if(faces[startSetIndex][pbface0].edges[e].first != pbedge) {
4208 multiFaces[currentMultiFace].edges.emplace_back(faces[startSetIndex][pbface0].edges[e].first,
4209 faces[startSetIndex][pbface0].edges[e].second);
4210 } else {
4211
4213 for(
MInt ee = 0; (unsigned)ee < faces[startSetIndex][pbface1].edges.size(); ee++) {
4214 if(faces[startSetIndex][pbface1].edges[ee].first == pbedge) {
4215 edgeIndex = ee;
4216 break;
4217 }
4218 }
4219 for(
MInt ee = 1; (unsigned)ee < faces[startSetIndex][pbface1].edges.size(); ee++) {
4220 MInt eee = (edgeIndex + ee) % faces[startSetIndex][pbface1].edges.size();
4221 multiFaces[currentMultiFace].edges.emplace_back(faces[startSetIndex][pbface1].edges[eee].first,
4222 faces[startSetIndex][pbface1].edges[eee].second);
4223 }
4224 }
4225 }
4226 somethingChanged = true;
4227
4228 while(somethingChanged) {
4229 somethingChanged = false;
4230 auto it = pureBodyEdges.begin();
4231 while(it != pureBodyEdges.end()) {
4232
4234 if(edgeCutCellPointer[edge] != c) {
4235 it++;
4236 continue;
4237 }
4238 MInt face0 = edges[startSetIndex][edge].face[0];
4239 MInt face1 = edges[startSetIndex][edge].face[1];
4240 if(multiFaceConnection[face0] == currentMultiFace && multiFaceConnection[face1] == currentMultiFace) {
4241 it = pureBodyEdges.erase(it);
4242 } else if(multiFaceConnection[face0] == currentMultiFace) {
4243 it = pureBodyEdges.erase(it);
4244 MBool mayBeConnected =
true;
4245 for(
MInt nf = 0; (unsigned)nf < multiFaces[currentMultiFace].faces.size(); nf++) {
4246 MInt nFace = multiFaces[currentMultiFace].faces[nf];
4248 mayBeConnected = false;
4249 break;
4250 }
4251 }
4252 if(!mayBeConnected) continue;
4253 multiFaceConnection[face1] = currentMultiFace;
4254 multiFaces[currentMultiFace].faces.push_back(face1);
4255 auto it2 = multiFaces[currentMultiFace].edges.begin();
4256 for(it2 = multiFaces[currentMultiFace].edges.begin();
4257 it2 != multiFaces[currentMultiFace].edges.end();
4258 it2++) {
4259 MInt edge2 = (*it2).first;
4260 if(edge2 == edge) {
4261 break;
4262 }
4263 }
4264
4266 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face1].edges.size(); e++) {
4267 if(faces[startSetIndex][face1].edges[e].first == edge) {
4268 edgeIndex = e;
4269 break;
4270 }
4271 }
4272 for(
MInt e = 1; (unsigned)e < faces[startSetIndex][face1].edges.size(); e++) {
4273 MInt ee = (edgeIndex + e) % faces[startSetIndex][face1].edges.size();
4274 multiFaces[currentMultiFace].edges.insert(
4275 it2, make_pair(faces[startSetIndex][face1].edges[ee].first,
4276 faces[startSetIndex][face1].edges[ee].second));
4277 }
4278 multiFaces[currentMultiFace].edges.erase(it2);
4279 somethingChanged = true;
4280 } else if(multiFaceConnection[face1] == currentMultiFace) {
4281 it = pureBodyEdges.erase(it);
4282 MBool mayBeConnected =
true;
4283 for(
MInt nf = 0; (unsigned)nf < multiFaces[currentMultiFace].faces.size(); nf++) {
4284 MInt nFace = multiFaces[currentMultiFace].faces[nf];
4286 mayBeConnected = false;
4287 break;
4288 }
4289 }
4290 if(!mayBeConnected) continue;
4291 multiFaceConnection[face0] = currentMultiFace;
4292 multiFaces[currentMultiFace].faces.push_back(face0);
4293 auto it2 = multiFaces[currentMultiFace].edges.begin();
4294 for(it2 = multiFaces[currentMultiFace].edges.begin();
4295 it2 != multiFaces[currentMultiFace].edges.end();
4296 it2++) {
4297 MInt edge2 = (*it2).first;
4298 if(edge2 == edge) {
4299 break;
4300 }
4301 }
4302
4304 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face0].edges.size(); e++) {
4305 if(faces[startSetIndex][face0].edges[e].first == edge) {
4306 edgeIndex = e;
4307 break;
4308 }
4309 }
4310 for(
MInt e = 1; (unsigned)e < faces[startSetIndex][face0].edges.size(); e++) {
4311 MInt ee = (edgeIndex + e) % faces[startSetIndex][face0].edges.size();
4312 multiFaces[currentMultiFace].edges.insert(
4313 it2, make_pair(faces[startSetIndex][face0].edges[ee].first,
4314 faces[startSetIndex][face0].edges[ee].second));
4315 }
4316 multiFaces[currentMultiFace].edges.erase(it2);
4317 somethingChanged = true;
4318 } else {
4319 it++;
4320 }
4321 }
4322 }
4323 currentMultiFace++;
4324 }
4325 } break;
4326
4327 default:
4328 mTerm(1, AT_,
"ERROR: invalid bodyFaceJoinMode specified. exiting...");
4329 break;
4330 }
4331
4332
4333 for(
MInt mf = 0; (unsigned)mf < multiFaces.size(); mf++) {
4334 MFloat normal[3] = {F0, F0, F0};
4335 MFloat center[3] = {F0, F0, F0};
4337 for(
MInt f = 0; (unsigned)f < multiFaces[mf].faces.size(); f++) {
4338 const MInt face = multiFaces[mf].faces[f];
4339 MFloat faceArea = faces[startSetIndex][face].area;
4340
4341 if(faceArea <
m_eps) {
4343 }
4344 area2 += faceArea;
4346 normal[i] += faces[startSetIndex][face].normal[i] * faceArea;
4347 center[i] += faces[startSetIndex][face].center[i] * faceArea;
4348 }
4349 }
4352 absNorm += normal[i] * normal[i];
4353 }
4354 absNorm = sqrt(absNorm);
4356 multiFaces[mf].center[i] = center[i] / area2;
4357 ASSERT(!std::isnan(multiFaces[mf].center[i]), "");
4358 multiFaces[mf].normal[i] = normal[i] / absNorm;
4359 ASSERT(!std::isnan(multiFaces[mf].normal[i]), "");
4360 }
4361 multiFaces[mf].area = absNorm;
4362 }
4363
4364
4365 if(!multiFaces.empty()) {
4366
4368 for(
MInt f = 0; (unsigned)f < cutCells[c].faces.size(); f++) {
4369 MInt face = cutCells[c].faces[f];
4370 if(multiFaceConnection[face] == -1) tmp_faces[noFaces++] = cutCells[c].faces[f];
4371 }
4372 cutCells[c].faces.clear();
4373 for(
MInt f = 0; f < noFaces; f++) {
4374 cutCells[c].faces.push_back(tmp_faces[f]);
4375 }
4376 }
4377
4378
4379 for(
MInt mf = 0; (unsigned)mf < multiFaces.size(); mf++) {
4380
4381
4382 const MInt newFace = faces[startSetIndex].size();
4383 const MInt faceId = faces[startSetIndex][multiFaces[mf].faces[0]].faceId;
4384 const MInt faceType = faces[startSetIndex][multiFaces[mf].faces[0]].faceType;
4385 const MInt bodyId = faces[startSetIndex][multiFaces[mf].faces[0]].bodyId;
4386
4387 faces[startSetIndex].emplace_back(faceId, faceType, bodyId);
4388 faces[startSetIndex][newFace].area = multiFaces[mf].area;
4390 faces[startSetIndex][newFace].center[i] = multiFaces[mf].center[i];
4391 faces[startSetIndex][newFace].normal[i] = multiFaces[mf].normal[i];
4392 }
4393 faces[startSetIndex][newFace].w = F0;
4394 faces[startSetIndex][newFace].cutCell = c;
4395
4396
4397
4398 MBool somethingChanged =
true;
4399 while(somethingChanged) {
4400 somethingChanged = false;
4401
4402
4403 auto it = multiFaces[mf].edges.begin();
4404 while(it != multiFaces[mf].edges.end()) {
4405 auto itLast = multiFaces[mf].edges.end();
4406 if(it == multiFaces[mf].edges.begin()) {
4407 itLast = multiFaces[mf].edges.end();
4408 } else {
4409 itLast = it;
4410 }
4411 itLast--;
4412 ASSERT(itLast != it, "");
4413 if((it->first == itLast->first) && (it->second == -(itLast->second))) {
4414 multiFaces[mf].edges.erase(itLast);
4415 it = multiFaces[mf].edges.erase(it);
4416 somethingChanged = true;
4417
4418 } else {
4419 it++;
4420 }
4421 }
4422 }
4423
4424
4425 for(auto& edge : multiFaces[mf].edges) {
4426 faces[startSetIndex][newFace].edges.emplace_back(edge.first, edge.second);
4427 MInt vertexId = edges[startSetIndex][edge.first].vertices[0];
4428 if(edge.second == -1) vertexId = edges[startSetIndex][edge.first].vertices[1];
4429 faces[startSetIndex][newFace].vertices.push_back(vertexId);
4430 }
4431 cutCells[c].faces.push_back(newFace);
4432 }
4433 }
4434 }
4435 RECORD_TIMER_STOP(tCutFace_5b);
4436
4437
4438
4439
4440
4441
4442
4443
4444 for(
MInt sc = 0; sc < (
MInt)cutCells.size(); sc++) {
4445 if(!cutCells[sc].faces.empty()) {
4448 MInt vertexRemap[maxNoVertices];
4449 fill_n(vertexRemap, maxNoVertices, -1);
4450 for(
MInt f = 0; (unsigned)f < cutCells[sc].faces.size(); f++) {
4451 MInt face = cutCells[sc].faces[f];
4452 noEdges += faces[startSetIndex][face].vertices.size();
4453 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face].vertices.size(); e++) {
4454 const MInt vertex = faces[startSetIndex][face].vertices[e];
4455 if(vertexRemap[vertex] < 0) {
4456 vertexRemap[vertex] = pcnt;
4457 pcnt++;
4458 }
4459 }
4460 }
4461
4462 if(pcnt + cutCells[sc].faces.size() != 2 + noEdges / 2 || noEdges % 2 == 1) {
4463 cerr <<
grid().domainId() <<
"Euler's polyhedral formula (1) not fulfilled ("
4464 << pcnt + (signed)cutCells[sc].faces.size() - 2 - noEdges / 2 << ") " << pcnt << " "
4465 << cutCells[sc].faces.size() <<
" " << noEdges / 2 <<
" " << noEdges <<
" for cell " <<
cellId <<
" "
4466 <<
globalTimeStep <<
" " << faces[startSetIndex].size() <<
" " << cutCells.size() <<
" "
4467 <<
grid().tree().globalId(cellId) <<
" " << sc << endl;
4469
4470 if(sc >= 1) {
4471 cerr << "Deleting splitchilds " << endl;
4472 cutCells.erase(cutCells.begin() + sc);
4473 }
4474 }
4475 }
4476 }
4477
4478
4479
4480
4481
4482
4483
4484
4485
4486
4487 MInt splitCellId = -1;
4488 MInt noSplitChildren = cutCells.size();
4489 if(noSplitChildren > CC::maxSplitCells) {
4490 mTerm(1, AT_,
"Too many split cells generated, see CutCell::maxSplitCells");
4491 }
4492 cutCellData[cutc].noSplitChilds = 0;
4493
4494
4495 if(noSplitChildren > 1) {
4496
4497
4498
4499
4500 for(
MInt sc = 0; sc < noSplitChildren; sc++) {
4501 MUint splitcutc = cutCellData.size();
4502 cutCellData[cutc].splitChildIds[cutCellData[cutc].noSplitChilds++] = splitcutc;
4503 cutCellData.emplace_back();
4504
4505
4506 cutCellData[splitcutc].cellId = -1;
4507
4508
4509 cutCellData[splitcutc].splitParentId = cutc;
4510 }
4511
4512
4513
4514
4515
4516
4517
4518
4519
4520
4521
4522
4523
4524
4525
4526
4527
4528
4529
4530
4531
4532
4533
4534
4535
4536
4537
4538 }
4539
4540 ASSERT(cutCellData[cutc].noSplitChilds == 0 || cutCellData[cutc].noSplitChilds > 1, "");
4541
4542 if(noSplitChildren > 1) {
4543 for(
MInt sc = 0; sc < noSplitChildren; sc++) {
4544 splitCellList[sc] = cutCellData[cutc].splitChildIds[sc];
4545 }
4546 } else {
4547 splitCellList[0] = cutc;
4548 }
4549
4550
4551
4552 MInt noFacesPerCartesianFaceAll[6] = {0, 0, 0, 0, 0, 0};
4553 MInt noBodySurfacesAll = 0;
4554
4555
4556
4557 for(
MInt sc = 0; sc < noSplitChildren; sc++) {
4558 RECORD_TIMER_START(tCutFace_6);
4559
4560
4561
4562
4563
4564
4565
4566
4567
4568
4569
4570 MInt cutCellId = splitCellList[sc];
4571
4572
4573
4574
4575
4576 if(cutCellData[cutCellId].splitParentId > -1) {
4578 cutCellData[cutCellId].cornerIsInsideGeometry[0][corner] = true;
4579 }
4580 for(
MInt f = 0; (unsigned)f < cutCells[sc].faces.size(); f++) {
4581 MInt face = cutCells[sc].faces[f];
4582 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face].vertices.size(); e++) {
4583 const MInt vStart = faces[startSetIndex][face].vertices[e];
4584 if(vertices[startSetIndex][vStart].
pointType == 0) {
4585 MInt corner = vertices[startSetIndex][vStart].pointId;
4586 cutCellData[cutCellId].cornerIsInsideGeometry[0][corner] = false;
4587 }
4588 }
4589 }
4590 }
4591
4592
4593
4594
4595 MInt noFacesPerCartesianFace[6] = {0, 0, 0, 0, 0, 0};
4596 MInt noBodySurfaces = 0;
4597 for(
MInt f = 0; (unsigned)f < cutCells[sc].faces.size(); f++) {
4598 MInt face = cutCells[sc].faces[f];
4599 if(faces[startSetIndex][face].faceType == 0) {
4600 noFacesPerCartesianFace[faces[startSetIndex][face].faceId]++;
4601 noFacesPerCartesianFaceAll[faces[startSetIndex][face].faceId]++;
4602 } else {
4603 noBodySurfaces++;
4604 noBodySurfacesAll++;
4605 }
4606 }
4607
4608
4609
4610
4611
4612
4613
4614
4615
4616
4617
4618
4619
4620
4621
4622
4623
4624
4625
4626
4627
4628
4629
4630
4631
4632
4633
4634
4635
4636
4637
4638
4639
4640
4641
4642
4643
4644
4645
4646
4647 if(noBodySurfaces > maxNoSurfaces || noBodySurfaces > CC::maxNoBoundarySurfaces) {
4648 cerr << "more than FvBndryCell<nDim>::m_maxNoSurfaces cut surfaces detected for a cell. This is not "
4649 "implemented in MB framework yet. cell "
4651 << noBodySurfaces << " maximum " << maxNoSurfaces << " " << CC::maxNoBoundarySurfaces << " "
4652 << noSplitChildren << " " << faces[startSetIndex].size() << endl;
4653 writeVTKFileOfCell(cellId, &faces[startSetIndex], &vertices[startSetIndex], startSetIndex);
4655 " more than 3 cut surfaces detected for a cell. This is not implemented in MB framework yet. "
4656 "exiting...");
4657 }
4658
4659
4660
4661
4662
4663
4664 if(cutCells[sc].faces.empty()) {
4665
4666
4667
4668
4669
4670
4671
4672 cutCellData[cutCellId].volume = F0;
4673
4674
4675 }
4678 }
4679
4680
4681
4682
4684 if(noFacesPerCartesianFace[f] == 0) {
4685
4686 cutCellData[cutCellId].externalFaces[f] = true;
4687
4688
4689
4690 if(splitCellId < 0) {
4691
4692
4693
4694
4695
4696
4697
4698
4699
4700
4701
4702
4703
4704
4705
4706
4707
4708
4709
4710
4711
4712
4713
4714
4715
4716
4717
4718
4719
4720
4721
4722
4723 }
4724 }
4725 }
4726
4727
4728
4729
4730 ASSERT(vertices[startSetIndex].size() <= (unsigned)maxNoVertices, "");
4731 for(
MInt cp = 0; cp < maxNoVertices; cp++)
4732 newCutPointId[cp] = -1;
4733 for(
MInt i = 0; (unsigned)i < vertices[startSetIndex].size(); i++)
4734 isPartOfBodySurface[i] = false;
4736
4737 pointEdgeId[e] = cutCellData[cutc].cutEdges[e];
4738 }
4739 noBodySurfaces = 0;
4740
4741 ASSERT(cutCells[sc].volume >= F0, "");
4742 ASSERT(!(std::isnan(cutCells[sc].volume)), "");
4743 ASSERT(!(std::isnan(cutCells[sc].center[sc])), "");
4744 ASSERT(!(std::isnan(cutCells[sc].center[1])), "");
4745 ASSERT(!(std::isnan(cutCells[sc].center[2])), "");
4746 ASSERT(!(std::isinf(cutCells[sc].volume)), "");
4747 ASSERT(!(std::isinf(cutCells[sc].center[sc])), "");
4748 ASSERT(!(std::isinf(cutCells[sc].center[1])), "");
4749 ASSERT(!(std::isinf(cutCells[sc].center[2])), "");
4750
4751
4752 cutCellData[cutCellId].volume = cutCells[sc].volume;
4753
4754
4755
4756
4757
4759
4760 cutCellData[cutCellId].volumetricCentroid[i] = cutCells[sc].center[i] -
a_coordinate(cellId, i);
4761 }
4762
4763
4764 cutCellData[cutCellId].noCartesianSurfaces = 0;
4765 cutCellData[cutCellId].noBoundarySurfaces = 0;
4766
4767
4768
4769 for(
MInt f = 0; (unsigned)f < cutCells[sc].faces.size(); f++) {
4770 MInt face = cutCells[sc].faces[f];
4771 ASSERT(faces[startSetIndex][face].area >= F0, "");
4772 ASSERT(!(std::isnan(faces[startSetIndex][face].area)), "");
4773 ASSERT(!(std::isnan(faces[startSetIndex][face].center[0])), "");
4774 ASSERT(!(std::isnan(faces[startSetIndex][face].center[1])), "");
4775 ASSERT(!(std::isnan(faces[startSetIndex][face].center[2])), "");
4776 ASSERT(!(std::isinf(faces[startSetIndex][face].area)), "");
4777 ASSERT(!(std::isinf(faces[startSetIndex][face].center[0])), "");
4778 ASSERT(!(std::isinf(faces[startSetIndex][face].center[1])), "");
4779 ASSERT(!(std::isinf(faces[startSetIndex][face].center[2])), "");
4780
4781
4782
4783
4784
4785 if(faces[startSetIndex][face].faceType == 1) {
4786
4787
4788 cutCellData[cutCellId].boundarySurfaceArea[noBodySurfaces] = faces[startSetIndex][face].area;
4789
4791
4792
4793 cutCellData[cutCellId].boundarySurfaceCentroid[noBodySurfaces][i] = faces[startSetIndex][face].center[i];
4794 cutCellData[cutCellId].boundarySurfaceNormal[noBodySurfaces][i] = faces[startSetIndex][face].normal[i];
4795 }
4796
4797
4798
4799 cutCellData[cutCellId].boundarySurfaceBndryCndId[noBodySurfaces] = -1;
4800 cutCellData[cutCellId].boundarySurfaceBodyId[noBodySurfaces] = faces[startSetIndex][face].bodyId;
4801
4802
4803 for(
MInt v = 0; (unsigned)v < faces[startSetIndex][face].vertices.size(); v++) {
4804 const MInt vertex = faces[startSetIndex][face].vertices[v];
4805
4806 isPartOfBodySurface[vertex] = true;
4807
4808 newCutPointId[vertex] = v;
4809 vertices[startSetIndex][vertex].cartSrfcId = noBodySurfaces;
4810
4811
4812
4813
4814
4815
4816
4817
4818
4819
4820
4821
4822
4823
4824
4825
4826
4827
4828
4829
4830
4831
4832
4833
4834 }
4835 noBodySurfaces++;
4836 cutCellData[cutCellId].noBoundarySurfaces++;
4837 ASSERT(cutCellData[cutCellId].noBoundarySurfaces == noBodySurfaces, "");
4838 } else {
4839
4840
4841
4842
4843
4844 const MInt dirId = faces[startSetIndex][face].faceId;
4845 const MInt spaceId = dirId / 2;
4846
4847
4848
4849
4850
4851
4852
4853
4854
4855
4856
4857
4858
4859
4860
4861
4862
4863
4864
4865
4866
4867
4868
4869
4870
4871
4872
4873
4874
4875
4876
4877
4878
4879
4880
4881
4882
4883
4884
4885
4886
4887
4888
4889
4890
4891
4892
4893
4894
4895
4896
4897
4898
4899
4900
4902 if(dirId % 2 == 0) sign = -F1;
4903
4904
4906
4907
4908 cutCellData[cutCellId].cartFaceCentroid[cutCellData[cutCellId].noCartesianSurfaces][i] =
4909 faces[startSetIndex][face].center[i];
4910 }
4911
4912
4913 cutCellData[cutCellId].cartFaceCentroid[cutCellData[cutCellId].noCartesianSurfaces][spaceId] =
4915
4916
4917 cutCellData[cutCellId].cartFaceArea[cutCellData[cutCellId].noCartesianSurfaces] =
4918 faces[startSetIndex][face].area;
4919 cutCellData[cutCellId].cartFaceDir[cutCellData[cutCellId].noCartesianSurfaces] =
4920 faces[startSetIndex][face].faceId;
4921
4922
4923
4924
4925
4926
4927
4928
4929
4930
4931
4932
4933
4934
4935
4936
4937
4938
4939
4940
4941
4942
4943
4944
4945
4946
4947
4948
4949
4950
4951
4952
4953
4954
4955
4956
4957
4958
4959
4960
4961
4962
4963
4964
4965
4966
4967
4968
4969
4970
4971
4972
4973
4974
4975
4976
4977
4978
4979
4980
4981 cutCellData[cutCellId].noCartesianSurfaces++;
4982 if(cutCellData[cutCellId].noCartesianSurfaces > CC::maxNoCartesianSurfaces) {
4983 mTerm(1, AT_,
"Increase CutCell::maxNoCartesianSurfaces");
4984 }
4985 }
4986 }
4987
4988 if(noSplitChildren > 1) {
4989
4990 cutCellData[cutc].noBoundarySurfaces = 0;
4991 }
4992
4993 cutCellData[cutCellId].noBoundarySurfaces = noBodySurfaces;
4994
4995 RECORD_TIMER_STOP(tCutFace_6);
4996
4997 for(
MInt k = 0; k < CC::noFaces; k++) {
4998 cutCellData[cutCellId].noFacesPerCartesianDir[k] = noFacesPerCartesianFace[k];
4999 }
5000
5001
5002
5003 RECORD_TIMER_START(tCutFace_7);
5004
5005
5006
5007 cutCellData[cutCellId].noTotalFaces = 0;
5008 cutCellData[cutCellId].noAdditionalVertices = 0;
5009 MInt vertexMap[maxNoVertices];
5010 fill_n(vertexMap, maxNoVertices, -1);
5011
5012 for(
MInt f = 0; (unsigned)f < cutCells[sc].faces.size(); f++) {
5013 MInt face = cutCells[sc].faces[f];
5014
5015
5016 const MInt facecnt = cutCellData[cutCellId].noTotalFaces;
5017 cutCellData[cutCellId].allFacesNoPoints[facecnt] = 0;
5018 cutCellData[cutCellId].allFacesBodyId[facecnt] = faces[startSetIndex][face].bodyId;
5019
5020 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face].vertices.size(); e++) {
5021 const MInt vertex = faces[startSetIndex][face].vertices[e];
5022 MInt pointId = vertices[startSetIndex][vertex].pointId;
5023 MInt storePointId = -1;
5024 if(vertices[startSetIndex][vertex].
pointType == 0) {
5025
5026
5027
5028 ASSERT(pointId < CC::noCorners, "");
5029 storePointId = pointId;
5030 } else if(pointId > -1) {
5031 ASSERT(vertices[startSetIndex][vertex].
pointType == 1,
"");
5032
5033
5034
5035
5036
5037
5038
5039
5040
5041
5042
5043 ASSERT(pointId < cutCellData[cutc].noCutPoints, "");
5044 storePointId = CC::noCorners + pointId;
5045 } else {
5046 ASSERT(vertices[startSetIndex][vertex].
pointType > 1,
"");
5047 if(vertexMap[vertex] < 0) {
5049 cutCellData[cutCellId].additionalVertices[cutCellData[cutCellId].noAdditionalVertices][i] =
5050 vertices[startSetIndex][vertex].coordinates[i];
5051 }
5052 vertexMap[vertex] = cutCellData[cutCellId].noAdditionalVertices;
5053 cutCellData[cutCellId].noAdditionalVertices++;
5054 if(cutCellData[cutCellId].noAdditionalVertices > CC::maxNoAdditionalVertices) {
5056 }
5057 ASSERT(cutCellData[cutCellId].noAdditionalVertices <= CC::maxNoAdditionalVertices, "");
5058 }
5059 storePointId = CC::noCorners + CC::maxNoCutPoints + vertexMap[vertex];
5060 }
5061 cutCellData[cutCellId].allFacesPointIds[facecnt][cutCellData[cutCellId].allFacesNoPoints[facecnt]] =
5062 storePointId;
5063 cutCellData[cutCellId].allFacesNoPoints[facecnt]++;
5064 ASSERT(cutCellData[cutCellId].allFacesNoPoints[facecnt] <= CC::maxNoFaceVertices,
5065 "has: " + to_string(cutCellData[cutCellId].allFacesNoPoints[facecnt]) + " max is "
5066 + to_string(CC::maxNoFaceVertices));
5067 }
5068
5069 cutCellData[cutCellId].noTotalFaces++;
5070 }
5071
5072
5073
5074
5075
5076
5077
5078
5079
5080
5081
5082
5083
5084
5085
5086
5087
5088
5089
5090
5091
5092
5093
5094
5095
5096
5097
5098
5099
5100
5101
5102
5103
5104
5105
5106
5107
5108
5109
5110
5111
5112
5113
5114
5115
5116
5117
5118
5119
5120
5121
5122
5123
5124
5125
5126
5127
5128
5129
5130
5131
5132
5133
5134
5135
5136
5137
5138
5139
5140
5141
5142
5143
5144
5145
5146
5147
5148
5149
5150
5151
5152
5153
5154
5155
5156
5157
5158
5159
5160
5161
5162
5163
5164
5165
5166
5167
5168
5169
5170
5171
5172
5173
5174
5175
5176
5177
5178
5179
5180
5181
5182
5183
5184
5185
5186
5187
5188
5189
5190
5191
5192
5193
5194
5195
5196
5197
5198 if(!cutCells[sc].faces.empty()) {
5201 MInt vertexRemap[maxNoVertices];
5202 fill_n(vertexRemap, maxNoVertices, -1);
5203 for(
MInt f = 0; (unsigned)f < cutCells[sc].faces.size(); f++) {
5204 MInt face = cutCells[sc].faces[f];
5205
5206 noEdges += faces[startSetIndex][face].vertices.size();
5207 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face].vertices.size(); e++) {
5208 const MInt vertex = faces[startSetIndex][face].vertices[e];
5209 if(vertexRemap[vertex] < 0) {
5210 vertexRemap[vertex] = pcnt;
5211 pcnt++;
5212 }
5213 }
5214 }
5215 if(pcnt + cutCells[sc].faces.size() != 2 + noEdges / 2 || noEdges % 2 == 1) {
5216 cerr <<
grid().domainId() <<
"Euler's polyhedral formula not fulfilled ("
5217 << pcnt + (signed)cutCells[sc].faces.size() - 2 - noEdges / 2 << ") " << pcnt << " "
5218 << cutCells[sc].faces.size() <<
" " << noEdges / 2 <<
" " << noEdges <<
" for cell " <<
cellId <<
" "
5219 <<
globalTimeStep <<
" " << faces[startSetIndex].size() <<
" " << cutCells.size() <<
" "
5220 <<
grid().tree().globalId(cellId) <<
" " << sc << endl;
5222 }
5223 }
5224 RECORD_TIMER_STOP(tCutFace_7);
5225 }
5226
5227
5228 if(noSplitChildren > 1) {
5229 cutCellData[cutc].volume = F0;
5230 for(
MInt sc = 0; sc < noSplitChildren; sc++) {
5231 MInt cutCellId = splitCellList[sc];
5232 cutCellData[cutc].volume += cutCellData[cutCellId].volume;
5233 }
5235 cutCellData[cutc].noFacesPerCartesianDir[dir] = noFacesPerCartesianFaceAll[dir];
5236
5237 if(cutCellData[cutc].noFacesPerCartesianDir[dir] == 0) {
5238 cutCellData[cutc].externalFaces[dir] = true;
5239 }
5240 }
5241 }
5242
5243
5244
5245
5246
5247
5248
5249
5250
5251
5252
5253
5254
5255
5256
5257
5258
5259
5260
5261
5262
5263
5264
5265
5266
5267
5268
5269
5270
5271
5272
5273
5274
5275
5276
5277
5278
5279
5280
5281
5282
5283
5284
5285
5286
5287
5288
5289
5290
5291
5292
5293
5294
5295
5296
5297
5298
5299
5300
5301
5302
5303
5304
5305
5306
5307
5308
5309
5310
5311
5312
5313
5314
5315
5316
5317
5318
5319
5320
5321
5322
5323
5324
5325
5326
5327
5328
5329
5330
5331
5332
5333
5334
5335
5336
5337
5338
5339
5340
5341
5342
5343
5344
5345
5346
5347
5348
5349
5350
5351
5352
5353
5354
5355
5356
5357
5358
5359
5360
5361
5362
5363
5364
5365
5366
5367
5368
5369
5370
5371
5372
5373
5374
5375
5376
5377
5378
5379
5380
5381
5382
5383
5384
5385
5386
5387
5388
5389
5390
5391
5392
5393
5394
5395
5396
5397
5398
5399
5400
5401
5402
5403
5404
5405
5406
5407
5408
5409
5410
5411
5412
5413
5414
5415
5416
5417
5418
5419
5420
5421
5422
5423
5424
5425
5426
5427
5428
5429
5430
5431
5432
5433
5434
5435
5436
5437
5438
5439
5440
5441
5442
5443
5444
5445
5446
5447
5448
5449
5450
5451
5452#ifdef CutCell_DEBUG
5453 if(
globalTimeStep == debugTimeStep && cellId == debugCellId &&
grid().domainId() == debugDomainId) {
5455
5456 const MChar* fileName =
"cell-Info_";
5457 stringstream fileName2;
5458 fileName2 << fileName <<
grid().tree().globalId(cellId) <<
".txt";
5459
5460 ofstream ofl;
5461 ofl.open((fileName2.str()).c_str(), ofstream::trunc);
5462
5463 if(ofl) {
5464 ofl.setf(ios::fixed);
5465 ofl.precision(12);
5466
5467 ofl <<
"Writing Cell " <<
grid().tree().globalId(cellId) <<
" " <<
cellId <<
" " << cutc <<
" noCuts "
5468 << noIndividualCuts << " noCutCells " << cutCells.size() << " vol " << cutCellData[cutc].volume
5469 << " SplitChilds " << cutCellData[cutc].noSplitChilds << " " << cutCellData[cutc].splitParentId << endl
5470 << endl;
5472 ofl << "Dir " << dir << " external " << cutCellData[cutc].externalFaces[dir] << " faces "
5473 << cutCellData[cutc].noFacesPerCartesianDir[dir] << endl
5474 << endl;
5475 }
5476
5477 ofl << " Cartesian-Surf: " << cutCellData[cutc].noCartesianSurfaces << endl << endl;
5478
5479 for(
MInt id = 0;
id < cutCellData[cutc].noCartesianSurfaces;
id++) {
5480 ofl <<
"Id " <<
id <<
" dir " << cutCellData[cutc].cartFaceDir[
id] <<
" area "
5481 << cutCellData[cutc].cartFaceArea[
id] <<
" coord " << cutCellData[cutc].cartFaceCentroid[
id][0]
5482 << cutCellData[cutc].cartFaceCentroid[
id][1] << cutCellData[cutc].cartFaceCentroid[
id][2] << endl
5483 << endl;
5484 }
5485
5486 ofl << " Bndry-Surfaces: " << cutCellData[cutc].noBoundarySurfaces << endl << endl;
5487
5488 for(
MInt id = 0;
id < cutCellData[cutc].noBoundarySurfaces;
id++) {
5489 ofl <<
"Id " <<
id <<
" normal " << cutCellData[cutc].boundarySurfaceNormal[
id][0]
5490 << cutCellData[cutc].boundarySurfaceNormal[
id][1] << cutCellData[cutc].boundarySurfaceNormal[
id][2]
5491 <<
" center " << cutCellData[cutc].boundarySurfaceCentroid[
id][0]
5492 << cutCellData[cutc].boundarySurfaceCentroid[
id][1] << cutCellData[cutc].boundarySurfaceCentroid[
id][2]
5493 <<
" area " << cutCellData[cutc].boundarySurfaceArea[
id] << endl
5494 << endl;
5495 }
5496 }
5497 }
5498#endif
5499 }
5500
5501
5502#ifdef CutCell_DEBUG
5504#endif
5505}
static constexpr MInt m_noCorners
static constexpr MInt m_noEdges
void compVolumeIntegrals_pyraBased3(std::vector< polyCutCell > *, std::vector< polyFace > *, const std::vector< polyVertex > *)
void writeVTKFileOfCell(MInt cellId, std::vector< polyFace > *faces, const std::vector< polyVertex > *vertices, MInt set)
MFloat & a_coordinate(const MInt cellId, const MInt dir)
Returns the coordinate of the cell cellId for dimension dir.
void addPoint(const std::vector< polyVertex > *, const MInt *, const MInt, std::array< MFloat, nDim >)
adds an additional MC point in the cell sums up all n points passed in indices. points have to be sto...
MInt m_noLevelSetsUsedForMb
void crossProduct(MFloat *c, const MFloat *a, const MFloat *b)
MFloat a_gridCellVolume(const MInt cellId) const
Returns the volume of the cell cellId.
static constexpr MInt m_noDirs
std::enable_if< nDim_==2, T >::type computeNormal(const std::array< MFloat, nDim > p0, const std::array< MFloat, nDim > p1, std::array< MFloat, nDim > res, MFloat &w)
returns the normal corresponding to the triangle abc and returns the result in res
MFloat a_cellLengthAtCell(const MInt cellId) const
Returns the cellLength of the cell cellId.
MBool computeCutFaceSimple(CutCell< nDim > &cutCell)
This class is a ScratchSpace.
void mTerm(const MInt errorCode, const MString &location, const MString &message)
constexpr Real POW2(const Real x)
constexpr T mMin(const T &x, const T &y)
constexpr T mMax(const T &x, const T &y)
constexpr MLong IPOW2(MInt x)
constexpr std::underlying_type< FcCell >::type p(const FcCell property)
Converts property name to underlying integer value.