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

void PrescribedGradientBCPeriodic :: computeTangent(FloatMatrix &E, TimeStep *tStep)
{
    EModelDefaultEquationNumbering fnum;
    DofIDEquationNumbering pnum(true, strain_id);
    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);
    Kfp->buildInternalStructure(rve, this->domain->giveNumber(), fnum, pnum);
    Kpp->buildInternalStructure(rve, this->domain->giveNumber(), pnum);
    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);

    int neq = Kfp->giveNumberOfRows();
    int nsd = this->domain->giveNumberOfSpatialDimensions();
    int ncomp = nsd * nsd;

    FloatMatrix grad_pert(ncomp, ncomp), rhs, sol(neq, ncomp);
    grad_pert.resize(ncomp, ncomp); // In fact, npeq should most likely equal ndev
    grad_pert.beUnitMatrix();

    // Compute the solution to each of the pertubation of eps
    Kfp->times(grad_pert, rhs);
    solver->solve(*Kff, rhs, sol);

    // Compute the solution to each of the pertubation of eps
    FloatMatrix E_tmp;
    Kfp->timesT(sol, E_tmp); // Assuming symmetry of stiffness matrix
    // This is probably always zero, but for generality
    FloatMatrix tmpMat;
    Kpp->times(grad_pert, tmpMat);
    E_tmp.subtract(tmpMat);
    E_tmp.times( - 1.0 / ( this->domainSize(this->giveDomain(), this->set) + this->domainSize(this->giveDomain(), this->masterSet) ));

    E.resize(E_tmp.giveNumberOfRows(), E_tmp.giveNumberOfColumns());
    if ( nsd == 3 ) {
        if ( E_tmp.giveNumberOfRows() == 6 ) {
            E.assemble(E_tmp, {1, 6, 5, 6, 2, 4, 5, 4, 3});
        } else {
            E.assemble(E_tmp, {1, 9, 8, 6, 2, 7, 5, 4, 3});
        }
    } else if ( nsd == 2 ) {
        E.assemble(E_tmp, {1, 4, 3, 2});
    } else {
        E = E_tmp;
    }
}

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 MixedGradientPressureWeakPeriodic :: computeTangents(FloatMatrix &Ed, FloatArray &Ep, FloatArray &Cd, double &Cp, TimeStep *tStep)
{
    //double size = this->domainSize();
    // Fetch some information from the engineering model
    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);
    EModelDefaultEquationNumbering fnum;
    Set *set = this->giveDomain()->giveSet(this->set);
    const IntArray &boundaries = set->giveBoundaryList();
    double rve_size = this->domainSize();

    // Set up and assemble tangent FE-matrix which will make up the sensitivity analysis for the macroscopic material tangent.
    std :: unique_ptr< SparseMtrx > Kff( classFactory.createSparseMtrx(stype) );
    if ( !Kff ) {
        OOFEM_ERROR("Couldn't create sparse matrix of type %d\n", stype);
    }
    Kff->buildInternalStructure(rve, this->domain->giveNumber(), fnum);
    rve->assemble(*Kff, tStep,TangentStiffnessMatrix, fnum, fnum, this->domain);

    // Setup up indices and locations
    int neq = Kff->giveNumberOfRows();

    // Indices and such of internal dofs
    int nsd = this->devGradient.giveNumberOfColumns();
    int nd = nsd * ( 1 + nsd ) / 2;

    IntArray t_loc, e_loc;

    // Fetch the positions;
    this->tractionsdman->giveLocationArray( t_id, t_loc, fnum );
    this->voldman->giveLocationArray( v_id, e_loc, fnum );

    // Matrices and arrays for sensitivities
    FloatMatrix ddev_pert(neq, nd);
    FloatArray p_pert(neq);

    // Solutions for the pertubations:
    FloatMatrix s_d(neq, nd);
    FloatArray s_p(neq);

    // Unit pertubations for d_dev
    for ( int i = 1; i <= nd; ++i ) {
        FloatMatrix ddev_tmp;
        FloatArray tmp(neq), fe, fetot, ddevVoigt(nd);
        ddevVoigt.at(i) = 1.0;
        this->constructFullMatrixForm(ddev_tmp, ddevVoigt);
        for ( int k = 1; k <= boundaries.giveSize() / 2; ++k ) {
            Element *e = this->giveDomain()->giveElement( boundaries.at(k * 2 - 1) );
            int boundary = boundaries.at(k * 2);

            this->integrateTractionDev(fe, e, boundary, ddev_tmp);
            fetot.subtract(fe);
        }
        tmp.assemble(fetot, t_loc);
        ddev_pert.setColumn(tmp, i);
    }

    // Unit pertubation for p
    p_pert.zero();
    p_pert.at( e_loc.at(1) ) = - 1.0 * rve_size;

    // Solve all sensitivities
    solver->solve(*Kff, ddev_pert, s_d);
    solver->solve(*Kff, p_pert, s_p);

    // Extract the tractions from the sensitivity solutions s_d and s_p:
    FloatArray tractions_p( t_loc.giveSize() );
    FloatMatrix tractions_d(t_loc.giveSize(), nd);
    tractions_p.beSubArrayOf(s_p, t_loc);
    for ( int i = 1; i <= nd; ++i ) {
        for ( int k = 1; k <= t_loc.giveSize(); ++k ) {
            tractions_d.at(k, i) = s_d.at(t_loc.at(k), i);
        }
    }

    // The de/dp tangent:
    Cp = - s_p.at( e_loc.at(1) );
    // The de/dd tangent:
    Cd.resize(nd);
    if ( nsd == 3 ) {
        Cd.at(1) = s_d.at(e_loc.at(1), 1);
        Cd.at(2) = s_d.at(e_loc.at(1), 2);
        Cd.at(3) = s_d.at(e_loc.at(1), 3);
        Cd.at(4) = s_d.at(e_loc.at(1), 4);
        Cd.at(5) = s_d.at(e_loc.at(1), 5);
        Cd.at(6) = s_d.at(e_loc.at(1), 6);
    } else if ( nsd == 2 ) {
        Cd.at(1) = s_d.at(e_loc.at(1), 1);
        Cd.at(2) = s_d.at(e_loc.at(1), 2);
        Cd.at(3) = s_d.at(e_loc.at(1), 3);
    }
    // The ds/de tangent:
    Ep = Cd; // This is much simpler, but it's "only" correct if there is a potential (i.e. symmetric problems). This could be generalized if needed.
    // The ds/dd tangent:
    Ed.resize(nd, nd);
    for ( int dpos = 1; dpos <= nd; ++dpos ) {
        FloatArray tractions, sigmaDev;
//.........这里部分代码省略.........