C++ EngngModel::giveCurrentStep方法代码示例

void
StaticFracture :: optimize(TimeStep *tStep)
{
    // Main optimization loop

    MinCompliance *objFunc = dynamic_cast< MinCompliance * >( this->objFuncList[0] ); 

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

        double cost = 0.0; // should lie in obj fnc
        double dce = 0.0;
        for (int i = 1; i <= d->giveNumberOfElements(); i++) {  
            Element *el = d->giveElement(i);
            cost += objFunc->evaluateYourself(el, dce, sp->giveCurrentStep() ); // add cost for each element
        }


        // Filter sensitivities
        this->filterSensitivities(objFunc);

        // Update design variables based on some method. For now use the 'standard' optimality criteria
        this->optimalityCriteria(objFunc);
     

        // Update material parameters
        for (int i = 1; i <= d->giveNumberOfMaterialModels(); i++) {
            DynamicInputRecord ir;
            Material *mat = d->giveMaterial(i);
            mat->giveInputRecord(ir);
            double E0 = 1.0;
            double fac = pow( objFunc->designVarList.at(i), objFunc->penalty);
            ir.setField( E0 * fac, _IFT_IsotropicLinearElasticMaterial_e);
            mat->initializeFrom(&ir);
        }

        printf("\n costfunction %e & sum sensitivity %e & sum x %e \n", cost, objFunc->sensitivityList.sum(), 
            objFunc->designVarList.sum() );

    }
}

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
                }
        }
    }
    

}

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

#if defined ( __PETSC_MODULE ) || defined ( __SLEPC_MODULE )
    int modulesArgc = modulesArgs.size();
    char **modulesArgv = const_cast< char ** >(& modulesArgs [ 0 ]);
#endif

#ifdef __PETSC_MODULE
    PetscInitialize(& modulesArgc, & modulesArgv, PETSC_NULL, PETSC_NULL);
#endif

#ifdef __SLEPC_MODULE
    SlepcInitialize(& modulesArgc, & modulesArgv, PETSC_NULL, PETSC_NULL);
#endif

#ifdef __PYTHON_MODULE
    Py_Initialize();
    // Adding . to the system path allows us to run Python functions stored in the working directory.
    PyRun_SimpleString("import sys");
    PyRun_SimpleString("sys.path.append(\".\")");
#endif

#ifdef __PARALLEL_MODE
    if ( parallelFlag ) {
        inputFileName << "." << rank;
        outputFileName << "." << rank;
        errOutputFileName << "." << rank;
    }
#endif
    if ( outputFileFlag ) {
        oofem_logger.appendLogTo( outputFileName.str() );
    }
    if ( errOutputFileFlag ) {
        oofem_logger.appendErrorTo( errOutputFileName.str() );
    }

    // print header to redirected output
    OOFEM_LOG_FORCED(PRG_HEADER_SM);

    OOFEMTXTDataReader dr( inputFileName.str() );
    EngngModel *problem = :: InstanciateProblem(dr, _processor, contextFlag, NULL, parallelFlag);
    dr.finish();
    if ( !problem ) {
        OOFEM_LOG_ERROR("Couldn't instanciate problem, exiting");
        exit(EXIT_FAILURE);
    }

    problem->checkProblemConsistency();
    problem->init();

    if ( renumberFlag ) {
        problem->setRenumberFlag();
    }

    if ( restartFlag ) {
        try {
            FileDataStream stream(problem->giveContextFileName(restartStep, 0), false);
            problem->restoreContext(stream, CM_State | CM_Definition);
        } catch ( const FileDataStream::CantOpen & e ) {
            printf("%s", e.what());
            exit(1);
        } catch ( ContextIOERR & c ) {
            c.print();
            exit(1);
        }
        problem->initStepIncrements();
    } else if ( adaptiveRestartFlag ) {
        problem->initializeAdaptive(adaptiveRestartFlag);
        problem->saveStepContext(problem->giveCurrentStep(),CM_State);
        // exit (1);
    }

    if ( debugFlag ) {
        oofem_debug(problem);
    }


    try {
        problem->solveYourself();
    } catch(OOFEM_Terminate & c) {
        delete problem;

        oofem_finalize_modules();

        return 1;
    }

    problem->terminateAnalysis();
#ifdef __PARALLEL_MODE
    if ( parallelFlag ) {
        DynamicCommunicationBuffer :: printInfo();
    }
#endif
    oofem_logger.printStatistics();
    delete problem;

    oofem_finalize_modules();

    return 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;
}