C++ EngngModel类代码示例

// Uses an input file for now, should eventually create the RVE itself.
bool FE2FluidMaterialStatus :: createRVE(int n, GaussPoint *gp, const std::string &inputfile)
{
    OOFEMTXTDataReader dr(inputfile.c_str());
    EngngModel *em = InstanciateProblem(&dr, _processor, 0); // Everything but nrsolver is updated.
    dr.finish();
    em->setProblemScale(microScale);
    em->checkProblemConsistency();
    em->initMetaStepAttributes( em->giveMetaStep( 1 ) );
    em->giveNextStep(); // Makes sure there is a timestep (which we will modify before solving a step)
    em->init();

    this->rve = dynamic_cast<StokesFlow*> (em);
    if (!this->rve) {
        return false;
    }
    std::ostringstream name;
    name << this->rve->giveOutputBaseFileName() << "-gp" << n;
#ifdef __PARALLEL_MODE
    if (this->domain->giveEngngModel()->isParallel() && this->domain->giveEngngModel()->giveNumberOfProcesses() > 1) {
        name << "." << this->domain->giveEngngModel()->giveRank();
    }
#endif
    this->rve->letOutputBaseFileNameBe(name.str());

    this->bc = dynamic_cast< MixedGradientPressureBC* >(this->rve->giveDomain(1)->giveBc(1));
    if (!this->bc) {
        OOFEM_ERROR("FE2FluidMaterialStatus :: createRVE - RVE doesn't have necessary boundary condition; should have MixedGradientPressure as first b.c. (in first domain)");
    }

    return true;
}

bool
StructuralFE2MaterialStatus :: createRVE(int n, GaussPoint *gp, const std :: string &inputfile)
{
    OOFEMTXTDataReader dr( inputfile.c_str() );
    EngngModel *em = InstanciateProblem(dr, _processor, 0); // Everything but nrsolver is updated.
    dr.finish();
    em->setProblemScale(microScale);
    em->checkProblemConsistency();
    em->initMetaStepAttributes( em->giveMetaStep(1) );
    em->giveNextStep(); // Makes sure there is a timestep (which we will modify before solving a step)
    em->init();

    this->rve.reset( em );

    std :: ostringstream name;
    name << this->rve->giveOutputBaseFileName() << "-gp" << n;
    if ( this->domain->giveEngngModel()->isParallel() && this->domain->giveEngngModel()->giveNumberOfProcesses() > 1 ) {
        name << "." << this->domain->giveEngngModel()->giveRank();
    }

    this->rve->letOutputBaseFileNameBe( name.str() );

    this->bc = dynamic_cast< PrescribedGradientHomogenization * >( this->rve->giveDomain(1)->giveBc(1) );
    if ( !this->bc ) {
        OOFEM_ERROR("RVE doesn't have necessary boundary condition; should have a type of PrescribedGradientHomogenization as first b.c.");
    }

    return true;
}

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
NodalAveragingRecoveryModel :: exchangeDofManValues(int ireg, FloatArray &lhs, IntArray &regionDofMansConnectivity,
                                                    IntArray &regionNodalNumbers, int regionValSize)
{
    EngngModel *emodel = domain->giveEngngModel();
    ProblemCommunicatorMode commMode = emodel->giveProblemCommMode();

    if ( commMode == ProblemCommMode__NODE_CUT ) {
        parallelStruct ls( &lhs, &regionDofMansConnectivity, &regionNodalNumbers, regionValSize);

        // exchange data for shared nodes
        communicator->packAllData(this, & ls, & NodalAveragingRecoveryModel :: packSharedDofManData);
        communicator->initExchange(789 + ireg);
        communicator->unpackAllData(this, & ls, & NodalAveragingRecoveryModel :: unpackSharedDofManData);
        communicator->finishExchange();
    } else {
        OOFEM_ERROR("NodalAveragingRecoveryModel :: exchangeDofManValues: Unsupported commMode");
    }
}

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

}

int ActiveDof :: askNewEquationNumber(TimeStep *tStep)
{
    if ( !this->isPrimaryDof() ) {
        return 0;
    }

    EngngModel *model = dofManager->giveDomain()->giveEngngModel();

    if ( dofManager->giveParallelMode() == DofManager_null ) {
        equationNumber = 0;
        return 0;
    }

    if ( this->hasBc(tStep) ) {
        equationNumber = -model->giveNewPrescribedEquationNumber(dofManager->giveDomain()->giveNumber(), this->dofID);
    } else {
        equationNumber = model->giveNewEquationNumber(dofManager->giveDomain()->giveNumber(), this->dofID);
    }

    return equationNumber;
}

int MasterDof :: askNewEquationNumber(TimeStep *tStep)
// Returns the newly obtained number of the equation in the governing system
// of equations that corres-
// ponds to the receiver. The equation number is 0 if the receiver is
// subjected to a boundary condition, else it is n+1, where n is the
// equation number of the most recently numbered degree of freedom.
{
    EngngModel *model = dofManager->giveDomain()->giveEngngModel();

    if ( dofManager->giveParallelMode() == DofManager_null ) {
        equationNumber = 0;
        return 0;
    }

    if ( this->hasBc(tStep) ) {
        equationNumber = -1 * model->giveNewPrescribedEquationNumber(dofManager->giveDomain()->giveNumber(), this->dofID);
    } else {
        equationNumber = model->giveNewEquationNumber(dofManager->giveDomain()->giveNumber(), this->dofID);
    }

    return equationNumber;
}

void
NodalAveragingRecoveryModel :: initCommMaps()
{
 #ifdef __PARALLEL_MODE
    if ( initCommMap ) {
        EngngModel *emodel = domain->giveEngngModel();
        ProblemCommunicatorMode commMode = emodel->giveProblemCommMode();
        if ( commMode == ProblemCommMode__NODE_CUT ) {
            commBuff = new CommunicatorBuff(emodel->giveNumberOfProcesses(), CBT_dynamic);
            communicator = new ProblemCommunicator(emodel, commBuff, emodel->giveRank(),
                                                   emodel->giveNumberOfProcesses(),
                                                   commMode);
            communicator->setUpCommunicationMaps(domain->giveEngngModel(), true, true);
            OOFEM_LOG_INFO("NodalAveragingRecoveryModel :: initCommMaps: initialized comm maps\n");
            initCommMap = false;
        } else {
            OOFEM_ERROR("NodalAveragingRecoveryModel :: initCommMaps: unsupported comm mode");
        }
    }

 #endif
}

EngngModel *InstanciateProblem(DataReader *dr, problemMode mode, int contextFlag, EngngModel *_master, bool parallelFlag)
{
    const char *__proc = "InstanciateProblem"; // Required by IR_GIVE_FIELD macro
    IRResultType result;                       // Required by IR_GIVE_FIELD macro
    EngngModel *problem;
    std::string problemName, dataOutputFileName, desc;

    dataOutputFileName = dr->giveOutputFileName();
    desc = dr->giveDescription();

    /* here we need copy of input record. The pointer returned by dr->giveInputRecord can (and will)
     * be updated as reading e-model components (nodes, etc). But we need this record being available
     * through the whole e-model instanciation
     */
    InputRecord *emodelir = dr->giveInputRecord(DataReader :: IR_emodelRec, 1)->GiveCopy();
    result = emodelir->giveRecordKeywordField(problemName); ///@todo Make this function robust, it can't be allowed to fail (the record keyword is not a normal field-id)
    if ( result != IRRT_OK ) {
        IR_IOERR("", __proc, "", emodelir, result);
    }

    problem = classFactory.createEngngModel(problemName.c_str(), 1, _master);
    if (!problem) {
        OOFEM_ERROR2("EngngModel::InstanciateProblem - Failed to construct engineering model if type \"%s\".\n", problemName.c_str());
    }
    problem->setProblemMode(mode);
    problem->setParallelMode(parallelFlag);

    if ( contextFlag ) {
        problem->setContextOutputMode(COM_Always);
    }

    problem->instanciateYourself(dr, emodelir, dataOutputFileName.c_str(), desc.c_str());

    delete emodelir;

    return problem;
}

void LSPrimaryVariableMapper :: mapPrimaryVariables(FloatArray &oU, Domain &iOldDom, Domain &iNewDom, ValueModeType iMode, TimeStep &iTStep)
{
    EngngModel *engngMod = iNewDom.giveEngngModel();
    EModelDefaultEquationNumbering num;


    const int dim = iNewDom.giveNumberOfSpatialDimensions();

    int numElNew = iNewDom.giveNumberOfElements();

    // Count dofs
    int numDofsNew = engngMod->giveNumberOfDomainEquations( 1, num );


    oU.resize(numDofsNew);
    oU.zero();

    FloatArray du(numDofsNew);
    du.zero();

    FloatArray res(numDofsNew);

#ifdef __PETSC_MODULE
    PetscSparseMtrx *K = dynamic_cast<PetscSparseMtrx*>( classFactory.createSparseMtrx(SMT_PetscMtrx) );
    SparseLinearSystemNM *solver = classFactory.createSparseLinSolver(ST_Petsc, & iOldDom, engngMod);
#else
    SparseMtrx *K = classFactory.createSparseMtrx(SMT_Skyline);
    SparseLinearSystemNM *solver = classFactory.createSparseLinSolver(ST_Direct, & iOldDom, engngMod);
#endif


    K->buildInternalStructure( engngMod, 1, num );

    int maxIter = 1;

    for ( int iter = 0; iter < maxIter; iter++ ) {
        K->zero();
        res.zero();


        // Contribution from elements
        for ( int elIndex = 1; elIndex <= numElNew; elIndex++ ) {
            StructuralElement *elNew = dynamic_cast< StructuralElement * >( iNewDom.giveElement(elIndex) );
            if ( elNew == NULL ) {
                OOFEM_ERROR("Failed to cast Element new to StructuralElement.");
            }

            ///////////////////////////////////
            // Compute residual

            // Count element dofs
            int numElNodes = elNew->giveNumberOfDofManagers();
            int numElDofs = 0;
            for ( int i = 1; i <= numElNodes; i++ ) {
                numElDofs += elNew->giveDofManager(i)->giveNumberOfDofs();
            }

            FloatArray elRes(numElDofs);
            elRes.zero();

            IntArray elDofsGlob;
            elNew->giveLocationArray( elDofsGlob, num );


            // Loop over Gauss points
            for ( int intRuleInd = 0; intRuleInd < elNew->giveNumberOfIntegrationRules(); intRuleInd++ ) {
                IntegrationRule *iRule = elNew->giveIntegrationRule(intRuleInd);

                for ( GaussPoint *gp: *iRule ) {

                    // New N-matrix
                    FloatMatrix NNew;
                    elNew->computeNmatrixAt(* ( gp->giveNaturalCoordinates() ), NNew);


                    //////////////
                    // Global coordinates of GP
                    const int nDofMan = elNew->giveNumberOfDofManagers();

                    FloatArray Nc;
                    FEInterpolation *interp = elNew->giveInterpolation();
                    const FloatArray &localCoord = * ( gp->giveNaturalCoordinates() );
                    interp->evalN( Nc, localCoord, FEIElementGeometryWrapper(elNew) );

                    const IntArray &elNodes = elNew->giveDofManArray();

                    FloatArray globalCoord(dim);
                    globalCoord.zero();

                    for ( int i = 1; i <= nDofMan; i++ ) {
                        DofManager *dMan = elNew->giveDofManager(i);

                        for ( int j = 1; j <= dim; j++ ) {
                            globalCoord.at(j) += Nc.at(i) * dMan->giveCoordinate(j);
                        }
                    }
                    //////////////


                    // Localize element and point in the old domain
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    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
SolutionbasedShapeFunction :: initializeSurfaceData(modeStruct *mode)
{
    EngngModel *m = mode->myEngngModel;
    double TOL2 = 1e-5;
    IntArray pNodes, mNodes, zNodes;

    Set *mySet = this->domain->giveSet( this->giveSetNumber() );
    IntArray BoundaryList = mySet->giveBoundaryList();

    // First add all nodes to pNodes or nNodes respectively depending on coordinate and normal.
    for ( int i = 0; i < BoundaryList.giveSize() / 2; i++ ) {
        int ElementID = BoundaryList(2 * i);
        int Boundary = BoundaryList(2 * i + 1);

        Element *e = m->giveDomain(1)->giveElement(ElementID);
        FEInterpolation *geoInterpolation = e->giveInterpolation();

        // Check all sides of element
        IntArray bnodes;

#define usePoints 1
#if usePoints == 1
        // Check if all nodes are on the boundary
        geoInterpolation->boundaryGiveNodes(bnodes, Boundary);
        for ( int k = 1; k <= bnodes.giveSize(); k++ ) {
            DofManager *dman = e->giveDofManager( bnodes.at(k) );
            for ( int l = 1; l <= dman->giveCoordinates()->giveSize(); l++ ) {
                if ( fabs( dman->giveCoordinates()->at(l) - maxCoord.at(l) ) < TOL2 ) {
                    pNodes.insertOnce( dman->giveNumber() );
                }
                if ( fabs( dman->giveCoordinates()->at(l) - minCoord.at(l) ) < TOL2 ) {
                    mNodes.insertOnce( dman->giveNumber() );
                }
            }
        }
#else
        // Check normal
        FloatArray lcoords;
        lcoords.resize(2);
        lcoords.at(1) = 0.33333;
        lcoords.at(2) = 0.33333;

        FloatArray normal;
        geoInterpolation->boundaryEvalNormal( normal, j, lcoords, FEIElementGeometryWrapper(e) );
        geoInterpolation->boundaryGiveNodes(bnodes, j);

        printf( "i=%u\tj=%u\t(%f\t%f\t%f)\n", i, j, normal.at(1), normal.at(2), normal.at(3) );
        for ( int k = 1; k <= normal.giveSize(); k++ ) {
            if ( fabs( ( fabs( normal.at(k) ) - 1 ) ) < 1e-4 ) { // Points in x, y or z direction
                addTo = NULL;
                if ( normal.at(k) > 0.5 ) {
                    addTo = & pNodes;
                }
                if ( normal.at(k) < -0.5 ) {
                    addTo = & mNodes;
                }
                if ( addTo != NULL ) {
                    for ( int l = 1; l <= bnodes.giveSize(); l++ ) {
                        bool isSurface = false;
                        DofManager *dman = e->giveDofManager( bnodes.at(l) );
                        dman->giveCoordinates()->printYourself();
                        for ( int m = 1; m <= dman->giveCoordinates()->giveSize(); m++ ) {
                            if ( ( fabs( dman->giveCoordinates()->at(m) - maxCoord.at(m) ) < TOL2 ) || ( fabs( dman->giveCoordinates()->at(m) - minCoord.at(m) ) < TOL2 ) ) {
                                isSurface = true;
                            }
                        }

                        if ( isSurface ) {
                            addTo->insertOnce( e->giveDofManagerNumber( bnodes.at(l) ) );
                        }
                    }
                }
            }
        }
#endif
    }

#if 0
    printf("p=[");
    for ( int i = 1; i < pNodes.giveSize(); i++ ) {
        printf( "%u, ", pNodes.at(i) );
    }
    printf("];\n");
    printf("m=[");
    for ( int i = 1; i < mNodes.giveSize(); i++ ) {
        printf( "%u, ", mNodes.at(i) );
    }
    printf("];\n");
#endif
    //The intersection of pNodes and mNodes constitutes zNodes
    {
        int i = 1, j = 1;
        while ( i <= pNodes.giveSize() ) {
            j = 1;
            while ( j <= mNodes.giveSize() && ( i <= pNodes.giveSize() ) ) {
                //printf("%u == %u?\n", pNodes.at(i), mNodes.at(j));
                if ( pNodes.at(i) == mNodes.at(j) ) {
                    zNodes.insertOnce( pNodes.at(i) );
                    pNodes.erase(i);
//.........这里部分代码省略.........

void TransportGradientPeriodic :: computeTangent(FloatMatrix &k, TimeStep *tStep)
{
    EModelDefaultEquationNumbering fnum;
    //DofIDEquationNumbering pnum(true, this->grad_ids);
    EModelDefaultPrescribedEquationNumbering pnum;
    int nsd = this->domain->giveNumberOfSpatialDimensions();

    EngngModel *rve = this->giveDomain()->giveEngngModel();
    ///@todo Get this from engineering model
    std :: unique_ptr< SparseLinearSystemNM > solver( classFactory.createSparseLinSolver( ST_Petsc, this->domain, this->domain->giveEngngModel() ) ); // = rve->giveLinearSolver();
    SparseMtrxType stype = solver->giveRecommendedMatrix(true);
    std :: unique_ptr< SparseMtrx > Kff( classFactory.createSparseMtrx( stype ) );
    std :: unique_ptr< SparseMtrx > Kfp( classFactory.createSparseMtrx( stype ) );
    std :: unique_ptr< SparseMtrx > Kpp( classFactory.createSparseMtrx( stype ) );

    Kff->buildInternalStructure(rve, this->domain->giveNumber(), fnum);
    int neq = Kff->giveNumberOfRows();
    Kfp->buildInternalStructure(rve, this->domain->giveNumber(), fnum, pnum);
    Kpp->buildInternalStructure(rve, this->domain->giveNumber(), pnum);
    //Kfp->buildInternalStructure(rve, neq, nsd, {}, {});
    //Kpp->buildInternalStructure(rve, nsd, nsd, {}, {});
#if 0
    rve->assemble(*Kff, tStep, TangentAssembler(TangentStiffness), fnum, this->domain);
    rve->assemble(*Kfp, tStep, TangentAssembler(TangentStiffness), fnum, pnum, this->domain);
    rve->assemble(*Kpp, tStep, TangentAssembler(TangentStiffness), pnum, this->domain);
#else
    auto ma = TangentAssembler(TangentStiffness);
    IntArray floc, ploc;
    FloatMatrix mat, R;

    int nelem = domain->giveNumberOfElements();
#ifdef _OPENMP
 #pragma omp parallel for shared(Kff, Kfp, Kpp) private(mat, R, floc, ploc)
#endif
    for ( int ielem = 1; ielem <= nelem; ielem++ ) {
        Element *element = domain->giveElement(ielem);
        // 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 || !element->isActivated(tStep) ) {
            continue;
        }

        ma.matrixFromElement(mat, *element, tStep);

        if ( mat.isNotEmpty() ) {
            ma.locationFromElement(floc, *element, fnum);
            ma.locationFromElement(ploc, *element, pnum);
            ///@todo This rotation matrix is not flexible enough.. it can only work with full size matrices and doesn't allow for flexibility in the matrixassembler.
            if ( element->giveRotationMatrix(R) ) {
                mat.rotatedWith(R);
            }

#ifdef _OPENMP
 #pragma omp critical
#endif
            {
                Kff->assemble(floc, mat);
                Kfp->assemble(floc, ploc, mat);
                Kpp->assemble(ploc, mat);
            }
        }
    }
    Kff->assembleBegin();
    Kfp->assembleBegin();
    Kpp->assembleBegin();

    Kff->assembleEnd();
    Kfp->assembleEnd();
    Kpp->assembleEnd();
#endif

    FloatMatrix grad_pert(nsd, nsd), rhs, sol(neq, nsd);
    grad_pert.resize(nsd, nsd);
    grad_pert.beUnitMatrix();
    // Workaround since the matrix size is inflexible with custom dof numbering (so far, planned to be fixed).
    IntArray grad_loc;
    this->grad->giveLocationArray(this->grad_ids, grad_loc, pnum);
    FloatMatrix pert(Kpp->giveNumberOfRows(), nsd);
    pert.assemble(grad_pert, grad_loc, {1,2,3});
    //pert.printYourself("pert");

    //printf("Kfp = %d x %d\n", Kfp->giveNumberOfRows(), Kfp->giveNumberOfColumns());
    //printf("Kff = %d x %d\n", Kff->giveNumberOfRows(), Kff->giveNumberOfColumns());
    //printf("Kpp = %d x %d\n", Kpp->giveNumberOfRows(), Kpp->giveNumberOfColumns());

    // Compute the solution to each of the pertubation of eps
    Kfp->times(pert, rhs);
    //rhs.printYourself("rhs");

    // Initial guess (Taylor assumption) helps KSP-iterations
    for ( auto &n : domain->giveDofManagers() ) {
        int k1 = n->giveDofWithID( this->dofs(0) )->__giveEquationNumber();
        if ( k1 ) {
            FloatArray *coords = n->giveCoordinates();
            for ( int i = 1; i <= nsd; ++i ) {
                sol.at(k1, i) = -(coords->at(i) - mCenterCoord.at(i));
            }
        }
    }
//.........这里部分代码省略.........

void
SolutionbasedShapeFunction :: loadProblem()
{
    for ( int i = 0; i < this->domain->giveNumberOfSpatialDimensions(); i++ ) {
        OOFEM_LOG_INFO("************************** Instanciating microproblem from file %s for dimension %u\n", filename.c_str(), i);

        // Set up and solve problem
        OOFEMTXTDataReader drMicro( filename.c_str() );
        EngngModel *myEngngModel = InstanciateProblem(& drMicro, _processor, 0, NULL, false);
        drMicro.finish();
        myEngngModel->checkProblemConsistency();
        myEngngModel->initMetaStepAttributes( myEngngModel->giveMetaStep(1) );
        thisTimestep = myEngngModel->giveNextStep();
        myEngngModel->init();
        this->setLoads(myEngngModel, i + 1);

        // Check
        for ( int j = 1; j <= myEngngModel->giveDomain(1)->giveNumberOfElements(); j++ ) {
            Element *e = myEngngModel->giveDomain(1)->giveElement(j);
            FloatArray centerCoord;
            int vlockCount = 0;
            centerCoord.resize(3);
            centerCoord.zero();

            for ( int k = 1; k <= e->giveNumberOfDofManagers(); k++ ) {
                DofManager *dman = e->giveDofManager(k);
                centerCoord.add( * dman->giveCoordinates() );
                for ( Dof *dof: *dman ) {
                    if ( dof->giveBcId() != 0 ) {
                        vlockCount++;
                    }
                }
            }
            if ( vlockCount == 30 ) {
                OOFEM_WARNING("Element over-constrained (%u)! Center coordinate: %f, %f, %f\n", e->giveNumber(), centerCoord.at(1) / 10, centerCoord.at(2) / 10, centerCoord.at(3) / 10);
            }
        }

        myEngngModel->solveYourselfAt(thisTimestep);
        isLoaded = true;

        // Set correct export filename
        std :: string originalFilename;
        originalFilename = myEngngModel->giveOutputBaseFileName();
        if ( i == 0 ) {
            originalFilename = originalFilename + "_X";
        }
        if ( i == 1 ) {
            originalFilename = originalFilename + "_Y";
        }
        if ( i == 2 ) {
            originalFilename = originalFilename + "_Z";
        }
        myEngngModel->letOutputBaseFileNameBe(originalFilename + "_1_Base");
        myEngngModel->doStepOutput(thisTimestep);

        modeStruct *mode = new(modeStruct);
        mode->myEngngModel = myEngngModel;

        // Check elements

        // Set unknowns to the mean value of opposite sides of the domain.
        // Loop thru all nodes and compute phi for all degrees of freedom on the boundary. Save phi in a list for later use.

        initializeSurfaceData(mode);
        // Update with factor
        double am = 1.0, ap = 1.0;
        computeCorrectionFactors(* mode, & dofs, & am, & ap);

        OOFEM_LOG_INFO("Correction factors: am=%f, ap=%f\n", am, ap);

        mode->ap = ap;
        mode->am = am;

        updateModelWithFactors(mode);

        //        myEngngModel->letOutputBaseFileNameBe(originalFilename + "_2_Updated");

        modes.push_back(mode);

        OOFEM_LOG_INFO("************************** Microproblem at %p instanciated \n", myEngngModel);
    }
}

void
SolutionbasedShapeFunction :: computeBaseFunctionValueAt(FloatArray &answer, FloatArray &coords, IntArray &dofIDs, EngngModel &myEngngModel)
{
    answer.resize( dofIDs.giveSize() );
    answer.zero();

    if ( useConstantBase ) {
        for ( int i = 1; i <= answer.giveSize(); i++ ) {
            answer.at(i) = 1;
        }
        ;
    } else {
        std :: vector< FloatArray * >checkcoords;
        std :: vector< int >permuteIndex;
        int n = 0;

        if ( !isLoaded ) {
            loadProblem();
        }

        myEngngModel.giveDomain(1)->giveSpatialLocalizer()->init(false);

        if ( this->giveDomain()->giveNumberOfSpatialDimensions() == 2 ) {
            coords.resize(2);
        }

        // Determine if current coordinate is at a max or min point and if so, which type (on a surface, edge or a corner?)
        // n is the number of dimensions at which the point is an extremum, permuteIndex tells in which dimension the coordinate is max/min

        int thisMask = 0;

        for ( int i = 1; i <= coords.giveSize(); i++ ) {
            if ( ( fabs( maxCoord.at(i) - coords.at(i) ) < TOL ) || ( fabs( minCoord.at(i) - coords.at(i) ) < TOL ) ) {
                permuteIndex.push_back(i);
                n++;
                //thisMask = thisMask + pow(2.0, i - 1);   // compiler warning on conversion from double to int
                thisMask = thisMask + ( 0x01 << ( i - 1 ) );
            }
        }
        int _s = 0x01 << n;
        for ( int i = 0; i < _s; i++ ) {
            int mask = i, counter = 1;
            FloatArray *newCoord = new(FloatArray) ( coords.giveSize() );
            * newCoord = coords;

            for ( int j = 1; j <= n; j++ ) {
                double d = 0.0; //TOL;
                if ( ( mask & 1 ) == 0 ) { // Max
                    newCoord->at( permuteIndex.at(counter - 1) ) = minCoord.at( permuteIndex.at(counter - 1) ) + d;
                } else { // Min
                    newCoord->at( permuteIndex.at(counter - 1) ) = maxCoord.at( permuteIndex.at(counter - 1) ) - d;
                }
                counter++;
                mask = mask >> 1;
            }
            checkcoords.push_back(newCoord);
        }

        // The followind define allows for use of weakly periodic bc to be copied. This does not comply with the theory but is used to check the validity of the code.
#define USEWPBC 0

#if USEWPBC == 1
        FloatArray *tempCoord = new FloatArray;
        * tempCoord = coords;
        checkcoords.clear();
        checkcoords.push_back(tempCoord);
#endif
        FloatArray values;
        for ( size_t i = 0; i < checkcoords.size(); i++ ) {
            giveValueAtPoint(values, * checkcoords.at(i), dofIDs, myEngngModel);
            //printf("Values at (%f, %f, %f) are [%f, %f, %f]\n", checkcoords.at(i)->at(1), checkcoords.at(i)->at(2), checkcoords.at(i)->at(3), values.at(1), values.at(2), values.at(3));
#if USEWPBC == 1
            for ( int j = 1; j <= values.giveSize(); j++ ) {
                answer.at(j) = values.at(j);
            }
#else
            for ( int j = 1; j <= values.giveSize(); j++ ) {
                answer.at(j) = answer.at(j) + values.at(j) / ( ( double ) pow(2.0, n) );
            }
#endif
            delete( checkcoords.at(i) );
        }
    }
}