void TattooMap::showCloseUp(int closeUpNum) {
	Events &events = *_vm->_events;
	Screen &screen = *_vm->_screen;

	// Reset scroll position
	screen._currentScroll = Common::Point(0, 0);

	// Get the closeup images
	Common::String fname = Common::String::format("res%02d.vgs", closeUpNum + 1);
	ImageFile pic(fname);

	Point32 closeUp(_data[closeUpNum].x * 100, _data[closeUpNum].y * 100);
	Point32 delta((SHERLOCK_SCREEN_WIDTH / 2 - closeUp.x / 100) * 100 / CLOSEUP_STEPS,
		(SHERLOCK_SCREEN_HEIGHT / 2 - closeUp.y / 100) * 100 / CLOSEUP_STEPS);
	Common::Rect oldBounds(closeUp.x / 100, closeUp.y / 100, closeUp.x / 100 + 1, closeUp.y / 100 + 1);
	int size = 64;
	int n = 256;
	int deltaVal = 512;
	bool minimize = false;
	int scaleVal, newSize;

	do {
		scaleVal = n;
		newSize = pic[0].sDrawXSize(n);

		if (newSize > size) {
			if (minimize)
				deltaVal /= 2;
			n += deltaVal;
		} else {
			minimize = true;
			deltaVal /= 2;
			n -= deltaVal;
			if (n < 1)
				n = 1;
		}
	} while (deltaVal && size != newSize);

	int deltaScale = (SCALE_THRESHOLD - scaleVal) / CLOSEUP_STEPS;

	for (int step = 0; step < CLOSEUP_STEPS; ++step) {
		Common::Point picSize(pic[0].sDrawXSize(scaleVal), pic[0].sDrawYSize(scaleVal));
		Common::Point pt(closeUp.x / 100 - picSize.x / 2, closeUp.y / 100 - picSize.y / 2);

		restoreArea(oldBounds);
		screen._backBuffer1.transBlitFrom(pic[0], pt, false, 0, scaleVal);

		screen.slamRect(oldBounds);
		screen.slamArea(pt.x, pt.y, picSize.x, picSize.y);

		oldBounds = Common::Rect(pt.x, pt.y, pt.x + picSize.x + 1, pt.y + picSize.y + 1);
		closeUp += delta;
		scaleVal += deltaScale;

		events.wait(1);
	}

	// Handle final drawing of closeup
	// TODO: Handle scrolling
	Common::Rect r(SHERLOCK_SCREEN_WIDTH / 2 - pic[0]._width / 2, SHERLOCK_SCREEN_HEIGHT / 2 - pic[0]._height / 2,
		SHERLOCK_SCREEN_WIDTH / 2 - pic[0]._width / 2 + pic[0]._width,
		SHERLOCK_SCREEN_HEIGHT / 2 - pic[0]._height / 2 + pic[0]._height);

	restoreArea(oldBounds);
	screen._backBuffer1.transBlitFrom(pic[0], Common::Point(r.left, r.top));
	screen.slamRect(oldBounds);
	screen.slamRect(r);
	events.wait(2);
}

//.........这里部分代码省略.........
	else if( (zmsgIs(type,Key) && zmsgIs(which,wheelforward)) || (zmsgIs(type,Key) && !strcmp(msg->getS("which"),",") ) ) {
		newzviewpointScale = zviewpointScale * 0.8f;
		scaleChange = 1;
		zMsgUsed();
	}
	else if( (zmsgIs(type,Key) && zmsgIs(which,wheelbackward)) || (zmsgIs(type,Key) && !strcmp(msg->getS("which"),".") ) ) {
		newzviewpointScale = zviewpointScale * 1.2f;
		scaleChange = 1;
		zMsgUsed();
	}

	if( scaleChange && zviewpointPermitScale ) {
		float x = (float)zMouseMsgX;
		float y = (float)zMouseMsgY;

		// I want the world coord at which the mouse is pointing before the zoom 
		// to be the same as the world coord at which the mouse is pointing after the zoom
		// @TODO: convert his mess into a single function
		DMat4 preScale = zviewpointReferenceModel;
        DMat4 preScaleMat( scale3D( DVec3(zviewpointScale, zviewpointScale, zviewpointScale) ) );
		preScale.cat( preScaleMat );

		DVec3 pre0, pre1;
		gluUnProject( x, y, 0.f, preScale, zviewpointReferenceProj, zviewpointReferenceViewport, &pre0.x, &pre0.y, &pre0.z );
		gluUnProject( x, y, 1.f, preScale, zviewpointReferenceProj, zviewpointReferenceViewport, &pre1.x, &pre1.y, &pre1.z );
		FVec3 wp0 = zviewpointLinePlaneIntersect( FVec3::Origin, FVec3::ZAxis, FVec3((float)pre0.x,(float)pre0.y,(float)pre0.z), FVec3((float)pre1.x,(float)pre1.y,(float)pre1.z) );

		DMat4 postScale = zviewpointReferenceModel;
        DMat4 postScaleMat = ( scale3D( DVec3(newzviewpointScale, newzviewpointScale, newzviewpointScale) ) );
		postScale.cat( postScaleMat );
		
		DVec3 post0, post1;
		gluUnProject( x, y, 0.f, postScale, zviewpointReferenceProj, zviewpointReferenceViewport, &post0.x, &post0.y, &post0.z );
		gluUnProject( x, y, 1.f, postScale, zviewpointReferenceProj, zviewpointReferenceViewport, &post1.x, &post1.y, &post1.z );
		FVec3 wp1 = zviewpointLinePlaneIntersect( FVec3::Origin, FVec3::ZAxis, FVec3((float)post0.x,(float)post0.y,(float)post0.z), FVec3((float)post1.x,(float)post1.y,(float)post1.z) );

		wp1.sub( wp0 );
		zviewpointTrans.add( wp1 );

		zviewpointScale = newzviewpointScale;
	}


	if( zmsgIs(type,Viewpoint_MouseDrag) ) {
		if( zmsgI(releaseDrag) ) {
			zMouseMsgCancelExclusiveDrag();
			zviewpointRotating = zviewpointTranslating = 0;
		}
		else {
			if( zviewpointRotating ) {
				float dx = +0.01f * ( viewpointMouseLast[0] - zmsgF(x) );
				float dy = -0.01f * ( viewpointMouseLast[1] - zmsgF(y) );
				float side = zmsgI( shift ) ? ( (zmsgF(localX) < zmsgF(w)/2) ? -1.f : 1.f ) : 0.f;
					// note: the side doesn't really work if the object is being rendered in a ZUI because the
					// w here refers to the screen, and not the ZUI.
				zviewpointRotateTrackball( dx, dy, side );
				viewpointMouseLast[0] = zmsgF(x);
				viewpointMouseLast[1] = zmsgF(y);
				zMsgUsed();
			}
//			else if( zmsgI(r) ) {
			else if( zviewpointTranslating ) {
				// COMPUTE how big is one pixel at the world z-plane of interest
				// For now, this plane normal to the camera and passes-though the origin
				// Do this by unprojecting rays at the mouse current and last and then
				// solving for the world position at the intersection of that plane

				// CAT in tranforms which take place before the translate
				DMat4 m = zviewpointReferenceModel;
                DMat4 _m( scale3D( DVec3(zviewpointScale, zviewpointScale, zviewpointScale) ) );
				m.cat( _m );
                DMat4 _m1( trans3D( DVec3(zviewpointTrans )) );
				m.cat( _m1 );

				DVec3 p0, p1, p2, p3;
				gluUnProject( zmsgF(x), zmsgF(y), 0.f, m, zviewpointReferenceProj, zviewpointReferenceViewport, &p0.x, &p0.y, &p0.z );
				gluUnProject( zmsgF(x), zmsgF(y), 1.f, m, zviewpointReferenceProj, zviewpointReferenceViewport, &p1.x, &p1.y, &p1.z );
				gluUnProject( viewpointMouseLast[0], viewpointMouseLast[1], 0.f, m, zviewpointReferenceProj, zviewpointReferenceViewport, &p2.x, &p2.y, &p2.z );
				gluUnProject( viewpointMouseLast[0], viewpointMouseLast[1], 1.f, m, zviewpointReferenceProj, zviewpointReferenceViewport, &p3.x, &p3.y, &p3.z );
				FVec3 wp0 = zviewpointLinePlaneIntersect( FVec3::Origin, FVec3::ZAxis, FVec3((float)p0.x,(float)p0.y,(float)p0.z), FVec3((float)p1.x,(float)p1.y,(float)p1.z) );
				FVec3 wp1 = zviewpointLinePlaneIntersect( FVec3::Origin, FVec3::ZAxis, FVec3((float)p2.x,(float)p2.y,(float)p2.z), FVec3((float)p3.x,(float)p3.y,(float)p3.z) );
				wp0.sub( wp1 );

				FVec3 delta( wp0.x, wp0.y, 0.f );

				if( zmsgI(shift) ) {
					delta.x = 0.f;
					delta.y = 0.f;
					delta.z = wp0.y * 4.f;
				}

				zviewpointTrans.add( delta );

				viewpointMouseLast[0] = zmsgF(x);
				viewpointMouseLast[1] = zmsgF(y);
				zMsgUsed();
			}
		}
	}
}

void Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
	ScGeom* geom = static_cast<ScGeom*>(ig.get()); 
	CFpmPhys* phys = static_cast<CFpmPhys*>(ip.get());
	const int &id1 = contact->getId1();
	const int &id2 = contact->getId2();
	Body* b1 = Body::byId(id1,scene).get();
	Body* b2 = Body::byId(id2,scene).get();
	
	Real displN = geom->penetrationDepth; // NOTE: the sign for penetrationdepth is different from ScGeom and Dem3DofGeom: geom->penetrationDepth>0 when spheres interpenetrate
	Real Dtensile=phys->FnMax/phys->kn;
	Real Dsoftening = phys->strengthSoftening*Dtensile; 
	
	/*to set the equilibrium distance between all cohesive elements when they first meet -> allows one to work with initial stress-free assembly*/
	if ( contact->isFresh(scene) ) { phys->initD = displN; phys->normalForce = Vector3r::Zero(); phys->shearForce = Vector3r::Zero();}
	Real D = displN - phys->initD; // interparticular distance is computed depending on the equilibrium distance

	/* Determination of interaction */
	if (D < 0){ //spheres do not touch 
	  if (!phys->isCohesive){ scene->interactions->requestErase(contact); return; } // destroy the interaction before calculation
	  if ((phys->isCohesive) && (abs(D) > (Dtensile + Dsoftening))) { // spheres are bonded and the interacting distance is greater than the one allowed ny the defined cohesion
	    phys->isCohesive=false; 
	    // update body state with the number of broken bonds
	    CFpmState* st1=dynamic_cast<CFpmState*>(b1->state.get());
	    CFpmState* st2=dynamic_cast<CFpmState*>(b2->state.get());
	    st1->numBrokenCohesive+=1;
	    st2->numBrokenCohesive+=1;
	    //// the same thing but from ConcretePM
	    //const shared_ptr<Body>& body1=Body::byId(contact->getId1(),scene), body2=Body::byId(contact->getId2(),scene); assert(body1); assert(body2);
	    //const shared_ptr<CFpmState>& st1=YADE_PTR_CAST<CFpmState>(body1->state), st2=YADE_PTR_CAST<CFpmState>(body2->state);
	    //{ boost::mutex::scoped_lock lock(st1->updateMutex); st1->numBrokenCohesive+=1; }
	    //{ boost::mutex::scoped_lock lock(st2->updateMutex); st2->numBrokenCohesive+=1; }
	    // end of update
	    scene->interactions->requestErase(contact); return;
	  }
	}	  
	
	/*NormalForce*/
	Real Fn=0, Dsoft=0;
		
	if ((D < 0) && (abs(D) > Dtensile)) { //to take into account strength softening
	  Dsoft = D+Dtensile; // Dsoft<0 for a negative value of Fn (attractive force)
	  Fn = -(phys->FnMax+(phys->kn/phys->strengthSoftening)*Dsoft); // computes FnMax - FnSoftening
	}
	else {
	  Fn = phys->kn*D;
	}
	phys->normalForce = Fn*geom->normal;  // NOTE normal is position2-position1 - It is directed from particle1 to particle2
	        
	/*ShearForce*/
	Vector3r& shearForce = phys->shearForce; 
		
	// using scGeom function rotateAndGetShear	
	State* st1 = Body::byId(id1,scene)->state.get();
	State* st2 = Body::byId(id2,scene)->state.get();

	geom->rotate(phys->shearForce);
	const Vector3r& dus = geom->shearIncrement();
	
	//Linear elasticity giving "trial" shear force
	shearForce -= phys->ks*dus;
	// needed for the next timestep
	phys->prevNormal = geom->normal;
		
	/* Morh-Coulomb criterion */
	Real maxFs = phys->FsMax + Fn*phys->tanFrictionAngle;

	if (shearForce.squaredNorm() > maxFs*maxFs){ 
	  shearForce*=maxFs/shearForce.norm(); // to fix the shear force to its yielding value
	}
	
	/* Apply forces */
	Vector3r f = phys->normalForce + shearForce;
	// these lines to adapt to periodic boundary conditions (NOTE applyForceAtContactPoint computes torque induced by normal and shear force too)
	if (!scene->isPeriodic)  
	applyForceAtContactPoint(f , geom->contactPoint , id2, st2->se3.position, id1, st1->se3.position);
	else { // in scg we do not wrap particles positions, hence "applyForceAtContactPoint" cannot be used when scene is periodic
		scene->forces.addForce(id1,-f);
		scene->forces.addForce(id2,f);
		scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(-f));
		scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(-f));
	}
	
	/* Moment Rotation Law */
	// NOTE this part could probably be computed in ScGeom to avoid copy/paste multiplication !!!
	Quaternionr delta( b1->state->ori * phys->initialOrientation1.conjugate() *phys->initialOrientation2 * b2->state->ori.conjugate()); delta.normalize(); //relative orientation
	AngleAxisr aa(delta); // axis of rotation - this is the Moment direction UNIT vector; angle represents the power of resistant ELASTIC moment
	if(aa.angle() > Mathr::PI) aa.angle() -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi 
	  
	phys->cumulativeRotation = aa.angle();
	  
	//Find angle*axis. That's all.  But first find angle about contact normal. Result is scalar. Axis is contact normal.
	Real angle_twist(aa.angle() * aa.axis().dot(geom->normal) ); //rotation about normal
	Vector3r axis_twist(angle_twist * geom->normal);
	Vector3r moment_twist(axis_twist * phys->kr);
	  
	Vector3r axis_bending(aa.angle()*aa.axis() - axis_twist); //total rotation minus rotation about normal
	Vector3r moment_bending(axis_bending * phys->kr);
	Vector3r moment = moment_twist + moment_bending;
	
	Real MomentMax = phys->maxBend*std::fabs(phys->normalForce.norm());
//.........这里部分代码省略.........

/* We resolve the second subproblem through sky-plane projection */
Sequence second_subproblem (Tensor5D& W_1, Tensor5D& W_2, Tensor5D& Y, double& mu, SequenceSet& allSeqs, vector<int> lenSeqs) {
    /*{{{*/
    int numSeq = allSeqs.size();
    int T2 = W_2[0][0].size();
    // reinitialize W_2 to all-zero matrix
    for (int n = 0; REINIT_W_ZERO_TOGGLE and n < numSeq; n ++) {
        int T1 = W_2[n].size();
        for (int i = 0; i < T1; i ++)  
            for (int j = 0; j < T2; j ++) 
                for (int d = 0; d < NUM_DNA_TYPE; d ++) 
                    for (int m = 0; m < NUM_MOVEMENT; m ++) 
                        W_2[n][i][j][d][m] = 0.0;
    }

    vector<Tensor4D> delta (numSeq, Tensor4D(0, Tensor(T2, Matrix(NUM_DNA_TYPE,
                        vector<double>(NUM_MOVEMENT, 0.0)))));  
    tensor5D_init (delta, allSeqs, lenSeqs, T2);
    Tensor tensor (T2, Matrix (NUM_DNA_TYPE, vector<double>(NUM_DNA_TYPE, 0.0)));
    Matrix mat_insertion (T2, vector<double>(NUM_DNA_TYPE, 0.0));

    Trace trace (0, Cell(2)); // 1d: j, 2d: ATCG
    int fw_iter = -1;
    while (fw_iter < MAX_2nd_FW_ITER) {
        fw_iter ++;
        // 1. compute delta
#ifdef PARRALLEL_COMPUTING
//#pragma omp parallel for
#endif
        for (int n = 0; n < numSeq; n ++) {
            int T1 = W_2[n].size();
            for (int i = 0; i < T1; i ++) { 
                for (int j = 0; j < T2; j ++) 
                    for (int d = 0; d < NUM_DNA_TYPE; d ++) 
                        for (int m = 0; m < NUM_MOVEMENT; m ++) {
                            delta[n][i][j][d][m] = -1.0* mu * (W_2[n][i][j][d][m] - W_1[n][i][j][d][m]) + Y[n][i][j][d][m];
#ifdef SECOND_SUBPROBLEM_DEBUG
                            if (delta[n][i][j][d][m] > 0)
                                cout <<"delta: " << n << "," << i << "," << j << "," << d  << "," << m << ": "
                                    << delta[n][i][j][d][m] << endl;
#endif
                            if (m == DELETION_A or m == MATCH_A)
                                tensor[j][d][dna2T3idx('A')] += max(0.0, delta[n][i][j][d][m]);
                            else if (m == DELETION_T or m == MATCH_T)
                                tensor[j][d][dna2T3idx('T')] += max(0.0, delta[n][i][j][d][m]);
                            else if (m == DELETION_C or m == MATCH_C)
                                tensor[j][d][dna2T3idx('C')] += max(0.0, delta[n][i][j][d][m]);
                            else if (m == DELETION_G or m == MATCH_G)
                                tensor[j][d][dna2T3idx('G')] += max(0.0, delta[n][i][j][d][m]);
                            else if (m == DELETION_START or m == MATCH_START)
                                tensor[j][d][dna2T3idx('*')] += max(0.0, delta[n][i][j][d][m]);
                            else if (m == DELETION_END or m == MATCH_END)
                                tensor[j][d][dna2T3idx('#')] += max(0.0, delta[n][i][j][d][m]);
                            else if (m == INSERTION) {
                                mat_insertion[j][d] += max(0.0, delta[n][i][j][d][m]);
                            }
                        }
            }
        }
#ifdef SECOND_SUBPROBLEM_DEBUG
        cout << "tensor transition input list:" << endl;
        for (int j = 0; j < T2; j ++) 
            for (int d = 0; d < NUM_DNA_TYPE; d ++) 
                for (int k = 0; k < NUM_DNA_TYPE; k ++) {
                    if (tensor[j][d][k] > 0)
                        cout << "(" << j << ", " << d << ", " << k << ")=" << tensor[j][d][k] << endl;
                }
#endif

        double delta_square = 0.0;
        for (int n = 0; n < numSeq; n ++) 
            delta_square += tensor4D_frob_prod (delta[n], delta[n]);
        //cout << "delta_square: " << delta_square << endl;
        if ( delta_square < 1e-12 ) {
            //cout << "small delta. early stop." << endl;
            break;
        }

        // 2. determine the trace: run viterbi algorithm
        trace.clear();
	refined_viterbi_algo (trace, tensor, mat_insertion);
        Tensor5D S (numSeq, Tensor4D(0, Tensor(T2, Matrix(NUM_DNA_TYPE, vector<double>(NUM_MOVEMENT, 0.0))))); 
        tensor5D_init (S, allSeqs, lenSeqs, T2);

        // 3. recover values for S 
        // 3b. set a number of selected elements to 1
        for (int t = 0; t < trace.size(); t++) {
            int sj = trace[t].location[0];
            int sd = trace[t].location[1];
            int sm = dna2T3idx(trace[t].acidB);
            // cout << trace[t].acidB;
            for (int n = 0; n < numSeq; n ++) {
                int T1 = S[n].size();
                for (int i = 0; i < T1; i ++) {
                    for (int m = 0; m < NUM_MOVEMENT; m ++)
                        if (delta[n][i][sj][sd][m] > 0.0) { 
                            if (m == DEL_BASE_IDX + sm or m == MTH_BASE_IDX + sm)
                                S[n][i][sj][sd][m] = 1.0;
                            else if (m == INSERTION and trace[t].action == INSERTION) {
                                S[n][i][sj][sd][m] = 1.0;
//.........这里部分代码省略.........

void Foam::inclinedFilmNusseltHeightFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const label patchI = patch().index();

    const scalar t = db().time().timeOutputValue();

    // retrieve the film region from the database

    const regionModels::regionModel& region =
        db().time().lookupObject<regionModels::regionModel>
        (
            "surfaceFilmProperties"
        );

    const regionModels::surfaceFilmModels::kinematicSingleLayer& film =
        dynamic_cast
        <
            const regionModels::surfaceFilmModels::kinematicSingleLayer&
        >(region);

    // calculate the vector tangential to the patch

    const vectorField n(patch().nf());

    const volVectorField& nHat = film.nHat();

    const vectorField nHatp(nHat.boundaryField()[patchI].patchInternalField());

    vectorField nTan(nHatp ^ n);
    nTan /= mag(nTan) + ROOTVSMALL;

    // calculate distance in patch tangential direction

    const vectorField& Cf = patch().Cf();
    scalarField d(nTan & Cf);

    // TODO: currently re-evaluating the entire gTan field to return this patch
    const scalarField gTan(film.gTan()().boundaryField()[patchI] & n);

    if (patch().size() && (max(mag(gTan)) < SMALL))
    {
       WarningIn
        (
            "void Foam::inclinedFilmNusseltHeightFvPatchScalarField::"
            "updateCoeffs()"
        )
            << "Tangential gravity component is zero.  This boundary condition "
            << "is designed to operate on patches inclined with respect to "
            << "gravity"
            << nl;
    }

    const volScalarField& mu = film.mu();
    // The use of the first internal cell for properties doesn't seem correct, kvm
    const scalarField mup(mu.boundaryField()[patchI].patchInternalField());
    // const scalarField mup(mu.boundaryField()[patchI]);

    const volScalarField& rho = film.rho();
    const scalarField rhop(rho.boundaryField()[patchI].patchInternalField());
    // const scalarField rhop(rho.boundaryField()[patchI]);
    
    // calculate the wavy film height

    const scalar GMean = GammaMean_->value(t);
    const scalar a = a_->value(t);
    const scalar omega = omega_->value(t);


    // solve for deltaMean via Bi-Section method
    scalar fxC = 10.0;
    scalar tol = 0.00001;
    label iter=0;
    scalar deltaMeanA = 2e-2;
    scalar deltaMeanB = 2e-6;
    scalar deltaMeanC;

    //  Info << " acoeff " << a << nl;
    scalarField C(pow(3.0*sqr(mup/rhop)/mup/(-gTan + ROOTVSMALL), 0.33333333));
    scalarField delta(deltaMeanA + a*sin(omega*constant::mathematical::twoPi*d));

    while(fabs(fxC)>tol){

        label deltaSize = delta.size();
        reduce(deltaSize, sumOp<label>());
        delta = (deltaMeanA + a*sin(omega*constant::mathematical::twoPi*d));
        scalar fxA = GMean - gSum(pow(delta/C,3))/deltaSize; 
        // DEBUG(delta.size());
        // DEBUG(deltaSize);
        if(fxA > 0.0){
            WarningIn
                (
                 "void Foam::inclinedFilmNusseltHeightFvPatchScalarField::"
                 "updateCoeffs()"
                )
                << "Initial guess for deltaMeanA too low:"
//.........这里部分代码省略.........

    void newton_raphson_constrained_system_method_test()
    {
      typedef eli::mutil::nls::newton_raphson_constrained_system_method<data__, 3, 1> nrcs_type;

      data__ delta(std::sqrt(std::numeric_limits<data__>::epsilon()));
      nrcs_type nrcm;
      typename nrcs_type::solution_matrix rhs, root, x0, x_exact;
      int stat;

      nrcm.set_absolute_tolerance(delta);
      nrcm.set_max_iteration(200);
      nrcm.set_norm_type(nrcs_type::max_norm);
      nrcm.set_lower_condition(0, 0, nrcs_type::NRC_EXCLUSIVE);
      nrcm.set_upper_condition(0, 4, nrcs_type::NRC_EXCLUSIVE);

      // decoupled system
      // set right hand side, initial guess & exact answer
      x_exact << 2, 3, 1;
      rhs = my_decoupled_system_function<data__>(x_exact);
      x0 << 1.5, 2.5, 1.5;
      nrcm.set_initial_guess(x0);

      // test using user defined functions
      stat = nrcm.find_root(root, std::ptr_fun(my_decoupled_system_function<data__>), std::ptr_fun(my_decoupled_system_function_derivative<data__>), rhs);
      TEST_ASSERT(stat==nrcs_type::converged);
      TEST_ASSERT(nrcm.get_iteration_count()<nrcm.get_max_iteration());
      TEST_ASSERT((x_exact-root).norm()<=2*delta);

      nrcm.set_max_iteration(2);
      stat = nrcm.find_root(root, std::ptr_fun(my_decoupled_system_function<data__>), std::ptr_fun(my_decoupled_system_function_derivative<data__>), rhs);
      TEST_ASSERT(stat==nrcs_type::max_iteration);

      // test using functor
      nrcm.set_absolute_tolerance(delta);
      nrcm.set_max_iteration(200);
      nrcm.set_initial_guess(x0);

      stat = nrcm.find_root(root, my_decoupled_system_functor<data__>(), my_decoupled_system_functor_derivative<data__>(), rhs);
      TEST_ASSERT(stat==nrcs_type::converged);
      TEST_ASSERT(nrcm.get_iteration_count()<nrcm.get_max_iteration());
      TEST_ASSERT((x_exact-root).norm()<=2*delta);

      // linear coupled system
      // set right hand side, initial guess & exact answer
      x_exact << 2, 3, 1;
      rhs = my_coupled_linear_system_function<data__>(x_exact);
      x0 << 1.5, 2.5, 1.5;
      nrcm.set_initial_guess(x0);

      // test using user defined functions
      stat = nrcm.find_root(root, std::ptr_fun(my_coupled_linear_system_function<data__>), std::ptr_fun(my_coupled_linear_system_function_derivative<data__>), rhs);
      TEST_ASSERT(stat==nrcs_type::converged);
      TEST_ASSERT(nrcm.get_iteration_count()<2);
      TEST_ASSERT((x_exact-root).norm()<delta);

      // test using functor
      nrcm.set_absolute_tolerance(delta);
      nrcm.set_max_iteration(200);
      nrcm.set_initial_guess(x0);

      stat = nrcm.find_root(root, my_coupled_linear_system_functor<data__>(), my_coupled_linear_system_functor_derivative<data__>(), rhs);
      TEST_ASSERT(stat==nrcs_type::converged);
      TEST_ASSERT(nrcm.get_iteration_count()<2);
      TEST_ASSERT((x_exact-root).norm()<delta);

      // nonlinear coupled system
      // set right hand side, initial guess & exact answer
      x_exact << 2, 3, 1;
      rhs = my_coupled_nonlinear_system_function<data__>(x_exact);
      x0 << 1.5, 2.5, 1.5;
      nrcm.set_initial_guess(x0);

      // test using user defined functions
      stat = nrcm.find_root(root, std::ptr_fun(my_coupled_nonlinear_system_function<data__>), std::ptr_fun(my_coupled_nonlinear_system_function_derivative<data__>), rhs);
      TEST_ASSERT(stat==nrcs_type::converged);
      TEST_ASSERT(nrcm.get_iteration_count()<nrcm.get_max_iteration());
      TEST_ASSERT((x_exact-root).norm()<delta);

      // test using functor
      nrcm.set_absolute_tolerance(delta);
      nrcm.set_max_iteration(200);
      nrcm.set_initial_guess(x0);

      stat = nrcm.find_root(root, my_coupled_nonlinear_system_functor<data__>(), my_coupled_nonlinear_system_functor_derivative<data__>(), rhs);
      TEST_ASSERT(stat==nrcs_type::converged);
      TEST_ASSERT(nrcm.get_iteration_count()<nrcm.get_max_iteration());
      TEST_ASSERT((x_exact-root).norm()<delta);
    }

//.........这里部分代码省略.........
				cProcessed++;
				bMaxLimitPerOneMustBeRecalculated = cProcessed != vProcessed.size ();
				continue;
			}

			if (!bMaxLimitPerOneMustBeRecalculated)
				vLimits [i] = uMaxLimitForThis;
		}
	}

	

	if (cProcessed != vProcessed.size ())
	{
		for (size_t i = 0; i < pItems->vItems.size (); i++)
		{
			if (vProcessed [i])
				continue;
			if (vLimits [i] == UINT64_MAX)
			{
				uEatenTotal = UINT64_MAX;
				break;
			}
			uEatenTotal += vLimits [i]; 
		}
	}

	if (uEatenTotal > uLimit)
		uEatenTotal = uLimit; 

	

	for (size_t i = 0; i < pItems->vItems.size (); i++)
	{
		LOG ("set limit of %I64d for %i item", vLimits [i], i);
		pItems->vItems [i].pDldr->setSpeedLimit (bForDownload, vLimits [i]);
	}

	

	NetworkStat& ns = refgetNetworkStat (bForDownload);

	
	UINT64 uLimit2 = uLimit;

	if (bIsHighestPriorityItems)
	{
		
		if (ns.uCurrentBandwidth < uLimit2)
		{
			uLimit2 = ns.uCurrentBandwidth;
			if (uEatenTotal != UINT64_MAX && uEatenTotal > uLimit2)
				uEatenTotal = uLimit2;
		}
	}

	if (uEatenTotal != UINT64_MAX)
	{
		
		uEatenTotal *= 1.05;
		if (uEatenTotal > uLimit2)
			uEatenTotal = uLimit2;

		if (uLimit2 != UINT64_MAX)
		{
			if (uEatenTotal != uLimit2)
			{
				
				if (uLimit2 / delta (uLimit2, uEatenTotal) >= 10)
					uEatenTotal = uLimit2; 
			}
		}
	}

	if (pItems->totalSpeedExplorer.isExploring ())
	{
		if (uEatenTotal != UINT64_MAX)
		{
			
			uEatenTotal += pItems->totalSpeedExplorer.getAdjustingValue ();
			if (uEatenTotal > uLimit2)
				uEatenTotal = uLimit2;
		}
	}

	*puUnusedLimit = uLimit2 > uEatenTotal ? uLimit2 - uEatenTotal : 0;

	if (*puUnusedLimit == 0)
		pItems->state |= ManageForSpeedItemsList::MSILS_USES_ALL_BANDWIDTH_PROVIDED;

	if (ns.currentBandwidthExplorer.isExploring ())
	{
		*puUnusedLimit += ns.currentBandwidthExplorer.getAdjustingValue ();
		if (*puUnusedLimit > uLimit)
		{
			assert (uEatenTotal <= uLimit);
			*puUnusedLimit = uLimit - uEatenTotal;
		}
	}
}

bool Weapon::update(Game* _game)
{
	int clip = 1;
	if(mClipSize > 0)
		clip = mClip;

	mCounter += _game->getDelta();
	if(mActive && mCounter >= mShootRate && clip > 0)
	{
		if(_game->mousePressed("left"))
		{
			// launch bullets
			EventManager* evm = EventManager::getSingleton();
			sf::Vector2i mousePosition = evm->getMousePos(_game);

			bool canShoot = false;
			if(mousePosition.y < mOwner->getPosition().y + mOwner->getAABB().height)
				canShoot = true;

			if(canShoot)
			{
				sf::View view = _game->getCamera();
				sf::Vector2f mousePosTrans = _game->getWindow()->mapPixelToCoords(mousePosition, view);
				sf::Vector2f shootPoint(mX + mShootPoint.x - 8, mY + mShootPoint.y);

				sf::Vector2f delta(mousePosTrans.x - shootPoint.x, mousePosTrans.y - shootPoint.y);
				float angle = atan2(delta.y, delta.x);
				float power = 500.f;

				int id = 7;
				if(mName == "FireballLauncher")
					id = 23;

				Bullet* bullet = new Bullet(id, shootPoint.x, shootPoint.y, sf::Vector2f(cos(angle) * power, sin(angle) * power), mDamage, EMITTERTYPE_Circle);
				bullet->setUniqueType("Bullet");
				bullet->setOwner(mOwner);
				ignore(bullet);
				mOwner->ignore(bullet);
				_game->pushObject(bullet);
				static_cast<Player*>(mOwner)->changeAnim("forward");

				// update weapon
				mClip--;
				mCounter = 0.f;
				std::stringstream ss;
				ss << "Clip: " << mClip << "/" << mClipCapacity;
				Text* text = static_cast<Text*>(_game->getInterface()->getInterface(2, "clipText"));
				text->setString(ss.str());
				
				if(mName == "FireballLauncher")
					_game->playSound("fireball");
				else if(mName == "Pistol")
					_game->playSound("pistol");
				else if(mName == "Shotgun")
					_game->playSound("shotgun");
			}
		}
	}

	mReloadCount += _game->getDelta();
	if(/*_game->keyPressed("r") &&*/ mClip <= 0 && !mReload /*&& mClip < mClipSize*/)
	{
		mReload   = true;
		mLastClip = mClip;
		mClip	  = 0;
		if(mReloadVisual)
			mReloadVisual->setActive(true);
	}

	if(mReload && mReloadCount > mReloadTime)
	{
		int capacity = mClipSize - mLastClip;
		if(capacity >= 0)
			mClipCapacity  -= mClipSize - mLastClip;
		else if(capacity < 0)
			mClipCapacity = 0;

		if(mClipCapacity > 0)
		{
			mClip			= mClipSize;
			mReloadCount	= .0f;
			mReload			= false;

			std::stringstream ss;
			ss << "Clip: " << mClip << "/" << mClipCapacity;
			Text* text = static_cast<Text*>(_game->getInterface()->getInterface(2, "clipText"));
			text->setString(ss.str());
			_game->playSound("reload");
			_game->playSound("reload");
		}

		if(mReloadVisual)
			mReloadVisual->setActive(false);
	}

	stringstream ss;
	ss << "Weapon Damage: " << mDamage;
	static_cast<Text*>(_game->getInterface()->getInterface(2, "damageText"))->setString(ss.str());
	
	Object::update(_game);
//.........这里部分代码省略.........

void QGraphicsWidgetPrivate::windowFrameMouseMoveEvent(QGraphicsSceneMouseEvent *event)
{
    Q_Q(QGraphicsWidget);
    ensureWindowData();
    if (!(event->buttons() & Qt::LeftButton) || windowData->hoveredSubControl != QStyle::SC_TitleBarLabel)
        return;

    QLineF delta(q->mapFromScene(event->buttonDownScenePos(Qt::LeftButton)), event->pos());
    QLineF parentDelta(q->mapToParent(delta.p1()), q->mapToParent(delta.p2()));
    QLineF parentXDelta(q->mapToParent(QPointF(delta.p1().x(), 0)), q->mapToParent(QPointF(delta.p2().x(), 0)));
    QLineF parentYDelta(q->mapToParent(QPointF(0, delta.p1().y())), q->mapToParent(QPointF(0, delta.p2().y())));

    QRectF newGeometry;
    switch (windowData->grabbedSection) {
    case Qt::LeftSection:
        newGeometry = QRectF(windowData->startGeometry.topLeft()
                             + QPointF(parentXDelta.dx(), parentXDelta.dy()),
                             windowData->startGeometry.size() - QSizeF(delta.dx(), delta.dy()));
        break;
    case Qt::TopLeftSection:
        newGeometry = QRectF(windowData->startGeometry.topLeft()
                             + QPointF(parentDelta.dx(), parentDelta.dy()),
                             windowData->startGeometry.size() - QSizeF(delta.dx(), delta.dy()));
        break;
    case Qt::TopSection:
        newGeometry = QRectF(windowData->startGeometry.topLeft()
                             + QPointF(parentYDelta.dx(), parentYDelta.dy()),
                             windowData->startGeometry.size() - QSizeF(0, delta.dy()));
        break;
    case Qt::TopRightSection:
        newGeometry = QRectF(windowData->startGeometry.topLeft()
                             + QPointF(parentYDelta.dx(), parentYDelta.dy()),
                             windowData->startGeometry.size() - QSizeF(-delta.dx(), delta.dy()));
        break;
    case Qt::RightSection:
        newGeometry = QRectF(windowData->startGeometry.topLeft(),
                             windowData->startGeometry.size() + QSizeF(delta.dx(), 0));
        break;
    case Qt::BottomRightSection:
        newGeometry = QRectF(windowData->startGeometry.topLeft(),
                             windowData->startGeometry.size() + QSizeF(delta.dx(), delta.dy()));
        break;
    case Qt::BottomSection:
        newGeometry = QRectF(windowData->startGeometry.topLeft(),
                             windowData->startGeometry.size() + QSizeF(0, delta.dy()));
        break;
    case Qt::BottomLeftSection:
        newGeometry = QRectF(windowData->startGeometry.topLeft()
                             + QPointF(parentXDelta.dx(), parentXDelta.dy()),
                             windowData->startGeometry.size() - QSizeF(delta.dx(), -delta.dy()));
        break;
    case Qt::TitleBarArea:
        newGeometry = QRectF(windowData->startGeometry.topLeft()
                             + QPointF(parentDelta.dx(), parentDelta.dy()),
                             windowData->startGeometry.size());
        break;
    case Qt::NoSection:
        break;
    }

    if (windowData->grabbedSection != Qt::NoSection) {
        _q_boundGeometryToSizeConstraints(windowData->startGeometry, &newGeometry,
                                          windowData->grabbedSection,
                                          q->effectiveSizeHint(Qt::MinimumSize),
                                          q->effectiveSizeHint(Qt::MaximumSize),
                                          q);
        q->setGeometry(newGeometry);
    }
}

osg::Node *createHUD(osgText::Text *updateText)
{
	osg::Camera *hudCamera = new osg::Camera;
	hudCamera->setReferenceFrame(osg::Transform::ABSOLUTE_RF);
	hudCamera->setProjectionMatrixAsOrtho2D(0,1280,0,1024);
	hudCamera->setViewMatrix(osg::Matrix::identity());
	hudCamera->setRenderOrder(osg::Camera::POST_RENDER);
	hudCamera->setClearMask(GL_DEPTH_BUFFER_BIT);

	//std::string timesFont(

	osg::Vec3 position(150.f,800.0f,0.0f);
	osg::Vec3 delta(0.0f,-60.0f,0.0f);

	{
		osg::Geode *geode = new osg::Geode;
		osg::StateSet *stateset = geode->getOrCreateStateSet();
		stateset->setMode(GL_LIGHTING,osg::StateAttribute::OFF);
		stateset->setMode(GL_DEPTH_TEST,osg::StateAttribute::OFF);
		geode->setName("simple");
		hudCamera->addChild(geode);

		osgText::Text *text = new osgText::Text;
		geode->addDrawable( text );

		text->setText("Picking in head up Displays is simple!");
		text->setCharacterSize(10.f);
		text->setPosition(position);

		position+=delta;
	}

	for ( int i=0; i<5; ++i ) {
		osg::Vec3 dy(0.0f,-30.0f,0.0f);
		osg::Vec3 dx(120.0f,0.0f,0.0f);
		osg::Geode *geode = new osg::Geode;
		osg::StateSet *stateset = geode->getOrCreateStateSet();
		const char *opts[] = {"one", "Two", "Three", "January", "Feb", "2003"};
		osg::Geometry *quad = new osg::Geometry;
		stateset->setMode(GL_LIGHTING,osg::StateAttribute::OFF);
		stateset->setMode(GL_DEPTH_TEST,osg::StateAttribute::OFF);
		std::string name = "subOption";
		name+=" ";
		name+=std::string(opts[i]);
		geode->setName(name);
		osg::Vec3Array *vertices = new osg::Vec3Array(4);
		osg::Vec4Array *colors = new osg::Vec4Array;
		colors->push_back(osg::Vec4(0.4-0.1*i,0.1*i,0.2*i,1.0));
		quad->setColorArray(colors);
		quad->setColorBinding(osg::Geometry::BIND_PER_PRIMITIVE);
		(*vertices)[0] = position;
		(*vertices)[1] = position+dx;
		(*vertices)[2] = position+dx+dy;
		(*vertices)[3] = position+dy;
		quad->setVertexArray(vertices);
		quad->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::QUADS,0,4));
		geode->addDrawable(quad);
		hudCamera->addChild(geode);

		position+=delta;
	}

	{
		osg::Geode *geode = new osg::Geode;
		osg::StateSet *stateset = geode->getOrCreateStateSet();
		stateset->setMode(GL_LIGHTING,osg::StateAttribute::OFF);
		stateset->setMode(GL_DEPTH_TEST,osg::StateAttribute::OFF);
		geode->setName("The text Label");
		geode->addDrawable(updateText);
		hudCamera->addChild(geode);

		updateText->setCharacterSize(10.f);
		updateText->setColor(osg::Vec4(1.0f,1.0f,0.0f,1.0f));
		updateText->setText("");
		updateText->setPosition(position);
		updateText->setDataVariance(osg::Object::DYNAMIC);

		position += delta;
	}

	return hudCamera;
}

        void onEmitCode(EmitArgs& args, GrGPArgs* gpArgs) override {
            const PLSQuadEdgeEffect& qe = args.fGP.cast<PLSQuadEdgeEffect>();
            GrGLSLVertexBuilder* vsBuilder = args.fVertBuilder;
            GrGLSLVaryingHandler* varyingHandler = args.fVaryingHandler;
            GrGLSLUniformHandler* uniformHandler = args.fUniformHandler;

            // emit attributes
            varyingHandler->emitAttributes(qe);

            GrGLSLVertToFrag uv(kVec2f_GrSLType);
            varyingHandler->addVarying("uv", &uv, kHigh_GrSLPrecision);
            vsBuilder->codeAppendf("%s = %s;", uv.vsOut(), qe.inUV()->fName);

            GrGLSLVertToFrag ep1(kVec2f_GrSLType);
            varyingHandler->addVarying("endpoint1", &ep1, kHigh_GrSLPrecision);
            vsBuilder->codeAppendf("%s = vec2(%s.x, %s.y);", ep1.vsOut(), 
                                  qe.inEndpoint1()->fName, qe.inEndpoint1()->fName);

            GrGLSLVertToFrag ep2(kVec2f_GrSLType);
            varyingHandler->addVarying("endpoint2", &ep2, kHigh_GrSLPrecision);
            vsBuilder->codeAppendf("%s = vec2(%s.x, %s.y);", ep2.vsOut(), 
                                  qe.inEndpoint2()->fName, qe.inEndpoint2()->fName);

            GrGLSLVertToFrag delta(kVec2f_GrSLType);
            varyingHandler->addVarying("delta", &delta, kHigh_GrSLPrecision);
            vsBuilder->codeAppendf("%s = vec2(%s.x - %s.x, %s.y - %s.y) * 0.5;", 
                                   delta.vsOut(), ep1.vsOut(), ep2.vsOut(), ep2.vsOut(), 
                                   ep1.vsOut());

            GrGLSLVertToFrag windings(kInt_GrSLType);
            varyingHandler->addFlatVarying("windings", &windings, kLow_GrSLPrecision);
            vsBuilder->codeAppendf("%s = %s;", 
                                   windings.vsOut(), qe.inWindings()->fName);

            // Setup position
            this->setupPosition(vsBuilder, gpArgs, qe.inPosition()->fName);

            // emit transforms
            this->emitTransforms(vsBuilder, varyingHandler, uniformHandler, gpArgs->fPositionVar, 
                                 qe.inPosition()->fName, qe.localMatrix(), args.fTransformsIn, 
                                 args.fTransformsOut);

            GrGLSLFragmentBuilder* fsBuilder = args.fFragBuilder;
            SkAssertResult(fsBuilder->enableFeature(
                           GrGLSLFragmentShaderBuilder::kPixelLocalStorage_GLSLFeature));
            SkAssertResult(fsBuilder->enableFeature(
                    GrGLSLFragmentShaderBuilder::kStandardDerivatives_GLSLFeature));
            static const int QUAD_ARGS = 2;
            GrGLSLShaderVar inQuadArgs[QUAD_ARGS] = {
                GrGLSLShaderVar("dot", kFloat_GrSLType, 0, kHigh_GrSLPrecision),
                GrGLSLShaderVar("uv", kVec2f_GrSLType, 0, kHigh_GrSLPrecision)
            };
            SkString inQuadName;

            const char* inQuadCode = "if (uv.x * uv.x <= uv.y) {"
                                     "return dot >= 0.0;"
                                     "} else {"
                                     "return false;"
                                     "}";
            fsBuilder->emitFunction(kBool_GrSLType, "in_quad", QUAD_ARGS, inQuadArgs, inQuadCode, 
                                    &inQuadName);
            fsBuilder->declAppendf(GR_GL_PLS_PATH_DATA_DECL);
            // keep the derivative instructions outside the conditional
            fsBuilder->codeAppendf("highp vec2 uvdX = dFdx(%s);", uv.fsIn());
            fsBuilder->codeAppendf("highp vec2 uvdY = dFdy(%s);", uv.fsIn());
            fsBuilder->codeAppend("highp vec2 uvIncX = uvdX * 0.45 + uvdY * -0.1;");
            fsBuilder->codeAppend("highp vec2 uvIncY = uvdX * 0.1 + uvdY * 0.55;");
            fsBuilder->codeAppendf("highp vec2 uv = %s.xy - uvdX * 0.35 - uvdY * 0.25;", 
                                   uv.fsIn());
            fsBuilder->codeAppendf("highp vec2 firstSample = %s.xy - vec2(0.25);",
                                   fsBuilder->fragmentPosition());
            fsBuilder->codeAppendf("highp float d = dot(%s, (firstSample - %s).yx) * 2.0;", 
                                   delta.fsIn(), ep1.fsIn());
            fsBuilder->codeAppendf("pls.windings[0] += %s(d, uv) ? %s : 0;", inQuadName.c_str(), 
                                   windings.fsIn());
            fsBuilder->codeAppend("uv += uvIncX;");
            fsBuilder->codeAppendf("d += %s.x;", delta.fsIn());
            fsBuilder->codeAppendf("pls.windings[1] += %s(d, uv) ? %s : 0;", inQuadName.c_str(), 
                                   windings.fsIn());
            fsBuilder->codeAppend("uv += uvIncY;");
            fsBuilder->codeAppendf("d += %s.y;", delta.fsIn());
            fsBuilder->codeAppendf("pls.windings[2] += %s(d, uv) ? %s : 0;", inQuadName.c_str(), 
                                   windings.fsIn());
            fsBuilder->codeAppend("uv -= uvIncX;");
            fsBuilder->codeAppendf("d -= %s.x;", delta.fsIn());
            fsBuilder->codeAppendf("pls.windings[3] += %s(d, uv) ? %s : 0;", inQuadName.c_str(), 
                                   windings.fsIn());
        }

Vector MultiColvarBase::getSeparation( const Vector& vec1, const Vector& vec2 ) const {
  if(usepbc){ return pbcDistance( vec1, vec2 ); }
  else{ return delta( vec1, vec2 ); }
}

 BBox BBox::expand(float f) {
     Vector3 delta(f, f, f);
     pMin -= delta;
     pMax += delta;
     return *this;
 }

//.........这里部分代码省略.........
            )   << "Mesh modifiers not read properly"
                << abort(FatalError);
        }

        setVirtualPistonPosition();
        checkAndCalculate();

        return;


    }

    Info << "checkAndCalculate()" << endl;
    checkAndCalculate();

    Info<< "Time = " << engTime().theta() << endl
        << "Adding zones to the engine mesh" << endl;


    //fz = 4: virtual piston, outSidePort, insidePort, cutFaceZone
    //pz = 2: piston points, cutPointZone
    //cz = 1: moving mask

    List<pointZone*> pz(3);
    List<faceZone*> fz(4);
    List<cellZone*> cz(1);

    label nPointZones = 0;
    label nFaceZones = 0;
    label nCellZones = 0;

    // Add the piston zone
    if (piston().patchID().active())
    {

        // Piston position

        Info << "Adding face zone for piston layer addition/removal" <&l