src/finiteArea/faMesh/faMeshDemandDrivenData.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "faMesh.H"
  30 #include "faMeshLduAddressing.H"
  31 #include "dimensionSet.H"
  32 #include "areaFields.H"
  33 #include "edgeFields.H"
  34 #include "primitiveFacePatch.H"
  35 #include "fac.H"
  36 #include "processorFaPatch.H"
  37 #include "wedgeFaPatch.H"
  38 #include "PstreamCombineReduceOps.H"
  39 #include "coordinateSystem.H"
  40 #include "scalarMatrices.H"
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 namespace Foam
  45 {
  46
  47 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  48
  49 void faMesh::calcLduAddressing() const
  50 {
  51     if (debug)
  52     {
  53         Info<< "void faMesh::calcLduAddressing() const : "
  54             << "Calculating addressing" << endl;
  55     }
  56
  57     if (lduPtr_)
  58     {
  59         FatalErrorIn
  60         (
  61             "void faMesh::calcLduAddressing() const"
  62         )   << "lduPtr_ already allocated"
  63             << abort(FatalError);
  64     }
  65
  66     lduPtr_ = new faMeshLduAddressing(*this);
  67 }
  68
  69
  70 void faMesh::calcPatchStarts() const
  71 {
  72     if (debug)
  73     {
  74         Info<< "void faMesh::calcPatchStarts() const : "
  75             << "Calculating patch starts" << endl;
  76     }
  77
  78     if (patchStartsPtr_)
  79     {
  80         FatalErrorIn
  81         (
  82             "void faMesh::calcPatchStarts() const"
  83         )   << "patchStartsPtr_ already allocated"
  84             << abort(FatalError);
  85     }
  86
  87     patchStartsPtr_ = new labelList(boundary().size(), -1);
  88     labelList& patchStarts = *patchStartsPtr_;
  89
  90     patchStarts[0] = nInternalEdges();
  91
  92     for (label i = 1; i < boundary().size(); i++)
  93     {
  94         patchStarts[i] =
  95             patchStarts[i - 1] + boundary()[i - 1].size();
  96     }
  97 }
  98
  99
 100 void faMesh::calcLe() const
 101 {
 102     if (debug)
 103     {
 104         Info<< "void faMesh::calcLe() const : "
 105             << "Calculating local edges" << endl;
 106     }
 107
 108     if (LePtr_)
 109     {
 110         FatalErrorIn
 111         (
 112             "void faMesh::calcLe() const"
 113         )   << "LePtr_ already allocated"
 114             << abort(FatalError);
 115     }
 116
 117     LePtr_ =
 118         new edgeVectorField
 119         (
 120             IOobject
 121             (
 122                 "Le",
 123                 mesh_.pointsInstance(),
 124                 meshSubDir,
 125                 mesh_
 126             ),
 127             *this,
 128             dimLength
 129         );
 130
 131     edgeVectorField& Le = *LePtr_;
 132
 133
 134     const pointField& pPoints = points();
 135     const edgeList& pEdges = edges();
 136
 137     const edgeVectorField& eCentres = edgeCentres();
 138     const areaVectorField& fCentres = areaCentres();
 139
 140     const edgeVectorField& edgeNormals = edgeAreaNormals();
 141
 142     vectorField& leInternal = Le.internalField();
 143     const vectorField& edgeNormalsInternal = edgeNormals.internalField();
 144     const vectorField& fCentresInternal = fCentres.internalField();
 145     const vectorField& eCentresInternal = eCentres.internalField();
 146     const scalarField& magLeInternal = magLe().internalField();
 147
 148     forAll (leInternal, edgeI)
 149     {
 150         leInternal[edgeI] =
 151             pEdges[edgeI].vec(pPoints) ^ edgeNormalsInternal[edgeI];
 152
 153         leInternal[edgeI] *=
 154           - sign
 155             (
 156                 leInternal[edgeI] &
 157                 (
 158                     fCentresInternal[owner()[edgeI]]
 159                   - eCentresInternal[edgeI]
 160                 )
 161             );
 162
 163         leInternal[edgeI] *=
 164             magLeInternal[edgeI]/mag(leInternal[edgeI]);
 165     }
 166
 167     forAll (boundary(), patchI)
 168     {
 169         const unallocLabelList& bndEdgeFaces =
 170             boundary()[patchI].edgeFaces();
 171
 172         const edgeList::subList bndEdges =
 173             boundary()[patchI].patchSlice(pEdges);
 174
 175         const vectorField& bndEdgeNormals =
 176             edgeNormals.boundaryField()[patchI];
 177
 178         vectorField& patchLe = Le.boundaryField()[patchI];
 179         const vectorField& patchECentres = eCentres.boundaryField()[patchI];
 180
 181         forAll (patchLe, edgeI)
 182         {
 183             patchLe[edgeI] =
 184                 bndEdges[edgeI].vec(pPoints) ^ bndEdgeNormals[edgeI];
 185
 186             patchLe[edgeI] *=
 187               - sign
 188                 (
 189                     patchLe[edgeI]&
 190                     (
 191                         fCentresInternal[bndEdgeFaces[edgeI]]
 192                       - patchECentres[edgeI]
 193                     )
 194                 );
 195
 196             patchLe[edgeI] *=
 197                 magLe().boundaryField()[patchI][edgeI]
 198                 /mag(patchLe[edgeI]);
 199         }
 200     }
 201 }
 202
 203
 204 void faMesh::calcMagLe() const
 205 {
 206     if (debug)
 207     {
 208         Info<< "void faMesh::calcMagLe() const : "
 209             << "Calculating local edge magnitudes" << endl;
 210     }
 211
 212     if (magLePtr_)
 213     {
 214         FatalErrorIn
 215         (
 216             "void faMesh::calcMagLe() const"
 217         )   << "magLePtr_ already allocated"
 218             << abort(FatalError);
 219     }
 220
 221     magLePtr_ =
 222         new edgeScalarField
 223         (
 224             IOobject
 225             (
 226                 "magLe",
 227                 mesh_.pointsInstance(),
 228                 meshSubDir,
 229                 mesh_
 230             ),
 231             *this,
 232             dimLength
 233         );
 234
 235     edgeScalarField& magLe = *magLePtr_;
 236
 237
 238     const pointField& localPoints = points();
 239
 240     const edgeList::subList internalEdges =
 241         edgeList::subList(edges(), nInternalEdges());
 242
 243
 244     forAll (internalEdges, edgeI)
 245     {
 246         magLe.internalField()[edgeI] =
 247             internalEdges[edgeI].mag(localPoints);
 248     }
 249
 250
 251     forAll (boundary(), patchI)
 252     {
 253         const edgeList::subList patchEdges =
 254             boundary()[patchI].patchSlice(edges());
 255
 256         forAll (patchEdges, edgeI)
 257         {
 258             magLe.boundaryField()[patchI][edgeI] =
 259                 patchEdges[edgeI].mag(localPoints);
 260         }
 261     }
 262 }
 263
 264 void faMesh::calcAreaCentres() const
 265 {
 266     if (debug)
 267     {
 268         Info<< "void faMesh::calcAreaCentres() const : "
 269             << "Calculating face centres" << endl;
 270     }
 271
 272     if (centresPtr_)
 273     {
 274         FatalErrorIn
 275         (
 276             "void faMesh::calcAreaCentres() const"
 277         )   << "centresPtr_ already allocated"
 278             << abort(FatalError);
 279     }
 280
 281     centresPtr_ =
 282         new areaVectorField
 283         (
 284             IOobject
 285             (
 286                 "centres",
 287                 mesh_.pointsInstance(),
 288                 meshSubDir,
 289                 mesh_
 290             ),
 291             *this,
 292             dimLength
 293         );
 294     areaVectorField& centres = *centresPtr_;
 295
 296     const pointField& localPoints = points();
 297     const faceList& localFaces = faces();
 298
 299     forAll (localFaces, faceI)
 300     {
 301         centres.internalField()[faceI] =
 302             localFaces[faceI].centre(localPoints);
 303     }
 304
 305     forAll (boundary(), patchI)
 306     {
 307         if (!boundary()[patchI].coupled())
 308         {
 309             const edgeList::subList patchEdges =
 310                 boundary()[patchI].patchSlice(edges());
 311
 312             forAll (patchEdges, edgeI)
 313             {
 314                 centres.boundaryField()[patchI][edgeI] =
 315                     patchEdges[edgeI].centre(localPoints);
 316             }
 317         }
 318     }
 319
 320     centres.correctBoundaryConditions();
 321 }
 322
 323
 324 void faMesh::calcEdgeCentres() const
 325 {
 326     if (debug)
 327     {
 328         Info<< "void faMesh::calcEdgeCentres() const : "
 329             << "Calculating edge centres" << endl;
 330     }
 331
 332     if (edgeCentresPtr_)
 333     {
 334         FatalErrorIn
 335         (
 336             "void faMesh::calcEdgeCentres() const"
 337         )   << "edgeCentresPtr_ already allocated"
 338             << abort(FatalError);
 339     }
 340
 341     edgeCentresPtr_ =
 342         new edgeVectorField
 343         (
 344             IOobject
 345             (
 346                 "edgeCentres",
 347                 mesh_.pointsInstance(),
 348                 meshSubDir,
 349                 mesh_
 350             ),
 351             *this,
 352             dimLength
 353         );
 354
 355     edgeVectorField& edgeCentres = *edgeCentresPtr_;
 356
 357     const pointField& localPoints = points();
 358
 359     const edgeList::subList internalEdges =
 360         edgeList::subList(edges(), nInternalEdges());
 361
 362
 363     forAll (internalEdges, edgeI)
 364     {
 365         edgeCentres.internalField()[edgeI] =
 366             internalEdges[edgeI].centre(localPoints);
 367     }
 368
 369
 370     forAll (boundary(), patchI)
 371     {
 372         const edgeList::subList patchEdges =
 373             boundary()[patchI].patchSlice(edges());
 374
 375         forAll (patchEdges, edgeI)
 376         {
 377             edgeCentres.boundaryField()[patchI][edgeI] =
 378                 patchEdges[edgeI].centre(localPoints);
 379         }
 380     }
 381 }
 382
 383
 384 void faMesh::calcS() const
 385 {
 386     if (debug)
 387     {
 388         Info<< "void faMesh::calcS() const : "
 389             << "Calculating areas" << endl;
 390     }
 391
 392     if (SPtr_)
 393     {
 394         FatalErrorIn
 395         (
 396             "void faMesh::calcS() const"
 397         )   << "SPtr_ already allocated"
 398             << abort(FatalError);
 399     }
 400
 401     SPtr_ = new DimensionedField<scalar, areaMesh>
 402     (
 403         IOobject
 404         (
 405             "S",
 406             time().timeName(),
 407             mesh_,
 408             IOobject::NO_READ,
 409             IOobject::NO_WRITE
 410         ),
 411         *this,
 412         dimVolume
 413     );
 414     DimensionedField<scalar, areaMesh>& S = *SPtr_;
 415
 416     const pointField& localPoints = points();
 417     const faceList& localFaces = faces();
 418
 419     forAll (S, faceI)
 420     {
 421         S[faceI] = localFaces[faceI].mag(localPoints);
 422     }
 423 }
 424
 425
 426 void faMesh::calcFaceAreaNormals() const
 427 {
 428     if (debug)
 429     {
 430         Info<< "void faMesh::calcFaceAreaNormals() const : "
 431             << "Calculating face area normals" << endl;
 432     }
 433
 434     if (faceAreaNormalsPtr_)
 435     {
 436         FatalErrorIn
 437         (
 438             "void faMesh::calcFaceAreaNormals() const"
 439         )   << "faceAreaNormalsPtr_ already allocated"
 440             << abort(FatalError);
 441     }
 442
 443     faceAreaNormalsPtr_ =
 444         new areaVectorField
 445         (
 446             IOobject
 447             (
 448                 "faceAreaNormals",
 449                 mesh_.pointsInstance(),
 450                 meshSubDir,
 451                 mesh_
 452             ),
 453             *this,
 454             dimless
 455         );
 456
 457     areaVectorField& faceAreaNormals = *faceAreaNormalsPtr_;
 458
 459     const pointField& localPoints = points();
 460     const faceList& localFaces = faces();
 461
 462     vectorField& nInternal = faceAreaNormals.internalField();
 463     forAll (localFaces, faceI)
 464     {
 465         nInternal[faceI] =
 466             localFaces[faceI].normal(localPoints)/
 467             localFaces[faceI].mag(localPoints);
 468     }
 469
 470     forAll (boundary(), patchI)
 471     {
 472         if (!boundary()[patchI].coupled())
 473         {
 474             faceAreaNormals.boundaryField()[patchI] =
 475                 edgeAreaNormals().boundaryField()[patchI];
 476         }
 477     }
 478
 479     faceAreaNormals.correctBoundaryConditions();
 480 }
 481
 482
 483 void faMesh::calcEdgeAreaNormals() const
 484 {
 485     if (debug)
 486     {
 487         Info<< "void faMesh::calcEdgeAreaNormals() const : "
 488             << "Calculating edge area normals" << endl;
 489     }
 490
 491     if (edgeAreaNormalsPtr_)
 492     {
 493         FatalErrorIn
 494         (
 495             "void faMesh::calcEdgeAreaNormals() const"
 496         )   << "edgeAreaNormalsPtr_ already allocated"
 497             << abort(FatalError);
 498     }
 499
 500     edgeAreaNormalsPtr_ =
 501         new edgeVectorField
 502         (
 503             IOobject
 504             (
 505                 "edgeAreaNormals",
 506                 mesh_.pointsInstance(),
 507                 meshSubDir,
 508                 mesh_
 509             ),
 510             *this,
 511             dimless
 512         );
 513
 514     edgeVectorField& edgeAreaNormals = *edgeAreaNormalsPtr_;
 515
 516
 517     // Point area normals
 518     const vectorField& pointNormals = pointAreaNormals();
 519
 520
 521 //     // Primitive patch edge normals
 522 //     const labelListList& patchPointEdges = patch().pointEdges();
 523
 524 //     vectorField patchEdgeNormals(nEdges(), vector::zero);
 525
 526 //     forAll (pointNormals, pointI)
 527 //     {
 528 //         const labelList& curPointEdges = patchPointEdges[pointI];
 529
 530 //         forAll (curPointEdges, edgeI)
 531 //         {
 532 //             label curEdge = curPointEdges[edgeI];
 533
 534 //             patchEdgeNormals[curEdge] += 0.5*pointNormals[pointI];
 535 //         }
 536 //     }
 537
 538 //     patchEdgeNormals /= mag(patchEdgeNormals);
 539
 540
 541 //     // Edge area normals
 542 //     label nIntEdges = patch().nInternalEdges();
 543
 544 //     for (label edgeI = 0; edgeI < nIntEdges; edgeI++)
 545 //     {
 546 //         edgeAreaNormals.internalField()[edgeI] =
 547 //             patchEdgeNormals[edgeI];
 548 //     }
 549
 550 //     forAll (boundary(), patchI)
 551 //     {
 552 //         const labelList& edgeLabels = boundary()[patchI];
 553
 554 //         forAll (edgeAreaNormals.boundaryField()[patchI], edgeI)
 555 //         {
 556 //             edgeAreaNormals.boundaryField()[patchI][edgeI] =
 557 //                 patchEdgeNormals[edgeLabels[edgeI]];
 558 //         }
 559 //     }
 560
 561
 562     forAll (edgeAreaNormals.internalField(), edgeI)
 563     {
 564         vector e = edges()[edgeI].vec(points());
 565         e /= mag(e);
 566
 567 //         scalar wStart =
 568 //             1.0 - sqr(mag(e^pointNormals[edges()[edgeI].end()]));
 569
 570 //         scalar wEnd =
 571 //             1.0 - sqr(mag(e^pointNormals[edges()[edgeI].start()]));
 572
 573 //         wStart = 1.0;
 574 //         wEnd = 1.0;
 575
 576 //         edgeAreaNormals.internalField()[edgeI] =
 577 //             wStart*pointNormals[edges()[edgeI].start()]
 578 //           + wEnd*pointNormals[edges()[edgeI].end()];
 579
 580 //         vector eC = 0.5*(points()[edges()[edgeI].start()] + points()[edges()[edgeI].end()]);
 581
 582 //         vector eCp = 0.5*
 583 //             (
 584 //                 points()[edges()[edgeI].start()] + pointNormals[edges()[edgeI].start()]
 585 //                 points()[edges()[edgeI].end()] +
 586 //             );
 587
 588         edgeAreaNormals.internalField()[edgeI] =
 589             pointNormals[edges()[edgeI].start()]
 590           + pointNormals[edges()[edgeI].end()];
 591
 592         edgeAreaNormals.internalField()[edgeI] -=
 593             e*(e&edgeAreaNormals.internalField()[edgeI]);
 594     }
 595
 596     edgeAreaNormals.internalField() /=
 597         mag(edgeAreaNormals.internalField());
 598
 599     forAll (boundary(), patchI)
 600     {
 601         const edgeList::subList patchEdges =
 602             boundary()[patchI].patchSlice(edges());
 603
 604         forAll (patchEdges, edgeI)
 605         {
 606             edgeAreaNormals.boundaryField()[patchI][edgeI] =
 607                 pointNormals[patchEdges[edgeI].start()]
 608               + pointNormals[patchEdges[edgeI].end()];
 609
 610             vector e = patchEdges[edgeI].vec(points());
 611             e /= mag(e);
 612
 613             edgeAreaNormals.boundaryField()[patchI][edgeI] -=
 614                 e*(e&edgeAreaNormals.boundaryField()[patchI][edgeI]);
 615         }
 616
 617         edgeAreaNormals.boundaryField()[patchI] /=
 618             mag(edgeAreaNormals.boundaryField()[patchI]);
 619     }
 620 }
 621
 622
 623 void faMesh::calcFaceCurvatures() const
 624 {
 625     if (debug)
 626     {
 627         Info<< "void faMesh::calcFaceCurvatures() const : "
 628             << "Calculating face curvatures" << endl;
 629     }
 630
 631     if (faceCurvaturesPtr_)
 632     {
 633         FatalErrorIn
 634         (
 635             "void faMesh::calcFaceCurvatures() const"
 636         )   << "faceCurvaturesPtr_ already allocated"
 637             << abort(FatalError);
 638     }
 639
 640     faceCurvaturesPtr_ =
 641         new areaScalarField
 642         (
 643             IOobject
 644             (
 645                 "faceCurvatures",
 646                 mesh_.pointsInstance(),
 647                 meshSubDir,
 648                 mesh_
 649             ),
 650             *this,
 651             dimless/dimLength
 652         );
 653
 654     areaScalarField& faceCurvatures = *faceCurvaturesPtr_;
 655
 656
 657 //     faceCurvatures =
 658 //         fac::edgeIntegrate(Le()*edgeLengthCorrection())
 659 //         &faceAreaNormals();
 660
 661     areaVectorField kN =
 662         fac::edgeIntegrate(Le()*edgeLengthCorrection());
 663
 664     faceCurvatures = sign(kN&faceAreaNormals())*mag(kN);
 665 }
 666
 667
 668 void faMesh::calcEdgeTransformTensors() const
 669 {
 670     if (debug)
 671     {
 672         Info<< "void faMesh::calcEdgeTransformTensors() const : "
 673             << "Calculating edge transformation tensors" << endl;
 674     }
 675
 676     if (edgeTransformTensorsPtr_)
 677     {
 678         FatalErrorIn
 679         (
 680             "void faMesh::calcEdgeTransformTensors() const"
 681         )   << "edgeTransformTensorsPtr_ already allocated"
 682             << abort(FatalError);
 683     }
 684
 685     edgeTransformTensorsPtr_ = new FieldField<Field, tensor>(nEdges());
 686     FieldField<Field, tensor>& edgeTransformTensors =
 687         *edgeTransformTensorsPtr_;
 688
 689     const areaVectorField& Nf = faceAreaNormals();
 690     const areaVectorField& Cf = areaCentres();
 691
 692     const edgeVectorField& Ne = edgeAreaNormals();
 693     const edgeVectorField& Ce = edgeCentres();
 694
 695     // Internal edges transformation tensors
 696     for (label edgeI=0; edgeI<nInternalEdges(); edgeI++)
 697     {
 698         edgeTransformTensors.set(edgeI, new Field<tensor>(3, I));
 699
 700         vector E = Ce.internalField()[edgeI];
 701
 702         if (skew())
 703         {
 704             E -= skewCorrectionVectors().internalField()[edgeI];
 705         }
 706
 707         // Edge transformation tensor
 708         vector il = E - Cf.internalField()[owner()[edgeI]];
 709
 710         il -= Ne.internalField()[edgeI]
 711             *(Ne.internalField()[edgeI]&il);
 712
 713         il /= mag(il);
 714
 715         vector kl = Ne.internalField()[edgeI];
 716         vector jl = kl^il;
 717
 718         edgeTransformTensors[edgeI][0] =
 719             tensor
 720             (
 721                 il.x(), il.y(), il.z(),
 722                 jl.x(), jl.y(), jl.z(),
 723                 kl.x(), kl.y(), kl.z()
 724             );
 725
 726         // Owner transformation tensor
 727         il = E - Cf.internalField()[owner()[edgeI]];
 728
 729         il -= Nf.internalField()[owner()[edgeI]]
 730             *(Nf.internalField()[owner()[edgeI]]&il);
 731
 732         il /= mag(il);
 733
 734         kl = Nf.internalField()[owner()[edgeI]];
 735         jl = kl^il;
 736
 737         edgeTransformTensors[edgeI][1] =
 738             tensor
 739             (
 740                 il.x(), il.y(), il.z(),
 741                 jl.x(), jl.y(), jl.z(),
 742                 kl.x(), kl.y(), kl.z()
 743             );
 744
 745         // Neighbour transformation tensor
 746         il = Cf.internalField()[neighbour()[edgeI]] - E;
 747
 748         il -= Nf.internalField()[neighbour()[edgeI]]
 749             *(Nf.internalField()[neighbour()[edgeI]]&il);
 750
 751         il /= mag(il);
 752
 753         kl = Nf.internalField()[neighbour()[edgeI]];
 754         jl = kl^il;
 755
 756         edgeTransformTensors[edgeI][2] =
 757             tensor
 758             (
 759                 il.x(), il.y(), il.z(),
 760                 jl.x(), jl.y(), jl.z(),
 761                 kl.x(), kl.y(), kl.z()
 762             );
 763     }
 764
 765     // Boundary edges transformation tensors
 766     forAll (boundary(), patchI)
 767     {
 768         if (boundary()[patchI].coupled())
 769         {
 770             const unallocLabelList& edgeFaces =
 771                 boundary()[patchI].edgeFaces();
 772
 773             vectorField ngbCf =
 774                 Cf.boundaryField()[patchI].patchNeighbourField();
 775
 776             vectorField ngbNf =
 777                 Nf.boundaryField()[patchI].patchNeighbourField();
 778
 779             forAll(edgeFaces, edgeI)
 780             {
 781                 edgeTransformTensors.set
 782                 (
 783                     boundary()[patchI].start() + edgeI,
 784                     new Field<tensor>(3, I)
 785                 );
 786
 787                 vector E = Ce.boundaryField()[patchI][edgeI];
 788
 789                 if (skew())
 790                 {
 791                     E -= skewCorrectionVectors()
 792                         .boundaryField()[patchI][edgeI];
 793                 }
 794
 795                 // Edge transformation tensor
 796                 vector il = E - Cf.internalField()[edgeFaces[edgeI]];
 797
 798                 il -= Ne.boundaryField()[patchI][edgeI]
 799                    *(Ne.boundaryField()[patchI][edgeI]&il);
 800
 801                 il /= mag(il);
 802
 803                 vector kl = Ne.boundaryField()[patchI][edgeI];
 804                 vector jl = kl^il;
 805
 806                 edgeTransformTensors[boundary()[patchI].start() + edgeI][0] =
 807                     tensor
 808                     (
 809                         il.x(), il.y(), il.z(),
 810                         jl.x(), jl.y(), jl.z(),
 811                         kl.x(), kl.y(), kl.z()
 812                     );
 813
 814                 // Owner transformation tensor
 815                 il = E - Cf.internalField()[edgeFaces[edgeI]];
 816
 817                 il -= Nf.internalField()[edgeFaces[edgeI]]
 818                    *(Nf.internalField()[edgeFaces[edgeI]]&il);
 819
 820                 il /= mag(il);
 821
 822                 kl = Nf.internalField()[edgeFaces[edgeI]];
 823                 jl = kl^il;
 824
 825                 edgeTransformTensors[boundary()[patchI].start() + edgeI][1] =
 826                     tensor
 827                     (
 828                         il.x(), il.y(), il.z(),
 829                         jl.x(), jl.y(), jl.z(),
 830                         kl.x(), kl.y(), kl.z()
 831                     );
 832
 833                 // Neighbour transformation tensor
 834                 il = ngbCf[edgeI] - E;
 835
 836                 il -= ngbNf[edgeI]*(ngbNf[edgeI]&il);
 837
 838                 il /= mag(il);
 839
 840                 kl = ngbNf[edgeI];
 841
 842                 jl = kl^il;
 843
 844                 edgeTransformTensors[boundary()[patchI].start() + edgeI][2] =
 845                     tensor
 846                     (
 847                         il.x(), il.y(), il.z(),
 848                         jl.x(), jl.y(), jl.z(),
 849                         kl.x(), kl.y(), kl.z()
 850                     );
 851             }
 852         }
 853         else
 854         {
 855             const unallocLabelList& edgeFaces = boundary()[patchI].edgeFaces();
 856
 857             forAll(edgeFaces, edgeI)
 858             {
 859                 edgeTransformTensors.set
 860                 (
 861                     boundary()[patchI].start() + edgeI,
 862                     new Field<tensor>(3, I)
 863                 );
 864
 865                 vector E = Ce.boundaryField()[patchI][edgeI];
 866
 867                 if (skew())
 868                 {
 869                     E -= skewCorrectionVectors()
 870                         .boundaryField()[patchI][edgeI];
 871                 }
 872
 873                 // Edge transformation tensor
 874                 vector il = E - Cf.internalField()[edgeFaces[edgeI]];
 875
 876                 il -= Ne.boundaryField()[patchI][edgeI]
 877                    *(Ne.boundaryField()[patchI][edgeI]&il);
 878
 879                 il /= mag(il);
 880
 881                 vector kl = Ne.boundaryField()[patchI][edgeI];
 882                 vector jl = kl^il;
 883
 884                 edgeTransformTensors[boundary()[patchI].start() + edgeI][0] =
 885                     tensor
 886                     (
 887                         il.x(), il.y(), il.z(),
 888                         jl.x(), jl.y(), jl.z(),
 889                         kl.x(), kl.y(), kl.z()
 890                     );
 891
 892                 // Owner transformation tensor
 893                 il = E - Cf.internalField()[edgeFaces[edgeI]];
 894
 895                 il -= Nf.internalField()[edgeFaces[edgeI]]
 896                    *(Nf.internalField()[edgeFaces[edgeI]]&il);
 897
 898                 il /= mag(il);
 899
 900                 kl = Nf.internalField()[edgeFaces[edgeI]];
 901                 jl = kl^il;
 902
 903                 edgeTransformTensors[boundary()[patchI].start() + edgeI][1] =
 904                     tensor
 905                     (
 906                         il.x(), il.y(), il.z(),
 907                         jl.x(), jl.y(), jl.z(),
 908                         kl.x(), kl.y(), kl.z()
 909                     );
 910             }
 911         }
 912     }
 913 }
 914
 915
 916 labelList faMesh::internalPoints() const
 917 {
 918     if (debug)
 919     {
 920         Info<< "labelList faMesh::internalPoints() const : "
 921             << "Calculating internal points" << endl;
 922     }
 923
 924     const edgeList& edges = patch().edges();
 925     label nIntEdges = patch().nInternalEdges();
 926
 927     List<bool> internal (nPoints(), true);
 928
 929     for (label curEdge = nIntEdges; curEdge < edges.size(); curEdge++)
 930     {
 931         internal[edges[curEdge].start()] = false;
 932
 933         internal[edges[curEdge].end()] = false;
 934     }
 935
 936     SLList<label> internalPoints;
 937
 938     forAll (internal, pointI)
 939     {
 940         if (internal[pointI])
 941         {
 942             internalPoints.append(pointI);
 943         }
 944     }
 945
 946     labelList result(internalPoints);
 947
 948     return result;
 949 }
 950
 951
 952 labelList faMesh::boundaryPoints() const
 953 {
 954     if (debug)
 955     {
 956         Info<< "labelList faMesh::boundaryPoints() const : "
 957             << "Calculating boundary points" << endl;
 958     }
 959
 960     const edgeList& edges = patch().edges();
 961     label nIntEdges = patch().nInternalEdges();
 962
 963     List<bool> internal (nPoints(), true);
 964
 965     for (label curEdge = nIntEdges; curEdge < edges.size(); curEdge++)
 966     {
 967         internal[edges[curEdge].start()] = false;
 968
 969         internal[edges[curEdge].end()] = false;
 970     }
 971
 972     SLList<label> boundaryPoints;
 973
 974     forAll (internal, pointI)
 975     {
 976         if (!internal[pointI])
 977         {
 978             boundaryPoints.append(pointI);
 979         }
 980     }
 981
 982     labelList result(boundaryPoints);
 983
 984     return result;
 985 }
 986
 987
 988 void faMesh::calcPointAreaNormals() const
 989 {
 990     if (pointAreaNormalsPtr_)
 991     {
 992         FatalErrorIn
 993         (
 994             "void faMesh::calcPointAreaNormals() const"
 995         )   << "pointAreaNormalsPtr_ already allocated"
 996             << abort(FatalError);
 997     }
 998
 999
1000     pointAreaNormalsPtr_ =
1001         new vectorField
1002         (
1003             nPoints(),
1004             vector::zero
1005         );
1006
1007     vectorField& result = *pointAreaNormalsPtr_;
1008
1009
1010     labelList intPoints = internalPoints();
1011     labelList bndPoints = boundaryPoints();
1012
1013     const pointField& points = patch().localPoints();
1014     const faceList& faces = patch().localFaces();
1015     const labelListList& pointFaces = patch().pointFaces();
1016
1017     forAll (intPoints, pointI)
1018     {
1019         label curPoint = intPoints[pointI];
1020
1021         faceList curFaceList(pointFaces[curPoint].size());
1022
1023         forAll (curFaceList, faceI)
1024         {
1025             curFaceList[faceI] = faces[pointFaces[curPoint][faceI]];
1026         }
1027
1028         primitiveFacePatch curPatch(curFaceList, points);
1029
1030         labelList curPointPoints = curPatch.edgeLoops()[0];
1031
1032         for (int i = 0; i < curPointPoints.size(); i++)
1033         {
1034             vector d1 =
1035                 points[curPatch.meshPoints()[curPointPoints[i]]]
1036               - points[curPoint];
1037
1038             label p = i + 1;
1039
1040             if (i == (curPointPoints.size() - 1))
1041             {
1042                 p = 0;
1043             }
1044
1045             vector d2 =
1046                 points[curPatch.meshPoints()[curPointPoints[p]]]
1047               - points[curPoint];
1048
1049             vector n = (d1^d2)/mag(d1^d2);
1050
1051             scalar sinAlpha = mag(d1^d2)/(mag(d1)*mag(d2));
1052
1053             scalar w = sinAlpha/(mag(d1)*mag(d2));
1054
1055             result[curPoint] += w*n;
1056         }
1057     }
1058
1059     forAll (bndPoints, pointI)
1060     {
1061         label curPoint = bndPoints[pointI];
1062
1063         faceList curFaceList(pointFaces[curPoint].size());
1064
1065         forAll (curFaceList, faceI)
1066         {
1067             curFaceList[faceI] = faces[pointFaces[curPoint][faceI]];
1068         }
1069
1070         primitiveFacePatch curPatch (curFaceList, points);
1071
1072         labelList agglomFacePoints = curPatch.edgeLoops()[0];
1073
1074         SLList<label> slList;
1075
1076         label curPointLabel = -1;
1077
1078         for (label i=0; i<agglomFacePoints.size(); i++)
1079         {
1080             if (curPatch.meshPoints()[agglomFacePoints[i]] == curPoint)
1081             {
1082                 curPointLabel = i;
1083             }
1084             else if ( curPointLabel != -1 )
1085             {
1086                 slList.append(curPatch.meshPoints()[agglomFacePoints[i]]);
1087             }
1088         }
1089
1090         for (label i=0; i<curPointLabel; i++)
1091         {
1092             slList.append(curPatch.meshPoints()[agglomFacePoints[i]]);
1093         }
1094
1095         labelList curPointPoints(slList);
1096
1097         for (label i=0; i < (curPointPoints.size() - 1); i++)
1098         {
1099             vector d1 = points[curPointPoints[i]] - points[curPoint];
1100
1101             vector d2 = points[curPointPoints[i + 1]] - points[curPoint];
1102
1103             vector n = (d1^d2)/mag(d1 ^ d2);
1104
1105             scalar sinAlpha = mag(d1 ^ d2)/(mag(d1)*mag(d2));
1106
1107             scalar w = sinAlpha/(mag(d1)*mag(d2));
1108
1109             result[curPoint] += w*n;
1110         }
1111     }
1112
1113     // Correcte wedge points
1114     forAll (boundary(), patchI)
1115     {
1116         if (boundary()[patchI].type() == wedgeFaPatch::typeName)
1117         {
1118             const wedgeFaPatch& wedgePatch =
1119                 refCast<const wedgeFaPatch>(boundary()[patchI]);
1120
1121             labelList patchPoints = wedgePatch.pointLabels();
1122
1123             vector N =
1124                 transform
1125                 (
1126                     wedgePatch.edgeT(),
1127                     wedgePatch.centreNormal()
1128                 );
1129
1130             N /= mag(N);
1131
1132             forAll (patchPoints, pointI)
1133             {
1134                 result[patchPoints[pointI]]
1135                     -= N*(N&result[patchPoints[pointI]]);
1136             }
1137         }
1138     }
1139
1140     // Axis point correction
1141     forAll (boundary(), patchI)
1142     {
1143         if (boundary()[patchI].type() == wedgeFaPatch::typeName)
1144         {
1145             const wedgeFaPatch& wedgePatch =
1146                 refCast<const wedgeFaPatch>(boundary()[patchI]);
1147
1148             if (wedgePatch.axisPoint() > -1)
1149             {
1150                 result[wedgePatch.axisPoint()] =
1151                     wedgePatch.axis()
1152                    *(
1153                        wedgePatch.axis()
1154                       &result[wedgePatch.axisPoint()]
1155                     );
1156             }
1157
1158             break;
1159         }
1160     }
1161
1162     // Boundary points correction
1163     forAll (boundary(), patchI)
1164     {
1165         if (correctPatchPointNormals(patchI) && !boundary()[patchI].coupled())
1166         {
1167             if (boundary()[patchI].ngbPolyPatchIndex() == -1)
1168             {
1169                 FatalErrorIn
1170                     (
1171                         "void faMesh::calcPointAreaNormals const"
1172                     )   << "Neighbour polyPatch index is not defined "
1173                         << "for faPatch " << boundary()[patchI].name()
1174                         << abort(FatalError);
1175             }
1176
1177             labelList patchPoints = boundary()[patchI].pointLabels();
1178             vectorField N = boundary()[patchI].ngbPolyPatchPointNormals();
1179
1180             forAll (patchPoints, pointI)
1181             {
1182                 result[patchPoints[pointI]]
1183                     -= N[pointI]*(N[pointI]&result[patchPoints[pointI]]);
1184             }
1185         }
1186     }
1187
1188     // Processor patch points correction
1189     forAll (boundary(), patchI)
1190     {
1191         if(boundary()[patchI].type() == processorFaPatch::typeName)
1192         {
1193             const processorFaPatch& procPatch =
1194                 refCast<const processorFaPatch>(boundary()[patchI]);
1195
1196             labelList patchPointLabels = procPatch.pointLabels();
1197
1198             vectorField patchPointNormals
1199             (
1200                 patchPointLabels.size(),
1201                 vector::zero
1202             );
1203
1204             forAll (patchPointNormals, pointI)
1205             {
1206                 patchPointNormals[pointI] =
1207                     result[patchPointLabels[pointI]];
1208             }
1209
1210             {
1211             OPstream::write
1212             (
1213                 Pstream::blocking,
1214                 procPatch.neighbProcNo(),
1215                 reinterpret_cast<const char*>(patchPointNormals.begin()),
1216                 patchPointNormals.byteSize()
1217             );
1218             }
1219
1220             vectorField ngbPatchPointNormals
1221             (
1222                 procPatch.neighbPoints().size(),
1223                 vector::zero
1224             );
1225
1226             {
1227             IPstream::read
1228             (
1229                 Pstream::blocking,
1230                 procPatch.neighbProcNo(),
1231                 reinterpret_cast<char*>(ngbPatchPointNormals.begin()),
1232                 ngbPatchPointNormals.byteSize()
1233             );
1234             }
1235
1236             const labelList& nonGlobalPatchPoints =
1237                 procPatch.nonGlobalPatchPoints();
1238
1239             forAll (nonGlobalPatchPoints, pointI)
1240             {
1241                 result[patchPointLabels[nonGlobalPatchPoints[pointI]]] +=
1242                     ngbPatchPointNormals
1243                     [
1244                         procPatch.neighbPoints()[nonGlobalPatchPoints[pointI]]
1245                     ];
1246             }
1247         }
1248     }
1249
1250
1251     // Correct global points
1252
1253     if (globalData().nGlobalPoints() > 0)
1254     {
1255         const labelList& spLabels =
1256             globalData().sharedPointLabels();
1257
1258         vectorField spNormals(spLabels.size(), vector::zero);
1259         forAll (spNormals, pointI)
1260         {
1261             spNormals[pointI] = result[spLabels[pointI]];
1262         }
1263
1264         const labelList& addr = globalData().sharedPointAddr();
1265
1266         vectorField gpNormals
1267         (
1268             globalData().nGlobalPoints(),
1269             vector::zero
1270         );
1271
1272         forAll (addr, i)
1273         {
1274             gpNormals[addr[i]] += spNormals[i];
1275         }
1276
1277         combineReduce(gpNormals, plusEqOp<vectorField>());
1278
1279         // Extract local data
1280         forAll (addr, i)
1281         {
1282             spNormals[i] = gpNormals[addr[i]];
1283         }
1284
1285         forAll (spNormals, pointI)
1286         {
1287             result[spLabels[pointI]] = spNormals[pointI];
1288         }
1289     }
1290
1291     result /= mag(result);
1292 }
1293
1294
1295 void faMesh::calcPointAreaNormalsByQuadricsFit() const
1296 {
1297     vectorField& result = *pointAreaNormalsPtr_;
1298
1299
1300     labelList intPoints = internalPoints();
1301     labelList bndPoints = boundaryPoints();
1302
1303     const pointField& points = patch().localPoints();
1304     const faceList& faces = patch().localFaces();
1305     const labelListList& pointFaces = patch().pointFaces();
1306
1307     forAll(intPoints, pointI)
1308     {
1309         label curPoint = intPoints[pointI];
1310
1311         labelHashSet faceSet;
1312         forAll(pointFaces[curPoint], faceI)
1313         {
1314             faceSet.insert(pointFaces[curPoint][faceI]);
1315         }
1316         labelList curFaces = faceSet.toc();
1317
1318         labelHashSet pointSet;
1319
1320         pointSet.insert(curPoint);
1321         for(label i=0; i<curFaces.size(); i++)
1322         {
1323             const labelList& facePoints = faces[curFaces[i]];
1324             for(label j=0; j<facePoints.size(); j++)
1325             {
1326                 if(!pointSet.found(facePoints[j]))
1327                 {
1328                     pointSet.insert(facePoints[j]);
1329                 }
1330             }
1331         }
1332         pointSet.erase(curPoint);
1333         labelList curPoints = pointSet.toc();
1334
1335         if (curPoints.size() < 5)
1336         {
1337             if (debug)
1338             {
1339                 Info << "WARNING: Extending point set for fitting." << endl;
1340             }
1341
1342             labelHashSet faceSet;
1343             forAll(pointFaces[curPoint], faceI)
1344             {
1345                 faceSet.insert(pointFaces[curPoint][faceI]);
1346             }
1347             labelList curFaces = faceSet.toc();
1348             forAll(curFaces, faceI)
1349             {
1350                 const labelList& curFaceFaces =
1351                     patch().faceFaces()[curFaces[faceI]];
1352
1353                 forAll(curFaceFaces, fI)
1354                 {
1355                     label curFaceFace = curFaceFaces[fI];
1356
1357                     if(!faceSet.found(curFaceFace))
1358                     {
1359                         faceSet.insert(curFaceFace);
1360                     }
1361                 }
1362             }
1363             curFaces = faceSet.toc();
1364
1365             labelHashSet pointSet;
1366
1367             pointSet.insert(curPoint);
1368             for(label i=0; i<curFaces.size(); i++)
1369             {
1370                 const labelList& facePoints = faces[curFaces[i]];
1371                 for(label j=0; j<facePoints.size(); j++)
1372                 {
1373                     if(!pointSet.found(facePoints[j]))
1374                     {
1375                         pointSet.insert(facePoints[j]);
1376                     }
1377                 }
1378             }
1379
1380             pointSet.erase(curPoint);
1381             curPoints = pointSet.toc();
1382         }
1383
1384         vectorField allPoints(curPoints.size());
1385         scalarField W(curPoints.size(), 1.0);
1386         for(label i=0; i<curPoints.size(); i++)
1387         {
1388             allPoints[i] = points[curPoints[i]];
1389             W[i] = 1.0/magSqr(allPoints[i] - points[curPoint]);
1390         }
1391
1392         // Transforme points
1393         vector origin = points[curPoint];
1394         vector axis = result[curPoint]/mag(result[curPoint]);
1395         vector dir = (allPoints[0] - points[curPoint]);
1396         dir -= axis*(axis&dir);
1397         dir /= mag(dir);
1398         coordinateSystem cs("cs", origin, axis, dir);
1399
1400         forAll(allPoints, pI)
1401         {
1402             allPoints[pI] = cs.localPosition(allPoints[pI]);
1403         }
1404
1405         scalarRectangularMatrix M
1406         (
1407             allPoints.size(),
1408             5,
1409             0.0
1410         );
1411
1412         for(label i = 0; i < allPoints.size(); i++)
1413         {
1414             M[i][0] = sqr(allPoints[i].x());
1415             M[i][1] = sqr(allPoints[i].y());
1416             M[i][2] = allPoints[i].x()*allPoints[i].y();
1417             M[i][3] = allPoints[i].x();
1418             M[i][4] = allPoints[i].y();
1419         }
1420
1421         scalarSquareMatrix MtM(5, 0.0);
1422
1423         for (label i = 0; i < MtM.n(); i++)
1424         {
1425             for (label j = 0; j < MtM.m(); j++)
1426             {
1427                 for (label k = 0; k < M.n(); k++)
1428                 {
1429                     MtM[i][j] += M[k][i]*M[k][j]*W[k];
1430                 }
1431             }
1432         }
1433
1434         scalarField MtR(5, 0);
1435
1436         for (label i=0; i<MtR.size(); i++)
1437         {
1438             for (label j=0; j<M.n(); j++)
1439             {
1440                 MtR[i] += M[j][i]*allPoints[j].z()*W[j];
1441             }
1442         }
1443
1444         scalarSquareMatrix::LUsolve(MtM, MtR);
1445
1446         vector curNormal = vector(MtR[3], MtR[4], -1);
1447
1448         curNormal = cs.globalVector(curNormal);
1449
1450         curNormal *= sign(curNormal&result[curPoint]);
1451
1452         result[curPoint] = curNormal;
1453     }
1454
1455
1456     forAll (boundary(), patchI)
1457     {
1458         if(boundary()[patchI].type() == processorFaPatch::typeName)
1459         {
1460             const processorFaPatch& procPatch =
1461                 refCast<const processorFaPatch>(boundary()[patchI]);
1462
1463             labelList patchPointLabels = procPatch.pointLabels();
1464
1465             labelList toNgbProcLsPointStarts(patchPointLabels.size(), 0);
1466             vectorField toNgbProcLsPoints
1467             (
1468                 10*patchPointLabels.size(),
1469                 vector::zero
1470             );
1471             label nPoints = 0;
1472
1473             for (label pointI=0; pointI<patchPointLabels.size(); pointI++)
1474             {
1475                 label curPoint = patchPointLabels[pointI];
1476
1477                 toNgbProcLsPointStarts[pointI] = nPoints;
1478
1479                 labelHashSet faceSet;
1480                 forAll(pointFaces[curPoint], faceI)
1481                 {
1482                     faceSet.insert(pointFaces[curPoint][faceI]);
1483                 }
1484                 labelList curFaces = faceSet.toc();
1485
1486                 labelHashSet pointSet;
1487
1488                 pointSet.insert(curPoint);
1489                 for (label i=0; i<curFaces.size(); i++)
1490                 {
1491                     const labelList& facePoints = faces[curFaces[i]];
1492                     for (label j=0; j<facePoints.size(); j++)
1493                     {
1494                         if(!pointSet.found(facePoints[j]))
1495                         {
1496                             pointSet.insert(facePoints[j]);
1497                         }
1498                     }
1499                 }
1500                 pointSet.erase(curPoint);
1501                 labelList curPoints = pointSet.toc();
1502
1503                 for (label i=0; i<curPoints.size(); i++)
1504                 {
1505                     toNgbProcLsPoints[nPoints++] =
1506                         points[curPoints[i]];
1507                 }
1508             }
1509
1510             toNgbProcLsPoints.setSize(nPoints);
1511
1512             {
1513                 OPstream toNeighbProc
1514                 (
1515                     Pstream::blocking,
1516                     procPatch.neighbProcNo(),
1517                     toNgbProcLsPoints.byteSize()
1518                   + toNgbProcLsPointStarts.byteSize()
1519                   + 10*sizeof(label)
1520                 );
1521
1522                 toNeighbProc << toNgbProcLsPoints
1523                     << toNgbProcLsPointStarts;
1524             }
1525         }
1526     }
1527
1528     forAll (boundary(), patchI)
1529     {
1530         if(boundary()[patchI].type() == processorFaPatch::typeName)
1531         {
1532             const processorFaPatch& procPatch =
1533                 refCast<const processorFaPatch>(boundary()[patchI]);
1534
1535             labelList patchPointLabels = procPatch.pointLabels();
1536
1537             labelList fromNgbProcLsPointStarts(patchPointLabels.size(), 0);
1538             vectorField fromNgbProcLsPoints;
1539
1540             {
1541                 IPstream fromNeighbProc
1542                 (
1543                     Pstream::blocking,
1544                     procPatch.neighbProcNo(),
1545                     10*patchPointLabels.size()*sizeof(vector)
1546                   + fromNgbProcLsPointStarts.byteSize()
1547                   + 10*sizeof(label)
1548                 );
1549
1550                 fromNeighbProc >> fromNgbProcLsPoints
1551                     >> fromNgbProcLsPointStarts;
1552             }
1553
1554             const labelList& nonGlobalPatchPoints =
1555                 procPatch.nonGlobalPatchPoints();
1556
1557             forAll(nonGlobalPatchPoints, pointI)
1558             {
1559                 label curPoint =
1560                     patchPointLabels[nonGlobalPatchPoints[pointI]];
1561                 label curNgbPoint =
1562                     procPatch.neighbPoints()[nonGlobalPatchPoints[pointI]];
1563
1564                 labelHashSet faceSet;
1565                 forAll(pointFaces[curPoint], faceI)
1566                 {
1567                     faceSet.insert(pointFaces[curPoint][faceI]);
1568                 }
1569                 labelList curFaces = faceSet.toc();
1570
1571                 labelHashSet pointSet;
1572
1573                 pointSet.insert(curPoint);
1574                 for(label i=0; i<curFaces.size(); i++)
1575                 {
1576                     const labelList& facePoints = faces[curFaces[i]];
1577                     for(label j=0; j<facePoints.size(); j++)
1578                     {
1579                         if(!pointSet.found(facePoints[j]))
1580                         {
1581                             pointSet.insert(facePoints[j]);
1582                         }
1583                     }
1584                 }
1585                 pointSet.erase(curPoint);
1586                 labelList curPoints = pointSet.toc();
1587
1588                 label nAllPoints = curPoints.size();
1589
1590                 if (curNgbPoint == fromNgbProcLsPointStarts.size() - 1)
1591                 {
1592                     nAllPoints +=
1593                         fromNgbProcLsPoints.size()
1594                       - fromNgbProcLsPointStarts[curNgbPoint];
1595                 }
1596                 else
1597                 {
1598                     nAllPoints +=
1599                         fromNgbProcLsPointStarts[curNgbPoint + 1]
1600                       - fromNgbProcLsPointStarts[curNgbPoint];
1601                 }
1602
1603                 vectorField allPointsExt(nAllPoints);
1604                 label counter = 0;
1605                 for (label i=0; i<curPoints.size(); i++)
1606                 {
1607                     allPointsExt[counter++] = points[curPoints[i]];
1608                 }
1609
1610                 if (curNgbPoint == fromNgbProcLsPointStarts.size() - 1)
1611                 {
1612                     for
1613                     (
1614                         label i=fromNgbProcLsPointStarts[curNgbPoint];
1615                         i<fromNgbProcLsPoints.size();
1616                         i++
1617                     )
1618                     {
1619                         allPointsExt[counter++] = fromNgbProcLsPoints[i];
1620                     }
1621                 }
1622                 else
1623                 {
1624                     for
1625                     (
1626                         label i=fromNgbProcLsPointStarts[curNgbPoint];
1627                         i<fromNgbProcLsPointStarts[curNgbPoint+1];
1628                         i++
1629                     )
1630                     {
1631                         allPointsExt[counter++] = fromNgbProcLsPoints[i];
1632                     }
1633                 }
1634
1635                 // Remove duplicate points
1636                 vectorField allPoints(nAllPoints, vector::zero);
1637                 boundBox bb(allPointsExt, false);
1638                 scalar tol = 0.001*mag(bb.max() - bb.min());
1639
1640                 nAllPoints = 0;
1641                 forAll(allPointsExt, pI)
1642                 {
1643                     bool duplicate = false;
1644                     for (label i=0; i<nAllPoints; i++)
1645                     {
1646                         if
1647                         (
1648                             mag
1649                             (
1650                                 allPoints[i]
1651                               - allPointsExt[pI]
1652                             )
1653                           < tol
1654                         )
1655                         {
1656                             duplicate = true;
1657                             break;
1658                         }
1659                     }
1660
1661                     if (!duplicate)
1662                     {
1663                         allPoints[nAllPoints++] =
1664                             allPointsExt[pI];
1665                     }
1666                 }
1667
1668                 allPoints.setSize(nAllPoints);
1669
1670                 if (nAllPoints < 5)
1671                 {
1672                     FatalErrorIn
1673                     (
1674                         "void faMesh::calcPointAreaNormals() const"
1675                     )   << "There are no enough points for quadrics "
1676                         << "fitting for a point at processor patch"
1677                         << abort(FatalError);
1678                 }
1679
1680                 // Transforme points
1681                 vector origin = points[curPoint];
1682                 vector axis = result[curPoint]/mag(result[curPoint]);
1683                 vector dir = (allPoints[0] - points[curPoint]);
1684                 dir -= axis*(axis&dir);
1685                 dir /= mag(dir);
1686                 coordinateSystem cs("cs", origin, axis, dir);
1687
1688                 scalarField W(allPoints.size(), 1.0);
1689
1690                 forAll(allPoints, pI)
1691                 {
1692                     W[pI] = 1.0/magSqr(allPoints[pI] - points[curPoint]);
1693
1694                     allPoints[pI] =
1695                         cs.localPosition(allPoints[pI]);
1696                 }
1697
1698                 scalarRectangularMatrix M
1699                 (
1700                     allPoints.size(),
1701                     5,
1702                     0.0
1703                 );
1704
1705                 for(label i=0; i<allPoints.size(); i++)
1706                 {
1707                     M[i][0] = sqr(allPoints[i].x());
1708                     M[i][1] = sqr(allPoints[i].y());
1709                     M[i][2] = allPoints[i].x()*allPoints[i].y();
1710                     M[i][3] = allPoints[i].x();
1711                     M[i][4] = allPoints[i].y();
1712                 }
1713
1714                 scalarSquareMatrix MtM(5, 0.0);
1715
1716                 for (label i = 0; i < MtM.n(); i++)
1717                 {
1718                     for (label j = 0; j < MtM.m(); j++)
1719                     {
1720                         for (label k = 0; k < M.n(); k++)
1721                         {
1722                             MtM[i][j] += M[k][i]*M[k][j]*W[k];
1723                         }
1724                     }
1725                 }
1726
1727                 scalarField MtR(5, 0);
1728
1729                 for (label i = 0; i < MtR.size(); i++)
1730                 {
1731                     for (label j = 0; j < M.n(); j++)
1732                     {
1733                         MtR[i] += M[j][i]*allPoints[j].z()*W[j];
1734                     }
1735                 }
1736
1737                 scalarSquareMatrix::LUsolve(MtM, MtR);
1738
1739                 vector curNormal = vector(MtR[3], MtR[4], -1);
1740
1741                 curNormal = cs.globalVector(curNormal);
1742
1743                 curNormal *= sign(curNormal&result[curPoint]);
1744
1745                 result[curPoint] = curNormal;
1746             }
1747         }
1748     }
1749
1750     // Correct global points
1751     if (globalData().nGlobalPoints() > 0)
1752     {
1753         const labelList& spLabels =
1754             globalData().sharedPointLabels();
1755
1756         const labelList& addr = globalData().sharedPointAddr();
1757
1758         for (label k=0; k<globalData().nGlobalPoints(); k++)
1759         {
1760             List<List<vector> > procLsPoints(Pstream::nProcs());
1761
1762             label curSharedPointIndex = findIndex(addr, k);
1763
1764             scalar tol = 0.0;
1765
1766             if (curSharedPointIndex != -1)
1767             {
1768                 label curPoint = spLabels[curSharedPointIndex];
1769
1770                 labelHashSet faceSet;
1771                 forAll(pointFaces[curPoint], faceI)
1772                 {
1773                     faceSet.insert(pointFaces[curPoint][faceI]);
1774                 }
1775                 labelList curFaces = faceSet.toc();
1776
1777                 labelHashSet pointSet;
1778                 pointSet.insert(curPoint);
1779                 for (label i=0; i<curFaces.size(); i++)
1780                 {
1781                     const labelList& facePoints = faces[curFaces[i]];
1782                     for (label j=0; j<facePoints.size(); j++)
1783                     {
1784                         if (!pointSet.found(facePoints[j]))
1785                         {
1786                             pointSet.insert(facePoints[j]);
1787                         }
1788                     }
1789                 }
1790                 pointSet.erase(curPoint);
1791                 labelList curPoints = pointSet.toc();
1792
1793                 vectorField locPoints(points, curPoints);
1794
1795                 procLsPoints[Pstream::myProcNo()] = locPoints;
1796
1797                 boundBox bb(locPoints, false);
1798                 tol = 0.001*mag(bb.max() - bb.min());
1799             }
1800
1801             Pstream::gatherList(procLsPoints);
1802             Pstream::scatterList(procLsPoints);
1803
1804             if (curSharedPointIndex != -1)
1805             {
1806                 label curPoint = spLabels[curSharedPointIndex];
1807
1808                 label nAllPoints = 0;
1809                 forAll(procLsPoints, procI)
1810                 {
1811                     nAllPoints += procLsPoints[procI].size();
1812                 }
1813
1814                 vectorField allPoints(nAllPoints, vector::zero);
1815
1816                 nAllPoints = 0;
1817                 forAll(procLsPoints, procI)
1818                 {
1819                     forAll(procLsPoints[procI], pointI)
1820                     {
1821                         bool duplicate = false;
1822                         for (label i=0; i<nAllPoints; i++)
1823                         {
1824                             if
1825                             (
1826                                 mag
1827                                 (
1828                                     allPoints[i]
1829                                   - procLsPoints[procI][pointI]
1830                                 )
1831                               < tol
1832                             )
1833                             {
1834                                 duplicate = true;
1835                                 break;
1836                             }
1837                         }
1838
1839                         if (!duplicate)
1840                         {
1841                             allPoints[nAllPoints++] =
1842                                 procLsPoints[procI][pointI];
1843                         }
1844                     }
1845                 }
1846
1847                 allPoints.setSize(nAllPoints);
1848
1849                 if (nAllPoints < 5)
1850                 {
1851                     FatalErrorIn
1852                     (
1853                         "void faMesh::calcPointAreaNormals() const"
1854                     )   << "There are no enough points for quadrics "
1855                         << "fitting for a global processor point "
1856                         << abort(FatalError);
1857                 }
1858
1859                 // Transforme points
1860                 vector origin = points[curPoint];
1861                 vector axis = result[curPoint]/mag(result[curPoint]);
1862                 vector dir = (allPoints[0] - points[curPoint]);
1863                 dir -= axis*(axis&dir);
1864                 dir /= mag(dir);
1865                 coordinateSystem cs("cs", origin, axis, dir);
1866
1867                 scalarField W(allPoints.size(), 1.0);
1868
1869                 forAll(allPoints, pointI)
1870                 {
1871                     W[pointI]=
1872                         1.0/magSqr(allPoints[pointI] - points[curPoint]);
1873
1874                     allPoints[pointI] =
1875                         cs.localPosition(allPoints[pointI]);
1876                 }
1877
1878                 scalarRectangularMatrix M
1879                 (
1880                     allPoints.size(),
1881                     5,
1882                     0.0
1883                 );
1884
1885                 for (label i=0; i<allPoints.size(); i++)
1886                 {
1887                     M[i][0] = sqr(allPoints[i].x());
1888                     M[i][1] = sqr(allPoints[i].y());
1889                     M[i][2] = allPoints[i].x()*allPoints[i].y();
1890                     M[i][3] = allPoints[i].x();
1891                     M[i][4] = allPoints[i].y();
1892                 }
1893
1894                 scalarSquareMatrix MtM(5, 0.0);
1895                 for (label i = 0; i < MtM.n(); i++)
1896                 {
1897                     for (label j = 0; j < MtM.m(); j++)
1898                     {
1899                         for (label k = 0; k < M.n(); k++)
1900                         {
1901                             MtM[i][j] += M[k][i]*M[k][j]*W[k];
1902                         }
1903                     }
1904                 }
1905
1906                 scalarField MtR(5, 0);
1907                 for (label i = 0; i < MtR.size(); i++)
1908                 {
1909                     for (label j = 0; j < M.n(); j++)
1910                     {
1911                         MtR[i] += M[j][i]*allPoints[j].z()*W[j];
1912                     }
1913                 }
1914
1915                 scalarSquareMatrix::LUsolve(MtM, MtR);
1916
1917                 vector curNormal = vector(MtR[3], MtR[4], -1);
1918
1919                 curNormal = cs.globalVector(curNormal);
1920
1921                 curNormal *= sign(curNormal&result[curPoint]);
1922
1923                 result[curPoint] = curNormal;
1924             }
1925         }
1926     }
1927
1928     result /= mag(result);
1929 }
1930
1931
1932 tmp<edgeScalarField> faMesh::edgeLengthCorrection() const
1933 {
1934     if (debug)
1935     {
1936         Info<< "tmp<edgeScalarField> faMesh::edgeLengthCorrection() const : "
1937             << "Calculating edge length correction" << endl;
1938     }
1939
1940     tmp<edgeScalarField> tcorrection
1941     (
1942         new edgeScalarField
1943         (
1944             IOobject
1945             (
1946                 "edgeLengthCorrection",
1947                 mesh_.pointsInstance(),
1948                 meshSubDir,
1949                 mesh_
1950             ),
1951             *this,
1952             dimless
1953         )
1954     );
1955     edgeScalarField& correction = tcorrection();
1956
1957
1958     const vectorField& pointNormals = pointAreaNormals();
1959
1960
1961     forAll (correction.internalField(), edgeI)
1962     {
1963         scalar sinAlpha = mag
1964         (
1965             pointNormals[edges()[edgeI].start()]^
1966             pointNormals[edges()[edgeI].end()]
1967         );
1968
1969         scalar alpha = asin(sinAlpha);
1970
1971         correction.internalField()[edgeI] = cos(alpha/2.0);
1972     }
1973
1974
1975     forAll (boundary(), patchI)
1976     {
1977         const edgeList::subList patchEdges =
1978             boundary()[patchI].patchSlice(edges());
1979
1980         forAll (patchEdges, edgeI)
1981         {
1982             scalar sinAlpha = mag
1983                 (
1984                     pointNormals[patchEdges[edgeI].start()]^
1985                     pointNormals[patchEdges[edgeI].end()]
1986                 );
1987
1988             scalar alpha = asin(sinAlpha);
1989
1990             correction.boundaryField()[patchI][edgeI] = cos(alpha/2.0);
1991         }
1992     }
1993
1994
1995         return tcorrection;
1996 }
1997
1998
1999 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2000
2001 } // End namespace Foam
2002
2003 // ************************************************************************* //