自学教程：C++ u函数代码示例

void dgApi dgRayToRayDistance (const dgVector& ray_p0, const dgVector& ray_p1, const dgVector& ray_q0, const dgVector& ray_q1, dgVector& pOut, dgVector& qOut){	dgFloat32 sN;	dgFloat32 tN;	dgVector u (ray_p1 - ray_p0);	dgVector v (ray_q1 - ray_q0);	dgVector w (ray_p0 - ray_q0);	dgFloat32 a = u.DotProduct3(u);        	dgFloat32 b = u.DotProduct3(v);	dgFloat32 c = v.DotProduct3(v);        	dgFloat32 d = u.DotProduct3(w);	dgFloat32 e = v.DotProduct3(w);	dgFloat32 D = a*c - b*b;   	dgFloat32 sD = D;				dgFloat32 tD = D;				// compute the line parameters of the two closest points	if (D < dgFloat32 (1.0e-8f)) { 		sN = dgFloat32 (0.0f);        		sD = dgFloat32 (1.0f);        		tN = e;		tD = c;	} else {                		// get the closest points on the infinite lines		sN = (b*e - c*d);		tN = (a*e - b*d);		if (sN < dgFloat32 (0.0f)) {       			// sc < 0 => the s=0 edge is visible			sN = dgFloat32 (0.0f);			tN = e;			tD = c;		} else if (sN > sD) {  			// sc > 1 => the s=1 edge is visible			sN = sD;			tN = e + b;			tD = c;		}	}	if (tN < dgFloat32 (0.0f)) {           		// tc < 0 => the t=0 edge is visible		tN = dgFloat32 (0.0f);		// recompute sc for this edge		if (-d < dgFloat32 (0.0f)) {			sN = dgFloat32 (0.0f);		} else if (-d > a) {			sN = sD;		} else {			sN = -d;			sD = a;		}	} else if (tN > tD) {      		// tc > 1 => the t=1 edge is visible		tN = tD;		// recompute sc for this edge		if ((-d + b) < dgFloat32 (0.0f)) {			sN = dgFloat32 (0.0f);		} else if ((-d + b) > a) {			sN = sD;		} else {			sN = (-d + b);			sD = a;		}	}	// finally do the division to get sc and tc	dgFloat32 sc = (dgAbs(sN) < dgFloat32(1.0e-8f) ? dgFloat32 (0.0f) : sN / sD);	dgFloat32 tc = (dgAbs(tN) < dgFloat32(1.0e-8f) ? dgFloat32 (0.0f) : tN / tD);	dgAssert (u.m_w == dgFloat32 (0.0f));	dgAssert (v.m_w == dgFloat32 (0.0f));	pOut = ray_p0 + u.Scale (sc);	qOut = ray_q0 + v.Scale (tc);}

//get new Ori and Pos arrays for the child after recombinationvoid AnimalClass::sampleMyPosOri(AnimalClass& father,                               AnimalClass& mother){    sireId  = father.myId;    damId   = mother.myId;		chromosome *currentFatherChrom;	chromosome *currentMotherChrom;		unsigned numChromosomePair = G.get_num_chrom();        //construct paternal chromosome for myId    ArrayXf tempPos;    ArrayXi tempOri;    tempPos.resize(100000);    tempOri.resize(100000);        for(unsigned i=0;i<numChromosomePair;i++){        double chrLength = G[i].chr_length;        unsigned binomialN = chrLength*3 + 1;        vector<float> rPos;        binomial_distribution<int> Binom(binomialN,chrLength/binomialN);        int numCrossover    =   Binom(randGen);                uniform_real_distribution<float> u(0,1);        for (unsigned k=0; k<numCrossover; k++) {            rPos.push_back(chrLength*u(randGen));        }        rPos.push_back(chrLength);        sort(rPos.begin(),rPos.end());                unsigned   startPosParent=0;        unsigned   startPosMe=0;                currentFatherChrom = (u(randGen)<0.5)?&father.GenomePat[i]:&father.GenomeMat[i];        for(unsigned j=0;j<rPos.size();j++){            unsigned numCopy=(currentFatherChrom->pos < rPos[j]).count()-startPosParent;            tempPos.block(startPosMe,0,numCopy,1)=currentFatherChrom->pos.block(startPosParent,0,numCopy,1);            tempOri.block(startPosMe,0,numCopy,1)=currentFatherChrom->ori.block(startPosParent,0,numCopy,1);                        currentFatherChrom= (currentFatherChrom==&father.GenomePat[i])?&father.GenomeMat[i]:&father.GenomePat[i];                        startPosParent=(currentFatherChrom->pos < rPos[j]).count();            startPosMe+=numCopy;            tempPos(startPosMe)=rPos[j];            tempOri(startPosMe)=currentFatherChrom->ori(startPosParent-1);            startPosMe++;        } //the length of tempPos should be startPosMe now. Tricky here.        //        unsigned posLengthMinusOne=startPosMe-1;//        unsigned nonzero=(tempPos.block(1,0,posLengthMinusOne,1)-tempPos.block(0,0,posLengthMinusOne,1)).count()+1;//        GenomePat[i].pos.resize(nonzero-1);//        GenomePat[i].ori.resize(nonzero-1);               unsigned keep=0;        GenomePat[i].pos.resize(startPosMe);        GenomePat[i].ori.resize(startPosMe);        GenomePat[i].pos[0]=tempPos[0];        GenomePat[i].ori[0]=tempOri[0];        for(unsigned m=1;m < startPosMe-1; m++){            if(tempOri[m]!=tempOri[m-1]){            keep=keep+1;            GenomePat[i].pos[keep]=tempPos[m];            GenomePat[i].ori[keep]=tempOri[m];            }        }                unsigned PosOriSize=keep+1;        GenomePat[i].pos.conservativeResize(PosOriSize);        GenomePat[i].ori.conservativeResize(PosOriSize);    }    ///    ///construct maternal chromosome for myId    for(unsigned i=0;i<numChromosomePair;i++){        double chrLength = G[i].chr_length;        unsigned binomialN = chrLength*3 + 1;        vector<float> rPos;        binomial_distribution<int> Binom(binomialN,chrLength/binomialN);        int numCrossover    =   Binom(randGen);                uniform_real_distribution<float> u(0,1);        for (unsigned k=0; k<numCrossover; k++) {            float rPos_Hao=chrLength*u(randGen);                rPos.push_back(rPos_Hao);        }        rPos.push_back(chrLength);        sort(rPos.begin(),rPos.end());                unsigned   startPosParent=0;        unsigned   startPosMe=0;           currentMotherChrom = (u(randGen)<0.5)?&mother.GenomePat[i]:&mother.GenomeMat[i];//.........这里部分代码省略.........

	void AddressBookSubscription::Request ()	{		// must be run in separate thread			LogPrint (eLogInfo, "Downloading hosts from ", m_Link, " ETag: ", m_Etag, " Last-Modified: ", m_LastModified);		bool success = false;			i2p::util::http::url u (m_Link);		i2p::data::IdentHash ident;		if (m_Book.GetIdentHash (u.host_, ident))		{			std::condition_variable newDataReceived;			std::mutex newDataReceivedMutex;			auto leaseSet = i2p::client::context.GetSharedLocalDestination ()->FindLeaseSet (ident);			if (!leaseSet)			{				std::unique_lock<std::mutex> l(newDataReceivedMutex);				i2p::client::context.GetSharedLocalDestination ()->RequestDestination (ident,					[&newDataReceived, &leaseSet](std::shared_ptr<i2p::data::LeaseSet> ls)				    {						leaseSet = ls;						newDataReceived.notify_all ();					});				if (newDataReceived.wait_for (l, std::chrono::seconds (SUBSCRIPTION_REQUEST_TIMEOUT)) == std::cv_status::timeout)					LogPrint (eLogError, "Subscription LeseseSet request timeout expired");			}			if (leaseSet)			{				std::stringstream request, response;				// standard header				request << "GET " << u.path_ << " HTTP/1.1/r/nHost: " << u.host_				<< "/r/nAccept: */*/r/n" << "User-Agent: Wget/1.11.4/r/n" << "Connection: close/r/n";				if (m_Etag.length () > 0) // etag					request << i2p::util::http::IF_NONE_MATCH << ": /"" << m_Etag << "/"/r/n";				if (m_LastModified.length () > 0) // if-modfief-since					request << i2p::util::http::IF_MODIFIED_SINCE << ": " << m_LastModified << "/r/n";				request << "/r/n"; // end of header				auto stream = i2p::client::context.GetSharedLocalDestination ()->CreateStream (leaseSet, u.port_);				stream->Send ((uint8_t *)request.str ().c_str (), request.str ().length ());								uint8_t buf[4096];				bool end = false;				while (!end)				{					stream->AsyncReceive (boost::asio::buffer (buf, 4096), 						[&](const boost::system::error_code& ecode, std::size_t bytes_transferred)						{							if (bytes_transferred)								response.write ((char *)buf, bytes_transferred);							if (ecode == boost::asio::error::timed_out || !stream->IsOpen ())								end = true;								newDataReceived.notify_all ();						},						30); // wait for 30 seconds					std::unique_lock<std::mutex> l(newDataReceivedMutex);					if (newDataReceived.wait_for (l, std::chrono::seconds (SUBSCRIPTION_REQUEST_TIMEOUT)) == std::cv_status::timeout)						LogPrint (eLogError, "Subscription timeout expired");				}				// process remaining buffer				while (size_t len = stream->ReadSome (buf, 4096))					response.write ((char *)buf, len);								// parse response				std::string version;				response >> version; // HTTP version				int status = 0;				response >> status; // status				if (status == 200) // OK				{					bool isChunked = false;					std::string header, statusMessage;					std::getline (response, statusMessage);					// read until new line meaning end of header					while (!response.eof () && header != "/r")					{						std::getline (response, header);						auto colon = header.find (':');						if (colon != std::string::npos)						{							std::string field = header.substr (0, colon);							header.resize (header.length () - 1); // delete /r								if (field == i2p::util::http::ETAG)								m_Etag = header.substr (colon + 1);							else if (field == i2p::util::http::LAST_MODIFIED)								m_LastModified = header.substr (colon + 1);							else if (field == i2p::util::http::TRANSFER_ENCODING)								isChunked = !header.compare (colon + 1, std::string::npos, "chunked");						}						}					LogPrint (eLogInfo, m_Link, " ETag: ", m_Etag, " Last-Modified: ", m_LastModified);					if (!response.eof ())						{						success = true;						if (!isChunked)							m_Book.LoadHostsFromStream (response);						else						{							// merge chunks							std::stringstream merged;							i2p::util::http::MergeChunkedResponse (response, merged);							m_Book.LoadHostsFromStream (merged);						}	//.........这里部分代码省略.........

//std::vector <osg::Vec3> LeafGrower::grow_planes(osg::Vec3 origin, double height, int angle, bool maya)std::vector <osg::Vec3> LeafGrower::grow_planes(){    std::vector <osg::Vec3> ret;    if(!_root)        return ret;    //copy to avoid modification    //BDLSkeletonNode *root = BDLSkeletonNode::copy_tree(_root);    //BDLSkeletonNode::rectify_skeleton_tree(root, origin, height, angle, maya);	BDLSkeletonNode *root = _root;    //set the expected leaf size be 0.35 of the trunk's radius    float leaf_size_hint = root->_radius * 0.35;    //store for billboard data    _all_pos.clear();    //bfs    std::queue <BDLSkeletonNode *> Queue;    Queue.push(root);    int ball_cnt = 0;    float outgrow = 1.3f;    //grow leaves only at tips or small radius nodes    while(!Queue.empty())    {        BDLSkeletonNode *front = Queue.front();        Queue.pop();        for(unsigned int i=0; i<front->_children.size(); i++)            Queue.push(front->_children[i]);        //skip the root        if(!front->_prev)            continue;        //graph's leafs only        if(front->_children.size() == 0)        {            osg::Vec3 pos(front->_sx, front->_sy, front->_sz);            osg::Vec3 pre(front->_prev->_sx, front->_prev->_sy, front->_prev->_sz);             osg::Vec3 mid = (pos+pre) / 2.0f;            osg::Vec3 tangent = pos - pre;            float cloud_width = tangent.length() * 0.50 * outgrow;            tangent.normalize();            //consider a plane with tangent as the normal and contains point pos            //find basis u, v of this plane            osg::Vec3 u(tangent.y(), -tangent.x(), 0.0);            osg::Vec3 v = tangent ^ u;            //consider another plane with u as normal and contains point pos            //the basis of this plane is chosen as tangent and v            //div_angle is the angle angle spanned by tangent and the root of the leaf            double div_angle = (rand()%30-15+70)/180.0*M_PI;            //number of ball of this node            int no_leaf = 1;            for(int i=1; i<=no_leaf; i++)            {                double angle = 360.0/no_leaf*i/180.0*M_PI;                osg::Vec3 div = u*cos(angle) + v*sin(angle);                div = tangent * cos(div_angle) + div * sin(div_angle);                //frame: fu, fv, div, at pos                osg::Vec3 fu(div.y(), -div.x(), 0.0f);                fu.normalize();                osg::Vec3 fv = div ^ fu;                //cloud around this node                std::vector <osg::Vec3> b_cloud;                std::vector <float> weights;                Transformer::sphere_points_transformed(b_cloud, weights, mid, cloud_width, tangent, 10-rand()%20, leaf_size_hint);                //Transformer::sphere_points_transformed(b_cloud, weights, mid, cloud_width, tangent, -89, leaf_size_hint);                //from points to leaf vertices                for(unsigned j=0; j<b_cloud.size(); j++)                {                    osg::Vec3 b_p = b_cloud[j];                    osg::Vec3 normal = b_p - mid;                    normal.normalize();                    osg::Vec3 normal2 = b_p - pos;                    normal2.normalize();                    //displace                    float d_ratio = weights[j] > 0.5 ? 0.40 : 1.0f/outgrow;                    b_p += normal2 * cloud_width * d_ratio * (weights[j]-0.5f);                    float leaf_s = leaf_size_hint * (1.0-0.8*std::max(0.0f, weights[j]-0.5f));                    //cloud pt inside segmentation?                    if(false || _pruner.is_inside_ortho(b_p))                    {                        Transformer::points_to_leafs(ret, b_p, normal, leaf_s, true);                        //for billboard data//.........这里部分代码省略.........

void FluidSim::solve_viscosity(float dt) {   int ni = liquid_phi.ni;   int nj = liquid_phi.nj;      //static obstacles for simplicity - for moving objects,    //use a spatially varying 2d array, and modify the linear system appropriately   float u_obj = 0;   float v_obj = 0;   Array2c u_state(ni+1,nj,(const char&)0);   Array2c v_state(ni,nj+1,(const char&)0);   const int SOLID = 1;   const int FLUID = 0;   printf("Determining states/n");   //just determine if the face position is inside the wall! That's it.   for(int j = 0; j < nj; ++j) {      for(int i = 0; i < ni+1; ++i) {         if(i - 1 < 0 || i >= ni || (nodal_solid_phi(i,j+1) + nodal_solid_phi(i,j))/2 <= 0)            u_state(i,j) = SOLID;         else            u_state(i,j) = FLUID;      }   }   for(int j = 0; j < nj+1; ++j)  {      for(int i = 0; i < ni; ++i) {         if(j - 1 < 0 || j >= nj || (nodal_solid_phi(i+1,j) + nodal_solid_phi(i,j))/2 <= 0)            v_state(i,j) = SOLID;         else             v_state(i,j) = FLUID;      }   }      printf("Building matrix/n");   int elts = (ni+1)*nj + ni*(nj+1);   if(vrhs.size() != elts) {      vrhs.resize(elts);      velocities.resize(elts);      vmatrix.resize(elts);   }   vmatrix.zero();      float factor = dt/sqr(dx);   for(int j = 1; j < nj-1; ++j) for(int i = 1; i < ni-1; ++i) {      if(u_state(i,j) == FLUID ) {         int index = u_ind(i,j);                        vrhs[index] = u_vol(i,j) * u(i,j);         vmatrix.set_element(index,index,u_vol(i,j));                  //uxx terms         float visc_right = viscosity(i,j);         float visc_left = viscosity(i-1,j);         float vol_right = c_vol(i,j);         float vol_left = c_vol(i-1,j);        //u_x_right         vmatrix.add_to_element(index,index, 2*factor*visc_right*vol_right);         if(u_state(i+1,j) == FLUID)            vmatrix.add_to_element(index,u_ind(i+1,j), -2*factor*visc_right*vol_right);         else if(u_state(i+1,j) == SOLID)            vrhs[index] -= -2*factor*visc_right*vol_right*u_obj;         //u_x_left         vmatrix.add_to_element(index,index, 2*factor*visc_left*vol_left);         if(u_state(i-1,j) == FLUID)            vmatrix.add_to_element(index,u_ind(i-1,j), -2*factor*visc_left*vol_left);         else if(u_state(i-1,j) == SOLID)            vrhs[index] -= -2*factor*visc_left*vol_left*u_obj;                  //uyy terms         float visc_top = 0.25f*(viscosity(i-1,j+1) + viscosity(i-1,j) + viscosity(i,j+1) + viscosity(i,j));         float visc_bottom = 0.25f*(viscosity(i-1,j) + viscosity(i-1,j-1) + viscosity(i,j) + viscosity(i,j-1));         float vol_top = n_vol(i,j+1);         float vol_bottom = n_vol(i,j);         //u_y_top         vmatrix.add_to_element(index,index, +factor*visc_top*vol_top);         if(u_state(i,j+1) == FLUID)            vmatrix.add_to_element(index,u_ind(i,j+1), -factor*visc_top*vol_top);         else if(u_state(i,j+1) == SOLID)            vrhs[index] -= -u_obj*factor*visc_top*vol_top;               //u_y_bottom         vmatrix.add_to_element(index,index, +factor*visc_bottom*vol_bottom);         if(u_state(i,j-1) == FLUID)            vmatrix.add_to_element(index,u_ind(i,j-1), -factor*visc_bottom*vol_bottom);         else if(u_state(i,j-1) == SOLID)            vrhs[index] -= -u_obj*factor*visc_bottom*vol_bottom;               //vxy terms         //v_x_top         if(v_state(i,j+1) == FLUID)            vmatrix.add_to_element(index,v_ind(i,j+1), -factor*visc_top*vol_top);         else if(v_state(i,j+1) == SOLID)            vrhs[index] -= -v_obj*factor*visc_top*vol_top;                  if(v_state(i-1,j+1) == FLUID)//.........这里部分代码省略.........

int main(int argc, char* argv[]) {  Teuchos::GlobalMPISession mpiSession(&argc, &argv);  Teuchos::RCP<const Teuchos::Comm<int> > comm    = Teuchos::DefaultComm<int>::getComm();  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.  int iprint = argc - 1;  Teuchos::RCP<std::ostream> outStream;  Teuchos::oblackholestream bhs; // outputs nothing  if (iprint > 0 && Teuchos::rank<int>(*comm)==0)    outStream = Teuchos::rcp(&std::cout, false);  else    outStream = Teuchos::rcp(&bhs, false);  int errorFlag  = 0;  try {    /**********************************************************************************************/    /************************* CONSTRUCT ROL ALGORITHM ********************************************/    /**********************************************************************************************/    // Get ROL parameterlist    std::string filename = "input.xml";    Teuchos::RCP<Teuchos::ParameterList> parlist = Teuchos::rcp( new Teuchos::ParameterList() );    Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );    // Build ROL algorithm    parlist->sublist("Status Test").set("Gradient Tolerance",1.e-8);    parlist->sublist("Status Test").set("Step Tolerance",1.e-14);    parlist->sublist("Status Test").set("Iteration Limit",100);    /**********************************************************************************************/    /************************* CONSTRUCT VECTORS **************************************************/    /**********************************************************************************************/    // Build control vectors    int nx = 256;    Teuchos::RCP<std::vector<RealT> > x_rcp  = Teuchos::rcp( new std::vector<RealT>(nx+2,0.0) );    ROL::StdVector<RealT> x(x_rcp);    Teuchos::RCP<std::vector<RealT> > d_rcp  = Teuchos::rcp( new std::vector<RealT>(nx+2,0.0) );    ROL::StdVector<RealT> d(d_rcp);    for ( int i = 0; i < nx+2; i++ ) {      (*x_rcp)[i] = random<RealT>(comm);      (*d_rcp)[i] = random<RealT>(comm);    }    Teuchos::RCP<ROL::Vector<RealT> > xp = Teuchos::rcp(&x,false);    // Build state and adjoint vectors    Teuchos::RCP<std::vector<RealT> > u_rcp  = Teuchos::rcp( new std::vector<RealT>(nx,1.0) );    ROL::StdVector<RealT> u(u_rcp);    Teuchos::RCP<std::vector<RealT> > p_rcp  = Teuchos::rcp( new std::vector<RealT>(nx,0.0) );    ROL::StdVector<RealT> p(p_rcp);    Teuchos::RCP<ROL::Vector<RealT> > up = Teuchos::rcp(&u,false);    Teuchos::RCP<ROL::Vector<RealT> > pp = Teuchos::rcp(&p,false);    /**********************************************************************************************/    /************************* CONSTRUCT SOL COMPONENTS *******************************************/    /**********************************************************************************************/    // Build samplers//    int dim = 4;//    int nSamp = 100;//    std::vector<RealT> tmp(2,0.0); tmp[0] = -1.0; tmp[1] = 1.0;//    std::vector<std::vector<RealT> > bounds(dim,tmp);//    Teuchos::RCP<ROL::BatchManager<RealT> > bman//      = Teuchos::rcp(new ROL::StdTeuchosBatchManager<RealT,int>(comm));//    Teuchos::RCP<ROL::SampleGenerator<RealT> > sampler//      = Teuchos::rcp(new ROL::MonteCarloGenerator<RealT>(nSamp,bounds,bman,false,false,100));    ROL::QuadratureInfo info;    info.dim = 4;    info.maxLevel = 7;    info.rule1D.clear(); info.rule1D.resize(info.dim,ROL::QUAD_CLENSHAWCURTIS);        info.growth1D.clear(); info.growth1D.resize(info.dim,ROL::GROWTH_DEFAULT);    info.normalized = true;    info.adaptive = false;    Teuchos::RCP<ROL::BatchManager<RealT> > bman      = Teuchos::rcp(new ROL::StdTeuchosBatchManager<RealT,int>(comm));    Teuchos::RCP<ROL::SampleGenerator<RealT> > sampler      = Teuchos::rcp(new ROL::SparseGridGenerator<RealT>(bman,info));    // Print samples    int width = 21;    std::stringstream name;    name << "samples_" << sampler->batchID() << ".txt";    std::ofstream file(name.str().c_str());    if (file.is_open()) {      file << std::scientific << std::setprecision(12);      for (int i = 0; i < sampler->numMySamples(); ++i) {        std::vector<RealT> pt = sampler->getMyPoint(i);        RealT wt = sampler->getMyWeight(i);        for (int j = 0; j < static_cast<int>(pt.size()); ++j) {          file << std::setw(width) << std::left << pt[j];        }        file << std::setw(width) << std::left << wt << std::endl;      }      file.close();    }    else {      TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,        ">>> (adapters/trikota/sol/test/test_01): Unable to open file!");    }    /**********************************************************************************************/    /************************* CONSTRUCT OBJECTIVE FUNCTION ***************************************/    /**********************************************************************************************/    // Build risk-averse objective function    RealT alpha = 1.e-3;//.........这里部分代码省略.........

bool QUrl::operator==( const QString& url ) const{    QUrl u( url );    return ( *this == u );}

    void tsp_vrplc::compute_constraints_impl(constraint_vector &c, const decision_vector& x) const     {        decision_vector::size_type n_cities = get_n_cities();        switch( get_encoding() )         {            case FULL:            {                // 1 - We set the equality constraints                for (size_t i = 0; i < n_cities; i++) {                    c[i] = 0;                    c[i+n_cities] = 0;                    for (size_t j = 0; j < n_cities; j++)                     {                        if(i==j) continue; // ignoring main diagonal                        decision_vector::size_type rows = compute_idx(i, j, n_cities);                        decision_vector::size_type cols = compute_idx(j, i, n_cities);                        c[i] += x[rows];                        c[i+n_cities] += x[cols];                    }                    c[i] = c[i]-1;                    c[i+n_cities] = c[i+n_cities]-1;                }                //2 - We set the inequality constraints                //2.1 - First we compute the uj (see http://en.wikipedia.org/wiki/Travelling_salesman_problem#Integer_linear_programming_formulation)                //      we start always out tour from the first city, without loosing generality                size_t next_city = 0,current_city = 0;                std::vector<int> u(n_cities);                for (size_t i = 0; i < n_cities; i++)                 {                    u[current_city] = i+1;                    for (size_t j = 0; j < n_cities; j++)                     {                        if (current_city==j) continue;                        if (x[compute_idx(current_city, j, n_cities)] == 1)                         {                            next_city = j;                            break;                        }                    }                    current_city = next_city;                }                int count=0;                for (size_t i = 1; i < n_cities; i++) {                    for (size_t j = 1; j < n_cities; j++)                     {                        if (i==j) continue;                        c[2*n_cities+count] = u[i]-u[j] + (n_cities+1) * x[compute_idx(i, j, n_cities)] - n_cities;                        count++;                    }                }                break;            }            case RANDOMKEYS:                break;            case CITIES:            {                std::vector<population::size_type> range(n_cities);                for (std::vector<population::size_type>::size_type i=0; i<range.size(); ++i)                 {                    range[i]=i;                }                c[0] = !std::is_permutation(x.begin(),x.end(),range.begin());                break;            }        }        return;    }

double m1(double t) { return u(0.0, t); }

double m2(double t) { return u(1.0, t); }

double fi(double x) { return u(x, 0.0) ; }

double U(double x) { return u(x, 1.0); }

static void tst9() {    params_ref      ps;    reslimit        rlim;    nlsat::solver s(rlim, ps);    anum_manager & am     = s.am();    nlsat::pmanager & pm  = s.pm();    nlsat::assignment as(am);    nlsat::explain& ex    = s.get_explain();    int num_lo = 4;    int num_hi = 5;    svector<nlsat::var> los, his;    for (int i = 0; i < num_lo; ++i) {        los.push_back(s.mk_var(false));        scoped_anum num(am);        am.set(num, - i - 1);        as.set(i, num);    }    for (int i = 0; i < num_hi; ++i) {        his.push_back(s.mk_var(false));        scoped_anum num(am);        am.set(num, i + 1);        as.set(num_lo + i, num);    }    nlsat::var _z = s.mk_var(false);    nlsat::var _x = s.mk_var(false);    polynomial_ref x(pm), z(pm);    x = pm.mk_polynomial(_x);    scoped_anum val(am);    am.set(val, 0);    as.set(num_lo + num_hi, val);    as.set(num_lo + num_hi + 1, val);    s.set_rvalues(as);    nlsat::scoped_literal_vector lits(s);    for (int i = 0; i < num_lo; ++i) {        polynomial_ref y(pm);        y = pm.mk_polynomial(los[i]);        lits.push_back(mk_gt(s, x - y));    }    for (int i = 0; i < num_hi; ++i) {        polynomial_ref y(pm);        y = pm.mk_polynomial(his[i]);        lits.push_back(mk_gt(s, y - x));    }    z = pm.mk_polynomial(_z);    lits.push_back(mk_eq(s, x - z));#define TEST_ON_OFF()                                   /    std::cout << "Off ";                                /    ex.set_signed_project(false);                       /    project(s, ex, _x, lits.size()-1, lits.c_ptr());    /    std::cout << "On ";                                 /    ex.set_signed_project(true);                        /    project(s, ex, _x, lits.size()-1, lits.c_ptr());    /    std::cout << "Off ";                                /    ex.set_signed_project(false);                       /    project(s, ex, _x, lits.size(), lits.c_ptr());      /    std::cout << "On ";                                 /    ex.set_signed_project(true);                        /    project(s, ex, _x, lits.size(), lits.c_ptr())       /    TEST_ON_OFF();    lits.reset();    polynomial_ref u(pm);    u = pm.mk_polynomial(his[1]);    for (int i = 0; i < num_lo; ++i) {        polynomial_ref y(pm);        y = pm.mk_polynomial(los[i]);        lits.push_back(mk_gt(s, u*x - y));    }    for (int i = 0; i < num_hi; ++i) {        polynomial_ref y(pm);        y = pm.mk_polynomial(his[i]);        lits.push_back(mk_gt(s, y - u*x));    }    z = pm.mk_polynomial(_z);    lits.push_back(mk_eq(s, u*x - z));    TEST_ON_OFF();    lits.reset();    u = pm.mk_polynomial(los[1]);    for (int i = 0; i < num_lo; ++i) {        polynomial_ref y(pm);        y = pm.mk_polynomial(los[i]);        lits.push_back(mk_gt(s, u*x - y));    }    for (int i = 0; i < num_hi; ++i) {        polynomial_ref y(pm);        y = pm.mk_polynomial(his[i]);        lits.push_back(mk_gt(s, y - u*x));    }    z = pm.mk_polynomial(_z);    lits.push_back(mk_eq(s, x - z));    TEST_ON_OFF();}

/*!  Compute and return the interaction matrix /f$ L /f$ from  a subset (/f$ /theta u_x, /theta u_y, /theta u_z/f$) of the possible  /f$ /theta u /f$ features that represent the 3D rotation  /f$^{c^*}R_c/f$ or /f$^{c}R_{c^*}/f$, with  /f[ L = [ 0_3 /; L_{/theta u}] /f]   See the vpFeatureThetaU class description for the equations of  /f$L_{/theta u}/f$.  /param select : Selection of a subset of the possible /f$ /theta u /f$  features.   - To compute the interaction matrix for all the three /f$ /theta u /f$     features use vpBasicFeature::FEATURE_ALL. In that case the dimension of the    interaction matrix is /f$ [3 /times 6] /f$  - To compute the interaction matrix for only one of the /f$ /theta u /f$     component feature (/f$/theta u_x, /theta u_y, /theta u_z/f$) use one of the    corresponding function selectTUx(), selectTUy() or selectTUz(). In    that case the returned interaction matrix is /f$ [1 /times 6] /f$    dimension.  /return The interaction matrix computed from the /f$ /theta u /f$  features that represent either the rotation /f$^{c^*}R_c/f$ or the  rotation /f$^{c}R_{c^*}/f$.  The code below shows how to compute the interaction matrix  associated to the visual feature /f$s = /theta u_x /f$.   /code  vpRotationMatrix cdMc;  // Creation of the current feature s  vpFeatureThetaU s(vpFeatureThetaU::cdRc);  s.buildFrom(cdMc);  vpMatrix L_x = s.interaction( vpFeatureThetaU::selectTUx() );  /endcode  The code below shows how to compute the interaction matrix  associated to the /f$s = (/theta u_x, /theta u_y) /f$  subset visual feature:  /code  vpMatrix L_xy = s.interaction( vpFeatureThetaU::selectTUx() | vpFeatureThetaU::selectTUy() );  /endcode  L_xy is here now a 2 by 6 matrix. The first line corresponds to  the /f$ /theta u_x /f$ visual feature while the second one to the /f$  /theta u_y /f$ visual feature.  It is also possible to build the interaction matrix from all the /f$  /theta u /f$ components by:  /code  vpMatrix L_xyz = s.interaction( vpBasicFeature::FEATURE_ALL );  /endcode  In that case, L_xyz is a 3 by 6 interaction matrix where the last  line corresponds to the /f$ /theta u_z /f$ visual feature.*/vpMatrixvpFeatureThetaU::interaction(const unsigned int select){  vpMatrix L ;  L.resize(0,6) ;  if (deallocate == vpBasicFeature::user)  {    for (unsigned int i = 0; i < nbParameters; i++)    {      if (flags[i] == false)      {        switch(i){        case 0:          vpTRACE("Warning !!!  The interaction matrix is computed but Tu_x was not set yet");        break;        case 1:          vpTRACE("Warning !!!  The interaction matrix is computed but Tu_y was not set yet");        break;        case 2:          vpTRACE("Warning !!!  The interaction matrix is computed but Tu_z was not set yet");        break;        default:          vpTRACE("Problem during the reading of the variable flags");        }      }    }    resetFlags();  }  // Lw computed using Lw = [theta/2 u]_x +/- (I + alpha [u]_x [u]_x)  vpColVector u(3)  ;  for (unsigned int i=0 ; i < 3 ; i++) {    u[i] = s[i]/2.0 ;   }  //.........这里部分代码省略.........

voidUrlWrapper::ArgvReceived(int32 argc, char** argv){    if (argc <= 1)        return;    const char* failc = " || read -p 'Press any key'";    const char* pausec = " ; read -p 'Press any key'";    char* args[] = { (char *)"/bin/sh", (char *)"-c", NULL, NULL};    status_t err;    BUrl url(argv[1]);    BString full = BUrl(url).SetProtocol(BString()).UrlString();    BString proto = url.Protocol();    BString host = url.Host();    BString port = BString() << url.Port();    BString user = url.UserInfo();    BString pass = url.Password();    BString path = url.Path();    if (!url.IsValid()) {        fprintf(stderr, "malformed url: '%s'/n", url.UrlString().String());        return;    }    // XXX: debug    PRINT(("PROTO='%s'/n", proto.String()));    PRINT(("HOST='%s'/n", host.String()));    PRINT(("PORT='%s'/n", port.String()));    PRINT(("USER='%s'/n", user.String()));    PRINT(("PASS='%s'/n", pass.String()));    PRINT(("PATH='%s'/n", path.String()));    if (proto == "about") {        app_info info;        BString sig;        // BUrl could get an accessor for the full - proto part...        sig = host << "/" << path;        BMessage msg(B_ABOUT_REQUESTED);        if (be_roster->GetAppInfo(sig.String(), &info) == B_OK) {            BMessenger msgr(sig.String());            msgr.SendMessage(&msg);            return;        }        if (be_roster->Launch(sig.String(), &msg) == B_OK)            return;        be_roster->Launch("application/x-vnd.Haiku-About");        return;    }    if (proto == "telnet") {        BString cmd("telnet ");        if (url.HasUserInfo())            cmd << "-l " << user << " ";        cmd << host;        if (url.HasPort())            cmd << " " << port;        PRINT(("CMD='%s'/n", cmd.String()));        cmd << failc;        args[2] = (char*)cmd.String();        be_roster->Launch(kTerminalSig, 3, args);        return;    }    // see draft:    // http://tools.ietf.org/wg/secsh/draft-ietf-secsh-scp-sftp-ssh-uri/    if (proto == "ssh") {        BString cmd("ssh ");        if (url.HasUserInfo())            cmd << "-l " << user << " ";        if (url.HasPort())            cmd << "-oPort=" << port << " ";        cmd << host;        PRINT(("CMD='%s'/n", cmd.String()));        cmd << failc;        args[2] = (char*)cmd.String();        be_roster->Launch(kTerminalSig, 3, args);        // TODO: handle errors        return;    }    if (proto == "ftp") {        BString cmd("ftp ");        cmd << proto << "://";        /*        if (user.Length())        	cmd << "-l " << user << " ";        cmd << host;        */        cmd << full;        PRINT(("CMD='%s'/n", cmd.String()));        cmd << failc;        args[2] = (char*)cmd.String();        be_roster->Launch(kTerminalSig, 3, args);        // TODO: handle errors        return;    }//.........这里部分代码省略.........

void ribi::tictactoe::Board::Test() noexcept{  {    static bool is_tested{false};    if (is_tested) return;    is_tested = true;  }  const TestTimer test_timer(__func__,__FILE__,1.0);  {    //Check empty board state    {      Board t;      const int s = t.GetSummarizedState();      Board u(s);      assert(u == t);    }    //Check one-move states    for (int i=0; i!=9; ++i)    {      Board t;      t.DoMove(i/3,i%3,Player::player1);      const int s = t.GetSummarizedState();      Board u(s);      assert(u == t);    }    //Check two-move states    for (int i=0; i!=8; ++i)    {      Board t;      t.DoMove(i/3,i%3,Player::player1);      t.DoMove(i/3,(i+1)%3,Player::player2);      const int s = t.GetSummarizedState();      Board u(s);      assert(u == t);    }    //Check draw detection    {      Board t;      t.DoMove(1,1,Player::player1);      t.DoMove(0,0,Player::player2);      t.DoMove(1,2,Player::player1);      t.DoMove(1,0,Player::player2);      t.DoMove(2,0,Player::player1);      t.DoMove(0,2,Player::player2);      t.DoMove(0,1,Player::player1);      t.DoMove(2,1,Player::player2);      t.DoMove(2,2,Player::player1);      assert(t.GetWinner() == Winner::draw);    }    //Check player1 wins horizontally detection    {      Board t;      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(0,0,Player::player1);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(1,0,Player::player2);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(0,1,Player::player1);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(1,1,Player::player2);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(0,2,Player::player1);      assert(t.GetWinner() == Winner::player1);    }    //Check player2 wins vertically detection    {      Board t;      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(0,0,Player::player1);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(1,0,Player::player2);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(0,1,Player::player1);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(1,1,Player::player2);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(2,2,Player::player1);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(1,2,Player::player2);      assert(t.GetWinner() == Winner::player2);    }    //Check player1 wins diagonally detection    {      Board t;      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(0,0,Player::player1);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(1,0,Player::player2);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(1,1,Player::player1);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(1,2,Player::player2);      assert(t.GetWinner() == Winner::no_winner);      t.DoMove(2,2,Player::player1);      assert(t.GetWinner() == Winner::player1);    }    //Check no-winner detection    {      Board t;      t.DoMove(1,1,Player::player1);//.........这里部分代码省略.........

voidUrlWrapper::RefsReceived(BMessage* msg){    char buff[B_PATH_NAME_LENGTH];    int32 index = 0;    entry_ref ref;    char* args[] = { const_cast<char*>("urlwrapper"), buff, NULL };    status_t err;    while (msg->FindRef("refs", index++, &ref) == B_OK) {        BFile f(&ref, B_READ_ONLY);        BNodeInfo ni(&f);        BString mimetype;        BString extension(ref.name);        extension.Remove(0, extension.FindLast('.') + 1);        if (f.InitCheck() == B_OK && ni.InitCheck() == B_OK) {            ni.GetType(mimetype.LockBuffer(B_MIME_TYPE_LENGTH));            mimetype.UnlockBuffer();            // Internet Explorer Shortcut            if (mimetype == "text/x-url" || extension == "url") {                // http://filext.com/file-extension/URL                // http://www.cyanwerks.com/file-format-url.html                off_t size;                if (f.GetSize(&size) < B_OK)                    continue;                BString contents;                BString url;                if (f.ReadAt(0LL, contents.LockBuffer(size), size) < B_OK)                    continue;                contents.UnlockBuffer();                while (contents.Length()) {                    BString line;                    int32 cr = contents.FindFirst('/n');                    if (cr < 0)                        cr = contents.Length();                    //contents.MoveInto(line, 0, cr);                    contents.CopyInto(line, 0, cr);                    contents.Remove(0, cr+1);                    line.RemoveAll("/r");                    if (!line.Length())                        continue;                    if (!line.ICompare("URL=", 4)) {                        line.MoveInto(url, 4, line.Length());                        break;                    }                }                if (url.Length()) {                    BUrl u(url.String());                    args[1] = (char*)u.UrlString().String();                    mimetype = kURLHandlerSigBase;                    mimetype += u.Protocol();                    err = be_roster->Launch(mimetype.String(), 1, args + 1);                    if (err != B_OK && err != B_ALREADY_RUNNING)                        err = be_roster->Launch(kAppSig, 1, args + 1);                    continue;                }            }            if (mimetype == "text/x-webloc" || extension == "webloc") {                // OSX url shortcuts                // XML file + resource fork                off_t size;                if (f.GetSize(&size) < B_OK)                    continue;                BString contents;                BString url;                if (f.ReadAt(0LL, contents.LockBuffer(size), size) < B_OK)                    continue;                contents.UnlockBuffer();                int state = 0;                while (contents.Length()) {                    BString line;                    int32 cr = contents.FindFirst('/n');                    if (cr < 0)                        cr = contents.Length();                    contents.CopyInto(line, 0, cr);                    contents.Remove(0, cr+1);                    line.RemoveAll("/r");                    if (!line.Length())                        continue;                    int32 s, e;                    switch (state) {                    case 0:                        if (!line.ICompare("<?xml", 5))                            state = 1;                        break;                    case 1:                        if (!line.ICompare("<plist", 6))                            state = 2;                        break;                    case 2:                        if (!line.ICompare("<dict>", 6))                            state = 3;                        break;                    case 3:                        if (line.IFindFirst("<key>URL</key>") > -1)                            state = 4;                        break;                    case 4:                        if ((s = line.IFindFirst("<string>")) > -1//.........这里部分代码省略.........

  bool is_kuratowski_subgraph(const Graph& g,                              ForwardIterator begin,                               ForwardIterator end,                               VertexIndexMap vm                              )  {    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;    typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator_t;    typedef typename graph_traits<Graph>::edge_descriptor edge_t;    typedef typename graph_traits<Graph>::edges_size_type e_size_t;    typedef typename graph_traits<Graph>::vertices_size_type v_size_t;    typedef typename std::vector<vertex_t> v_list_t;    typedef typename v_list_t::iterator v_list_iterator_t;    typedef iterator_property_map      <typename std::vector<v_list_t>::iterator, VertexIndexMap>       vertex_to_v_list_map_t;    typedef adjacency_list<vecS, vecS, undirectedS> small_graph_t;    detail::target_graph_t target_graph = detail::tg_k_3_3; //unless we decide otherwise later    static small_graph_t K_5(detail::make_K_5<small_graph_t>());    static small_graph_t K_3_3(detail::make_K_3_3<small_graph_t>());    v_size_t n_vertices(num_vertices(g));    v_size_t max_num_edges(3*n_vertices - 5);    std::vector<v_list_t> neighbors_vector(n_vertices);    vertex_to_v_list_map_t neighbors(neighbors_vector.begin(), vm);    e_size_t count = 0;    for(ForwardIterator itr = begin; itr != end; ++itr)      {        if (count++ > max_num_edges)          return false;        edge_t e(*itr);        vertex_t u(source(e,g));        vertex_t v(target(e,g));        neighbors[u].push_back(v);        neighbors[v].push_back(u);      }    for(v_size_t max_size = 2; max_size < 5; ++max_size)      {        vertex_iterator_t vi, vi_end;        for(tie(vi,vi_end) = vertices(g); vi != vi_end; ++vi)          {            vertex_t v(*vi);            //a hack to make sure we don't contract the middle edge of a path            //of four degree-3 vertices            if (max_size == 4 && neighbors[v].size() == 3)              {                if (neighbors[neighbors[v][0]].size() +                    neighbors[neighbors[v][1]].size() +                    neighbors[neighbors[v][2]].size()                    < 11 // so, it has two degree-3 neighbors                    )                  continue;              }            while (neighbors[v].size() > 0 && neighbors[v].size() < max_size)              {                // Find one of v's neighbors u such that that v and u                // have no neighbors in common. We'll look for such a                 // neighbor with a naive cubic-time algorithm since the                 // max size of any of the neighbor sets we'll consider                 // merging is 3                                bool neighbor_sets_intersect = false;                                vertex_t min_u = graph_traits<Graph>::null_vertex();                vertex_t u;                v_list_iterator_t v_neighbor_end = neighbors[v].end();                for(v_list_iterator_t v_neighbor_itr = neighbors[v].begin();                    v_neighbor_itr != v_neighbor_end;                     ++v_neighbor_itr                    )                  {                    neighbor_sets_intersect = false;                    u = *v_neighbor_itr;                    v_list_iterator_t u_neighbor_end = neighbors[u].end();                    for(v_list_iterator_t u_neighbor_itr =                           neighbors[u].begin();                        u_neighbor_itr != u_neighbor_end &&                           !neighbor_sets_intersect;                         ++u_neighbor_itr                        )                      {                        for(v_list_iterator_t inner_v_neighbor_itr =                               neighbors[v].begin();                            inner_v_neighbor_itr != v_neighbor_end; //.........这里部分代码省略.........

bool PelletTransport::evaluate(NOX::Epetra::Interface::Required::FillType fType, 				    const Epetra_Vector* soln, 				    Epetra_Vector* tmp_rhs, 				    Epetra_RowMatrix* tmp_matrix){  if( fType == NOX::Epetra::Interface::Required::Jac )   {    cout << "This problem only works for Finite-Difference or Matrix-Free Based Jacobians."          << " No analytic Jacobian fill available !!" << endl;    throw "Problem ERROR";  }  else     rhs = new Epetra_Vector(*OverlapMap);  // Create the overlapped solution and position vectors  Epetra_Vector u(*OverlapMap);  Epetra_Vector uold(*OverlapMap);  Epetra_Vector xvec(*OverlapNodeMap);  // Export Solution to Overlap vector  uold.Import(*oldSolution, *Importer, Insert);  xvec.Import(*xptr, *nodeImporter, Insert);  u.Import(*soln, *Importer, Insert);  // Declare required variables  int i;  int OverlapNumMyNodes = OverlapNodeMap->NumMyElements();  MaterialPropBase::PropData materialProps;  // Hard-code use of UO2 for now - RWH 5/14/2007  string fixedType = "UO2";    int OverlapMinMyNodeGID;  if (MyPID==0) OverlapMinMyNodeGID = StandardNodeMap->MinMyGID();  else OverlapMinMyNodeGID = StandardNodeMap->MinMyGID()-1;  int row1, row2;  double term1, term2;  double flux1, flux2;  double xx[2];  double uu[2*NUMSPECIES]; // Use of the anonymous enum is needed for SGI builds  double uuold[2*NUMSPECIES];  Basis basis(NumSpecies);    // Zero out the objects that will be filled  rhs->PutScalar(0.0);  ACTIVE_REGIONS propType;  // Loop Over # of Finite Elements on Processor  for( int ne = 0; ne < OverlapNumMyNodes - 1; ++ne )  {    propType = (ACTIVE_REGIONS) (*elemTypes)[ne];    // Loop Over Gauss Points    for( int gp = 0; gp < 2; ++gp )    {      // Get the solution and coordinates at the nodes       xx[0] = xvec[ne];      xx[1] = xvec[ne+1];      for (int k=0; k<NumSpecies; k++) {        uu[NumSpecies * 0 + k] = u[NumSpecies * ne + k];        uu[NumSpecies * 1 + k] = u[NumSpecies * (ne+1) + k];        uuold[NumSpecies * 0 + k] = uold[NumSpecies * ne + k];        uuold[NumSpecies * 1 + k] = uold[NumSpecies * (ne+1) + k];      }      // Calculate the basis function at the gauss point      basis.getBasis(gp, xx, uu, uuold);      MaterialPropFactory::computeProps( propType, basis.uu[0], basis.uu[1], materialProps );      double & rho1    = materialProps.density  ;      double & k1      = materialProps.k_thermal;      double & Cp1     = materialProps.Cp       ;      double & Qstar1  = materialProps.Qstar    ;      double & Qdot1   = materialProps.Qdot     ;      double & Ffunc   = materialProps.thermoF  ;      double & D_diff1 = materialProps.D_diff   ;	                  // Loop over Nodes in Element      for( i = 0; i < 2; ++i )       {	row1=OverlapMap->GID(NumSpecies * (ne+i));	row2=OverlapMap->GID(NumSpecies * (ne+i) + 1);        flux1 = -k1*basis.duu[0]/basis.dx;        flux2 = -0.5*D_diff1*(basis.duu[1]/basis.dx               + basis.uu[1]*Qstar1/(Ffunc*8.3142*basis.uu[0]*basis.uu[0])*basis.duu[0]/basis.dx);        term1 = basis.wt*basis.dx*basis.xx*(                     rho1*Cp1*(basis.uu[0] - basis.uuold[0])/dt * basis.phi[i]                   - flux1*basis.dphide[i]/basis.dx                  - Qdot1*basis.phi[i]                );        term2 = basis.wt*basis.dx*basis.xx*(                    0.5*(basis.uu[1] - basis.uuold[1])/dt * basis.phi[i]                   - flux2*basis.dphide[i]/basis.dx                 );        (*rhs)[NumSpecies*(ne+i)]   += term1;        (*rhs)[NumSpecies*(ne+i)+1] += term2;//.........这里部分代码省略.........

void test_comp() {	// signed gt	assert(i(10)>i(5)==1);	assert(i(5)>i(10)==0);	assert(i(5)>i(5)==0);	assert(i(-10)>i(5)==0);	assert(i(-5)>i(10)==0);	assert(i(-5)>i(5)==0);	assert(i(10)>i(-5)==1);	assert(i(5)>i(-10)==1);	assert(i(5)>i(-5)==1);	assert(i(-10)>i(-5)==0);	assert(i(-5)>i(-10)==1);	assert(i(-5)>i(-5)==0);	// signed ge	assert(i(10)>=i(5)==1);	assert(i(5)>=i(10)==0);	assert(i(5)>=i(5)==1);	assert(i(-10)>=i(5)==0);	assert(i(-5)>=i(10)==0);	assert(i(-5)>=i(5)==0);	assert(i(10)>=i(-5)==1);	assert(i(5)>=i(-10)==1);	assert(i(5)>=i(-5)==1);	assert(i(-10)>=i(-5)==0);	assert(i(-5)>=i(-10)==1);	assert(i(-5)>=i(-5)==1);	// signed lt	assert(i(10)<i(5)==0);	assert(i(5)<i(10)==1);	assert(i(5)<i(5)==0);	assert(i(-10)<i(5)==1);	assert(i(-5)<i(10)==1);	assert(i(-5)<i(5)==1);	assert(i(10)<i(-5)==0);	assert(i(5)<i(-10)==0);	assert(i(5)<i(-5)==0);	assert(i(-10)<i(-5)==1);	assert(i(-5)<i(-10)==0);	assert(i(-5)<i(-5)==0);	// signed le	assert(i(10)<=i(5)==0);	assert(i(5)<=i(10)==1);	assert(i(5)<=i(5)==1);	assert(i(-10)<=i(5)==1);	assert(i(-5)<=i(10)==1);	assert(i(-5)<=i(5)==1);	assert(i(10)<=i(-5)==0);	assert(i(5)<=i(-10)==0);	assert(i(5)<=i(-5)==0);	assert(i(-10)<=i(-5)==1);	assert(i(-5)<=i(-10)==0);	assert(i(-5)<=i(-5)==1);	// signed eq	assert(i(10)==i(5)==0);	assert(i(5)==i(10)==0);	assert(i(5)==i(5)==1);	assert(i(-10)==i(5)==0);	assert(i(-5)==i(10)==0);	assert(i(-5)==i(5)==0);	assert(i(10)==i(-5)==0);	assert(i(5)==i(-10)==0);	assert(i(5)==i(-5)==0);	assert(i(-10)==i(-5)==0);	assert(i(-5)==i(-10)==0);	assert(i(-5)==i(-5)==1);	// signed ne	assert(i(10)!=i(5)==1);	assert(i(5)!=i(10)==1);	assert(i(5)!=i(5)==0);	assert(i(-10)!=i(5)==1);	assert(i(-5)!=i(10)==1);	assert(i(-5)!=i(5)==1);	assert(i(10)!=i(-5)==1);	assert(i(5)!=i(-10)==1);	assert(i(5)!=i(-5)==1);	assert(i(-10)!=i(-5)==1);	assert(i(-5)!=i(-10)==1);	assert(i(-5)!=i(-5)==0);	// unsigned gt	assert(u(10)>u(5)==1);	assert(u(5)>u(10)==0);	assert(u(5)>u(5)==0);	// unsigned ge	assert(u(10)>=u(5)==1);	assert(u(5)>=u(10)==0);	assert(u(5)>=u(5)==1);	// unsigned lt	assert(u(10)<u(5)==0);	assert(u(5)<u(10)==1);	assert(u(5)<u(5)==0);//.........这里部分代码省略.........

//An implementation of the variational pressure projection solve for static geometryvoid FluidSim::solve_pressure(float dt) {      //This linear system could be simplified, but I've left it as is for clarity    //and consistency with the standard naive discretization      int ni = v.ni;   int nj = u.nj;   int system_size = ni*nj;   if(rhs.size() != system_size) {      rhs.resize(system_size);      pressure.resize(system_size);      matrix.resize(system_size);   }   matrix.zero();      //Build the linear system for pressure   for(int j = 1; j < nj-1; ++j) {      for(int i = 1; i < ni-1; ++i) {         int index = i + ni*j;         rhs[index] = 0;         pressure[index] = 0;         float centre_phi = liquid_phi(i,j);         if(centre_phi < 0) {            //right neighbour            float term = u_weights(i+1,j) * dt / sqr(dx);            float right_phi = liquid_phi(i+1,j);            if(right_phi < 0) {               matrix.add_to_element(index, index, term);               matrix.add_to_element(index, index + 1, -term);            }            else {               float theta = fraction_inside(centre_phi, right_phi);               if(theta < 0.01f) theta = 0.01f;               matrix.add_to_element(index, index, term/theta);            }            rhs[index] -= u_weights(i+1,j)*u(i+1,j) / dx;                        //left neighbour            term = u_weights(i,j) * dt / sqr(dx);            float left_phi = liquid_phi(i-1,j);            if(left_phi < 0) {               matrix.add_to_element(index, index, term);               matrix.add_to_element(index, index - 1, -term);            }            else {               float theta = fraction_inside(centre_phi, left_phi);               if(theta < 0.01f) theta = 0.01f;               matrix.add_to_element(index, index, term/theta);            }            rhs[index] += u_weights(i,j)*u(i,j) / dx;                        //top neighbour            term = v_weights(i,j+1) * dt / sqr(dx);            float top_phi = liquid_phi(i,j+1);            if(top_phi < 0) {               matrix.add_to_element(index, index, term);               matrix.add_to_element(index, index + ni, -term);            }            else {               float theta = fraction_inside(centre_phi, top_phi);               if(theta < 0.01f) theta = 0.01f;               matrix.add_to_element(index, index, term/theta);            }            rhs[index] -= v_weights(i,j+1)*v(i,j+1) / dx;                        //bottom neighbour            term = v_weights(i,j) * dt / sqr(dx);            float bot_phi = liquid_phi(i,j-1);            if(bot_phi < 0) {               matrix.add_to_element(index, index, term);               matrix.add_to_element(index, index - ni, -term);            }            else {               float theta = fraction_inside(centre_phi, bot_phi);               if(theta < 0.01f) theta = 0.01f;               matrix.add_to_element(index, index, term/theta);            }            rhs[index] += v_weights(i,j)*v(i,j) / dx;         }      }   }   //Solve the system using Robert Bridson's incomplete Cholesky PCG solver      double tolerance;   int iterations;   bool success = solver.solve(matrix, rhs, pressure, tolerance, iterations);   if(!success) {      printf("WARNING: Pressure solve failed!************************************************/n");   }      //Apply the velocity update   u_valid.assign(0);   for(int j = 0; j < u.nj; ++j) for(int i = 1; i < u.ni-1; ++i) {      int index = i + j*ni;      if(u_weights(i,j) > 0 && (liquid_phi(i,j) < 0 || liquid_phi(i-1,j) < 0)) {         float theta = 1;         if(liquid_phi(i,j) >= 0 || liquid_phi(i-1,j) >= 0)//.........这里部分代码省略.........

inline double initialPosition(Node& p) {	return u(p, 0.);}

Vetor produto_vetorial(Vetor v, Vetor w){	Vetor u(v[1]*w[2] - v[2]*w[1], v[0]*w[2] - v[2]*w[0], v[0]*w[1] - v[1]*w[0]);	return u;}

inline double G1_D(Node& p, double t) {	return u(p, t);}

//-----------------------------------------------------------------------------------------//int main()	{	hash_table_t *t = CreateHashTable();	{	Set( t, u("abc"), 3, u("ABC"), 3 );	uint8_t tmp[4];	size_t rs = 0;	Get( t, u("abc"), 3, tmp, sizeof( tmp ), &rs );	tmp[rs] = '/0';	assert( rs == 3 );	assert( memcmp( tmp, "ABC", rs ) == 0 );	}		{	for( size_t i = 0 ; i < 0xffff; ++i )		{		assert( Set( t, u( &i ), sizeof( i ), u( &i ), sizeof( i ) ) );		}		for( size_t i = 0 ; i < 0xffff; ++i )		{		size_t n = 0;		size_t s = 0;		assert( Get( t, u( &i ), sizeof( i ), u( &n ), sizeof( n ), &s ) == true );		assert( s == sizeof( i ) );		assert( memcmp( &i, &n, s ) == 0 );		}	}	{	for( size_t i = 0 ; i < 0xffff; ++i )		{		Set( t, u( &i ), sizeof( i ), u( &i ), sizeof( i ) );		}	for( size_t i = 0 ; i < 0xffff; ++i )		{		assert( Remove( t, u( &i ), sizeof( i ) ) );		}	for( size_t i = 0 ; i < 0xffff; ++i )		{		size_t n = 0;		assert( Get( t, u( &i ), sizeof( i ), u( &n ), sizeof( n ), NULL ) == false );		}	}	{	assert( Set( t, u( "abc" ), 3, u( "ABC" ), 3 ) );	assert( Append( t, u( "abc" ), 3, u( "DEF" ), 3 ) );	uint8_t tmp[6];	size_t n = 0;	assert( Get( t, u( "abc" ), 3, tmp, 6, &n ) == true );	assert( n == 6 );	assert( memcmp( "ABCDEF", tmp, 6 ) == 0 );	}	DeleteHashTable( t );		return 0;	}

voidrollback( const unsigned g, PrivGlobs& globs ) {    unsigned numX = globs.myX.size(),             numY = globs.myY.size();    unsigned numZ = max(numX,numY);    unsigned i, j;    REAL dtInv = 1.0/(globs.myTimeline[g+1]-globs.myTimeline[g]);    vector<vector<REAL> > u(numY, vector<REAL>(numX));   // [numY][numX]    vector<vector<REAL> > v(numX, vector<REAL>(numY));   // [numX][numY]    vector<REAL> a(numZ), b(numZ), c(numZ), y(numZ);     // [max(numX,numY)]    vector<REAL> yy(numZ);  // temporary used in tridag  // [max(numX,numY)]    //	explicit x    // parallelizable directly since all reads and writes are independent.    // Degree of parallelism: numX*numY.    // TODO: Examine how tiling/shared memory can be used on globs (.myResult).    for(i=0;i<numX;i++) { //par        for(j=0;j<numY;j++) { //par            u[j][i] = dtInv*globs.myResult[i][j];            if(i > 0) {              u[j][i] += 0.5*( 0.5*globs.myVarX[i][j]*globs.myDxx[i][0] )                            * globs.myResult[i-1][j];            }            u[j][i]  +=  0.5*( 0.5*globs.myVarX[i][j]*globs.myDxx[i][1] )                            * globs.myResult[i][j];            if(i < numX-1) {              u[j][i] += 0.5*( 0.5*globs.myVarX[i][j]*globs.myDxx[i][2] )                            * globs.myResult[i+1][j];            }        }    }    //	explicit y    // parallelizable directly since all reads and writes are independent.    // Degree of parallelism: numY*numX.    // TODO: Examine how tiling/shared memory can be used on globs (.myResult).    // and u.?    for(j=0;j<numY;j++) { //par        for(i=0;i<numX;i++) { //par            v[i][j] = 0.0;            if(j > 0) {              v[i][j] +=  ( 0.5*globs.myVarY[i][j]*globs.myDyy[j][0] )                         *  globs.myResult[i][j-1];            }            v[i][j]  +=   ( 0.5*globs.myVarY[i][j]*globs.myDyy[j][1] )                         *  globs.myResult[i][j];            if(j < numY-1) {              v[i][j] +=  ( 0.5*globs.myVarY[i][j]*globs.myDyy[j][2] )                         *  globs.myResult[i][j+1];            }            u[j][i] += v[i][j];        }    }    //	implicit x    // ASSUMING tridag is independent.    // parallelizable directly since all reads and writes are independent.    // Degree of parallelism: numY*numX.    // TODO: Examine tridag    for(j=0;j<numY;j++) { // par        // parallelizable via loop distribution / array expansion.        for(i=0;i<numX;i++) {  // par // here a, b,c should have size [numX]            a[i] =		 - 0.5*(0.5*globs.myVarX[i][j]*globs.myDxx[i][0]);            b[i] = dtInv - 0.5*(0.5*globs.myVarX[i][j]*globs.myDxx[i][1]);            c[i] =		 - 0.5*(0.5*globs.myVarX[i][j]*globs.myDxx[i][2]);        }        // here yy should have size [numX]        tridag(a,b,c,u[j],numX,u[j],yy);    }    //	implicit y    // ASSUMING tridag is independent.    // parallelizable directly since all reads and writes are independent.    // Degree of parallelism: numY*numX.    // TODO: Examine tridag    for(i=0;i<numX;i++) { // par        // parallelizable via loop distribution / array expansion.        for(j=0;j<numY;j++) { // par  // here a, b, c should have size [numY]            a[j] =		 - 0.5*(0.5*globs.myVarY[i][j]*globs.myDyy[j][0]);            b[j] = dtInv - 0.5*(0.5*globs.myVarY[i][j]*globs.myDyy[j][1]);            c[j] =		 - 0.5*(0.5*globs.myVarY[i][j]*globs.myDyy[j][2]);        }        for(j=0;j<numY;j++)            y[j] = dtInv*u[j][i] - 0.5*v[i][j];        // here yy should have size [numY]        tridag(a,b,c,y,numY,globs.myResult[i],yy);    }}

unsigned int cubicSolver(double * ce, double *roots)//cubic equation solver// x^3 + ce[0] x^2 + ce[1] x + ce[2] = 0{    // depressed cubic, Tschirnhaus transformation, x= t - b/(3a)    // t^3 + p t +q =0    unsigned int ret=0;    double shift=(1./3)*ce[0];    double p=ce[1] -shift*ce[0];    double q=ce[0]*( (2./27)*ce[0]*ce[0]-(1./3)*ce[1])+ce[2];    //Cardano's method,    //	t=u+v    //	u^3 + v^3 + ( 3 uv + p ) (u+v) + q =0    //	select 3uv + p =0, then,    //	u^3 + v^3 = -q    //	u^3 v^3 = - p^3/27    //	so, u^3 and v^3 are roots of equation,    //	z^2 + q z - p^3/27 = 0    //	and u^3,v^3 are,    //		-q/2 /pm sqrt(q^2/4 + p^3/27)    //	discriminant= q^2/4 + p^3/27    std::cout<<"cubicSolver:: p="<<p<<"/tq="<<q<<std::endl;    double discriminant= (1./27)*p*p*p+(1./4)*q*q;    if ( fabs(p)< 1.0e-75) {        ret=1;        *roots=(q>0)?-pow(q,(1./3)):pow(-q,(1./3));        *roots -= shift;        return ret;    }    std::cout<<"cubicSolver:: discriminant="<<discriminant<<std::endl;    if(discriminant>0) {        double ce2[2]= {q, -1./27*p*p*p},u3[2];        ret=quadraticSolver(ce2,u3);        if (! ret ) { //should not happen            std::cerr<<"cubicSolver::Error cubicSolver("<<ce[0]<<' '<<ce[1]<<' '<<ce[2]<<")/n";        }        ret=1;        double u,v;        u= (q<=0) ? pow(u3[0], 1./3): -pow(-u3[1],1./3);        //u=(q<=0)?pow(-0.5*q+sqrt(discriminant),1./3):-pow(0.5*q+sqrt(discriminant),1./3);        v=(-1./3)*p/u;        std::cout<<"cubicSolver:: u="<<u<<"/tv="<<v<<std::endl;        std::cout<<"cubicSolver:: u^3="<<u*u*u<<"/tv^3="<<v*v*v<<std::endl;        *roots=u+v - shift;        return ret;    }    ret=3;    std::complex<double> u(q,0),rt[3];    u=pow(-0.5*u-sqrt(0.25*u*u+p*p*p/27),1./3);    rt[0]=u-p/(3.*u)-shift;    std::complex<double> w(-0.5,sqrt(3.)/2);    rt[1]=u*w-p/(3.*u*w)-shift;    rt[2]=u/w-p*w/(3.*u)-shift;//	std::cout<<"Roots:/n";//	std::cout<<rt[0]<<std::endl;//	std::cout<<rt[1]<<std::endl;//	std::cout<<rt[2]<<std::endl;    roots[0]=rt[0].real();    roots[1]=rt[1].real();    roots[2]=rt[2].real();    return ret;}