C++ Dof类代码示例

double PrescribedGradientBCPeriodic :: giveUnknown(double val, ValueModeType mode, TimeStep *tStep, ActiveDof *dof)
{
    DofManager *master = this->domain->giveDofManager(this->slavemap[dof->giveDofManager()->giveNumber()]);
    DofIDItem id = dof->giveDofID();
    FloatArray *coords = dof->giveDofManager()->giveCoordinates();
    FloatArray *masterCoords = master->giveCoordinates();
    FloatArray dx, uM;
    dx.beDifferenceOf(* coords, * masterCoords );

    int ind;
    if ( id == D_u || id == V_u || id == P_f || id == T_f ) {
        ind = 1;
    } else if ( id == D_v || id == V_v ) {
        ind = 2;
    } else { /*if ( id == D_w || id == V_w )*/   // 3D only:
        ind = 3;
    }

    FloatMatrix grad(3, 3);
    for ( int i = 0; i < this->strain_id.giveSize(); ++i ) {
        Dof *dof = this->strain->giveDofWithID(strain_id[i]);
        grad(i % 3, i / 3) = dof->giveUnknown(mode, tStep);
    }
    uM.beProductOf(grad, dx); // The "jump" part of the unknown ( u^+ = [[u^M]] + u^- )

    return val + uM.at(ind);
}

bool
NodeErrorCheckingRule :: check(Domain *domain, TimeStep *tStep)
{
    // Rule doesn't apply yet.
    if ( tStep->giveNumber() != tstep ) {
        return true;
    }

    DofManager *dman = domain->giveGlobalDofManager(number);
    if ( !dman ) {
        if ( domain->giveEngngModel()->isParallel() ) {
            return true;
        } else {
            OOFEM_WARNING("Dof manager %d not found.", number);
            return false;
        }
    }

    if ( dman->giveParallelMode() == DofManager_remote || dman->giveParallelMode() == DofManager_null ) {
        return true;
    }

    Dof *dof = dman->giveDofWithID(dofid);

    double dmanValue = dof->giveUnknown(mode, tStep);
    bool check = checkValue(dmanValue);
    if ( !check ) {
        OOFEM_WARNING("Check failed in: tstep %d, node %d, dof %d, mode %d:\n"
                      "value is %.8e, but should be %.8e ( error is %e but tolerance is %e )",
                      tstep, number, dofid, mode,
                      dmanValue, value, fabs(dmanValue-value), tolerance );
    }
    return check;
}

void Tr21Stokes :: giveIntegratedVelocity(FloatMatrix &answer, TimeStep *tStep )
{
    /*
    * Integrate velocity over element
    */

    IntegrationRule *iRule = integrationRulesArray [ 0 ];
    FloatMatrix v, v_gamma, ThisAnswer, boundaryV, Nmatrix;
    double detJ;
    FloatArray *lcoords, N;
    int i, j, k=0;
    Dof *d;
    GaussPoint *gp;

    v.resize(12,1);
    v.zero();
    boundaryV.resize(2,1);


    for (i=1; i<=this->giveNumberOfDofManagers(); i++) {
        for (j=1; j<=this->giveDofManager(i)->giveNumberOfDofs(); j++) {
            d = this->giveDofManager(i)->giveDof(j);
            if ((d->giveDofID()==V_u) || (d->giveDofID()==V_v)) {
                k=k+1;
                v.at(k,1)=d->giveUnknown(EID_ConservationEquation, VM_Total, tStep);
            /*} else if (d->giveDofID()==A_x) {
                boundaryV.at(1,1)=d->giveUnknown(EID_ConservationEquation, VM_Total, tStep);
            } else if (d->giveDofID()==A_y) {
                boundaryV.at(2,1)=d->giveUnknown(EID_ConservationEquation, VM_Total, tStep);*/
            }
        }
    }

    answer.resize(2,1);
    answer.zero();

    Nmatrix.resize(2,12);

    for (i=0; i<iRule->getNumberOfIntegrationPoints(); i++) {

        gp = iRule->getIntegrationPoint(i);

        lcoords = gp->giveCoordinates();

        this->interpolation_quad.evalN(N, *lcoords, FEIElementGeometryWrapper(this));
        detJ = this->interpolation_quad.giveTransformationJacobian(*lcoords, FEIElementGeometryWrapper(this));

        N.times(detJ*gp->giveWeight());

        for (j=1; j<=6;j++) {
            Nmatrix.at(1,j*2-1)=N.at(j);
            Nmatrix.at(2,j*2)=N.at(j);
        }

        ThisAnswer.beProductOf(Nmatrix,v);
        answer.add(ThisAnswer);

    }

}

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
SolutionbasedShapeFunction :: copyDofManagersToSurfaceData(modeStruct *mode, IntArray nodeList, bool isPlus, bool isMinus, bool isZeroBoundary)
{
    for ( int i = 1; i <= nodeList.giveSize(); i++ ) {
        FloatArray values;
        DofManager *dman = mode->myEngngModel->giveDomain(1)->giveDofManager( nodeList.at(i) );

        computeBaseFunctionValueAt(values, * dman->giveCoordinates(), this->dofs, * mode->myEngngModel);

        for ( int j = 1; j <= this->dofs.giveSize(); j++ ) {
            SurfaceDataStruct *surfaceData = new(SurfaceDataStruct);
            Dof *d = dman->giveDofWithID( dofs.at(j) );

            surfaceData->DofID = ( DofIDItem ) this->dofs.at(j);
            surfaceData->DofMan = dman;
            surfaceData->isPlus = isPlus;
            surfaceData->isMinus = isMinus;
            surfaceData->isZeroBoundary = isZeroBoundary;
            surfaceData->isFree = d->giveBcId() == 0;
            surfaceData->value = values.at(j);

            mode->SurfaceData.push_back(surfaceData);
        }
    }
}

int
EIPrimaryUnknownMapper :: mapAndUpdate(FloatArray &answer, ValueModeType mode,
                                       Domain *oldd, Domain *newd,  TimeStep *tStep)
{
    int inode, nd_nnodes = newd->giveNumberOfDofManagers();
    int nsize = newd->giveEngngModel()->giveNumberOfDomainEquations( newd->giveNumber(), EModelDefaultEquationNumbering() );
    FloatArray unknownValues;
    IntArray dofidMask, locationArray;
    IntArray reglist;
#ifdef OOFEM_MAPPING_CHECK_REGIONS
    ConnectivityTable *conTable = newd->giveConnectivityTable();
    const IntArray *nodeConnectivity;
#endif

    answer.resize(nsize);
    answer.zero();

    for ( inode = 1; inode <= nd_nnodes; inode++ ) {
        DofManager *node = newd->giveNode(inode);
        /* HUHU CHEATING */
#ifdef __PARALLEL_MODE
        if ( ( node->giveParallelMode() == DofManager_null ) ||
            ( node->giveParallelMode() == DofManager_remote ) ) {
            continue;
        }

#endif

#ifdef OOFEM_MAPPING_CHECK_REGIONS
        // build up region list for node
        nodeConnectivity = conTable->giveDofManConnectivityArray(inode);
        reglist.resize( nodeConnectivity->giveSize() );
        reglist.clear();
        for ( int indx = 1; indx <= nodeConnectivity->giveSize(); indx++ ) {
            reglist.insertSortedOnce( newd->giveElement( nodeConnectivity->at(indx) )->giveRegionNumber() );
        }

#endif
        ///@todo Shouldn't we pass a primary field or something to this function?
        if ( this->evaluateAt(unknownValues, dofidMask, mode, oldd, * node->giveCoordinates(), reglist, tStep) ) {
            ///@todo This doesn't respect local coordinate systems in nodes. Supporting that would require major reworking.
            for ( int ii = 1; ii <= dofidMask.giveSize(); ii++ ) {
                // exclude slaves; they are determined from masters
                auto it = node->findDofWithDofId((DofIDItem)dofidMask.at(ii));
                if ( it != node->end() ) {
                    Dof *dof = *it;
                    if ( dof->isPrimaryDof() ) {
                        int eq = dof->giveEquationNumber(EModelDefaultEquationNumbering());
                        answer.at( eq ) += unknownValues.at(ii);
                    }
                }
            }
        } else {
            OOFEM_ERROR("evaluateAt service failed for node %d", inode);
        }
    }

    return 1;
}

    bool FileInfoSkel<SkeletonType>::saveFile( const char* _fName) const {
        using namespace std;
        ofstream output(_fName);
        if(output.fail()){
            cout<<"Unable to save file "<<_fName<<endl;
            return false;
        }
        output<<fixed<<setprecision(6);

        int _numLinks = mSkel->mNumNodes;

        output<<"dofs {\n";
        for(int i=0; i<mSkel->getNumDofs(); i++){
            Dof *d = mSkel->getDof(i);
            if(d->getJoint()->getChildNode()->getSkelIndex() >= _numLinks) break;
            output<<'\t'<<d->getName()<<" { "<<d->getValue()<<", "<<d->getMin()<<", "<<d->getMax()<<" }\n";
        }
        output<<unitlength<<" { 1.0, 1.0, 1.0 }\n";
        output<<"}\n";

        // masses
        output<<"\nmass {\n";
        for(int i=0; i<mSkel->mNumNodes; i++){
            if(i>=_numLinks) break;
            string massname = mSkel->getNode(i)->getName();
            massname += "_mass";
            output<<massname<<" { "<<mSkel->getNode(i)->getMass()<<" }\n";

        }
        output<<"}\n";

        // nodes
        saveBodyNodeTree(mSkel->mRoot, output, _numLinks);

        // markers
        output<<"\nmarkers {\n";
        for(int i=0; i<mSkel->getNumMarkers(); i++){
            if(mSkel->getMarker(i)->getNode()->getSkelIndex()>=_numLinks) break;
            string hname = mSkel->getMarker(i)->getName();
            if(hname.empty()){
                ostringstream ss;
                ss<<"marker"<<i;
                hname = ss.str();
            }
            Eigen::Vector3d lc = mSkel->getMarker(i)->getLocalCoords();
            output<<hname<<" { "<<"<"<<lc[0]<<", "<<lc[1]<<", "<<lc[2]<<">, "<<i<<", "<<mSkel->getMarker(i)->getNode()->getName()<<" } "<<endl;

        }
        output<<"}\n";

        cout<<"\n";

        return true;
    }

void
Node2NodeContactL :: computeContactTractionAt(GaussPoint *gp, FloatArray &t, FloatArray &gap, TimeStep *tStep)
{
    // should be replaced with a call to constitutive model
    // gap should be in a local system
    if ( gap.at(1) < 0.0 ) {
        
        Dof *dof = masterNode->giveDofWithID( this->giveDofIdArray().at(1) );
        double lambda = dof->giveUnknown(VM_Total, tStep);
        t = {lambda, 0.0, 0.0};
        //printf("lambda %e \n\n", lambda);
    } else {
        t = {0.0, 0.0, 0.0};
    }
  
} 

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

void
MatlabExportModule :: doOutputData(TimeStep *tStep, FILE *FID)
{
    Domain *domain  = emodel->giveDomain(1);
    std :: vector< int >DofIDList;
    std :: vector< int > :: iterator it;
    std :: vector< std :: vector< double > * >valuesList;
    std :: vector< double > *values;

    for ( int i = 1; i <= domain->giveNumberOfDofManagers(); i++ ) {
        for ( int j = 1; j <= domain->giveDofManager(i)->giveNumberOfDofs(); j++ ) {
            Dof *thisDof;
            thisDof = domain->giveDofManager(i)->giveDof(j);
            it = std :: find( DofIDList.begin(), DofIDList.end(), thisDof->giveDofID() );

            if ( it == DofIDList.end() ) {
                DofIDList.push_back( thisDof->giveDofID() );
                values = new( std :: vector< double > );
                valuesList.push_back(values);
            } else {
                int pos = it - DofIDList.begin();
                values = valuesList.at(pos);
            }

            double value = thisDof->giveUnknown(EID_MomentumBalance, VM_Total, tStep);
            values->push_back(value);
        }
    }

    fprintf(FID, "\tdata.DofIDs=[");
    for ( size_t i = 0; i < DofIDList.size(); i++ ) {
        fprintf( FID, "%u, ", DofIDList.at(i) );
    }

    fprintf(FID, "];\n");

    for ( size_t i = 0; i < valuesList.size(); i++ ) {
        fprintf(FID, "\tdata.a{%lu}=[", static_cast<long unsigned int>(i) + 1);
        for ( size_t j = 0; j < valuesList.at(i)->size(); j++ ) {
            fprintf( FID, "%f,", valuesList.at(i)->at(j) );
        }

        fprintf(FID, "];\n");
    }
}

void
PhaseFieldElement :: computeLocationArrayOfDofIDs( const IntArray &dofIdArray, IntArray &answer )
{
    // Routine to extract compute the location array an element given an dofid array.
    answer.clear();
    NLStructuralElement *el = this->giveElement();
    int k = 0;
    for ( int i = 1; i <= el->giveNumberOfDofManagers(); i++ ) {
        DofManager *dMan = el->giveDofManager( i );
        for ( int j = 1; j <= dofIdArray.giveSize( ); j++ ) {

            if ( dMan->hasDofID( (DofIDItem) dofIdArray.at( j ) ) ) {
                Dof *d = dMan->giveDofWithID( dofIdArray.at( j ) );
                answer.followedBy( k + d->giveNumber( ) );
            }
        }
        k += dMan->giveNumberOfDofs( );
    }
}

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

void
CoupledFieldsElement :: computeLocationArrayOfDofIDs(const IntArray &dofIdArray, IntArray &answer)
{
    // Routine to extract compute the location array an element given an dofid array.
    answer.resize(0);
    
    int k = 0;
    for ( int i = 1; i <= numberOfDofMans; i++ ) {
        DofManager *dMan = this->giveDofManager(i);        
        for (int j = 1; j <= dofIdArray.giveSize(); j++ ) {   
            
            if ( dMan->hasDofID( (DofIDItem) dofIdArray.at(j) ) ) {
                Dof *d = dMan->giveDofWithID( dofIdArray.at(j) );
                answer.followedBy( k + d->giveNumber() );
                //answer.followedBy( k + j );
            }
        }
        k += dMan->giveNumberOfDofs( );
    }
}

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
CoupledFieldsElement :: computeVectorOfDofIDs(const IntArray &dofIdArray, ValueModeType valueMode, TimeStep *stepN, FloatArray &answer)
{
    // Routine to extract the solution vector for an element given an dofid array.
    // Size will be numberOfDofs and if a certain dofId does not exist a zero is used as value. 
    
    answer.resize( numberOfDofMans * dofIdArray.giveSize() ); // equal number of nodes for all fields
    answer.zero();
    int k = 1;
    for ( int i = 1; i <= numberOfDofMans; i++ ) {
        DofManager *dMan = this->giveDofManager(i);        
        for (int j = 1; j <= dofIdArray.giveSize(); j++ ) {   
            
            if ( dMan->hasDofID( (DofIDItem) dofIdArray.at(j) ) ) {
                Dof *d = dMan->giveDofWithID( dofIdArray.at(j) );
                answer.at(k) = d->giveUnknown(valueMode, stepN);
            }
            k++;
        }
    }
}