void parcelport::add_received_parcel(parcel p, std::size_t num_thread)
    {
        // do some work (notify event handlers)
        if(applier_)
        {
            while (threads::threadmanager_is(state_starting))
            {
                boost::this_thread::sleep(boost::get_system_time() +
                    boost::posix_time::milliseconds(HPX_NETWORK_RETRIES_SLEEP));
            }

            // Give up if we're shutting down.
            if (threads::threadmanager_is(state_stopping))
            {
    //             LPT_(debug) << "parcelport: add_received_parcel: dropping late "
    //                             "parcel " << p;
                return;
            }

            // write this parcel to the log
    //         LPT_(debug) << "parcelport: add_received_parcel: " << p;

            applier_->schedule_action(std::move(p));
        }
        // If the applier has not been set yet, we are in bootstrapping and
        // need to execute the action directly
        else
        {
            // TODO: Make assertions exceptions
            // decode the action-type in the parcel
            actions::base_action * act = p.get_action();

            // early parcels should only be plain actions
            HPX_ASSERT(actions::base_action::plain_action == act->get_action_type());

            // early parcels can't have continuations
            HPX_ASSERT(!p.get_continuation());

            // We should not allow any exceptions to escape the execution of the
            // action as this would bring down the ASIO thread we execute in.
            try {
                act->get_thread_function(0)
                    (threads::thread_state_ex(threads::wait_signaled));
            }
            catch (...) {
                hpx::report_error(boost::current_exception());
            }
        }
    }

void Foam::noAtomization::atomizeParcel
(
    parcel& p,
    const scalar deltaT,
    const vector& vel,
    const liquidMixtureProperties& fuels
) const
{
    p.liquidCore() = 0.0;
}

    void parcelport::put_parcel(parcel const& p, write_handler_type const& f)
    {
        typedef pending_parcels_map::iterator iterator;
        typedef pending_parcels_map::mapped_type mapped_type;

        naming::locality locality_id = p.get_destination_locality();
        naming::gid_type parcel_id = p.get_parcel_id();

        // enqueue the incoming parcel ...
        {
            lcos::local::spinlock::scoped_lock l(mtx_);

            mapped_type& e = pending_parcels_[locality_id];
            e.first.push_back(p);
            e.second.push_back(f);
        }

        get_connection_and_send_parcels(locality_id, parcel_id);
    }

void parcel_sent_handler(parcelhandler::write_handler_type & f,
                         boost::system::error_code const & ec, parcel const & p)
{
    // inform termination detection of a sent message
    if (!p.does_termination_detection())
    {
        hpx::detail::dijkstra_make_black();
    }

    // invoke the original handler
    f(ec, p);
}

    void parcelhandler::put_parcel(parcel& p, write_handler_type f)
    {
        // properly initialize parcel
        init_parcel(p);

        std::vector<naming::gid_type> const& gids = p.get_destinations();
        std::vector<naming::address>& addrs = p.get_destination_addrs();

        if (1 == gids.size()) {
            if (!addrs[0])
                resolver_.resolve(gids[0], addrs[0]);
        }
        else {
            boost::dynamic_bitset<> locals;
            resolver_.resolve(gids, addrs, locals);
        }

        if (!p.get_parcel_id())
            p.set_parcel_id(parcel::generate_unique_id());

        pp_.put_parcel(p, f);
    }

    ///////////////////////////////////////////////////////////////////////////
    // the code below is needed to bootstrap the parcel layer
    void parcelport::early_pending_parcel_handler(
        boost::system::error_code const& ec, parcel const & p)
    {
        if (ec) {
            // all errors during early parcel handling are fatal
            std::exception_ptr exception =
                HPX_GET_EXCEPTION(ec,
                    "early_pending_parcel_handler",
                    "error while handling early parcel: " +
                        ec.message() + "(" +
                        std::to_string(ec.value()) +
                        ")" + parcelset::dump_parcel(p));

            hpx::report_error(exception);
            return;
        }

#if defined(HPX_HAVE_APEX) && defined(HPX_HAVE_PARCEL_PROFILING)
        // tell APEX about the sent parcel
        apex::send(p.parcel_id().get_lsb(), p.size(),
            p.destination_locality_id());
#endif
    }

    void serialize_certificate(Archive& archive, Connection & connection,
        std::set<boost::uint32_t>& localities, parcel const& p)
    {
        // We send the certificate corresponding to the originating locality
        // of the parcel if this is the first message over this connection
        // or if the originating locality is not the current one.
        boost::uint32_t locality_id =
            naming::get_locality_id_from_gid(p.get_parcel_id());
        error_code ec(lightweight);
        boost::uint32_t this_locality_id = get_locality_id(ec);
        if (ec) {
            // this should only happen during bootstrap
            HPX_ASSERT(hpx::is_starting());
            this_locality_id = locality_id;
        }

        bool has_certificate = false;
        if ((connection.first_message_ || locality_id != this_locality_id) &&
            localities.find(locality_id) == localities.end())
        {
            // the first message must originate from this locality
            HPX_ASSERT(!connection.first_message_ || locality_id == this_locality_id);

            components::security::signed_certificate const& certificate =
                hpx::get_locality_certificate(locality_id, ec);

            if (!ec) {
                has_certificate = true;
                if (locality_id == this_locality_id)
                    connection.first_message_ = false;
                archive << has_certificate << certificate;

                // keep track of all certificates already prepended for this message
                localities.insert(locality_id);
            }
            else {
                // if the certificate is not available we have to still be on
                // the 'first' message (it's too early for a certificate)
                HPX_ASSERT(connection.first_message_);
                archive << has_certificate;
            }
        }
        else {
            archive << has_certificate;
        }
    }

 boost::shared_ptr<parcel_buffer_type> get_buffer(parcel const & p, std::size_t arg_size)
 {
     // generate the name for this data_buffer
     std::string data_buffer_name(p.get_parcel_id().to_string());
     if(!buffer_)
     {
         // clear and preallocate out_buffer_ (or fetch from cache)
         buffer_ = boost::make_shared<parcel_buffer_type>(
             get_data_buffer((arg_size * 12) / 10 + 1024,
                 data_buffer_name)
         );
     }
     else
     {
         buffer_->data_ =
             get_data_buffer((arg_size * 12) / 10 + 1024,
                 data_buffer_name);
     }
     return buffer_;
 }

        bool add_parcel(parcel const& p)
        {
            naming::gid_type id(p.get_parcel_id());

            // Add parcel to queue.
            {
                mutex_type::scoped_lock l(mtx_);
                std::pair<parcel_map_type::iterator, bool> ret =
                    parcels_.insert(parcel_map_type::value_type(id, p));

                if (!ret.second) {
                    HPX_THROW_EXCEPTION(bad_parameter,
                        "global_parcelhandler_queue::add_parcel",
                        "Could not add received parcel to the parcelhandler "
                        "queue");
                    return false;
                }
            }

            // do some work (notify event handlers)
            HPX_ASSERT(ph_ != 0);
            notify_(*ph_, id);
            return true;
        }

    void parcelport::put_parcel(parcel const& p, write_handler_type f)
    {
        typedef pending_parcels_map::iterator iterator;
        naming::locality locality_id = p.get_destination_locality();

        parcelport_connection_ptr client_connection;

        // enqueue the incoming parcel ...
        {
            util::spinlock::scoped_lock l(mtx_);
            pending_parcels_[locality_id].first.push_back(p);
            pending_parcels_[locality_id].second.push_back(f);
        }

        // Get a connection or reserve space for a new connection. 
        if (!connection_cache_.get_or_reserve(locality_id, client_connection))
            return;

        // Check if we need to create the new connection.
        if (!client_connection)
        {
            client_connection.reset(new parcelport_connection(
                io_service_pool_.get_io_service(), locality_id,
                connection_cache_, timer_, parcels_sent_));

            // Connect to the target locality, retry if needed.
            boost::system::error_code error = boost::asio::error::try_again;
            for (int i = 0; i < HPX_MAX_NETWORK_RETRIES; ++i)
            {
                try {
                    naming::locality::iterator_type end = locality_id.connect_end();
                    for (naming::locality::iterator_type it =
                            locality_id.connect_begin(io_service_pool_.get_io_service());
                         it != end; ++it)
                    {
                        client_connection->socket().close();
                        client_connection->socket().connect(*it, error);
                        if (!error)
                            break;
                    }
                    if (!error)
                        break;

                    // we wait for a really short amount of time
                    // TODO: Should this be an hpx::threads::suspend?
                    boost::this_thread::sleep(boost::get_system_time() +
                        boost::posix_time::milliseconds(HPX_NETWORK_RETRIES_SLEEP));
                }
                catch (boost::system::error_code const& e) {
                    HPX_THROW_EXCEPTION(network_error,
                        "parcelport::send_parcel", e.message());
                }
            }

            if (error) {
                client_connection->socket().close();

                hpx::util::osstream strm;
                strm << error.message() << " (while trying to connect to: "
                     << locality_id << ")";
                HPX_THROW_EXCEPTION(network_error,
                    "parcelport::send_parcel",
                    hpx::util::osstream_get_string(strm));
            }
#if defined(HPX_DEBUG)
            else {
                std::string connection_addr = 
                    client_connection->socket().remote_endpoint().address().to_string();
                boost::uint16_t connection_port = 
                    client_connection->socket().remote_endpoint().port();
                BOOST_ASSERT(locality_id.get_address() == connection_addr);
                BOOST_ASSERT(locality_id.get_port() == connection_port);
            }
#endif
        }
#if defined(HPX_DEBUG)
        else {
            //LPT_(info) << "parcelport: reusing existing connection to: "
            //           << addr.locality_;
            BOOST_ASSERT(locality_id == client_connection->destination());

            std::string connection_addr = client_connection->socket().remote_endpoint().address().to_string();
            boost::uint16_t connection_port = client_connection->socket().remote_endpoint().port();
            BOOST_ASSERT(locality_id.get_address() == connection_addr);
            BOOST_ASSERT(locality_id.get_port() == connection_port);
        }
#endif

        std::vector<parcel> parcels;
        std::vector<write_handler_type> handlers;

        util::spinlock::scoped_lock l(mtx_);
        iterator it = pending_parcels_.find(locality_id);

        if (it != pending_parcels_.end())
        {
            BOOST_ASSERT(it->first == locality_id);
            std::swap(parcels, it->second.first);
            std::swap(handlers, it->second.second);
        }
//.........这里部分代码省略.........

void ETAB::breakupParcel
(
    parcel& p,
    const scalar deltaT,
    const vector& Ug,
    const liquidMixture& fuels
) const
{

    scalar T  = p.T();
    scalar pc  = spray_.p()[p.cell()];
    scalar r  = 0.5*p.d();
    scalar r2 = r*r;
    scalar r3 = r*r2;

    scalar rho   = fuels.rho(pc, T, p.X());
    scalar sigma = fuels.sigma(pc, T, p.X());
    scalar mu    = fuels.mu(pc, T, p.X());

    // inverse of characteristic viscous damping time
    scalar rtd = 0.5*Cmu_*mu/(rho*r2);

    // oscillation frequency (squared)
    scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;

    if (omega2 > 0)
    {
        scalar omega = sqrt(omega2);
        scalar romega = 1.0/omega;

        scalar rhog = spray_.rho()[p.cell()];
        scalar We = p.We(Ug, rhog, sigma);
        scalar Wetmp = We/WeCrit_;

        scalar y1 = p.dev() - Wetmp;
        scalar y2 = p.ddev()*romega;

        scalar a = sqrt(y1*y1 + y2*y2);

        // scotty we may have break-up
        if (a+Wetmp > 1.0)
        {
            scalar phic = y1/a;

            // constrain phic within -1 to 1
            phic = max(min(phic, 1), -1);

            scalar phit = acos(phic);
            scalar phi = phit;
            scalar quad = -y2/a;
            if (quad < 0)
            {
                phi = 2*mathematicalConstant::pi - phit;
            }

            scalar tb = 0;

            if (mag(p.dev()) < 1.0)
            {
                scalar theta = acos((1.0 - Wetmp)/a);

                if (theta < phi)
                {
                    if (2*mathematicalConstant::pi-theta >= phi)
                    {
                        theta = -theta;
                    }
                    theta += 2*mathematicalConstant::pi;
                }
                tb = (theta-phi)*romega;

                // breakup occurs
                if (deltaT > tb)
                {
                    p.dev() = 1.0;
                    p.ddev() = -a*omega*sin(omega*tb + phi);
                }
            }

            // update droplet size
            if (deltaT > tb)
            {
                scalar sqrtWe = AWe_*pow(We, 4.0) + 1.0;
                scalar Kbr = k1_*omega*sqrtWe;

                if (We > WeTransition_)
                {
                    sqrtWe = sqrt(We);
                    Kbr =k2_*omega*sqrtWe;
                }

                scalar rWetmp = 1.0/Wetmp;
                scalar cosdtbu = max(-1.0, min(1.0, 1.0-rWetmp));
                scalar dtbu = romega*acos(cosdtbu);
                scalar decay = exp(-Kbr*dtbu);

                scalar rNew = decay*r;
                if (rNew < r)
                {
                    p.d() = 2*rNew;
//.........这里部分代码省略.........

void parcelhandler::put_parcel(parcel p, write_handler_type f)
{
    HPX_ASSERT(resolver_);

    naming::id_type const* ids = p.destinations();
    naming::address* addrs = p.addrs();

    // During bootstrap this is handled separately (see
    // addressing_service::resolve_locality.
    if (0 == hpx::threads::get_self_ptr() && !hpx::is_starting())
    {
        HPX_ASSERT(resolver_);
        naming::gid_type locality =
            naming::get_locality_from_gid(ids[0].get_gid());
        if (!resolver_->has_resolved_locality(locality))
        {
            // reschedule request as an HPX thread to avoid hangs
            void (parcelhandler::*put_parcel_ptr) (
                parcel p, write_handler_type f
            ) = &parcelhandler::put_parcel;

            threads::register_thread_nullary(
                util::bind(
                    util::one_shot(put_parcel_ptr), this,
                    std::move(p), std::move(f)),
                "parcelhandler::put_parcel", threads::pending, true,
                threads::thread_priority_boost);
            return;
        }
    }

    // properly initialize parcel
    init_parcel(p);

    bool resolved_locally = true;

#if !defined(HPX_SUPPORT_MULTIPLE_PARCEL_DESTINATIONS)
    if (!addrs[0])
    {
        resolved_locally = resolver_->resolve_local(ids[0], addrs[0]);
    }
#else
    std::size_t size = p.size();

    if (0 == size) {
        HPX_THROW_EXCEPTION(network_error, "parcelhandler::put_parcel",
                            "no destination address given");
        return;
    }

    if (1 == size) {
        if (!addrs[0])
            resolved_locally = resolver_->resolve_local(ids[0], addrs[0]);
    }
    else {
        boost::dynamic_bitset<> locals;
        resolved_locally = resolver_->resolve_local(ids, addrs, size, locals);
    }
#endif

    if (!p.parcel_id())
        p.parcel_id() = parcel::generate_unique_id();

    using util::placeholders::_1;
    using util::placeholders::_2;
    write_handler_type wrapped_f =
        util::bind(&detail::parcel_sent_handler, std::move(f), _1, _2);

    // If we were able to resolve the address(es) locally we send the
    // parcel directly to the destination.
    if (resolved_locally)
    {
        // dispatch to the message handler which is associated with the
        // encapsulated action
        typedef std::pair<boost::shared_ptr<parcelport>, locality> destination_pair;
        destination_pair dest = find_appropriate_destination(addrs[0].locality_);

        if (load_message_handlers_)
        {
            policies::message_handler* mh =
                p.get_message_handler(this, dest.second);

            if (mh) {
                mh->put_parcel(dest.second, std::move(p), std::move(wrapped_f));
                return;
            }
        }

        dest.first->put_parcel(dest.second, std::move(p), std::move(wrapped_f));
        return;
    }

    // At least one of the addresses is locally unknown, route the parcel
    // to the AGAS managing the destination.
    ++count_routed_;

    resolver_->route(std::move(p), std::move(wrapped_f));
}

void reitzDiwakar::breakupParcel
(
    parcel& p,
    const scalar deltaT,
    const vector& vel,
    const liquidMixture& fuels
) const
{
    /*
        These are the default values for this model...
        static const scalar Cbag   = 6.0;
        static const scalar Cb     = 0.785;
        static const scalar Cstrip = 0.5;
        static const scalar Cs     = 10.0;
    */

    const PtrList<volScalarField>& Y = spray_.composition().Y();

    label Ns = Y.size();
    label cellI = p.cell();
    scalar pressure = spray_.p()[cellI];
    scalar temperature = spray_.T()[cellI];
    scalar Taverage = p.T() + (temperature - p.T())/3.0;

    scalar muAverage = 0.0;
    scalar Winv = 0.0;
    for(label i=0; i<Ns; i++)
    {
        Winv += Y[i][cellI]/spray_.gasProperties()[i].W();
        muAverage += Y[i][cellI]*spray_.gasProperties()[i].mu(Taverage);
    }
    scalar R = specie::RR()*Winv;

    // ideal gas law to evaluate density
    scalar rhoAverage = pressure/R/Taverage;
    scalar nuAverage = muAverage/rhoAverage;
    scalar sigma = fuels.sigma(pressure, p.T(), p.X());

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    //     The We and Re numbers are to be evaluated using the 1/3 rule.
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    scalar WeberNumber = p.We(vel, rhoAverage, sigma);
    scalar ReynoldsNumber = p.Re(vel, nuAverage);

    scalar sqRey = sqrt(ReynoldsNumber);

    if (WeberNumber > Cbag_)
    {
        if (WeberNumber > Cstrip_*sqRey)
        {
            scalar dStrip =
                pow(2.0*Cstrip_*sigma, 2.0)/
                (
                    rhoAverage
                  * pow(mag(p.Urel(vel)), 3.0)
                  * muAverage
                );

            scalar tauStrip =
                Cs_ * p.d()
              * sqrt
                (
                    fuels.rho(pressure, p.T(), p.X())
                    / rhoAverage
                )
              / mag(p.Urel(vel));

            scalar fraction = deltaT/tauStrip;

            // new droplet diameter, implicit calculation
            p.d() = (fraction*dStrip + p.d())/(1.0 + fraction);
        }
        else
        {
            scalar dBag =
                2.0 * Cbag_ * sigma
              / (
                  rhoAverage
                * pow(mag(p.Urel(vel)), 2.0)
                );

            scalar tauBag =
                Cb_ * p.d()
                * sqrt
                  (
                      fuels.rho(pressure, p.T(), p.X())
                    * p.d()
                    / sigma
                  );

            scalar fraction = deltaT/tauBag;

            // new droplet diameter, implicit calculation
            p.d() = (fraction*dBag + p.d())/(1.0 + fraction);
        }

    }

}

void blobsSheetAtomization::atomizeParcel
(
    parcel& p,
    const scalar deltaT,
    const vector& vel,
    const liquidMixture& fuels
) const
{

    const PtrList<volScalarField>& Y = spray_.composition().Y();

    label Ns = Y.size();
    label cellI = p.cell();
    scalar pressure = spray_.p()[cellI];
    scalar temperature = spray_.T()[cellI];
    scalar Taverage = p.T() + (temperature - p.T())/3.0;

    scalar Winv = 0.0;
    for(label i=0; i<Ns; i++)
    {
        Winv += Y[i][cellI]/spray_.gasProperties()[i].W();
    }
    scalar R = specie::RR*Winv;

    // ideal gas law to evaluate density
    scalar rhoAverage = pressure/R/Taverage;
    scalar sigma = fuels.sigma(pressure, p.T(), p.X());

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    //     The We and Re numbers are to be evaluated using the 1/3 rule.
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    scalar rhoFuel = fuels.rho(1.0e+5, p.T(), p.X());

    scalar U = mag(p.Urel(vel));

    const injectorType& it =
        spray_.injectors()[label(p.injector())].properties();

    vector itPosition(vector::zero);
    label nHoles = it.nHoles();
    if (nHoles > 1)
    {
        for(label i=0; i<nHoles;i++)
        {
            itPosition += it.position(i);
        }
        itPosition /= nHoles;
    }
    else
    {
        itPosition = it.position(0);
    }
//    const vector itPosition = it.position();


    scalar lBU = B_ * sqrt
    (
        rhoFuel * sigma * p.d() * cos(angle_*mathematicalConstant::pi/360.0)
      / sqr(rhoAverage*U)
    );

    scalar pWalk = mag(p.position() - itPosition);

    if(pWalk > lBU && p.liquidCore() == 1.0)
    {
        p.liquidCore() = 0.0;
    }
}

void Foam::SHF::breakupParcel
(
    parcel& p,
    const scalar deltaT,
    const vector& vel,
    const liquidMixtureProperties& fuels
) const
{
    label cellI = p.cell();
    scalar T = p.T();
    scalar pc = spray_.p()[cellI];

    scalar sigma = fuels.sigma(pc, T, p.X());
    scalar rhoLiquid = fuels.rho(pc, T, p.X());
    scalar muLiquid = fuels.mu(pc, T, p.X());
    scalar rhoGas = spray_.rho()[cellI];

    scalar weGas = p.We(vel, rhoGas, sigma);
    scalar weLiquid = p.We(vel, rhoLiquid, sigma);

    // correct the Reynolds number. Reitz is using radius instead of diameter

    scalar reLiquid   = p.Re(rhoLiquid, vel, muLiquid);
    scalar ohnesorge  = sqrt(weLiquid)/(reLiquid + VSMALL);

    vector vRel = p.Urel(vel);

    scalar weGasCorr = weGas/(1.0 + weCorrCoeff_*ohnesorge);

    // droplet deformation characteristic time

    scalar tChar = p.d()/mag(vRel)*sqrt(rhoLiquid/rhoGas);

    scalar tFirst = cInit_*tChar;

    scalar tSecond = 0;
    scalar tCharSecond = 0;


    //  updating the droplet characteristic time
    p.ct() += deltaT;

    if (weGas > weConst_)
    {
        if (weGas < weCrit1_)
        {
            tCharSecond = c1_*pow((weGas - weConst_),cExp1_);
        }
        else if (weGas >= weCrit1_ && weGas <= weCrit2_)
        {
            tCharSecond = c2_*pow((weGas - weConst_),cExp2_);
        }
        else
        {
            tCharSecond = c3_*pow((weGas - weConst_),cExp3_);
        }
    }

    scalar weC = weBuCrit_*(1.0+ohnCoeffCrit_*pow(ohnesorge, ohnExpCrit_));
    scalar weB = weBuBag_*(1.0+ohnCoeffBag_*pow(ohnesorge, ohnExpBag_));
    scalar weMM = weBuMM_*(1.0+ohnCoeffMM_*pow(ohnesorge, ohnExpMM_));

    bool bag = (weGas > weC && weGas < weB);

    bool multimode = (weGas >= weB && weGas <= weMM);

    bool shear = (weGas > weMM);

    tSecond = tCharSecond*tChar;

    scalar tBreakUP = tFirst + tSecond;
    if (p.ct() > tBreakUP)
    {

        scalar d32 =
            coeffD_*p.d()*pow(ohnesorge, onExpD_)*pow(weGasCorr, weExpD_);

        if (bag || multimode)
        {

            scalar d05 = d32Coeff_*d32;

            scalar x = 0.0;
            scalar y = 0.0;
            scalar d = 0.0;
            scalar px = 0.0;

            do
            {
                x = cDmaxBM_*rndGen_.sample01<scalar>();
                d = sqr(x)*d05;
                y = rndGen_.sample01<scalar>();

                px =
                    x
                   /(2.0*sqrt(constant::mathematical::twoPi)*sigma_)
                   *exp(-0.5*sqr((x-mu_)/sigma_));

            } while (y >= px);

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

// Return 'keepParcel'
bool reflectParcel::wallTreatment
(
    parcel& p,
    const label globalFacei
) const
{
    label patchi = p.patch(globalFacei);
    label facei = p.patchFace(patchi, globalFacei);

    const polyMesh& mesh = spray_.mesh();

    if (mesh_.boundaryMesh()[patchi].isWall())
    {
        // wallNormal defined to point outwards of domain
        vector Sf = mesh_.Sf().boundaryField()[patchi][facei];
        Sf /= mag(Sf);

        if (!mesh.moving())
        {
            // static mesh
            scalar Un = p.U() & Sf;

            if (Un > 0)
            {
                p.U() -= (1.0 + elasticity_)*Un*Sf;
            }

        }
        else
        {
            // moving mesh
            vector Ub1 = U_.boundaryField()[patchi][facei];
            vector Ub0 = U_.oldTime().boundaryField()[patchi][facei];

            scalar dt = spray_.runTime().deltaT().value();
            const vectorField& oldPoints = mesh.oldPoints();

            const vector& Cf1 = mesh.faceCentres()[globalFacei];

            vector Cf0 = mesh.faces()[globalFacei].centre(oldPoints);
            vector Cf = Cf0 + p.stepFraction()*(Cf1 - Cf0);
            vector Sf0 = mesh.faces()[globalFacei].normal(oldPoints);

            // for layer addition Sf0 = vector::zero and we use Sf
            if (mag(Sf0) > SMALL)
            {
                Sf0 /= mag(Sf0);
            }
            else
            {
                Sf0 = Sf;
            }

            scalar magSfDiff = mag(Sf - Sf0);

            vector Ub = Ub0 + p.stepFraction()*(Ub1 - Ub0);

            if (magSfDiff > SMALL)
            {
                // rotation + translation
                vector Sfp = Sf0 + p.stepFraction()*(Sf - Sf0);

                vector omega = Sf0 ^ Sf;
                scalar magOmega = mag(omega);
                omega /= magOmega+SMALL;

                scalar phiVel = ::asin(magOmega)/dt;

                scalar dist = (p.position() - Cf) & Sfp;
                vector pos = p.position() - dist*Sfp;
                vector vrot = phiVel*(omega ^ (pos - Cf));

                vector v = Ub + vrot;

                scalar Un = ((p.U() - v) & Sfp);

                if (Un > 0.0)
                {
                    p.U() -= (1.0 + elasticity_)*Un*Sfp;
                }
            }
            else
            {
                // translation
                vector Ur = p.U() - Ub;
                scalar Urn = Ur & Sf;
                /*
                if (mag(Ub-v) > SMALL)
                {
                    Info << "reflectParcel:: v = " << v
                        << ", Ub = " << Ub
                        << ", facei = " << facei
                        << ", patchi = " << patchi
                        << ", globalFacei = " << globalFacei
                        << ", name = " << mesh_.boundaryMesh()[patchi].name()
                        << endl;
                }
                    */
                if (Urn > 0.0)
//.........这里部分代码省略.........

void myLISA_3_InjPos::atomizeParcel
(
    parcel& p,
    const scalar deltaT,
    const vector& vel,
    const liquidMixture& fuels
) const
{


    const PtrList<volScalarField>& Y = spray_.composition().Y();

    label Ns = Y.size();
    label cellI = p.cell();
    scalar pressure = spray_.p()[cellI];
    scalar temperature = spray_.T()[cellI];
//--------------------------------------AL____101015--------------------------------//
//  scalar Taverage = p.T() + (temperature - p.T())/3.0;
    scalar Taverage = temperature;
//-----------------------------------------END--------------------------------------//
    scalar Winv = 0.0;

    for(label i=0; i<Ns; i++)
    {
        Winv += Y[i][cellI]/spray_.gasProperties()[i].W();
    }

    scalar R = specie::RR*Winv;

    // ideal gas law to evaluate density
    scalar rhoAverage = pressure/R/Taverage;
    //scalar nuAverage = muAverage/rhoAverage;
    scalar sigma = fuels.sigma(pressure, p.T(), p.X());

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    //     The We and Re numbers are to be evaluated using the 1/3 rule.
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    scalar WeberNumber = p.We(vel, rhoAverage, sigma);

    scalar tau = 0.0;
    scalar dL = 0.0;
    scalar k = 0.0;
    scalar muFuel = fuels.mu(pressure, p.T(), p.X());
    scalar rhoFuel = fuels.rho(1.0e+5, p.T(), p.X());
    scalar nuFuel = muFuel/rhoFuel;

    vector uDir = p.U()/mag(p.U());

    scalar uGas = mag(vel & uDir);
    vector Ug = uGas*uDir;

    /*
        TL
        It might be the relative velocity between Liquid and Gas, but I use the
        absolute velocity of the parcel as suggested by the authors
    */

//    scalar U = mag(p.Urel(vel));
    scalar U = mag(p.U());
    p.ct() += deltaT;

    scalar Q = rhoAverage/rhoFuel;

    const injectorType& it =
        spray_.injectors()[label(p.injector())].properties();

    if (it.nHoles() > 1)
    {
        Info << "Warning: This atomization model is not suitable for multihole injector." << endl
             << "Only the first hole will be used." << endl;
    }

    const vector direction = it.direction(0, spray_.runTime().value());
//--------------------------------CH 101108--------------------------------------------------//
//    const vector itPosition = it.position(0);
    const injectorModel& im = spray_.injection();
    const vector itPosition = it.position(0) + im.injDist(0)*direction/mag(direction);
//------------------------------------END----------------------------------------------------//
    scalar pWalk = mag(p.position() - itPosition);

//  Updating liquid sheet tickness... that is the droplet diameter

//    const vector direction = it.direction(0, spray_.runTime().value());

    scalar h = (p.position() - itPosition) & direction;

    scalar d = sqrt(sqr(pWalk)-sqr(h));

    scalar time = pWalk/mag(p.U());

    scalar elapsedTime = spray_.runTime().value();

    scalar massFlow = it.massFlowRate(max(0.0,elapsedTime-time));

    scalar hSheet = massFlow/(mathematicalConstant::pi*d*rhoFuel*mag(p.U()));

    p.d() = min(hSheet,p.d());

    if(WeberNumber > 27.0/16.0)
//.........这里部分代码省略.........