45 void Foam::faMesh::calcLduAddressing()
const
48 <<
"Calculating addressing" <<
endl;
53 <<
"lduPtr_ already allocated"
57 lduPtr_ =
new faMeshLduAddressing(*
this);
61 void Foam::faMesh::calcPatchStarts()
const
64 <<
"Calculating patch starts" <<
endl;
69 <<
"patchStartsPtr_ already allocated"
74 labelList& patchStarts = *patchStartsPtr_;
76 patchStarts[0] = nInternalEdges();
81 patchStarts[i - 1] +
boundary()[i - 1].faPatch::size();
86 void Foam::faMesh::calcLe()
const
89 <<
"Calculating local edges" <<
endl;
94 <<
"LePtr_ already allocated"
104 mesh().pointsInstance(),
124 const vectorField& edgeNormalsInternal = edgeNormals.internalField();
125 const vectorField& fCentresInternal = fCentres.internalField();
126 const vectorField& eCentresInternal = eCentres.internalField();
127 const scalarField& magLeInternal = magLe().internalField();
132 pEdges[edgeI].vec(pPoints) ^ edgeNormalsInternal[edgeI];
139 fCentresInternal[owner()[edgeI]]
140 - eCentresInternal[edgeI]
145 magLeInternal[edgeI]/
mag(leInternal[edgeI]);
153 boundary()[patchI].patchSlice(pEdges);
156 edgeNormals.boundaryField()[patchI];
158 vectorField& patchLe = Le.boundaryFieldRef()[patchI];
159 const vectorField& patchECentres = eCentres.boundaryField()[patchI];
164 bndEdges[edgeI].vec(pPoints) ^ bndEdgeNormals[edgeI];
171 fCentresInternal[bndEdgeFaces[edgeI]]
172 - patchECentres[edgeI]
177 magLe().boundaryField()[patchI][edgeI]
178 /
mag(patchLe[edgeI]);
184 void Foam::faMesh::calcMagLe()
const
187 <<
"Calculating local edge magnitudes" <<
endl;
192 <<
"magLePtr_ already allocated"
202 mesh().pointsInstance(),
218 forAll(internalEdges, edgeI)
220 magLe.ref()[edgeI] = internalEdges[edgeI].mag(localPoints);
227 boundary()[patchI].patchSlice(edges());
231 magLe.boundaryFieldRef()[patchI][edgeI] =
232 patchEdges[edgeI].mag(localPoints);
238 void Foam::faMesh::calcAreaCentres()
const
241 <<
"Calculating face centres" <<
endl;
246 <<
"centresPtr_ already allocated"
256 mesh().pointsInstance(),
267 const faceList& localFaces = faces();
271 centres.ref()[faceI] = localFaces[faceI].centre(localPoints);
277 boundary()[patchI].patchSlice(edges());
281 centres.boundaryFieldRef()[patchI][edgeI] =
282 patchEdges[edgeI].centre(localPoints);
286 forAll(centres.boundaryField(), patchI)
291 isA<processorFaPatchVectorField>
293 centres.boundaryField()[patchI]
297 centres.boundaryFieldRef()[patchI].initEvaluate();
298 centres.boundaryFieldRef()[patchI].evaluate();
304 void Foam::faMesh::calcEdgeCentres()
const
307 <<
"Calculating edge centres" <<
endl;
312 <<
"edgeCentresPtr_ already allocated"
322 mesh().pointsInstance(),
338 forAll(internalEdges, edgeI)
340 edgeCentres.ref()[edgeI] = internalEdges[edgeI].centre(localPoints);
347 boundary()[patchI].patchSlice(edges());
351 edgeCentres.boundaryFieldRef()[patchI][edgeI] =
352 patchEdges[edgeI].centre(localPoints);
358 void Foam::faMesh::calcS()
const
361 <<
"Calculating areas" <<
endl;
366 <<
"SPtr_ already allocated"
370 SPtr_ =
new DimensionedField<scalar, areaMesh>
383 DimensionedField<scalar, areaMesh>& S = *SPtr_;
386 const faceList& localFaces = faces();
390 S[faceI] = localFaces[faceI].mag(localPoints);
395 void Foam::faMesh::calcFaceAreaNormals()
const
398 <<
"Calculating face area normals" <<
endl;
400 if (faceAreaNormalsPtr_)
403 <<
"faceAreaNormalsPtr_ already allocated"
407 faceAreaNormalsPtr_ =
413 mesh().pointsInstance(),
424 const faceList& localFaces = faces();
429 nInternal[faceI] = localFaces[faceI].unitNormal(localPoints);
434 faceAreaNormals.boundaryFieldRef()[patchI] =
435 edgeAreaNormals().boundaryField()[patchI];
438 forAll(faceAreaNormals.boundaryField(), patchI)
442 isA<processorFaPatchVectorField>
444 faceAreaNormals.boundaryField()[patchI]
448 faceAreaNormals.boundaryFieldRef()[patchI].initEvaluate();
449 faceAreaNormals.boundaryFieldRef()[patchI].evaluate();
455 void Foam::faMesh::calcEdgeAreaNormals()
const
458 <<
"Calculating edge area normals" <<
endl;
460 if (edgeAreaNormalsPtr_)
463 <<
"edgeAreaNormalsPtr_ already allocated"
467 edgeAreaNormalsPtr_ =
473 mesh().pointsInstance(),
485 const vectorField& pointNormals = pointAreaNormals();
529 forAll(edgeAreaNormals.internalField(), edgeI)
560 edgeAreaNormals.ref()[edgeI] =
561 pointNormals[edges()[edgeI].start()]
562 + pointNormals[edges()[edgeI].end()];
564 edgeAreaNormals.ref()[edgeI] -=
565 e*(
e&edgeAreaNormals.internalField()[edgeI]);
568 edgeAreaNormals.ref() /=
mag(edgeAreaNormals.internalField());
573 boundary()[patchI].patchSlice(edges());
577 edgeAreaNormals.boundaryFieldRef()[patchI][edgeI] =
578 pointNormals[patchEdges[edgeI].start()]
579 + pointNormals[patchEdges[edgeI].end()];
583 edgeAreaNormals.boundaryFieldRef()[patchI][edgeI] -=
584 e*(
e&edgeAreaNormals.boundaryField()[patchI][edgeI]);
587 edgeAreaNormals.boundaryFieldRef()[patchI] /=
588 mag(edgeAreaNormals.boundaryField()[patchI]);
593 void Foam::faMesh::calcFaceCurvatures()
const
596 <<
"Calculating face curvatures" <<
endl;
598 if (faceCurvaturesPtr_)
601 <<
"faceCurvaturesPtr_ already allocated"
611 mesh().pointsInstance(),
628 faceCurvatures =
sign(kN&faceAreaNormals())*
mag(kN);
632 void Foam::faMesh::calcEdgeTransformTensors()
const
635 <<
"Calculating edge transformation tensors" <<
endl;
637 if (edgeTransformTensorsPtr_)
640 <<
"edgeTransformTensorsPtr_ already allocated"
644 edgeTransformTensorsPtr_ =
new FieldField<Field, tensor>(nEdges());
645 FieldField<Field, tensor>& edgeTransformTensors =
646 *edgeTransformTensorsPtr_;
655 for (
label edgeI=0; edgeI<nInternalEdges(); ++edgeI)
657 edgeTransformTensors.set(edgeI,
new Field<tensor>(3,
I));
659 vector E = Ce.internalField()[edgeI];
663 E -= skewCorrectionVectors().internalField()[edgeI];
667 vector il = E - Cf.internalField()[owner()[edgeI]];
669 il -= Ne.internalField()[edgeI]
670 *(Ne.internalField()[edgeI]&il);
674 vector kl = Ne.internalField()[edgeI];
677 edgeTransformTensors[edgeI][0] =
680 il.x(), il.y(), il.z(),
681 jl.x(), jl.y(), jl.z(),
682 kl.x(), kl.y(), kl.z()
686 il = E - Cf.internalField()[owner()[edgeI]];
688 il -= Nf.internalField()[owner()[edgeI]]
689 *(Nf.internalField()[owner()[edgeI]]&il);
693 kl = Nf.internalField()[owner()[edgeI]];
696 edgeTransformTensors[edgeI][1] =
699 il.x(), il.y(), il.z(),
700 jl.x(), jl.y(), jl.z(),
701 kl.x(), kl.y(), kl.z()
705 il = Cf.internalField()[neighbour()[edgeI]] - E;
707 il -= Nf.internalField()[neighbour()[edgeI]]
708 *(Nf.internalField()[neighbour()[edgeI]]&il);
712 kl = Nf.internalField()[neighbour()[edgeI]];
715 edgeTransformTensors[edgeI][2] =
718 il.x(), il.y(), il.z(),
719 jl.x(), jl.y(), jl.z(),
720 kl.x(), kl.y(), kl.z()
732 vectorField ngbCf(Cf.boundaryField()[patchI].patchNeighbourField());
734 vectorField ngbNf(Nf.boundaryField()[patchI].patchNeighbourField());
738 edgeTransformTensors.set
741 new Field<tensor>(3,
I)
744 vector E = Ce.boundaryField()[patchI][edgeI];
748 E -= skewCorrectionVectors()
749 .boundaryField()[patchI][edgeI];
753 vector il = E - Cf.internalField()[edgeFaces[edgeI]];
755 il -= Ne.boundaryField()[patchI][edgeI]
756 *(Ne.boundaryField()[patchI][edgeI]&il);
760 vector kl = Ne.boundaryField()[patchI][edgeI];
763 edgeTransformTensors[
boundary()[patchI].start() + edgeI][0] =
767 il = E - Cf.internalField()[edgeFaces[edgeI]];
769 il -= Nf.internalField()[edgeFaces[edgeI]]
770 *(Nf.internalField()[edgeFaces[edgeI]]&il);
774 kl = Nf.internalField()[edgeFaces[edgeI]];
777 edgeTransformTensors[
boundary()[patchI].start() + edgeI][1] =
781 il = ngbCf[edgeI] - E;
783 il -= ngbNf[edgeI]*(ngbNf[edgeI]&il);
791 edgeTransformTensors[
boundary()[patchI].start() + edgeI][2] =
801 edgeTransformTensors.set
804 new Field<tensor>(3,
I)
807 vector E = Ce.boundaryField()[patchI][edgeI];
811 E -= skewCorrectionVectors()
812 .boundaryField()[patchI][edgeI];
816 vector il = E - Cf.internalField()[edgeFaces[edgeI]];
818 il -= Ne.boundaryField()[patchI][edgeI]
819 *(Ne.boundaryField()[patchI][edgeI]&il);
823 vector kl = Ne.boundaryField()[patchI][edgeI];
826 edgeTransformTensors[
boundary()[patchI].start() + edgeI][0] =
830 il = E - Cf.internalField()[edgeFaces[edgeI]];
832 il -= Nf.internalField()[edgeFaces[edgeI]]
833 *(Nf.internalField()[edgeFaces[edgeI]]&il);
837 kl = Nf.internalField()[edgeFaces[edgeI]];
840 edgeTransformTensors[
boundary()[patchI].start() + edgeI][1] =
851 <<
"Calculating internal points" <<
endl;
858 for (
label curEdge = nIntEdges; curEdge < edges.size(); ++curEdge)
860 internal[edges[curEdge].start()] =
false;
862 internal[edges[curEdge].end()] =
false;
869 if (
internal[pointI])
871 internalPoints.
append(pointI);
884 <<
"Calculating boundary points" <<
endl;
891 for (
label curEdge = nIntEdges; curEdge < edges.size(); ++curEdge)
893 internal[edges[curEdge].start()] =
false;
895 internal[edges[curEdge].end()] =
false;
902 if (!
internal[pointI])
904 boundaryPoints.
append(pointI);
914 void Foam::faMesh::calcPointAreaNormals()
const
916 if (pointAreaNormalsPtr_)
919 <<
"pointAreaNormalsPtr_ already allocated"
935 for (
const label curPoint : intPoints)
937 faceList curFaceList(pointFaces[curPoint].size());
939 forAll(curFaceList, faceI)
941 curFaceList[faceI] = faces[pointFaces[curPoint][faceI]];
946 labelList curPointPoints = curPatch.edgeLoops()[0];
948 for (
int i = 0; i < curPointPoints.size(); ++i)
951 points[curPatch.meshPoints()[curPointPoints[i]]]
956 if (i == (curPointPoints.size() - 1))
962 points[curPatch.meshPoints()[curPointPoints[
p]]]
967 scalar sinAlpha =
mag(d1^d2)/(
mag(d1)*
mag(d2));
969 scalar w = sinAlpha/(
mag(d1)*
mag(d2));
971 result[curPoint] += w*
n;
975 for (
const label curPoint : bndPoints)
977 faceList curFaceList(pointFaces[curPoint].size());
979 forAll(curFaceList, faceI)
981 curFaceList[faceI] = faces[pointFaces[curPoint][faceI]];
986 labelList agglomFacePoints = curPatch.edgeLoops()[0];
988 SLList<label> slList;
990 label curPointLabel = -1;
992 for (
label i=0; i<agglomFacePoints.size(); ++i)
994 if (curPatch.meshPoints()[agglomFacePoints[i]] == curPoint)
998 else if ( curPointLabel != -1 )
1000 slList.append(curPatch.meshPoints()[agglomFacePoints[i]]);
1004 for (
label i=0; i<curPointLabel; ++i)
1006 slList.append(curPatch.meshPoints()[agglomFacePoints[i]]);
1011 for (
label i=0; i < (curPointPoints.size() - 1); ++i)
1019 scalar sinAlpha =
mag(d1 ^ d2)/(
mag(d1)*
mag(d2));
1021 scalar w = sinAlpha/(
mag(d1)*
mag(d2));
1023 result[curPoint] += w*
n;
1030 const faPatch& fap =
boundary()[patchI];
1032 if (isA<wedgeFaPatch>(fap))
1034 const wedgeFaPatch& wedgePatch = refCast<const wedgeFaPatch>(fap);
1036 const labelList& patchPoints = wedgePatch.pointLabels();
1042 wedgePatch.centreNormal()
1047 for (
const label pti : patchPoints)
1049 result[pti] -=
N*(
N&result[pti]);
1053 if (wedgePatch.axisPoint() > -1)
1055 result[wedgePatch.axisPoint()] =
1059 &result[wedgePatch.axisPoint()]
1067 for (
const faPatch& fap :
boundary())
1071 const processorFaPatch& procPatch =
1072 refCast<const processorFaPatch>(fap);
1074 const labelList& patchPointLabels = procPatch.pointLabels();
1078 patchPointLabels.size(),
1082 forAll(patchPointNormals, pointI)
1084 patchPointNormals[pointI] = result[patchPointLabels[pointI]];
1091 procPatch.neighbProcNo(),
1092 reinterpret_cast<const char*
>(patchPointNormals.begin()),
1093 patchPointNormals.byteSize()
1099 procPatch.neighbPoints().size(),
1107 procPatch.neighbProcNo(),
1108 reinterpret_cast<char*
>(ngbPatchPointNormals.begin()),
1109 ngbPatchPointNormals.byteSize()
1114 procPatch.nonGlobalPatchPoints();
1116 for (
const label pti : nonGlobalPatchPoints)
1118 result[patchPointLabels[pti]] +=
1119 ngbPatchPointNormals[procPatch.neighbPoints()[pti]];
1126 if (globalData().nGlobalPoints() > 0)
1128 const labelList& spLabels(globalData().sharedPointLabels());
1131 forAll(spNormals, pointI)
1133 spNormals[pointI] = result[spLabels[pointI]];
1136 const labelList& addr = globalData().sharedPointAddr();
1140 globalData().nGlobalPoints(),
1146 gpNormals[addr[i]] += spNormals[i];
1154 spNormals[i] = gpNormals[addr[i]];
1157 forAll(spNormals, pointI)
1159 result[spLabels[pointI]] = spNormals[pointI];
1167 const faPatch& fap =
boundary()[patchI];
1169 if (correctPatchPointNormals(patchI) && !fap.coupled())
1171 if (fap.ngbPolyPatchIndex() == -1)
1174 <<
"Neighbour polyPatch index is not defined "
1175 <<
"for faPatch " << fap.name()
1179 const labelList& patchPoints = fap.pointLabels();
1183 forAll(patchPoints, pointI)
1185 result[patchPoints[pointI]]
1186 -=
N[pointI]*(
N[pointI]&result[patchPoints[pointI]]);
1191 result /=
mag(result);
1195 void Foam::faMesh::calcPointAreaNormalsByQuadricsFit()
const
1199 const labelList intPoints(internalPoints());
1200 const labelList bndPoints(boundaryPoints());
1206 forAll(intPoints, pointI)
1208 label curPoint = intPoints[pointI];
1211 const labelList curFaces(faceSet.toc());
1215 pointSet.
insert(curPoint);
1216 for (
label i=0; i<curFaces.size(); ++i)
1218 const labelList& facePoints = faces[curFaces[i]];
1219 for (
label j=0; j<facePoints.size(); ++j)
1221 if (!pointSet.found(facePoints[j]))
1223 pointSet.insert(facePoints[j]);
1227 pointSet.erase(curPoint);
1230 if (curPoints.size() < 5)
1234 Info <<
"WARNING: Extending point set for fitting." <<
endl;
1242 patch().faceFaces()[curFaces[faceI]];
1246 label curFaceFace = curFaceFaces[fI];
1248 if (!faceSet.found(curFaceFace))
1250 faceSet.insert(curFaceFace);
1254 curFaces = faceSet.toc();
1258 pointSet.
insert(curPoint);
1259 for (
label i=0; i<curFaces.size(); ++i)
1261 const labelList& facePoints = faces[curFaces[i]];
1262 for (
label j=0; j<facePoints.size(); ++j)
1264 if (!pointSet.found(facePoints[j]))
1266 pointSet.insert(facePoints[j]);
1271 pointSet.erase(curPoint);
1272 curPoints = pointSet.toc();
1277 for (
label i=0; i<curPoints.size(); ++i)
1284 coordSystem::cartesian cs
1293 p = cs.localPosition(
p);
1314 for (
label i = 0; i < MtM.n(); ++i)
1316 for (
label j = 0; j < MtM.m(); ++j)
1320 MtM[i][j] +=
M[
k][i]*
M[
k][j]*W[
k];
1327 for (
label i=0; i<MtR.size(); ++i)
1329 for (
label j=0; j<
M.n(); ++j)
1339 curNormal = cs.globalVector(curNormal);
1341 curNormal *=
sign(curNormal&result[curPoint]);
1343 result[curPoint] = curNormal;
1349 const faPatch& fap =
boundary()[patchI];
1353 const processorFaPatch& procPatch =
1354 refCast<const processorFaPatch>(
boundary()[patchI]);
1356 const labelList& patchPointLabels = procPatch.pointLabels();
1358 labelList toNgbProcLsPointStarts(patchPointLabels.size(),
Zero);
1361 10*patchPointLabels.size(),
1366 for (
label pointI=0; pointI<patchPointLabels.size(); ++pointI)
1368 label curPoint = patchPointLabels[pointI];
1370 toNgbProcLsPointStarts[pointI] =
nPoints;
1373 const labelList curFaces(faceSet.toc());
1377 pointSet.
insert(curPoint);
1378 for (
label i=0; i<curFaces.size(); ++i)
1380 const labelList& facePoints = faces[curFaces[i]];
1381 for (
label j=0; j<facePoints.size(); ++j)
1383 if (!pointSet.found(facePoints[j]))
1385 pointSet.insert(facePoints[j]);
1389 pointSet.erase(curPoint);
1392 for (
label i=0; i<curPoints.size(); ++i)
1398 toNgbProcLsPoints.setSize(
nPoints);
1401 OPstream toNeighbProc
1404 procPatch.neighbProcNo(),
1405 toNgbProcLsPoints.byteSize()
1406 + toNgbProcLsPointStarts.byteSize()
1411 << toNgbProcLsPoints
1412 << toNgbProcLsPointStarts;
1417 for (
const faPatch& fap :
boundary())
1421 const processorFaPatch& procPatch =
1422 refCast<const processorFaPatch>(fap);
1424 const labelList& patchPointLabels = procPatch.pointLabels();
1426 labelList fromNgbProcLsPointStarts(patchPointLabels.size(),
Zero);
1430 IPstream fromNeighbProc
1433 procPatch.neighbProcNo(),
1434 10*patchPointLabels.size()*
sizeof(
vector)
1435 + fromNgbProcLsPointStarts.byteSize()
1440 >> fromNgbProcLsPoints
1441 >> fromNgbProcLsPointStarts;
1445 procPatch.nonGlobalPatchPoints();
1447 forAll(nonGlobalPatchPoints, pointI)
1450 patchPointLabels[nonGlobalPatchPoints[pointI]];
1452 procPatch.neighbPoints()[nonGlobalPatchPoints[pointI]];
1455 faceSet.
insert(pointFaces[curPoint]);
1461 pointSet.
insert(curPoint);
1462 for (
label i=0; i<curFaces.size(); ++i)
1464 const labelList& facePoints = faces[curFaces[i]];
1465 for (
label j=0; j<facePoints.size(); ++j)
1467 if (!pointSet.found(facePoints[j]))
1469 pointSet.insert(facePoints[j]);
1473 pointSet.erase(curPoint);
1476 label nAllPoints = curPoints.size();
1478 if (curNgbPoint == fromNgbProcLsPointStarts.size() - 1)
1481 fromNgbProcLsPoints.size()
1482 - fromNgbProcLsPointStarts[curNgbPoint];
1487 fromNgbProcLsPointStarts[curNgbPoint + 1]
1488 - fromNgbProcLsPointStarts[curNgbPoint];
1493 for (
label i=0; i<curPoints.size(); ++i)
1495 allPointsExt[counter++] =
points[curPoints[i]];
1498 if (curNgbPoint == fromNgbProcLsPointStarts.size() - 1)
1502 label i=fromNgbProcLsPointStarts[curNgbPoint];
1503 i<fromNgbProcLsPoints.size();
1507 allPointsExt[counter++] = fromNgbProcLsPoints[i];
1514 label i=fromNgbProcLsPointStarts[curNgbPoint];
1515 i<fromNgbProcLsPointStarts[curNgbPoint+1];
1519 allPointsExt[counter++] = fromNgbProcLsPoints[i];
1525 boundBox bb(allPointsExt,
false);
1526 scalar tol = 0.001*
mag(bb.max() - bb.min());
1531 bool duplicate =
false;
1532 for (
label i=0; i<nAllPoints; ++i)
1551 allPoints[nAllPoints++] = allPointsExt[pI];
1560 <<
"There are no enough points for quadrics "
1561 <<
"fitting for a point at processor patch"
1569 dir -= axis*(axis&dir);
1572 const coordinateSystem cs(
"cs", origin, axis, dir);
1601 for (
label i = 0; i < MtM.n(); ++i)
1603 for (
label j = 0; j < MtM.m(); ++j)
1607 MtM[i][j] +=
M[
k][i]*
M[
k][j]*W[
k];
1614 for (
label i = 0; i < MtR.size(); ++i)
1616 for (
label j = 0; j <
M.n(); ++j)
1626 curNormal = cs.globalVector(curNormal);
1628 curNormal *=
sign(curNormal&result[curPoint]);
1630 result[curPoint] = curNormal;
1636 if (globalData().nGlobalPoints() > 0)
1638 const labelList& spLabels = globalData().sharedPointLabels();
1640 const labelList& addr = globalData().sharedPointAddr();
1642 for (
label k=0;
k<globalData().nGlobalPoints(); ++
k)
1646 const label curSharedPointIndex = addr.find(
k);
1650 if (curSharedPointIndex != -1)
1652 label curPoint = spLabels[curSharedPointIndex];
1655 const labelList curFaces(faceSet.toc());
1658 pointSet.
insert(curPoint);
1659 for (
label i=0; i<curFaces.size(); ++i)
1661 const labelList& facePoints = faces[curFaces[i]];
1662 for (
label j=0; j<facePoints.size(); ++j)
1664 if (!pointSet.found(facePoints[j]))
1666 pointSet.insert(facePoints[j]);
1670 pointSet.erase(curPoint);
1677 const boundBox bb(locPoints,
false);
1678 tol = 0.001*
mag(bb.max() - bb.min());
1684 if (curSharedPointIndex != -1)
1686 label curPoint = spLabels[curSharedPointIndex];
1688 label nAllPoints = 0;
1689 forAll(procLsPoints, procI)
1691 nAllPoints += procLsPoints[procI].size();
1697 forAll(procLsPoints, procI)
1699 forAll(procLsPoints[procI], pointI)
1701 bool duplicate =
false;
1702 for (
label i=0; i<nAllPoints; ++i)
1709 - procLsPoints[procI][pointI]
1722 procLsPoints[procI][pointI];
1732 <<
"There are no enough points for quadratic "
1733 <<
"fitting for a global processor point "
1741 dir -= axis*(axis&dir);
1744 coordinateSystem cs(
"cs", origin, axis, dir);
1773 for (
label i = 0; i < MtM.n(); ++i)
1775 for (
label j = 0; j < MtM.m(); ++j)
1779 MtM[i][j] +=
M[
k][i]*
M[
k][j]*W[
k];
1785 for (
label i = 0; i < MtR.size(); ++i)
1787 for (
label j = 0; j <
M.n(); ++j)
1797 curNormal = cs.globalVector(curNormal);
1799 curNormal *=
sign(curNormal&result[curPoint]);
1801 result[curPoint] = curNormal;
1806 result /=
mag(result);
1813 <<
"Calculating edge length correction" <<
endl;
1821 "edgeLengthCorrection",
1822 mesh().pointsInstance(),
1832 const vectorField& pointNormals = pointAreaNormals();
1836 scalar sinAlpha =
mag
1838 pointNormals[edges()[edgeI].
start()]^
1839 pointNormals[edges()[edgeI].
end()]
1852 boundary()[patchI].patchSlice(edges())
1855 forAll(patchEdges, edgeI)
1860 pointNormals[patchEdges[edgeI].
start()]
1861 ^ pointNormals[patchEdges[edgeI].
end()]