Foam::solverPerformance Foam::paralution_PBiCG::solve
(
    scalarField& psi,
    const scalarField& source,
    const direction cmpt
) const
{

    word precond_name = lduMatrix::preconditioner::getName(controlDict_);
    double div   = controlDict_.lookupOrDefault<double>("div", 1e+08);
    bool accel   = controlDict_.lookupOrDefault<bool>("useAccelerator", true);
    word mformat = controlDict_.lookupOrDefault<word>("MatrixFormat", "CSR");
    word pformat = controlDict_.lookupOrDefault<word>("PrecondFormat", "CSR");
    int ILUp     = controlDict_.lookupOrDefault<int>("ILUp", 0);
    int ILUq     = controlDict_.lookupOrDefault<int>("ILUq", 1);
    int MEp      = controlDict_.lookupOrDefault<int>("MEp", 1);
    word LBPre   = controlDict_.lookupOrDefault<word>("LastBlockPrecond", "paralution_Jacobi");

    solverPerformance solverPerf(typeName + '(' + precond_name + ')', fieldName_);

    register label nCells = psi.size();

    scalarField pA(nCells);
    scalarField wA(nCells);

    // --- Calculate A.psi
    matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);

    // --- Calculate initial residual field
    scalarField rA(source - wA);

    // --- Calculate normalisation factor
    scalar normFactor = this->normFactor(psi, source, wA, pA);

    // --- Calculate normalised residual norm
    solverPerf.initialResidual() = gSumMag(rA)/normFactor;
    solverPerf.finalResidual() = solverPerf.initialResidual();

    if (!solverPerf.checkConvergence(tolerance_, relTol_)) {

      paralution::_matrix_format mf = paralution::CSR;
      if      (mformat == "CSR")   mf = paralution::CSR;
      else if (mformat == "DIA")   mf = paralution::DIA;
      else if (mformat == "HYB")   mf = paralution::HYB;
      else if (mformat == "ELL")   mf = paralution::ELL;
      else if (mformat == "MCSR")  mf = paralution::MCSR;
      else if (mformat == "BCSR")  mf = paralution::BCSR;
      else if (mformat == "COO")   mf = paralution::COO;
      else if (mformat == "DENSE") mf = paralution::DENSE;

      paralution::LocalVector<double> x;
      paralution::LocalVector<double> rhs;
      paralution::LocalMatrix<double> mat;

      paralution::BiCGStab<paralution::LocalMatrix<double>,
                     paralution::LocalVector<double>,
                     double> ls;

      paralution::import_openfoam_matrix(matrix(), &mat);
      paralution::import_openfoam_vector(source, &rhs);
      paralution::import_openfoam_vector(psi, &x);

      ls.Clear();

      paralution::Preconditioner<paralution::LocalMatrix<double>,
                                 paralution::LocalVector<double>,
                                 double > *precond = NULL;

      precond = paralution::GetPreconditioner<double>(precond_name, LBPre, pformat, ILUp, ILUq, MEp);
      if (precond != NULL) ls.SetPreconditioner(*precond);

      ls.SetOperator(mat);

      // Switch to L1 norm to be consistent with OpenFOAM solvers
      ls.SetResidualNorm(1);

      ls.Init(tolerance_*normFactor, // abs
              relTol_,    // rel
              div,        // div
              maxIter_);  // max iter

      if (accel) {
        mat.MoveToAccelerator();
        rhs.MoveToAccelerator();
        x.MoveToAccelerator();
      }

      ls.Build();

      switch(mf) {
        case paralution::DENSE:
          mat.ConvertToDENSE();
          break;
        case paralution::CSR:
          mat.ConvertToCSR();
          break;
        case paralution::MCSR:
          mat.ConvertToMCSR();
          break;
        case paralution::BCSR:
//.........这里部分代码省略.........

Foam::lduSolverPerformance Foam::PCG::solve
(
    scalarField& x,
    const scalarField& b,
    const direction cmpt
) const
{
    // --- Setup class containing solver performance data
    lduSolverPerformance solverPerf
    (
        lduMatrix::preconditioner::getName(dict()) + typeName,
        fieldName()
    );

    register label nCells = x.size();

    scalar* __restrict__ xPtr = x.begin();

    scalarField pA(nCells);
    scalar* __restrict__ pAPtr = pA.begin();

    scalarField wA(nCells);
    scalar* __restrict__ wAPtr = wA.begin();

    // Calculate A.x
    matrix_.Amul(wA, x, coupleBouCoeffs_, interfaces_, cmpt);

    // Calculate initial residual field
    scalarField rA(b - wA);
    scalar* __restrict__ rAPtr = rA.begin();

    // Calculate normalisation factor
    scalar normFactor = this->normFactor(x, b, wA, pA, cmpt);

    if (lduMatrix::debug >= 2)
    {
        Info<< "   Normalisation factor = " << normFactor << endl;
    }

    // Calculate normalised residual norm
    solverPerf.initialResidual() = gSumMag(rA)/normFactor;
    solverPerf.finalResidual() = solverPerf.initialResidual();

    // Check convergence, solve if not converged
    if (!stop(solverPerf))
    {
        scalar wArA = matrix_.great_;
        scalar wArAold = wArA;

        // Select and construct the preconditioner
        autoPtr<lduPreconditioner> preconPtr;

        preconPtr =
            lduPreconditioner::New
            (
                matrix_,
                coupleBouCoeffs_,
                coupleIntCoeffs_,
                interfaces_,
                dict()
            );

        // Solver iteration
        do
        {
            // Store previous wArA
            wArAold = wArA;

            // Precondition residual
            preconPtr->precondition(wA, rA, cmpt);

            // Update search directions:
            wArA = gSumProd(wA, rA);

            if (solverPerf.nIterations() == 0)
            {
                for (register label cell=0; cell<nCells; cell++)
                {
                    pAPtr[cell] = wAPtr[cell];
                }
            }
            else
            {
                scalar beta = wArA/wArAold;

                for (register label cell=0; cell<nCells; cell++)
                {
                    pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
                }
            }


            // Update preconditioned residual
            matrix_.Amul(wA, pA, coupleBouCoeffs_, interfaces_, cmpt);

            scalar wApA = gSumProd(wA, pA);


            // Test for singularity
            if (solverPerf.checkSingularity(mag(wApA)/normFactor)) break;
//.........这里部分代码省略.........

void Foam::GaussSeidelSmoother::smooth
(
    const word& fieldName_,
    scalarField& psi,
    const lduMatrix& matrix_,
    const scalarField& source,
    const FieldField<Field, scalar>& interfaceBouCoeffs_,
    const lduInterfaceFieldPtrsList& interfaces_,
    const direction cmpt,
    const label nSweeps
)
{

	//add by Xiaow:begin
	Foam::Time::enterSec("GaussSeidelSmoother");
	//add by Xiaow:end

    register scalar* __restrict__ psiPtr = psi.begin();

    register const label nCells = psi.size();

    scalarField bPrime(nCells);
    register scalar* __restrict__ bPrimePtr = bPrime.begin();

    register const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
    register const scalar* const __restrict__ upperPtr =
        matrix_.upper().begin();
    register const scalar* const __restrict__ lowerPtr =
        matrix_.lower().begin();

    register const label* const __restrict__ uPtr =
        matrix_.lduAddr().upperAddr().begin();

    register const label* const __restrict__ ownStartPtr =
        matrix_.lduAddr().ownerStartAddr().begin();


    // Parallel boundary initialisation.  The parallel boundary is treated
    // as an effective jacobi interface in the boundary.
    // Note: there is a change of sign in the coupled
    // interface update.  The reason for this is that the
    // internal coefficients are all located at the l.h.s. of
    // the matrix whereas the "implicit" coefficients on the
    // coupled boundaries are all created as if the
    // coefficient contribution is of a source-kind (i.e. they
    // have a sign as if they are on the r.h.s. of the matrix.
    // To compensate for this, it is necessary to turn the
    // sign of the contribution.

    FieldField<Field, scalar> mBouCoeffs(interfaceBouCoeffs_.size());

    forAll(mBouCoeffs, patchi)
    {
        if (interfaces_.set(patchi))
        {
            mBouCoeffs.set(patchi, -interfaceBouCoeffs_[patchi]);
        }
    }

	//add by Xiaow:begin
	Foam::label interid =Foam::Time::commProfiler_.enterIterSec();
	//add by Xiaow:end

    for (label sweep=0; sweep<nSweeps; sweep++)
    {
        bPrime = source;

        matrix_.initMatrixInterfaces
        (
            mBouCoeffs,
            interfaces_,
            psi,
            bPrime,
            cmpt
        );

        matrix_.updateMatrixInterfaces
        (
            mBouCoeffs,
            interfaces_,
            psi,
            bPrime,
            cmpt
        );

        register scalar curPsi;
        register label fStart;
        register label fEnd = ownStartPtr[0];

        for (register label cellI=0; cellI<nCells; cellI++)
        {
            // Start and end of this row
            fStart = fEnd;
            fEnd = ownStartPtr[cellI + 1];

            // Get the accumulated neighbour side
            curPsi = bPrimePtr[cellI];

            // Accumulate the owner product side
            for (register label curFace=fStart; curFace<fEnd; curFace++)
//.........这里部分代码省略.........

void Foam::meshRefinement::snapToSurface
(
    labelList& pointSurfaceRegion,
    labelList& edgeSurfaceRegion,
    scalarField& edgeWeight
)
{
    const edgeList& edges = mesh_.edges();
    const pointField& points = mesh_.points();

    pointSurfaceRegion.setSize(points.size());
    pointSurfaceRegion = -1;

    edgeSurfaceRegion.setSize(edges.size());
    edgeSurfaceRegion = -1;

    edgeWeight.setSize(edges.size());
    edgeWeight = -GREAT;


    // Do test for intersections
    // ~~~~~~~~~~~~~~~~~~~~~~~~~

    labelList surface1;
    List<pointIndexHit> hit1;
    labelList region1;
    //vectorField normal1;

    labelList surface2;
    List<pointIndexHit> hit2;
    labelList region2;
    //vectorField normal2;
    {
        vectorField start(edges.size());
        vectorField end(edges.size());
        forAll(edges, edgei)
        {
            const edge& e = edges[edgei];
            start[edgei] = points[e[0]];
            end[edgei] = points[e[1]];
        }
        
        surfaces_.findNearestIntersection
        (
            //labelList(1, 0),    //identity(surfaces_.surfaces().size()),
            identity(surfaces_.surfaces().size()),
            start,
            end,

            surface1,
            hit1,
            region1,
            //normal1,

            surface2,
            hit2,
            region2
            //normal2
        );
    }

    // Adjust location
    // ~~~~~~~~~~~~~~~

    pointField newPoints(points);
    label nAdjusted = 0;

    const labelListList& pointEdges = mesh_.pointEdges();
    forAll(pointEdges, pointi)
    {
        const point& pt = points[pointi];
        const labelList& pEdges = pointEdges[pointi];

        // Get the nearest intersection
        label minEdgei = -1;
        scalar minFraction = 0.5;   // Harpoon 0.25; // Samm?
        forAll(pEdges, pEdgei)
        {
            label edgei = pEdges[pEdgei];
            if (hit1[edgei].hit())
            {
                const point& hitPt = hit1[edgei].hitPoint();

                const edge& e = edges[edgei];
                label otherPointi = e.otherVertex(pointi);
                const point& otherPt = points[otherPointi];

                vector eVec(otherPt-pt);
                scalar f = eVec&(hitPt-pt)/magSqr(eVec);

                if (f < minFraction)
                {
                    minEdgei = edgei;
                    minFraction = f;
                }
            }
        }
        if (minEdgei != -1 && minFraction >= 0.01)
        {
            // Move point to intersection with minEdgei
//.........这里部分代码省略.........

// Call scotch with options from dictionary.
Foam::label Foam::ptscotchDecomp::decompose
(
    const fileName& meshPath,
    const label adjncySize,
    const label adjncy[],
    const label xadjSize,
    const label xadj[],
    const scalarField& cWeights,

    List<label>& finalDecomp
) const
{
    if (debug)
    {
        Pout<< "ptscotchDecomp : entering with xadj:" << xadjSize << endl;
    }

    // Dump graph
    if (decompositionDict_.found("scotchCoeffs"))
    {
        const dictionary& scotchCoeffs =
            decompositionDict_.subDict("scotchCoeffs");

        if (scotchCoeffs.lookupOrDefault("writeGraph", false))
        {
            OFstream str
            (
               meshPath + "_" + Foam::name(Pstream::myProcNo()) + ".dgr"
            );

            Pout<< "Dumping Scotch graph file to " << str.name() << endl
                << "Use this in combination with dgpart." << endl;

            globalIndex globalCells(xadjSize-1);

            // Distributed graph file (.grf)
            label version = 2;
            str << version << nl;
            // Number of files (procglbnbr)
            str << Pstream::nProcs();
            // My file number (procloc)
            str << ' ' << Pstream::myProcNo() << nl;

            // Total number of vertices (vertglbnbr)
            str << globalCells.size();
            // Total number of connections (edgeglbnbr)
            str << ' ' << returnReduce(xadj[xadjSize-1], sumOp<label>())
                << nl;
            // Local number of vertices (vertlocnbr)
            str << xadjSize-1;
            // Local number of connections (edgelocnbr)
            str << ' ' << xadj[xadjSize-1] << nl;
            // Numbering starts from 0
            label baseval = 0;
            // 100*hasVertlabels+10*hasEdgeWeights+1*hasVertWeighs
            str << baseval << ' ' << "000" << nl;
            for (label cellI = 0; cellI < xadjSize-1; cellI++)
            {
                label start = xadj[cellI];
                label end = xadj[cellI+1];
                str << end-start;

                for (label i = start; i < end; i++)
                {
                    str << ' ' << adjncy[i];
                }
                str << nl;
            }
        }
    }

    // Strategy
    // ~~~~~~~~

    // Default.
    SCOTCH_Strat stradat;
    check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");

    if (decompositionDict_.found("scotchCoeffs"))
    {
        const dictionary& scotchCoeffs =
            decompositionDict_.subDict("scotchCoeffs");


        string strategy;
        if (scotchCoeffs.readIfPresent("strategy", strategy))
        {
            if (debug)
            {
                Info<< "ptscotchDecomp : Using strategy " << strategy << endl;
            }
            SCOTCH_stratDgraphMap(&stradat, strategy.c_str());
            //fprintf(stdout, "S\tStrat=");
            //SCOTCH_stratSave(&stradat, stdout);
            //fprintf(stdout, "\n");
        }
    }


//.........这里部分代码省略.........