74 return s < ROOTVSMALL ?
Zero : (a ^
b) /
s;
84 const point& basePoint,
85 const point& rightPoint,
86 const point& leftPoint,
87 const point& centrePoint
90 const vector mid(centrePoint - basePoint);
96 0.5*(rightPoint - basePoint),
102 0.5*(leftPoint - basePoint)
124 const point& faceCentre
131 (edgeLine.
a() - faceCentre),
132 (edgeLine.
b() - faceCentre)
143 const point& faceCentre
172 dirn.removeCollinear(unitAxis1);
196 bitSet markPoints(
p.nPoints());
197 for (
const edge&
e :
p.boundaryEdges())
200 markPoints.
set(
e.first());
201 markPoints.set(
e.second());
213 void Foam::faMesh::calcLduAddressing()
const 216 <<
"Calculating addressing" <<
endl;
221 <<
"lduPtr_ already allocated" 225 lduPtr_ =
new faMeshLduAddressing(*
this);
229 void Foam::faMesh::calcPatchStarts()
const 232 <<
"Calculating patch starts" <<
endl;
237 <<
"patchStarts already allocated" 245 void Foam::faMesh::calcWhichPatchFaces()
const 248 if (polyPatchFacesPtr_ || polyPatchIdsPtr_)
251 <<
"Already allocated polyPatchFaces/polyPatchIds" 257 polyPatchFacesPtr_.reset
265 for (
const labelPair& tup : *polyPatchFacesPtr_)
276 bitOrOp<labelHashSet>(),
281 polyPatchIdsPtr_.reset(
new labelList(ids.sortedToc()));
296 const point& faceCentre,
301 const vector centreToEdge(edgeLine.
centre() - faceCentre);
303 vector leVector(edgeLine.
vec() ^ edgeNormal);
305 scalar
s(
mag(leVector));
313 leVector = centreToEdge;
314 leVector.removeCollinear(edgeLine.
unitVec());
323 leVector *= edgeLine.
mag()/
s;
328 leVector *= edgeLine.
mag()/
s *
sign(leVector & centreToEdge);
348 const vector minVector(vector::uniform(0.57735*SMALL));
349 const scalar minLenSqr(SMALL*SMALL);
353 if (v.magSqr() < minLenSqr)
369 auto& edgeNormals = tedgeNormals.ref();
387 <<
"Using geometryOrder < 0 : " 388 "simplified edge area-normals, without processor connectivity" 392 for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
394 const linePointRef edgeLine(edges_[edgei].line(localPoints));
401 fCentres[edgeOwner()[edgei]]
407 fCentres[edgeNeighbour()[edgei]]
413 for (label edgei = nInternalEdges_; edgei < nEdges_; ++edgei)
415 const linePointRef edgeLine(edges_[edgei].line(localPoints));
422 fCentres[edgeOwner()[edgei]]
435 <<
"Using geometryOrder == 0 : " 436 "weighted edge normals, without processor connectivity" 442 for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
444 const linePointRef edgeLine(edges_[edgei].line(localPoints));
451 fCentres[edgeOwner()[edgei]]
457 fCentres[edgeNeighbour()[edgei]]
463 for (label edgei = nInternalEdges_; edgei < nEdges_; ++edgei)
465 const linePointRef edgeLine(edges_[edgei].line(localPoints));
472 fCentres[edgeOwner()[edgei]]
480 bitSet nbrBoundaryAdjust;
485 DynamicList<UPstream::Request> requests
490 PtrList<vectorField> nbrEdgeNormals(
boundary().size());
497 label patchOffset = nInternalEdges_;
501 const faPatch& fap =
boundary()[patchi];
503 auto edgeNorms = edgeNormals.
slice(patchOffset, fap.nEdges());
504 patchOffset += edgeNorms.size();
506 const auto* fapp = isA<processorFaPatch>(fap);
510 const label nbrProci = fapp->neighbProcNo();
514 nbrEdgeNormals.emplace(patchi, edgeNorms.size());
519 requests.emplace_back(),
521 nbrEdgeNormals[patchi]
527 requests.emplace_back(),
533 else if (isA<wedgeFaPatch>(fap))
536 const auto& wedgePatch = refCast<const wedgeFaPatch>(fap);
543 wedgePatch.centreNormal()
549 e.removeCollinear(wedgeNorm);
557 else if (correctPatchPointNormals(patchi) && !fap.coupled())
560 nbrBoundaryAdjust.set(patchi);
567 if (!requests.empty())
571 patchOffset = nInternalEdges_;
575 const faPatch& fap =
boundary()[patchi];
577 auto edgeNorms = edgeNormals.
slice(patchOffset, fap.nEdges());
578 patchOffset += edgeNorms.size();
580 if (nbrEdgeNormals.set(patchi))
582 const vectorField& nbrNorms = nbrEdgeNormals[patchi];
584 forAll(edgeNorms, patchEdgei)
586 edgeNorms[patchEdgei] += nbrNorms[patchEdgei];
598 && hasPatchPointNormalsCorrection()
603 <<
"Apply " << nbrBoundaryAdjust.count()
604 <<
" boundary neighbour corrections" <<
nl;
608 (void)this->haloFaceNormals();
611 label patchOffset = nInternalEdges_;
613 for (
const label patchi : nbrBoundaryAdjust)
615 const faPatch& fap =
boundary()[patchi];
617 auto edgeNorms = edgeNormals.
slice(patchOffset, fap.nEdges());
618 patchOffset += edgeNorms.size();
620 if (fap.ngbPolyPatchIndex() < 0)
623 <<
"Neighbour polyPatch index is not defined " 624 <<
"for faPatch " << fap.name()
638 vectorField nbrNorms(this->haloFaceNormals(patchi));
640 const label nPatchEdges =
min(edgeNorms.size(), nbrNorms.size());
642 for (label patchEdgei = 0; patchEdgei < nPatchEdges; ++patchEdgei)
644 edgeNorms[patchEdgei].removeCollinear(nbrNorms[patchEdgei]);
652 forAll(edgeNormals, edgei)
654 const linePointRef edgeLine(edges_[edgei].line(localPoints));
656 edgeNormals[edgei].removeCollinear(edgeLine.
unitVec());
657 edgeNormals[edgei].normalise();
666 void Foam::faMesh::calcLe()
const 669 <<
"Calculating local edges" <<
endl;
674 <<
"LePtr_ already allocated" 684 mesh().pointsInstance(),
712 for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
716 fCentres[edgeOwner()[edgei]],
717 edges_[edgei].line(localPoints),
726 const faPatch& fap =
boundary()[patchi];
730 edgeNormals.boundaryField()[patchi];
732 label edgei = fap.start();
738 fCentres[edgeOwner()[edgei]],
739 edges_[edgei].line(localPoints),
740 bndEdgeNormals[patchEdgei]
764 void Foam::faMesh::calcMagLe()
const 767 <<
"Calculating local edge magnitudes" <<
endl;
772 <<
"magLe() already allocated" 782 mesh().pointsInstance(),
799 auto iter = magLe.primitiveFieldRef().begin();
801 for (
const edge&
e : internalEdges())
803 *iter =
e.mag(localPoints);
806 if (
mag(*iter) < SMALL)
817 auto& bfld = magLe.boundaryFieldRef();
821 auto iter = bfld[patchi].begin();
823 for (
const edge&
e :
boundary()[patchi].patchSlice(edges_))
825 *iter =
e.mag(localPoints);
828 if (
mag(*iter) < SMALL)
840 void Foam::faMesh::calcFaceCentres()
const 843 <<
"Calculating face centres" <<
endl;
848 <<
"areaCentres already allocated" 858 mesh().pointsInstance(),
878 if (
mesh().hasFaceCentres())
882 centres.primitiveFieldRef()
888 auto iter = centres.primitiveFieldRef().begin();
890 for (
const face&
f : faces())
892 *iter =
f.centre(localPoints);
900 auto& bfld = centres.boundaryFieldRef();
904 auto iter = bfld[patchi].begin();
906 for (
const edge&
e :
boundary()[patchi].patchSlice(edges_))
908 *iter =
e.centre(localPoints);
917 centres.boundaryFieldRef().evaluateCoupled<processorFaPatch>();
922 void Foam::faMesh::calcEdgeCentres()
const 925 <<
"Calculating edge centres" <<
endl;
930 <<
"edgeCentres already allocated" 940 mesh().pointsInstance(),
960 auto iter = centres.primitiveFieldRef().begin();
962 for (
const edge&
e : internalEdges())
964 *iter =
e.centre(localPoints);
971 auto& bfld = centres.boundaryFieldRef();
975 auto iter = bfld[patchi].begin();
977 for (
const edge&
e :
boundary()[patchi].patchSlice(edges_))
979 *iter =
e.centre(localPoints);
987 void Foam::faMesh::calcS()
const 990 <<
"Calculating areas" <<
endl;
995 <<
"S() already allocated" 999 SPtr_ =
new DimensionedField<scalar, areaMesh>
1013 auto& areas = *SPtr_;
1020 UIndirectList<vector> meshFaceAreas(
mesh().faceAreas(),
faceLabels());
1022 auto&
fld = areas.field();
1029 if (
mag(
fld[facei]) < SMALL)
1041 auto iter = areas.field().begin();
1043 for (
const face&
f : faces())
1045 *iter =
f.mag(localPoints);
1048 if (
mag(*iter) < SMALL)
1059 void Foam::faMesh::calcFaceAreaNormals()
const 1062 <<
"Calculating face area normals" <<
endl;
1064 if (faceAreaNormalsPtr_)
1067 <<
"faceAreaNormals already allocated" 1071 faceAreaNormalsPtr_ =
1077 mesh().pointsInstance(),
1094 auto&
fld = faceNormals.primitiveFieldRef();
1096 if (
mesh().hasFaceAreas())
1105 auto iter =
fld.begin();
1107 for (
const face&
f : faces())
1109 *iter =
f.areaNormal(localPoints);
1123 const auto& edgeNormalsBoundary = edgeAreaNormals().boundaryField();
1127 faceNormals.boundaryFieldRef()[patchi]
1128 = edgeNormalsBoundary[patchi];
1135 faceNormals.boundaryFieldRef().evaluateCoupled<processorFaPatch>();
1140 void Foam::faMesh::calcEdgeAreaNormals()
const 1143 <<
"Calculating edge area normals" <<
endl;
1145 if (edgeAreaNormalsPtr_)
1148 <<
"edgeAreaNormalsPtr_ already allocated" 1152 edgeAreaNormalsPtr_ =
1158 mesh().pointsInstance(),
1182 edgeAreaNormals.primitiveFieldRef()
1187 auto& bfld = edgeAreaNormals.boundaryFieldRef();
1189 label patchOffset = nInternalEdges_;
1193 const faPatch& fap =
boundary()[patchi];
1198 bfld[patchi].size(),
1202 patchOffset += fap.nEdges();
1215 const vectorField& pointNormals = pointAreaNormals();
1222 auto&
fld = edgeAreaNormals.primitiveFieldRef();
1226 const edge&
e = edges_[edgei];
1230 fld[edgei] = (pointNormals[
e.first()] + pointNormals[
e.second()]);
1232 fld[edgei].removeCollinear(edgeLine.
unitVec());
1233 fld[edgei].normalise();
1241 auto& bfld = edgeAreaNormals.boundaryFieldRef();
1245 const faPatch& fap =
boundary()[patchi];
1247 auto& pfld = bfld[patchi];
1251 forAll(bndEdges, patchEdgei)
1253 const edge&
e = bndEdges[patchEdgei];
1258 (pointNormals[
e.first()] + pointNormals[
e.second()]);
1260 pfld[patchEdgei].removeCollinear(edgeLine.
unitVec());
1261 pfld[patchEdgei].normalise();
1270 void Foam::faMesh::calcFaceCurvatures()
const 1273 <<
"Calculating face curvatures" <<
endl;
1275 if (faceCurvaturesPtr_)
1278 <<
"faceCurvaturesPtr_ already allocated" 1282 faceCurvaturesPtr_ =
1288 mesh().pointsInstance(),
1308 faceCurvatures =
sign(kN&faceAreaNormals())*
mag(kN);
1312 void Foam::faMesh::calcEdgeTransformTensors()
const 1315 <<
"Calculating edge transformation tensors" <<
endl;
1317 if (edgeTransformTensorsPtr_)
1320 <<
"edgeTransformTensorsPtr_ already allocated" 1324 edgeTransformTensorsPtr_ =
new FieldField<Field, tensor>(nEdges_);
1325 auto& edgeTransformTensors = *edgeTransformTensorsPtr_;
1328 for (label edgei = 0; edgei < nEdges_; ++edgei)
1330 edgeTransformTensors.set(edgei,
new Field<tensor>(3,
tensor::I));
1340 for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
1342 const label ownFacei = owner()[edgei];
1343 const label neiFacei = neighbour()[edgei];
1344 auto& tensors = edgeTransformTensors[edgei];
1346 vector edgeCtr = Ce.internalField()[edgei];
1350 edgeCtr -= skewCorrectionVectors().internalField()[edgei];
1357 Ne.internalField()[edgei],
1358 (edgeCtr - Cf.internalField()[ownFacei])
1365 Nf.internalField()[ownFacei],
1366 (edgeCtr - Cf.internalField()[ownFacei])
1373 Nf.internalField()[neiFacei],
1374 (Cf.internalField()[neiFacei] - edgeCtr)
1385 vectorField nbrCf(Cf.boundaryField()[patchi].patchNeighbourField());
1386 vectorField nbrNf(Nf.boundaryField()[patchi].patchNeighbourField());
1388 forAll(edgeFaces, bndEdgei)
1390 const label ownFacei = edgeFaces[bndEdgei];
1391 const label meshEdgei(
boundary()[patchi].start() + bndEdgei);
1393 auto& tensors = edgeTransformTensors[meshEdgei];
1395 vector edgeCtr = Ce.boundaryField()[patchi][bndEdgei];
1399 edgeCtr -= skewCorrectionVectors()
1400 .boundaryField()[patchi][bndEdgei];
1407 Ne.boundaryField()[patchi][bndEdgei],
1408 (edgeCtr - Cf.internalField()[ownFacei])
1415 Nf.internalField()[ownFacei],
1416 (edgeCtr - Cf.internalField()[ownFacei])
1424 (nbrCf[bndEdgei] - edgeCtr)
1430 forAll(edgeFaces, bndEdgei)
1432 const label ownFacei = edgeFaces[bndEdgei];
1433 const label meshEdgei(
boundary()[patchi].start() + bndEdgei);
1435 auto& tensors = edgeTransformTensors[meshEdgei];
1437 vector edgeCtr = Ce.boundaryField()[patchi][bndEdgei];
1441 edgeCtr -= skewCorrectionVectors()
1442 .boundaryField()[patchi][bndEdgei];
1449 Ne.boundaryField()[patchi][bndEdgei],
1450 (edgeCtr - Cf.internalField()[ownFacei])
1457 Nf.internalField()[ownFacei],
1458 (edgeCtr - Cf.internalField()[ownFacei])
1472 <<
"Calculating internal points" <<
endl;
1477 return markPoints.sortedToc();
1484 <<
"Calculating boundary points" <<
endl;
1518 void Foam::faMesh::calcPointAreaNormals(
vectorField& result)
const 1521 <<
"Calculating pointAreaNormals : face dual method" <<
endl;
1523 result.resize_nocopy(
nPoints());
1532 for (
const face&
f : faces)
1534 const label nVerts(
f.
size());
1537 for (
const label fp :
f)
1539 centrePoint +=
points[fp];
1541 centrePoint /= nVerts;
1543 for (label i = 0; i < nVerts; ++i)
1545 const label pt0 =
f.thisLabel(i);
1560 bitSet nbrBoundaryAdjust;
1564 const faPatch& fap =
boundary()[patchi];
1566 if (isA<wedgeFaPatch>(fap))
1570 const auto& wedgePatch = refCast<const wedgeFaPatch>(fap);
1572 const labelList& patchPoints = wedgePatch.pointLabels();
1579 wedgePatch.centreNormal()
1583 for (
const label pti : patchPoints)
1589 if (wedgePatch.axisPoint() > -1)
1591 result[wedgePatch.axisPoint()] =
1595 &result[wedgePatch.axisPoint()]
1602 const auto& procPatch = refCast<const processorFaPatch>(fap);
1604 const labelList& patchPoints = procPatch.pointLabels();
1605 const labelList& nbrPatchPoints = procPatch.neighbPoints();
1608 procPatch.nonGlobalPatchPoints();
1614 UIndirectList<vector>(result, patchPoints)
1621 procPatch.neighbProcNo(),
1627 patchPointNormals.resize_nocopy(nbrPatchPoints.size());
1633 procPatch.neighbProcNo(),
1638 for (
const label pti : nonGlobalPatchPoints)
1640 result[patchPoints[pti]] +=
1641 patchPointNormals[nbrPatchPoints[pti]];
1649 else if (correctPatchPointNormals(patchi) && !fap.coupled())
1652 nbrBoundaryAdjust.set(patchi);
1658 if (globalData().nGlobalPoints())
1660 const labelList& spLabels = globalData().sharedPointLabels();
1661 const labelList& addr = globalData().sharedPointAddr();
1665 UIndirectList<vector>(result, spLabels)
1670 globalData().nGlobalPoints(),
1676 gpNormals[addr[i]] += spNormals[i];
1684 spNormals[i] = gpNormals[addr[i]];
1687 forAll(spNormals, pointI)
1689 result[spLabels[pointI]] = spNormals[pointI];
1696 hasPatchPointNormalsCorrection()
1701 <<
"Apply " << nbrBoundaryAdjust.count()
1702 <<
" boundary neighbour corrections" <<
nl;
1708 const auto& haloNormals = this->haloFaceNormals();
1710 Map<vector> fpNormals(4*nBoundaryEdges());
1712 for (
const label patchi : nbrBoundaryAdjust)
1714 const faPatch& fap =
boundary()[patchi];
1716 if (fap.ngbPolyPatchIndex() < 0)
1719 <<
"Neighbour polyPatch index is not defined " 1720 <<
"for faPatch " << fap.name()
1725 for (
const label edgei : fap.edgeLabels())
1727 const edge&
e =
patch().edges()[edgei];
1730 const vector& fnorm = haloNormals[edgei - nInternalEdges_];
1732 fpNormals(
e.first()) += fnorm;
1733 fpNormals(
e.second()) += fnorm;
1749 const label pointi = iter.key();
1750 result[pointi].removeCollinear(
normalised(iter.val()));
1758 void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(
vectorField& result)
const 1760 const labelList intPoints(internalPoints());
1761 const labelList bndPoints(boundaryPoints());
1767 forAll(intPoints, pointI)
1769 label curPoint = intPoints[pointI];
1772 const labelList curFaces(faceSet.toc());
1776 for (
const label facei : curFaces)
1778 const labelList& facePoints = faces[facei];
1779 pointSet.insert(facePoints);
1781 pointSet.erase(curPoint);
1784 if (pointSet.size() < 5)
1787 <<
"WARNING: Extending point set for fitting." <<
endl;
1791 for (
const label facei : curFaces)
1794 faceSet.insert(curFaceFaces);
1796 curFaces = faceSet.toc();
1800 for (
const label facei : curFaces)
1802 const labelList& facePoints = faces[facei];
1803 pointSet.insert(facePoints);
1805 pointSet.erase(curPoint);
1806 curPoints = pointSet.toc();
1811 for (label i=0; i<curPoints.size(); ++i)
1818 coordSystem::cartesian cs
1827 p = cs.localPosition(
p);
1837 for (label i = 0; i <
allPoints.size(); ++i)
1848 for (label i = 0; i < MtM.n(); ++i)
1850 for (label j = 0; j < MtM.m(); ++j)
1852 for (label
k = 0;
k <
M.n(); ++
k)
1854 MtM[i][j] +=
M[
k][i]*
M[
k][j]*W[
k];
1861 for (label i=0; i<MtR.size(); ++i)
1863 for (label j=0; j<
M.n(); ++j)
1873 curNormal = cs.globalVector(curNormal);
1875 curNormal *=
sign(curNormal&result[curPoint]);
1877 result[curPoint] = curNormal;
1883 const faPatch& fap =
boundary()[patchI];
1887 const processorFaPatch& procPatch =
1888 refCast<const processorFaPatch>(
boundary()[patchI]);
1890 const labelList& patchPointLabels = procPatch.pointLabels();
1892 labelList toNgbProcLsPointStarts(patchPointLabels.size(),
Zero);
1895 10*patchPointLabels.size(),
1900 for (label pointI=0; pointI<patchPointLabels.size(); ++pointI)
1902 label curPoint = patchPointLabels[pointI];
1904 toNgbProcLsPointStarts[pointI] =
nPoints;
1907 const labelList curFaces(faceSet.toc());
1911 for (
const label facei : curFaces)
1913 const labelList& facePoints = faces[facei];
1914 pointSet.insert(facePoints);
1916 pointSet.erase(curPoint);
1919 for (label i=0; i<curPoints.size(); ++i)
1928 OPstream toNeighbProc
1931 procPatch.neighbProcNo(),
1932 toNgbProcLsPoints.size_bytes()
1933 + toNgbProcLsPointStarts.size_bytes()
1938 << toNgbProcLsPoints
1939 << toNgbProcLsPointStarts;
1944 for (
const faPatch& fap :
boundary())
1948 const processorFaPatch& procPatch =
1949 refCast<const processorFaPatch>(fap);
1951 const labelList& patchPointLabels = procPatch.pointLabels();
1953 labelList fromNgbProcLsPointStarts(patchPointLabels.size(),
Zero);
1957 IPstream fromNeighbProc
1960 procPatch.neighbProcNo(),
1961 10*patchPointLabels.size()*
sizeof(
vector)
1962 + fromNgbProcLsPointStarts.size_bytes()
1967 >> fromNgbProcLsPoints
1968 >> fromNgbProcLsPointStarts;
1972 procPatch.nonGlobalPatchPoints();
1974 forAll(nonGlobalPatchPoints, pointI)
1977 patchPointLabels[nonGlobalPatchPoints[pointI]];
1979 procPatch.neighbPoints()[nonGlobalPatchPoints[pointI]];
1982 faceSet.
insert(pointFaces[curPoint]);
1988 for (
const label facei : curFaces)
1990 const labelList& facePoints = faces[facei];
1991 pointSet.insert(facePoints);
1993 pointSet.erase(curPoint);
1996 label nAllPoints = curPoints.
size();
1998 if (curNgbPoint == fromNgbProcLsPointStarts.size() - 1)
2001 fromNgbProcLsPoints.size()
2002 - fromNgbProcLsPointStarts[curNgbPoint];
2007 fromNgbProcLsPointStarts[curNgbPoint + 1]
2008 - fromNgbProcLsPointStarts[curNgbPoint];
2013 for (label i=0; i<curPoints.size(); ++i)
2015 allPointsExt[counter++] =
points[curPoints[i]];
2018 if (curNgbPoint == fromNgbProcLsPointStarts.size() - 1)
2022 label i=fromNgbProcLsPointStarts[curNgbPoint];
2023 i<fromNgbProcLsPoints.size();
2027 allPointsExt[counter++] = fromNgbProcLsPoints[i];
2034 label i=fromNgbProcLsPointStarts[curNgbPoint];
2035 i<fromNgbProcLsPointStarts[curNgbPoint+1];
2039 allPointsExt[counter++] = fromNgbProcLsPoints[i];
2045 const scalar tol = 0.001*treeBoundBox(allPointsExt).mag();
2048 for (
const point& pt : allPointsExt)
2050 bool duplicate =
false;
2051 for (label i = 0; i < nAllPoints; ++i)
2071 <<
"There are no enough points for quadrics " 2072 <<
"fitting for a point at processor patch" 2080 dir.removeCollinear(axis);
2083 const coordinateSystem cs(
"cs", origin, axis, dir);
2101 for (label i=0; i<
allPoints.size(); ++i)
2112 for (label i = 0; i < MtM.n(); ++i)
2114 for (label j = 0; j < MtM.m(); ++j)
2116 for (label
k = 0;
k <
M.n(); ++
k)
2118 MtM[i][j] +=
M[
k][i]*
M[
k][j]*W[
k];
2125 for (label i = 0; i < MtR.size(); ++i)
2127 for (label j = 0; j <
M.n(); ++j)
2137 curNormal = cs.globalVector(curNormal);
2139 curNormal *=
sign(curNormal&result[curPoint]);
2141 result[curPoint] = curNormal;
2147 if (globalData().nGlobalPoints() > 0)
2149 const labelList& spLabels = globalData().sharedPointLabels();
2151 const labelList& addr = globalData().sharedPointAddr();
2153 for (label
k=0;
k<globalData().nGlobalPoints(); ++
k)
2157 const label curSharedPointIndex = addr.find(
k);
2161 if (curSharedPointIndex != -1)
2163 label curPoint = spLabels[curSharedPointIndex];
2166 const labelList curFaces(faceSet.toc());
2169 for (
const label facei : curFaces)
2171 const labelList& facePoints = faces[facei];
2172 pointSet.insert(facePoints);
2174 pointSet.erase(curPoint);
2181 tol = 0.001*treeBoundBox(locPoints).mag();
2186 if (curSharedPointIndex != -1)
2188 label curPoint = spLabels[curSharedPointIndex];
2190 label nAllPoints = 0;
2191 forAll(procLsPoints, procI)
2193 nAllPoints += procLsPoints[procI].size();
2199 forAll(procLsPoints, procI)
2201 forAll(procLsPoints[procI], pointI)
2203 bool duplicate =
false;
2204 for (label i=0; i<nAllPoints; ++i)
2211 - procLsPoints[procI][pointI]
2224 procLsPoints[procI][pointI];
2234 <<
"There are no enough points for quadratic " 2235 <<
"fitting for a global processor point " 2243 dir.removeCollinear(axis);
2246 coordinateSystem cs(
"cs", origin, axis, dir);
2265 for (label i=0; i<
allPoints.size(); ++i)
2275 for (label i = 0; i < MtM.n(); ++i)
2277 for (label j = 0; j < MtM.m(); ++j)
2279 for (label
k = 0;
k <
M.n(); ++
k)
2281 MtM[i][j] +=
M[
k][i]*
M[
k][j]*W[
k];
2287 for (label i = 0; i < MtR.size(); ++i)
2289 for (label j = 0; j <
M.n(); ++j)
2299 curNormal = cs.globalVector(curNormal);
2301 curNormal *=
sign(curNormal&result[curPoint]);
2303 result[curPoint] = curNormal;
2318 <<
"Calculating edge length correction" <<
endl;
2324 "edgeLengthCorrection",
2325 mesh().pointsInstance(),
2337 const vectorField& pointNormals = pointAreaNormals();
2339 const auto angleCorrection =
2352 fld[edgei] = angleCorrection
2354 pointNormals[edges_[edgei].start()],
2355 pointNormals[edges_[edgei].
end()]
2366 const faPatch& fap =
boundary()[patchi];
2369 label edgei = fap.start();
2373 pfld[patchEdgei] = angleCorrection
2375 pointNormals[edges_[edgei].start()],
2376 pointNormals[edges_[edgei].
end()]
2395 mesh().pointsInstance(),
2404 (this->Le() / this->magLe())
2411 return tunitVectors;
dimensionedScalar sign(const dimensionedScalar &ds)
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix by subtracting the matrix multiplied by the current fi...
labelList internalPoints() const
Return internal point labels.
Point unitVec() const
Return the unit vector (start-to-end)
const polyBoundaryMesh & pbm
void size(const label n)
Older name for setAddressableSize.
"blocking" : (MPI_Bsend, MPI_Recv)
void set(const bitSet &bitset)
Set specified bits from another bitset.
tmp< edgeScalarField > edgeLengthCorrection() const
Return edge length correction.
static label read(const UPstream::commsTypes commsType, const int fromProcNo, char *buf, const std::streamsize bufSize, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm, UPstream::Request *req=nullptr)
Read buffer contents from given processor.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
scalar mag() const
The magnitude (length) of the line.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedTensor skew(const dimensionedTensor &dt)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr char nl
The newline '\n' character (0x0a)
labelList sortedToc() const
The indices of the on bits as a sorted labelList.
Ostream & endl(Ostream &os)
Add newline and flush stream.
labelPair whichPatchFace(const label meshFacei) const
Lookup mesh face index and return (patchi, patchFacei) tuple or (-1, meshFacei) for internal faces...
static bool & parRun() noexcept
Test if this a parallel run.
static int & msgType() noexcept
Message tag of standard messages.
SubList< vector > subList
Declare type of subList.
static vector areaInvDistSqrWeightedNormalFaceTriangle(const linePointRef &edgeLine, const point &faceCentre)
label k
Boltzmann constant.
Ignore writing from objectRegistry::writeObject()
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
const dimensionSet dimless
Dimensionless.
static int myProcNo(const label communicator=worldComm)
Rank of this process in the communicator (starting from masterNo()). Can be negative if the process i...
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
labelList faceLabels(nFaceLabels)
List< labelList > labelListList
List of labelList.
Namespace of functions to calculate explicit derivatives.
static vector areaNormalFaceTriangle(const linePointRef &edgeLine, const point &faceCentre)
virtual const objectRegistry & thisDb() const
Return reference to the mesh database.
dimensionedScalar asin(const dimensionedScalar &ds)
static void waitRequests()
Wait for all requests to finish.
UList< label > labelUList
A UList of labels.
static void imposeMinVectorLength(vectorField &fld)
#define forAll(list, i)
Loop across all elements in list.
static int geometryOrder() noexcept
Return the current geometry treatment.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
static void allGatherList(List< T > &values, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Gather data, but keep individual values separate. Uses linear/tree communication. ...
List< face > faceList
List of faces.
A list of faces which address into the list of points.
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator). It is 1 for serial run. ...
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
dimensionedScalar cos(const dimensionedScalar &ds)
line< point, const point & > linePointRef
A line using referred points.
SubList< T > slice(const label pos, label len=-1)
Return SubList slice (non-const access) - no range checking.
static const Identity< scalar > I
An edge is a list of two vertex labels. This can correspond to a directed graph edge or an edge on a ...
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ const Field< vector > edgeCentres(faMeshTools::flattenEdgeField(aMesh.edgeCentres(), true))
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
#define DebugInFunction
Report an information message using Foam::Info.
GeometricField< vector, faePatchField, edgeMesh > edgeVectorField
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
PointRef a() const noexcept
The first point.
GeometricField< scalar, faePatchField, edgeMesh > edgeScalarField
PointRef b() const noexcept
The second point.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
static void combineReduce(const List< commsStruct > &comms, T &value, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) applying cop to inplace combine value from different processors...
#define DebugInfo
Report an information message using Foam::Info.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
Pair< label > labelPair
A pair of labels.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
static Foam::bitSet markupBoundaryPoints(const uindirectPrimitivePatch &p)
bool hasFaceAreas() const noexcept
tmp< GeometricField< Type, faPatchField, areaMesh > > edgeIntegrate(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
const Vector< label > N(dict.get< Vector< label >>("N"))
Point centre() const
Return centre (centroid)
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< ' ';}gmvFile<< nl;for(const word &name :lagrangianScalarNames){ IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Point vec() const
Return start-to-end vector.
const vectorField & faceCentres() const
labelList boundaryPoints() const
Return boundary point labels.
static word meshSubDir
The mesh sub-directory name (usually "faMesh")
vector point
Point is a vector.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
static vector calcLeVector(const point &faceCentre, const linePointRef &edgeLine, const vector &edgeNormal)
Calculate the 'Le' vector from faceCentre to edge centre.
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
const vectorField & faceAreas() const
const std::string patch
OpenFOAM patch number as a std::string.
static bool write(const UPstream::commsTypes commsType, const int toProcNo, const char *buf, const std::streamsize bufSize, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm, UPstream::Request *req=nullptr, const UPstream::sendModes sendMode=UPstream::sendModes::normal)
Write buffer contents to given processor.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< label > labelList
A List of labels.
A class for managing temporary objects.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Vector< Cmpt > & removeCollinear(const Vector< Cmpt > &unitVec)
Inplace removal of components that are collinear to the given unit vector.
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
Tensor of scalars, i.e. Tensor<scalar>.
Defines the attributes of an object for which implicit objectRegistry management is supported...
tmp< edgeVectorField > unitLe() const
Return normalised edge length vectors.
SquareMatrix< scalar > scalarSquareMatrix
GeometricField< vector, faPatchField, areaMesh > areaVectorField
static vector areaInvDistSqrWeightedNormalDualEdge(const point &basePoint, const point &rightPoint, const point &leftPoint, const point ¢rePoint)
Do not request registration (bool: false)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
static tensor rotation_e3e1(const vector &unitAxis1, vector dirn)
const dimensionSet dimArea(sqr(dimLength))
forAllConstIters(mixture.phases(), phase)
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
Solve the matrix using LU decomposition with pivoting returning the LU form of the matrix and the sol...
RectangularMatrix< scalar > scalarRectangularMatrix
vector areaNormal() const
The area normal - with magnitude equal to area of triangle.
GeometricField< scalar, faPatchField, areaMesh > areaScalarField
static constexpr const zero Zero
Global zero (0)
static vector areaInvDistSqrWeightedNormal(const vector &a, const vector &b)