C++ Dof::__giveEquationNumber方法代码示例

int
StructuralEngngModel :: packDofManagers(FloatArray *src, ProcessCommunicator &processComm, bool prescribedEquations)
{
    int result = 1;
    int i, size;
    int j, ndofs, eqNum;
    Domain *domain = this->giveDomain(1);
    IntArray const *toSendMap = processComm.giveToSendMap();
    ProcessCommunicatorBuff *pcbuff = processComm.giveProcessCommunicatorBuff();
    DofManager *dman;
    Dof *jdof;

    size = toSendMap->giveSize();
    for ( i = 1; i <= size; i++ ) {
        dman = domain->giveDofManager( toSendMap->at(i) );
        ndofs = dman->giveNumberOfDofs();
        for ( j = 1; j <= ndofs; j++ ) {
            jdof = dman->giveDof(j);
            if ( prescribedEquations ) {
                eqNum = jdof->__givePrescribedEquationNumber();
            } else {
                eqNum = jdof->__giveEquationNumber();
            }
            if ( jdof->isPrimaryDof() && eqNum ) {
                    result &= pcbuff->packDouble( src->at(eqNum) );
            }
        }
    }

    return result;
}

void IncrementalLinearStatic :: updateDofUnknownsDictionary(DofManager *inode, TimeStep *tStep)
{
    // update DOF unknowns dictionary, where
    // unknowns are hold instead of keeping them in global unknowns
    // vectors in engng instances
    // this is necessary, because during solution equation numbers for
    // particular DOFs may changed, and it is necessary to keep them
    // in DOF level.

    int ndofs = inode->giveNumberOfDofs();
    Dof *iDof;
    double val;
    for ( int i = 1; i <= ndofs; i++ ) {
        iDof = inode->giveDof(i);
        // skip slave DOFs (only master (primary) DOFs have to be updated).
        if (!iDof->isPrimaryDof()) continue;
        val = iDof->giveUnknown(VM_Total, tStep);
        if ( !iDof->hasBc(tStep) ) {
            val += this->incrementOfDisplacementVector.at( iDof->__giveEquationNumber() );
        }

        iDof->updateUnknownsDictionary(tStep, VM_Total_Old, val);
        iDof->updateUnknownsDictionary(tStep, VM_Total, val);
    }
}

int
StructuralEngngModel :: unpackDofManagers(FloatArray *dest, ProcessCommunicator &processComm, bool prescribedEquations)
{
    int result = 1;
    int i, size;
    int j, ndofs, eqNum;
    Domain *domain = this->giveDomain(1);
    dofManagerParallelMode dofmanmode;
    IntArray const *toRecvMap = processComm.giveToRecvMap();
    ProcessCommunicatorBuff *pcbuff = processComm.giveProcessCommunicatorBuff();
    DofManager *dman;
    Dof *jdof;
    double value;


    size = toRecvMap->giveSize();
    for ( i = 1; i <= size; i++ ) {
        dman = domain->giveDofManager( toRecvMap->at(i) );
        ndofs = dman->giveNumberOfDofs();
        dofmanmode = dman->giveParallelMode();
        for ( j = 1; j <= ndofs; j++ ) {
            jdof = dman->giveDof(j);
            if ( prescribedEquations ) {
                eqNum = jdof->__givePrescribedEquationNumber();
            } else {
                eqNum = jdof->__giveEquationNumber();
            }
            if ( jdof->isPrimaryDof() && eqNum ) {
                result &= pcbuff->unpackDouble(value);
                if ( dofmanmode == DofManager_shared ) {
                    dest->at(eqNum) += value;
                } else if ( dofmanmode == DofManager_remote ) {
                    dest->at(eqNum)  = value;
                } else {
                    _error("unpackReactions: unknown dof namager parallel mode");
                }
            }
        }
    }

    return result;
}

int
NonLinearDynamic :: estimateMaxPackSize(IntArray &commMap, CommunicationBuffer &buff, int packUnpackType)
{
    int mapSize = commMap.giveSize();
    int i, j, ndofs, count = 0, pcount = 0;
    IntArray locationArray;
    Domain *domain = this->giveDomain(1);
    DofManager *dman;
    Dof *jdof;

    if ( packUnpackType == ProblemCommMode__ELEMENT_CUT ) {
        for ( i = 1; i <= mapSize; i++ ) {
            count += domain->giveDofManager( commMap.at(i) )->giveNumberOfDofs();
        }

        return ( buff.givePackSize(MPI_DOUBLE, 1) * count );
    } else if ( packUnpackType == ProblemCommMode__NODE_CUT ) {
        for ( i = 1; i <= mapSize; i++ ) {
            ndofs = ( dman = domain->giveDofManager( commMap.at(i) ) )->giveNumberOfDofs();
            for ( j = 1; j <= ndofs; j++ ) {
                jdof = dman->giveDof(j);
                if ( jdof->isPrimaryDof() && ( jdof->__giveEquationNumber() ) ) {
                    count++;
                } else {
                    pcount++;
                }
            }
        }

        //printf ("\nestimated count is %d\n",count);
        return ( buff.givePackSize(MPI_DOUBLE, 1) * max(count, pcount) );
    } else if ( packUnpackType == ProblemCommMode__REMOTE_ELEMENT_MODE ) {
        for ( i = 1; i <= mapSize; i++ ) {
            count += domain->giveElement( commMap.at(i) )->estimatePackSize(buff);
        }

        return count;
    }

    return 0;
}

void
NonStationaryTransportProblem :: applyIC(TimeStep *stepWhenIcApply)
{
    Domain *domain = this->giveDomain(1);
    int neq =  this->giveNumberOfEquations(EID_ConservationEquation);
    FloatArray *solutionVector;
    double val;

#ifdef VERBOSE
    OOFEM_LOG_INFO("Applying initial conditions\n");
#endif
    int nDofs, j, k, jj;
    int nman  = domain->giveNumberOfDofManagers();
    DofManager *node;
    Dof *iDof;

    UnknownsField->advanceSolution(stepWhenIcApply);
    solutionVector = UnknownsField->giveSolutionVector(stepWhenIcApply);
    solutionVector->resize(neq);
    solutionVector->zero();

    for ( j = 1; j <= nman; j++ ) {
        node = domain->giveDofManager(j);
        nDofs = node->giveNumberOfDofs();

        for ( k = 1; k <= nDofs; k++ ) {
            // ask for initial values obtained from
            // bc (boundary conditions) and ic (initial conditions)
            iDof  =  node->giveDof(k);
            if ( !iDof->isPrimaryDof() ) {
                continue;
            }

            jj = iDof->__giveEquationNumber();
            if ( jj ) {
                val = iDof->giveUnknown(EID_ConservationEquation, VM_Total, stepWhenIcApply);
                solutionVector->at(jj) = val;
                //update in dictionary, if the problem is growing/decreasing
                if ( this->changingProblemSize ) {
                    iDof->updateUnknownsDictionary(stepWhenIcApply, EID_MomentumBalance, VM_Total, val);
                }
            }
        }
    }

    int nelem = domain->giveNumberOfElements();
    
    //project initial temperature to integration points

//     for ( j = 1; j <= nelem; j++ ) {
//         domain->giveElement(j)->updateInternalState(stepWhenIcApply);
//     }

#ifdef __CEMHYD_MODULE
    // Not relevant in linear case, but needed for CemhydMat for temperature averaging before solving balance equations
    // Update element state according to given ic
    TransportElement *element;
    CemhydMat *cem;
    for ( j = 1; j <= nelem; j++ ) {
        element = ( TransportElement * ) domain->giveElement(j);
        //assign status to each integration point on each element
        if ( element->giveMaterial()->giveClassID() == CemhydMatClass ) {
            element->giveMaterial()->initMaterial(element); //create microstructures and statuses on specific GPs
            element->updateInternalState(stepWhenIcApply);   //store temporary unequilibrated temperature
            element->updateYourself(stepWhenIcApply);   //store equilibrated temperature
            cem = ( CemhydMat * ) element->giveMaterial();
            cem->clearWeightTemperatureProductVolume(element);
            cem->storeWeightTemperatureProductVolume(element, stepWhenIcApply);
        }
    }

    //perform averaging on each material instance of CemhydMatClass
    int nmat = domain->giveNumberOfMaterialModels();
    for ( j = 1; j <= nmat; j++ ) {
        if ( domain->giveMaterial(j)->giveClassID() == CemhydMatClass ) {
            cem = ( CemhydMat * ) domain->giveMaterial(j);
            cem->averageTemperature();
        }
    }
#endif //__CEMHYD_MODULE
}

//.........这里部分代码省略.........
                OOFEM_LOG_RELEVANT("DEIDynamic: deltaT reduced to %e\n", maxDt);
                deltaT = maxDt;
                tStep->setTimeIncrement(deltaT);
            }
        }


        //
        // special init step - compute displacements at tstep 0
        //
        displacementVector.resize(neq);
        displacementVector.zero();
        nextDisplacementVector.resize(neq);
        nextDisplacementVector.zero();
        velocityVector.resize(neq);
        velocityVector.zero();
        accelerationVector.resize(neq);
        accelerationVector.zero();


        for ( j = 1; j <= nman; j++ ) {
            node = domain->giveDofManager(j);
            nDofs = node->giveNumberOfDofs();

            for ( k = 1; k <= nDofs; k++ ) {
                // ask for initial values obtained from
                // bc (boundary conditions) and ic (initial conditions)
                // now we are setting initial cond. for step -1.
                iDof  =  node->giveDof(k);
                if ( !iDof->isPrimaryDof() ) {
                    continue;
                }

                jj = iDof->__giveEquationNumber();
                if ( jj ) {
                    nextDisplacementVector.at(jj) = iDof->giveUnknown(EID_MomentumBalance, VM_Total, tStep);
                    // become displacementVector after init
                    velocityVector.at(jj)     = iDof->giveUnknown(EID_MomentumBalance, VM_Velocity, tStep);
                    // accelerationVector = iDof->giveUnknown(AccelerartionVector,tStep) ;
                }
            }
        }

        for ( j = 1; j <= neq; j++ ) {
            nextDisplacementVector.at(j) -= velocityVector.at(j) * ( deltaT );
        }

        return;
    } // end of init step

#ifdef VERBOSE
    OOFEM_LOG_INFO("Assembling right hand side\n");
#endif


    c1 = ( 1. / ( deltaT * deltaT ) );
    c2 = ( 1. / ( 2. * deltaT ) );
    c3 = ( 2. / ( deltaT * deltaT ) );

    previousDisplacementVector = displacementVector;
    displacementVector         = nextDisplacementVector;

    //
    // assembling the element part of load vector
    //
    loadVector.resize( this->giveNumberOfEquations(EID_MomentumBalance) );

void
NonLinearDynamic :: proceedStep(int di, TimeStep *tStep)
{
    // creates system of governing eq's and solves them at given time step
    // first assemble problem at current time step

    int neq = this->giveNumberOfDomainEquations(1, EModelDefaultEquationNumbering());

    // Time-stepping constants
    this->determineConstants(tStep);

    if ( ( tStep->giveNumber() == giveNumberOfFirstStep() ) && initFlag ) {
        // Initialization
        incrementOfDisplacement.resize(neq);
        incrementOfDisplacement.zero();
        totalDisplacement.resize(neq);
        totalDisplacement.zero();
        velocityVector.resize(neq);
        velocityVector.zero();
        accelerationVector.resize(neq);
        accelerationVector.zero();
        internalForces.resize(neq);
        internalForces.zero();
        previousIncrementOfDisplacement.resize(neq);
        previousIncrementOfDisplacement.zero();
        previousTotalDisplacement.resize(neq);
        previousTotalDisplacement.zero();
        previousVelocityVector.resize(neq);
        previousVelocityVector.zero();
        previousAccelerationVector.resize(neq);
        previousAccelerationVector.zero();
        previousInternalForces.resize(neq);
        previousInternalForces.zero();

        TimeStep *stepWhenIcApply = new TimeStep(giveNumberOfTimeStepWhenIcApply(), this, 0,
                                                 -deltaT, deltaT, 0);

        int nDofs, j, k, jj;
        int nman = this->giveDomain(di)->giveNumberOfDofManagers();
        DofManager *node;
        Dof *iDof;

        // Considering initial conditions.
        for ( j = 1; j <= nman; j++ ) {
            node = this->giveDomain(di)->giveDofManager(j);
            nDofs = node->giveNumberOfDofs();

            for ( k = 1; k <= nDofs; k++ ) {
                // Ask for initial values obtained from
                // bc (boundary conditions) and ic (initial conditions).
                iDof  =  node->giveDof(k);
                if ( !iDof->isPrimaryDof() ) {
                    continue;
                }

                jj = iDof->__giveEquationNumber();
                if ( jj ) {
                    totalDisplacement.at(jj)  = iDof->giveUnknown(VM_Total, stepWhenIcApply);
                    velocityVector.at(jj)     = iDof->giveUnknown(VM_Velocity, stepWhenIcApply);
                    accelerationVector.at(jj) = iDof->giveUnknown(VM_Acceleration, stepWhenIcApply);
                }
            }
        }

        this->giveInternalForces(internalForces, true, di, tStep);
    }

    if ( initFlag ) {
        // First assemble problem at current time step.
        // Option to take into account initial conditions.
        if ( !effectiveStiffnessMatrix ) {
            effectiveStiffnessMatrix = classFactory.createSparseMtrx(sparseMtrxType);
            massMatrix = classFactory.createSparseMtrx(sparseMtrxType);
        }

        if ( effectiveStiffnessMatrix == NULL || massMatrix == NULL ) {
            _error("proceedStep: sparse matrix creation failed");
        }

        if ( nonlocalStiffnessFlag ) {
            if ( !effectiveStiffnessMatrix->isAsymmetric() ) {
                _error("proceedStep: effectiveStiffnessMatrix does not support asymmetric storage");
            }
        }

        effectiveStiffnessMatrix->buildInternalStructure( this, di, EID_MomentumBalance, EModelDefaultEquationNumbering() );
        massMatrix->buildInternalStructure( this, di, EID_MomentumBalance, EModelDefaultEquationNumbering() );

        // Assemble mass matrix
        this->assemble(massMatrix, tStep, EID_MomentumBalance, MassMatrix,
                       EModelDefaultEquationNumbering(), this->giveDomain(di));

        // Initialize vectors
        help.resize(neq);
        help.zero();
        rhs.resize(neq);
        rhs.zero();
        rhs2.resize(neq);
        rhs2.zero();

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

void NlDEIDynamic :: solveYourselfAt(TimeStep *tStep)
{
    //
    // Creates system of governing eq's and solves them at given time step.
    //

    Domain *domain = this->giveDomain(1);
    int neq = this->giveNumberOfEquations(EID_MomentumBalance);
    int nman  = domain->giveNumberOfDofManagers();

    DofManager *node;
    Dof *iDof;

    int nDofs;
    int i, k, j, jj;
    double coeff, maxDt, maxOm = 0.;
    double prevIncrOfDisplacement, incrOfDisplacement;

    if ( initFlag ) {
#ifdef VERBOSE
        OOFEM_LOG_DEBUG("Assembling mass matrix\n");
#endif

        //
        // Assemble mass matrix.
        //
        this->computeMassMtrx(massMatrix, maxOm, tStep);

        if ( drFlag ) {
            // If dynamic relaxation: Assemble amplitude load vector.
            loadRefVector.resize(neq);
            loadRefVector.zero();

            this->computeLoadVector(loadRefVector, VM_Total, tStep);

#ifdef __PARALLEL_MODE
            // Compute the processor part of load vector norm pMp
            this->pMp = 0.0;
            double my_pMp = 0.0, coeff = 1.0;
            int eqNum, ndofs, ndofman = domain->giveNumberOfDofManagers();
            dofManagerParallelMode dofmanmode;
            DofManager *dman;
            Dof *jdof;
            for ( int dm = 1; dm <= ndofman; dm++ ) {
                dman = domain->giveDofManager(dm);
                ndofs = dman->giveNumberOfDofs();
                dofmanmode = dman->giveParallelMode();

                // Skip all remote and null dofmanagers
                coeff = 1.0;
                if ( ( dofmanmode == DofManager_remote ) || ( ( dofmanmode == DofManager_null ) ) ) {
                    continue;
                } else if ( dofmanmode == DofManager_shared ) {
                    coeff = 1. / dman->givePartitionsConnectivitySize();
                }

                // For shared nodes we add locally an average = 1/givePartitionsConnectivitySize()*contribution,
                for ( j = 1; j <= ndofs; j++ ) {
                    jdof = dman->giveDof(j);
                    if ( jdof->isPrimaryDof() && ( eqNum = jdof->__giveEquationNumber() ) ) {
                        my_pMp += coeff * loadRefVector.at(eqNum) * loadRefVector.at(eqNum) / massMatrix.at(eqNum);
                    }
                }
            }

            // Sum up the contributions from processors.
            MPI_Allreduce(& my_pMp, & pMp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
#else
            this->pMp = 0.0;
            for ( i = 1; i <= neq; i++ ) {
                pMp += loadRefVector.at(i) * loadRefVector.at(i) / massMatrix.at(i);
            }
#endif
            // Solve for rate of loading process (parameter "c") (undamped system assumed),
            if ( dumpingCoef < 1.e-3 ) {
                c = 3.0 * this->pyEstimate / pMp / Tau / Tau;
            } else {
                c = this->pyEstimate * Tau * dumpingCoef * dumpingCoef * dumpingCoef / pMp /
                    ( -3.0 / 2.0 + dumpingCoef * Tau + 2.0 * exp(-dumpingCoef * Tau) - 0.5 * exp(-2.0 * dumpingCoef * Tau) );
            }
        }

        initFlag = 0;
    }


    if ( tStep->giveNumber() == giveNumberOfFirstStep() ) {
        //
        // Special init step - Compute displacements at tstep 0.
        //
        displacementVector.resize(neq);
        displacementVector.zero();
        previousIncrementOfDisplacementVector.resize(neq);
        previousIncrementOfDisplacementVector.zero();
        velocityVector.resize(neq);
        velocityVector.zero();
        accelerationVector.resize(neq);
        accelerationVector.zero();

        for ( j = 1; j <= nman; j++ ) {
//.........这里部分代码省略.........

void
NonLinearDynamic :: proceedStep(int di, TimeStep *tStep)
{
    // creates system of governing eq's and solves them at given time step
    // first assemble problem at current time step

    int neq = this->giveNumberOfEquations(EID_MomentumBalance);

    // Time-stepping constants
    double dt2 = deltaT * deltaT;

    if ( tStep->giveTimeDiscretization() == TD_Newmark ) {
        OOFEM_LOG_DEBUG("Solving using Newmark-beta method\n");
        a0 = 1 / ( beta * dt2 );
        a1 = gamma / ( beta * deltaT );
        a2 = 1 / ( beta * deltaT );
        a3 = 1 / ( 2 *  beta ) - 1;
        a4 = ( gamma / beta ) - 1;
        a5 = deltaT / 2 * ( gamma / beta - 2 );
        a6 = 0;
    } else if ( ( tStep->giveTimeDiscretization() == TD_TwoPointBackward ) || ( tStep->giveNumber() == giveNumberOfFirstStep() ) ) {
        if ( tStep->giveTimeDiscretization() != TD_ThreePointBackward ) {
            OOFEM_LOG_DEBUG("Solving using Backward Euler method\n");
        } else {
            OOFEM_LOG_DEBUG("Solving initial step using Three-point Backward Euler method\n");
        }
        a0 = 1 / dt2;
        a1 = 1 / deltaT;
        a2 = 1 / deltaT;
        a3 = 0;
        a4 = 0;
        a5 = 0;
        a6 = 0;
    } else if ( tStep->giveTimeDiscretization() == TD_ThreePointBackward ) {
        OOFEM_LOG_DEBUG("Solving using Three-point Backward Euler method\n");
        a0 = 2 / dt2;
        a1 = 3 / ( 2 * deltaT );
        a2 = 2 / deltaT;
        a3 = 0;
        a4 = 0;
        a5 = 0;
        a6 = 1 / ( 2 * deltaT );
    } else {
        _error("NonLinearDynamic: Time-stepping scheme not found!\n")
    }

    if ( tStep->giveNumber() == giveNumberOfFirstStep() ) {
        // Initialization
        previousIncrementOfDisplacement.resize(neq);
        previousIncrementOfDisplacement.zero();
        previousTotalDisplacement.resize(neq);
        previousTotalDisplacement.zero();
        totalDisplacement.resize(neq);
        totalDisplacement.zero();
        previousInternalForces.resize(neq);
        previousInternalForces.zero();
        incrementOfDisplacement.resize(neq);
        incrementOfDisplacement.zero();
        velocityVector.resize(neq);
        velocityVector.zero();
        accelerationVector.resize(neq);
        accelerationVector.zero();

        TimeStep *stepWhenIcApply = new TimeStep(giveNumberOfTimeStepWhenIcApply(), this, 0,
                                                 -deltaT, deltaT, 0);

        int nDofs, j, k, jj;
        int nman  = this->giveDomain(di)->giveNumberOfDofManagers();
        DofManager *node;
        Dof *iDof;

        // Considering initial conditions.
        for ( j = 1; j <= nman; j++ ) {
            node = this->giveDomain(di)->giveDofManager(j);
            nDofs = node->giveNumberOfDofs();

            for ( k = 1; k <= nDofs; k++ ) {
                // Ask for initial values obtained from
                // bc (boundary conditions) and ic (initial conditions).
                iDof  =  node->giveDof(k);
                if ( !iDof->isPrimaryDof() ) {
                    continue;
                }

                jj = iDof->__giveEquationNumber();
                if ( jj ) {
                    incrementOfDisplacement.at(jj) = iDof->giveUnknown(EID_MomentumBalance, VM_Total, stepWhenIcApply);
                    velocityVector.at(jj) = iDof->giveUnknown(EID_MomentumBalance, VM_Velocity, stepWhenIcApply);
                    accelerationVector.at(jj) = iDof->giveUnknown(EID_MomentumBalance, VM_Acceleration, stepWhenIcApply);
                }
            }
        }
    } else {
        incrementOfDisplacement.resize(neq);
        incrementOfDisplacement.zero();
    }

    if ( initFlag ) {
        // First assemble problem at current time step.
        // Option to take into account initial conditions.
//.........这里部分代码省略.........