void GnuplotExportModule::outputBoundaryCondition(PrescribedGradientBCNeumann &iBC, TimeStep *tStep)
{
    FloatArray stress;
    iBC.computeField(stress, tStep);

    printf("Mean stress computed in Gnuplot export module: "); stress.printYourself();

    double time = 0.0;

    TimeStep *ts = emodel->giveCurrentStep();
    if ( ts != NULL ) {
        time = ts->giveTargetTime();
    }

    int bcIndex = iBC.giveNumber();

    std :: stringstream strMeanStress;
    strMeanStress << "PrescribedGradientGnuplotMeanStress" << bcIndex << "Time" << time << ".dat";
    std :: string nameMeanStress = strMeanStress.str();
    std::vector<double> componentArray, stressArray;

    for(int i = 1; i <= stress.giveSize(); i++) {
        componentArray.push_back(i);
        stressArray.push_back(stress.at(i));
    }

    XFEMDebugTools::WriteArrayToGnuplot(nameMeanStress, componentArray, stressArray);

    // Homogenized strain
    
    FloatArray grad;
    iBC.giveGradientVoigt(grad);
    outputGradient(iBC.giveNumber(), *iBC.giveDomain(), grad, tStep);
}

void FE2FluidMaterialStatus :: setTimeStep(TimeStep *tStep)
{
    TimeStep *rveTStep = this->rve->giveCurrentStep(); // Should i create a new one if it is empty?
    rveTStep->setNumber( tStep->giveNumber() );
    rveTStep->setTime( tStep->giveTargetTime() );
    rveTStep->setTimeIncrement( tStep->giveTimeIncrement() );
}

void main()
{
	LOG("Burger Bash running on the Sifteo SDK\n");
	
	EventHandler eventHandler;
	eventHandler.init();
	
	// Plate cube initialize
	pCube.init(gNumCubes - 1);
	
	// Init all other cubes
	for (int i = 0; i < gNumCubes - 1; i++)
	{
		cubes[i].init(i);
	}
	
	// Game loop
	TimeStep ts;
    while (1) {
		for (int i = 0; i < gNumCubes; i++)
		{
			cubes[i].update();
		}
		
        System::paint();
        ts.next();
    }
    
    //AudioTracker::play(Music);
}

void GnuplotExportModule::outputXFEMGeometry(const std::vector< std::vector<FloatArray> > &iEnrItemPoints)
{
    double time = 0.0;

    TimeStep *ts = emodel->giveCurrentStep();
    if ( ts != NULL ) {
        time = ts->giveTargetTime();
    }

    std :: stringstream strCracks;
    strCracks << "CracksTime" << time << ".dat";
    std :: string nameCracks = strCracks.str();
    WritePointsToGnuplot(nameCracks, iEnrItemPoints);
}

contextIOResultType
PrimaryField :: restoreContext(DataStream *stream, ContextMode mode)
{
    int i, class_id;
    contextIOResultType iores;
    // read class header
    if ( !stream->read(& class_id, 1) ) {
        return CIO_IOERR;
    }

    if ( class_id != PrimaryFieldClass ) {
        return CIO_BADVERSION;
    }

    if ( !stream->read(& actualStepNumber, 1) ) {
        THROW_CIOERR(CIO_IOERR);
    }

    if ( !stream->read(& actualStepIndx, 1) ) {
        THROW_CIOERR(CIO_IOERR);
    }

    for ( i = 0; i <= nHistVectors; i++ ) {
        if ( ( iores = solutionVectors.at(i + 1)->restoreYourself(stream, mode) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
    }

    int flag;
    TimeStep *iStep;
    for ( i = 0; i <= nHistVectors; i++ ) {
        if ( !stream->read(& flag, 1) ) {
            THROW_CIOERR(CIO_IOERR);
        }

        if ( flag ) {
            iStep = new TimeStep(emodel);
            if ( ( iores = iStep->restoreContext(stream, mode) ) != CIO_OK ) {
                THROW_CIOERR(iores);
            }
        } else {
            iStep = NULL;
        }

        solStepList.put(i + 1, iStep);
    }

    return CIO_OK;
}

double NLTransientTransportProblem :: giveUnknownComponent(ValueModeType mode, TimeStep *tStep, Domain *d, Dof *dof)
// returns unknown quantity like displacement, velocity of equation
// This function translates this request to numerical method language
{
    if ( this->requiresUnknownsDictionaryUpdate() ) {
        if ( mode == VM_Incremental ) { //get difference between current and previous time variable
            return dof->giveUnknowns()->at(0) - dof->giveUnknowns()->at(1);
        }
        int hash = this->giveUnknownDictHashIndx(mode, tStep);
        if ( dof->giveUnknowns()->includes(hash) ) {
            return dof->giveUnknowns()->at(hash);
        } else {
            OOFEM_ERROR("Dof unknowns dictionary does not contain unknown of value mode (%s)", __ValueModeTypeToString(mode) );
        }
    }

    double t = tStep->giveTargetTime();
    TimeStep *previousStep = this->givePreviousStep(), *currentStep = this->giveCurrentStep();

    if ( dof->__giveEquationNumber() == 0 ) {
        OOFEM_ERROR("invalid equation number on DoF %d", dof->giveDofID() );
    }

    if ( ( t >= previousStep->giveTargetTime() ) && ( t <= currentStep->giveTargetTime() ) ) {
        ///@todo Shouldn't it be enough to just run this?
        //UnknownsField->giveUnknownValue(dof, mode, currentStep);
        double rtdt = UnknownsField->giveUnknownValue(dof, VM_Total, currentStep);
        double rt   = UnknownsField->giveUnknownValue(dof, VM_Total, previousStep);
        double psi = ( t - previousStep->giveTargetTime() ) / currentStep->giveTimeIncrement();
        if ( mode == VM_Velocity ) {
            return ( rtdt - rt ) / currentStep->giveTimeIncrement();
        } else if ( mode == VM_Total ) {
            return psi * rtdt + ( 1. - psi ) * rt;
        } else if ( mode == VM_Incremental ) {
            if ( previousStep->isIcApply() ) {
                return 0;
            } else {
                return ( rtdt - rt );
            }
        } else {
            OOFEM_ERROR("Unknown mode %s is undefined for this problem", __ValueModeTypeToString(mode) );
        }
    } else {
        OOFEM_ERROR("time value %f not within bounds %f and %f", t, previousStep->giveTargetTime(), currentStep->giveTargetTime() );
    }

    return 0.; // to make compiler happy;
}

void PlayerComponent::OnUpdate(TimeStep timeStep)
{
	Application& app = Application::GetApplication();
	float speed = (float)(0.1 * timeStep.GetMills());
	vec3 pos = m_Transform->m_Transform.GetPosition();

	if (app.IsKeyPressed(GLFW_KEY_UP))
	{
		pos.y += speed;
	}
	else if (app.IsKeyPressed(GLFW_KEY_DOWN))
	{
		pos.y -= speed;
	}

	if (app.IsKeyPressed(GLFW_KEY_LEFT))
	{
		pos.x -= speed;
	}
	else if (app.IsKeyPressed(GLFW_KEY_RIGHT))
	{
		pos.x += speed;
	}

	m_Transform->m_Transform.SetPosition(pos);
}

void main()
{
    assetConfig.append(MainSlot, BetterflowAssets);
    loader.init();

    for (CubeID cube : CubeSet::connected())
    {
        onConnect(NULL, cube);
    }
    Events::cubeConnect.set(&onConnect);

    gVideo[0].setMode(BG0);

    MainMenu gameMenu(&gVideo[0]);

    int colour;
    Colormap colourList;
    colourList.setEGA();

    

    while(1)
    {
        colour = gameMenu.runMenu();

        gVideo[0].setMode(SOLID_MODE);

        gVideo[0].colormap[0].set(colourList[colour].get());
        System::paint();

        TimeStep ts;
        TimeTicker ticker = TimeTicker(1);
        int tickIncr;
        int tickCount  = 0;
        do
        {
            tickIncr = ticker.tick( ts.delta() );
            ts.next();
            tickCount += tickIncr;
        } while (tickCount < 3);

        gVideo[0].setMode(BG0);
        gVideo[0].bg0.erase(StripeTile);
    }

}

int  TimeIntegrationScheme :: requiresNewLhsAt (TimeStep* stepN)
   // Returns True if the time increment deltaT of stepN differs from that
   // of the previous step, else returns False. This is the default imple-
   // mentation of the method.
{
   TimeStep* pastStep ;
   int       n,answer ;

   n = stepN->giveNumber() ;
   if (n <= 1)
      answer = TRUE ;
   else {
      pastStep = new TimeStep(n-1,this) ;
      answer = (pastStep->giveTimeIncrement() != stepN->giveTimeIncrement());
      delete pastStep ;}

   return answer ;
}

void GnuplotExportModule::outputMesh(Domain &iDomain)
{
    std::vector< std::vector<FloatArray> > pointArray;

    if(iDomain.giveNumberOfSpatialDimensions() == 2) {
        // Write all element edges to gnuplot
        for ( auto &el : iDomain.giveElements() ) {
            int numEdges = el->giveNumberOfNodes();


            for ( int edgeIndex = 1; edgeIndex <= numEdges; edgeIndex++ ) {
                std::vector<FloatArray> points;

                IntArray bNodes;
                el->giveInterpolation()->boundaryGiveNodes(bNodes, edgeIndex);

                int niLoc = bNodes.at(1);
                const FloatArray &xS = *(el->giveNode(niLoc)->giveCoordinates() );
                points.push_back(xS);

                int njLoc = bNodes.at( bNodes.giveSize() );
                const FloatArray &xE = *(el->giveNode(njLoc)->giveCoordinates() );
                points.push_back(xE);

                pointArray.push_back(points);
            }

        }


        double time = 0.0;

        TimeStep *ts = emodel->giveCurrentStep();
        if ( ts != NULL ) {
            time = ts->giveTargetTime();
        }

        std :: stringstream strMesh;
        strMesh << "MeshGnuplotTime" << time << ".dat";
        std :: string nameMesh = strMesh.str();

        WritePointsToGnuplot(nameMesh, pointArray);
    }
}

void
StaticFracture :: solveYourself()
{
    MetaStep *activeMStep;
    FILE *out = this->giveOutputStream();
    this->timer.startTimer(EngngModelTimer :: EMTT_AnalysisTimer);
    this->giveNumberOfSlaveProblems();
    
    this->timer.startTimer(EngngModelTimer :: EMTT_SolutionStepTimer);
    this->timer.initTimer(EngngModelTimer :: EMTT_NetComputationalStepTimer);


    int numMetaSteps = this->giveNumberOfMetaSteps();
    for (int imstep = 1; imstep <= numMetaSteps; imstep++) { // don't know what will happen if we have several meta steps?
        activeMStep = this->giveMetaStep(imstep);

        int nTimeSteps = activeMStep->giveNumberOfSteps();
        for ( int tStepNum = 1; tStepNum <= nTimeSteps; tStepNum++ ) { //loop over time steps in opt analysis

            for ( int subProb = 1; subProb <= this->giveNumberOfSlaveProblems(); subProb++ ) {
        
                EngngModel *sp = this->giveSlaveProblem(subProb); 
                sp->solveYourself();

                //this->updateYourself( this->giveCurrentStep()); // not neccessary
            
                // optimization
                this->optimize( this->giveCurrentStep() );    

                // Resetting the time step number for each sp after each optimization time step
                TimeStep *tStep =  sp->giveCurrentStep();
                tStep->setNumber(tStepNum); 
                sp->giveExportModuleManager()->doOutput(tStep); // turn off export during regular analysis
            
                tStep->setNumber(0); // otherwise the anlysis wont restart at time 0
                }
        }
    }
    

}

void TimeSteps::clear()
{
    Tools::ClassManager *cm = Tools::ClassManager::getInstance();
    vector<TimeStep *>::iterator tsIterator = timeSteps.begin();

    while (tsIterator != timeSteps.end())
    {
        TimeStep *timeStep = *tsIterator;

        if (timeStep != NULL)
        {
            timeStep->clear();
            cm->deleteObject(timeStep->getID());

            timeStep = NULL;
        }

        ++tsIterator;
    }

    timeSteps.clear();
}

void GnuplotExportModule::outputXFEM(Crack &iCrack, TimeStep *tStep)
{
    const std::vector<GaussPoint*> &czGaussPoints = iCrack.giveCohesiveZoneGaussPoints();

    std::vector<double> arcLengthPositions, normalJumps, tangJumps, normalTractions;

    const BasicGeometry *bg = iCrack.giveGeometry();

    for( GaussPoint *gp: czGaussPoints ) {

        StructuralInterfaceMaterialStatus *matStat = dynamic_cast<StructuralInterfaceMaterialStatus*> ( gp->giveMaterialStatus() );
        if(matStat != NULL) {

            // Compute arc length position of the Gauss point
            const FloatArray &coord = (gp->giveGlobalCoordinates());
            double tangDist = 0.0, arcPos = 0.0;
            bg->computeTangentialSignDist(tangDist, coord, arcPos);
            arcLengthPositions.push_back(arcPos);

            // Compute displacement jump in normal and tangential direction
            // Local numbering: (tang_z, tang, normal)
            const FloatArray &jumpLoc = matStat->giveJump();

            double normalJump = jumpLoc.at(3);
            normalJumps.push_back(normalJump);


            tangJumps.push_back( jumpLoc.at(2) );


            const FloatArray &trac = matStat->giveFirstPKTraction();
            normalTractions.push_back(trac.at(3));
        }
    }



    Domain *domain = emodel->giveDomain(1);
    XfemManager *xMan = domain->giveXfemManager();
    if ( xMan != NULL ) {
        double time = 0.0;

        TimeStep *ts = emodel->giveCurrentStep();
        if ( ts != NULL ) {
            time = ts->giveTargetTime();
        }

        int eiIndex = iCrack.giveNumber();

        std :: stringstream strNormalJump;
        strNormalJump << "NormalJumpGnuplotEI" << eiIndex << "Time" << time << ".dat";
        std :: string nameNormalJump = strNormalJump.str();
        XFEMDebugTools::WriteArrayToGnuplot(nameNormalJump, arcLengthPositions, normalJumps);

        std :: stringstream strTangJump;
        strTangJump << "TangJumpGnuplotEI" << eiIndex << "Time" << time << ".dat";
        std :: string nameTangJump = strTangJump.str();
        XFEMDebugTools::WriteArrayToGnuplot(nameTangJump, arcLengthPositions, tangJumps);

        std :: stringstream strNormalTrac;
        strNormalTrac << "NormalTracGnuplotEI" << eiIndex << "Time" << time << ".dat";
        std :: string nameNormalTrac = strNormalTrac.str();
        XFEMDebugTools::WriteArrayToGnuplot(nameNormalTrac, arcLengthPositions, normalTractions);


        std::vector<FloatArray> matForcesStart, matForcesEnd;
        std::vector<double> radii;

        // Material forces
        for(double matForceRadius : mMatForceRadii) {

            radii.push_back(matForceRadius);

            EnrichmentFront *efStart = iCrack.giveEnrichmentFrontStart();
            const TipInfo &tipInfoStart = efStart->giveTipInfo();

            FloatArray matForceStart;
            mpMatForceEvaluator->computeMaterialForce(matForceStart, *domain, tipInfoStart, tStep, matForceRadius);

            if(matForceStart.giveSize() > 0) {
                matForcesStart.push_back(matForceStart);
            }
            else {
                matForcesStart.push_back({0.0,0.0});
            }


            EnrichmentFront *efEnd = iCrack.giveEnrichmentFrontEnd();
            const TipInfo &tipInfoEnd = efEnd->giveTipInfo();

            FloatArray matForceEnd;
            mpMatForceEvaluator->computeMaterialForce(matForceEnd, *domain, tipInfoEnd, tStep, matForceRadius);

            if(matForceEnd.giveSize() > 0) {
                matForcesEnd.push_back(matForceEnd);
            }
            else {
                matForcesEnd.push_back({0.0,0.0});
            }

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

//.........这里部分代码省略.........
    for(size_t i = 0; i < mTriangles.size(); i++) {
    	if( mTriangles[i].getArea() > triTol ) {
    		triToKeep.push_back(i);
    	}
    }

    int nPointsTot = nPoints * triToKeep.size();
    FloatArray coords_xi1, coords_xi2, weights;
    this->giveTriCoordsAndWeights(nPoints, coords_xi1, coords_xi2, weights);
    this->gaussPointArray = new GaussPoint * [ nPointsTot ];
    ////////////////////////////////////////////


    std :: vector< FloatArray >newGPCoord;

    double parentArea = this->elem->computeArea();

    // Loop over triangles
    for ( int i = 0; i < int( triToKeep.size() ); i++ ) {
        // TODO: Probably unnecessary to allocate here
        const FloatArray **coords = new const FloatArray * [ mTriangles [ triToKeep[i] ].giveNrVertices() ];
        // this we should put into the function before
        for ( int k = 0; k < mTriangles [ triToKeep[i] ].giveNrVertices(); k++ ) {
            coords [ k ] = new FloatArray( ( mTriangles [ triToKeep[i] ].giveVertex(k + 1) ) );
        }

        // Can not be used because it writes to the start of the array instead of appending.
        //		int nPointsTri = GaussIntegrationRule :: SetUpPointsOnTriangle(nPoints, mode);

        for ( int j = 0; j < nPoints; j++ ) {
            FloatArray global;
            GaussPoint * &gp = this->gaussPointArray [ pointsPassed ];

            FloatArray *coord = new FloatArray(2);
            coord->at(1) = coords_xi1.at(j + 1);
            coord->at(2) = coords_xi2.at(j + 1);
            gp = new GaussPoint(this, pointsPassed + 1, coord, weights.at(j + 1), mode);



            mTriInterp.local2global( global, * gp->giveCoordinates(),
                                     FEIVertexListGeometryWrapper(mTriangles [ triToKeep[i] ].giveNrVertices(), coords) );

            newGPCoord.push_back(global);


            FloatArray local;
            this->elem->computeLocalCoordinates(local, global);

            gp->setCoordinates(local);




            double refElArea = this->elem->giveParentElSize();

            gp->setWeight(2.0 * refElArea * gp->giveWeight() * mTriangles [ triToKeep[i] ].getArea() / parentArea); // update integration weight


            pointsPassed++;
        }


        for ( int k = 0; k < mTriangles [ triToKeep[i] ].giveNrVertices(); k++ ) {
            delete coords [ k ];
        }

        delete [] coords;
    }

#if PATCH_INT_DEBUG > 0

    double time = 0.0;

    Element *el = this->elem;
    if(el != NULL) {
    	Domain *dom = el->giveDomain();
    	if(dom != NULL) {
    		EngngModel *em = dom->giveEngngModel();
    		if(em != NULL) {
    			TimeStep *ts = em->giveCurrentStep();
    			if(ts != NULL) {
    				time = ts->giveTargetTime();
    			}
    		}
    	}
    }
    int elIndex = this->elem->giveGlobalNumber();
    std :: stringstream str;
    str << "GaussPointsTime" << time << "El" << elIndex << ".vtk";
    std :: string name = str.str();

    XFEMDebugTools :: WritePointsToVTK(name, newGPCoord);
#endif

    numberOfIntegrationPoints = pointsPassed;


    return numberOfIntegrationPoints;
}

void main()
{
    GameState state = StateResults;
	Random randomGen;
  	randomGen.seed();

    for (unsigned i = 0; i < arraysize(cubeVideo); i++)
    {
        InitCube(i, i);
    }

    EventHandler handler;
    handler.init();

    TimeStep ts;
    float Delay;
    bool isDone;

    while (1)
    {
        switch (state)
        {
            case StateInit:
				CurrentGame->init();
				Delay = 3;
				playSfx (CountSound);
				LOG ("Init\n");
				state = StateStart;
            break;

            case StateStart:
                Delay -= float(ts.delta());

                if (Delay <= 0)
                {
                    CurrentGame->start();
                    for (unsigned i = 0; i < gNumCubes; i++)
                    {
                        CurrentGameCube[i]->start();
                    }
                    LOG ("Start Done\n");
		            state = StatePlay;
                }
            break;

            case StatePlay:
				isDone = CurrentGame->update(ts.delta());

				if (isDone)
				{
					LOG ("Game Done\n");
					state = StateResults;
					Delay = 3;
					playSfx (CheersSound);
				}
				else
				{
					for (unsigned i = 0; i < gNumCubes; i++)
					{
						CurrentGameCube[i]->update(ts.delta());
					}
				}
            break;

            case StateResults:
                Delay -= float(ts.delta());

                if (Delay <= 0)
                {
					state = StateInit;
					LOG ("Results Done\n");
					
					int rand = randomGen.randint(0, 3);
					switch (rand)
					{
					    case 0:
						case 2:
						case 3:
						LOG ("Playing Flip It\n");
						CurrentGame = &flipItGame;
					    break;
						
					    case 1:
						LOG ("Playing Shake\n");
						CurrentGame = &shakeGame;
						break;
					} 
                }
            break;

        }

        System::paint();
        ts.next();
    }
}

void
NLTransientTransportProblem :: assembleAlgorithmicPartOfRhs(FloatArray &answer,
                                                            const UnknownNumberingScheme &ns, TimeStep *tStep)
{
    //
    // Computes right hand side on all nodes
    //
    double t = tStep->giveTargetTime();
    IntArray loc;
    FloatMatrix charMtrxCond, charMtrxCap;
    FloatArray r, drdt, contrib, help;
    Element *element;
    TimeStep *previousStep = this->givePreviousStep(); //r_t
    TimeStep *currentStep = this->giveCurrentStep(); //r_{t+\Delta t}. Note that *tStep is a Tau step between r_t and r_{t+\Delta t}

    Domain *domain = this->giveDomain(1);
    int nelem = domain->giveNumberOfElements();

    for ( int i = 1; i <= nelem; i++ ) {
        element = domain->giveElement(i);
        // skip remote elements (these are used as mirrors of remote elements on other domains
        // when nonlocal constitutive models are used. They introduction is necessary to
        // allow local averaging on domains without fine grain communication between domains).
        if ( element->giveParallelMode() == Element_remote ) {
            continue;
        }

        if ( !element->isActivated(tStep) ) {
            continue;
        }

        element->giveLocationArray(loc, ns);

        element->giveCharacteristicMatrix(charMtrxCond, TangentStiffnessMatrix, tStep);
        element->giveCharacteristicMatrix(charMtrxCap, CapacityMatrix, tStep);


        /*
         *  element -> computeVectorOf (VM_Total, tStep, r);
         *  element -> computeVectorOf (VM_Velocity, tStep, drdt);
         */

        if ( ( t >= previousStep->giveTargetTime() ) && ( t <= currentStep->giveTargetTime() ) ) {
            FloatArray rp, rc;
            element->computeVectorOf(VM_Total, currentStep, rc);
            element->computeVectorOf(VM_Total, previousStep, rp);

            //approximate derivative with a difference
            drdt.beDifferenceOf(rc, rp);
            drdt.times( 1. / currentStep->giveTimeIncrement() );
            //approximate current solution from linear interpolation
            rp.times(1 - alpha);
            rc.times(alpha);
            r = rc;
            r.add(rp);
        } else {
            OOFEM_ERROR("unsupported time value");
        }


        if ( lumpedCapacityStab ) {
            int size = charMtrxCap.giveNumberOfRows();
            double s;
            for ( int j = 1; j <= size; j++ ) {
                s = 0.0;
                for ( int k = 1; k <= size; k++ ) {
                    s += charMtrxCap.at(j, k);
                    charMtrxCap.at(j, k) = 0.0;
                }

                charMtrxCap.at(j, j) = s;
            }
        }

        help.beProductOf(charMtrxCap, drdt);

        contrib.beProductOf(charMtrxCond, r);
        contrib.add(help);
        contrib.negated();

        answer.assemble(contrib, loc);
    }
}

void main()
{
	//Create our audio channels
	AudioChannel a1(0);
	AudioChannel a2(1);
	AudioChannel a3(2);
	AudioChannel a4(3);
	AudioChannel a5(4);
	AudioChannel a6(5);
	channels[0] = &a1;
	channels[1] = &a2;
	channels[2] = &a3;
	channels[3] = &a4;
	channels[4] = &a5;
	channels[5] = &a6;
	r.seed();
	
	channels[BALL_ROLL_CHANNEL]->play(sRollLoop, REPEAT);	//Start playing rolling marble noise
	channels[BALL_ROLL_CHANNEL]->setVolume(0);
	for(int i = 0; i < NUM_CUBES; i++)
		g_starsCollected[i] = 0;
	for(int i = 0; i < NUM_COLORS; i++)
		g_bColorsUsed[i] = false;
	g_iBoardReset = -1;
	g_iScore = -1;
    TimeStep ts;
	float fTapPromptDelay = 0.0;
	int iBoardDied;

	//Initialize our boards
	for(int i = 0; i < NUM_CUBES; i++)
		boards[i].init(i);
	Events::neighborAdd.set(onNeighborAdd);	//Function for when two cubes touch each other
	Events::cubeTouch.set(onTouch);			//Function for when a cube is tapped
	
	//Add the marble to one of them
	Float2 fVel;
	fVel.set(0,0);
	Float2 fPos = LCD_center;
	fPos.x += TILE_WIDTH/2.0;
	fPos.y += TILE_HEIGHT/2.0;
	boards[0].addMarble(fPos, fVel);
	
	TextDraw td;
	bool bFirstSound = true;
	
	//Main loop
    while (1) 
	{
        ts.next();
		for(int i = 0; i < NUM_CUBES; i++)
		{
			int iMode;
			switch(g_iCurMode)
			{
				case BOARD_NOTHING:
					iMode = boards[i].update(float(ts.delta()));
					//Update our rolling sound to the right volume
					if(boards[i].hasMarble())
					{
						float fVol = boards[i].getMarbleVelocity() * 0.9;
						if(fVol > MAX_VOLUME)
							fVol = MAX_VOLUME;
						channels[BALL_ROLL_CHANNEL]->setVolume(fVol);
					}
					if(iMode & BOARD_GOTPOINT)
					{
						iMode ^= BOARD_GOTPOINT;
						g_iScore++;
						if(++g_starsCollected[i] == NUM_STARS_CUBE)
							g_iBoardReset = i;
						//Play sound for getting a star, but not right on reset
						if(bFirstSound)
							bFirstSound = false;
						else
							channels[BALL_SFX_CHANNEL]->play(sGetStar);
					}
					if(iMode & BOARD_DIED)
					{
						//Show game over screen
						iMode ^= BOARD_DIED;
						iMode |= MODE_GAMEOVER;
						td.draw(boards[i].getVid(), "Game over", 6);
						String<64> s;
						s << "Score: " << g_iScore;
						td.draw(boards[i].getVid(), s.c_str(), 8);
						fTapPromptDelay = 0.0;
						iBoardDied = i;
						bFirstSound = true;
						channels[BALL_ROLL_CHANNEL]->setVolume(0);
						channels[BALL_SFX_CHANNEL]->play(sDie);
					}
					if(iMode & BOARD_LEFT)
					{
						iMode ^= BOARD_LEFT;
						if(g_iBoardReset == i)
						{
							boards[i].reset(g_bColorsUsed);
							g_starsCollected[i] = 0;
							g_iBoardReset = -1;
//.........这里部分代码省略.........

void StructuralMaterialEvaluator :: solveYourself()
{
    Domain *d = this->giveDomain(1);

    MaterialMode mode = _3dMat;
    FloatArray initialStrain(6);
    gps.clear();
    gps.reserve(d->giveNumberOfMaterialModels());
    for ( int i = 1; i <= d->giveNumberOfMaterialModels(); i++ ) {
        std :: unique_ptr< GaussPoint > gp = std::make_unique<GaussPoint>(nullptr, i, FloatArray(0), 1, mode);
        gps.emplace_back( std :: move(gp) );
        // Initialize the strain vector;
        StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( d->giveMaterial(i)->giveStatus( gps[i-1].get() ) );
        status->letStrainVectorBe(initialStrain);
    }

    std :: string outname = this->giveOutputBaseFileName() + ".matdata";
    this->outfile.open( outname.c_str() );

    this->timer.startTimer(EngngModelTimer :: EMTT_AnalysisTimer);

    TimeStep *tStep = giveNextStep();

    // Note, strain == strain-rate (kept as strain for brevity)
    int maxiter = 100; // User input?
    FloatArray stressC, deltaStrain, strain, stress, res;
    stressC.resize( sControl.giveSize() );
    res.resize( sControl.giveSize() );

    FloatMatrix tangent, reducedTangent;
    for ( int istep = 1; istep <= this->numberOfSteps; ++istep ) {
        this->timer.startTimer(EngngModelTimer :: EMTT_SolutionStepTimer);
        for ( int imat = 1; imat <= d->giveNumberOfMaterialModels(); ++imat ) {
            GaussPoint *gp = gps[imat-1].get();
            StructuralMaterial *mat = static_cast< StructuralMaterial * >( d->giveMaterial(imat) );
            StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( mat->giveStatus(gp) );

            strain = status->giveStrainVector();
            // Update the controlled parts
            for ( int j = 1; j <= eControl.giveSize(); ++j ) {
                int p = eControl.at(j);
                strain.at(p) = d->giveFunction( cmpntFunctions.at(p) )->evaluateAtTime( tStep->giveIntrinsicTime() );
            }

            for ( int j = 1; j <= sControl.giveSize(); ++j ) {
                int p = sControl.at(j);
                stressC.at(j) = d->giveFunction( cmpntFunctions.at(p) )->evaluateAtTime( tStep->giveIntrinsicTime() );
            }

            //strain.add(-100, {6.27e-06,  6.27e-06, 6.27e-06, 0, 0, 0});
            for ( int iter = 1; iter < maxiter; iter++ ) {
#if 0
                // Debugging:
                mat->give3dMaterialStiffnessMatrix(tangent, TangentStiffness, gp, tStep);
                tangent.printYourself("# tangent");
                
                strain.zero();
                mat->giveRealStressVector_3d(stress, gp, strain, tStep);
                FloatArray strain2;
                tangent.solveForRhs(stress, strain2);
                strain2.printYourself("# thermal expansion");
                break;
#endif

                strain.printYourself("Macro strain guess");
                mat->giveRealStressVector_3d(stress, gp, strain, tStep);
                for ( int j = 1; j <= sControl.giveSize(); ++j ) {
                    res.at(j) = stressC.at(j) - stress.at( sControl.at(j) );
                }

                OOFEM_LOG_INFO("*** Time step: %d (t = %.2e), Material %d, Iteration: %d,  Residual = %e (tolerance %.2e)\n", 
                              istep, tStep->giveIntrinsicTime(), imat, iter, res.computeNorm(), tolerance);

                if ( res.computeNorm() <= tolerance ) {
                    break;
                } else {
                    if ( tangent.giveNumberOfRows() == 0 || !keepTangent ) {
                        mat->give3dMaterialStiffnessMatrix(tangent, TangentStiffness, gp, tStep);
                    }

                    // Pick out the stress-controlled part;
                    reducedTangent.beSubMatrixOf(tangent, sControl, sControl);

                    // Update stress-controlled part of the strain
                    reducedTangent.solveForRhs(res, deltaStrain);
                    //deltaStrain.printYourself("deltaStrain");
                    for ( int j = 1; j <= sControl.giveSize(); ++j ) {
                        strain.at( sControl.at(j) ) += deltaStrain.at(j);
                    }
                }
            }

            if ( res.computeNorm() > tolerance ) {
                OOFEM_WARNING("Residual did not converge!");
            }

            // This material model has converged, so we update it and go on to the next.
            gp->updateYourself(tStep);
        }

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

//.........这里部分代码省略.........
#ifdef __PARALLEL_MODE
        if ( this->giveDomain(1)->giveElement(ielem)->giveParallelMode() == Element_remote ) {
            continue;
        }

#endif

        result &= this->giveDomain(1)->giveElement(ielem)->adaptiveUpdate( this->giveCurrentStep() );
    }

    // finish mapping process
    for ( ielem = 1; ielem <= nelem; ielem++ ) {
        /* HUHU CHEATING */
#ifdef __PARALLEL_MODE
        if ( this->giveDomain(1)->giveElement(ielem)->giveParallelMode() == Element_remote ) {
            continue;
        }

#endif

        result &= this->giveDomain(1)->giveElement(ielem)->adaptiveFinish( this->giveCurrentStep() );
    }

    nMethod->reinitialize();


    // increment time step if mapped state will be considered as new solution stepL
    // this->giveNextStep();
    if ( equilibrateMappedConfigurationFlag ) {
        // we need to  assemble the load vector in same time as the restarted step,
        // so new time step is generated with same intrincic time as has the
        // previous step if equilibrateMappedConfigurationFlag is set.
        // this allows to equlibrate the previously reached state
        TimeStep *cts = this->giveCurrentStep();
        // increment version of solution step
        cts->incrementVersion();

        //cts->setTime(cts->giveTime()-cts->giveTimeIncrement());
        //cts = this->givePreviousStep();
        //cts->setTime(cts->giveTime()-cts->giveTimeIncrement());
    }

    if ( this->giveCurrentStep()->giveNumber() == this->giveCurrentMetaStep()->giveFirstStepNumber() ) {
        this->updateAttributes( this->giveCurrentMetaStep() );
    }

    // increment solution state counter - not needed, IPs are updated by adaptiveUpdate previously called
    // and there is no change in primary vars.
    // this->giveCurrentStep()->incrementStateCounter();

    // assemble new initial load for new discretization
    this->assembleInitialLoadVector( initialLoadVector, initialLoadVectorOfPrescribed,
                                    this, 1, this->giveCurrentStep() );
    this->assembleIncrementalReferenceLoadVectors( incrementalLoadVector, incrementalLoadVectorOfPrescribed,
                                                  refLoadInputMode, this->giveDomain(1), EID_MomentumBalance,
                                                  this->giveCurrentStep() );

    // assemble new total load for new discretization
    // this->assembleCurrentTotalLoadVector (totalLoadVector, totalLoadVectorOfPrescribed, this->giveCurrentStep());
    // set bcloadVector to zero (no increment within same step)

    timer.stopTimer();
    mc3 = timer.getUtime();

    // compute processor time used by the program
    OOFEM_LOG_INFO("user time consumed by primary mapping: %.2fs\n", mc1);