C++ DofManager类代码示例

void
VTKXMLExportModule :: exportPrimVarAs(UnknownType valID, IntArray &mapG2L, IntArray &mapL2G,
                                      int regionDofMans, int ireg,
#ifdef __VTK_MODULE
                                      vtkSmartPointer<vtkUnstructuredGrid> &stream,
#else
                                      FILE *stream,
#endif
                                      TimeStep *tStep)
{
    Domain *d = emodel->giveDomain(1);
    InternalStateValueType type = ISVT_UNDEFINED;

    if ( ( valID == DisplacementVector ) || ( valID == EigenVector ) || ( valID == VelocityVector ) ) {
        type = ISVT_VECTOR;
    } else if ( ( valID == FluxVector ) || ( valID == PressureVector ) || ( valID == Temperature ) ) {
        type = ISVT_SCALAR;
    } else {
        OOFEM_ERROR2( "VTKXMLExportModule::exportPrimVarAs: unsupported UnknownType %s", __UnknownTypeToString(valID) );
    }

#ifdef __VTK_MODULE
    vtkSmartPointer<vtkDoubleArray> primVarArray = vtkSmartPointer<vtkDoubleArray>::New();
    primVarArray->SetName(__UnknownTypeToString(valID));

    if ( type == ISVT_SCALAR ) {
        primVarArray->SetNumberOfComponents(1);
    } else if ( type == ISVT_VECTOR ) {
        primVarArray->SetNumberOfComponents(3);
    } else {
        fprintf( stderr, "VTKXMLExportModule::exportPrimVarAs: unsupported variable type %s\n", __UnknownTypeToString(valID) );
    }
    primVarArray->SetNumberOfTuples(regionDofMans);
#else
    if ( type == ISVT_SCALAR ) {
        fprintf( stream, "<DataArray type=\"Float64\" Name=\"%s\" format=\"ascii\"> ", __UnknownTypeToString(valID) );
    } else if ( type == ISVT_VECTOR ) {
        fprintf( stream, "<DataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"3\" format=\"ascii\"> ", __UnknownTypeToString(valID) );
    } else {
        fprintf( stderr, "VTKXMLExportModule::exportPrimVarAs: unsupported variable type %s\n", __UnknownTypeToString(valID) );
    }
#endif

    DofManager *dman;
    FloatArray iVal, iValLCS;

    for ( int inode = 1; inode <= regionDofMans; inode++ ) {
        dman = d->giveNode( mapL2G.at(inode) );

        this->getPrimaryVariable(iVal, dman, tStep, valID, ireg);

        if ( type == ISVT_SCALAR ) {
#ifdef __VTK_MODULE
            primVarArray->SetTuple1(inode-1, iVal.at(1));
#else
            fprintf(stream, "%e ", iVal.at(1) );
#endif
        } else if ( type == ISVT_VECTOR ) {

            //rotate back from nodal CS to global CS if applies
            if ( (dman->giveClassID() == NodeClass) && d->giveNode( dman->giveNumber() )->hasLocalCS() ) {
                iVal.resize(3);
                iValLCS = iVal;
                iVal.beTProductOf(* d->giveNode( dman->giveNumber() )->giveLocalCoordinateTriplet(), iValLCS);
            }

#ifdef __VTK_MODULE
            primVarArray->SetTuple3(inode-1, iVal.at(1), iVal.at(2), iVal.at(3));
#else
            fprintf( stream, "%e %e %e ", iVal.at(1), iVal.at(2), iVal.at(3) );
#endif
        }
    } // end loop over nodes
#ifdef __VTK_MODULE
    if ( type == ISVT_SCALAR ) {
        stream->GetPointData()->SetActiveScalars(__UnknownTypeToString(valID));
        stream->GetPointData()->SetScalars(primVarArray);
    } else if ( type == ISVT_VECTOR ) {
        stream->GetPointData()->SetActiveVectors(__UnknownTypeToString(valID));
        stream->GetPointData()->SetVectors(primVarArray);
    }
#else
    fprintf(stream, "</DataArray>\n");
#endif
}

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
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
TrPlaneStress2dXFEM :: giveCompositeExportData(std::vector< VTKPiece > &vtkPieces, IntArray &primaryVarsToExport, IntArray &internalVarsToExport, IntArray cellVarsToExport, TimeStep *tStep)
{
    vtkPieces.resize(1);

    const int numCells = mSubTri.size();

    if(numCells == 0) {
        // Enriched but uncut element
        // Visualize as a quad
        vtkPieces[0].setNumberOfCells(1);

        int numTotalNodes = 3;
        vtkPieces[0].setNumberOfNodes(numTotalNodes);

        // Node coordinates
        std :: vector< FloatArray >nodeCoords;
        for(int i = 1; i <= 3; i++) {
            FloatArray &x = *(giveDofManager(i)->giveCoordinates());
            nodeCoords.push_back(x);

            vtkPieces[0].setNodeCoords(i, x);
        }

        // Connectivity
        IntArray nodes1 = {1, 2, 3};
        vtkPieces[0].setConnectivity(1, nodes1);

        // Offset
        int offset = 3;
        vtkPieces[0].setOffset(1, offset);

        // Cell types
        vtkPieces[0].setCellType(1, 5); // Linear triangle




        // Export nodal variables from primary fields
        vtkPieces[0].setNumberOfPrimaryVarsToExport(primaryVarsToExport.giveSize(), numTotalNodes);

        for ( int fieldNum = 1; fieldNum <= primaryVarsToExport.giveSize(); fieldNum++ ) {
            UnknownType type = ( UnknownType ) primaryVarsToExport.at(fieldNum);

            for ( int nodeInd = 1; nodeInd <= numTotalNodes; nodeInd++ ) {

                if ( type == DisplacementVector ) { // compute displacement

                        FloatArray u = {0.0, 0.0, 0.0};

                        // Fetch global coordinates (in undeformed configuration)
                        const FloatArray &x = nodeCoords[nodeInd-1];

                        // Compute local coordinates
                        FloatArray locCoord;
                        computeLocalCoordinates(locCoord, x);

                        // Compute displacement in point
                        FloatMatrix NMatrix;
                        computeNmatrixAt(locCoord, NMatrix);
                        FloatArray solVec;
                        computeVectorOf(VM_Total, tStep, solVec);
                        FloatArray uTemp;
                        uTemp.beProductOf(NMatrix, solVec);

                        if(uTemp.giveSize() == 3) {
                            u = uTemp;
                        }
                        else {
                            u = {uTemp[0], uTemp[1], 0.0};
                        }

                        vtkPieces[0].setPrimaryVarInNode(fieldNum, nodeInd, u);
                } else {
                    printf("fieldNum: %d\n", fieldNum);
                    // TODO: Implement
//                    ZZNodalRecoveryMI_recoverValues(values, layer, ( InternalStateType ) 1, tStep); // does not work well - fix
//                    for ( int j = 1; j <= numCellNodes; j++ ) {
//                        vtkPiece.setPrimaryVarInNode(fieldNum, nodeNum, values [ j - 1 ]);
//                        nodeNum += 1;
//                    }
                }
            }
        }


        // Export nodal variables from internal fields
        vtkPieces[0].setNumberOfInternalVarsToExport(0, numTotalNodes);


        // Export cell variables
        vtkPieces[0].setNumberOfCellVarsToExport(cellVarsToExport.giveSize(), 1);
        for ( int i = 1; i <= cellVarsToExport.giveSize(); i++ ) {
            InternalStateType type = ( InternalStateType ) cellVarsToExport.at(i);
            FloatArray average;
            std :: unique_ptr< IntegrationRule > &iRule = integrationRulesArray [ 0 ];
            VTKXMLExportModule :: computeIPAverage(average, iRule.get(), this, type, tStep);

            FloatArray averageV9(9);
            averageV9.at(1) = average.at(1);
//.........这里部分代码省略.........

void PLHoopStressCirc :: propagateInterfaces(Domain &iDomain, EnrichmentDomain &ioEnrDom)
{
    // Fetch crack tip data
    TipInfo tipInfoStart, tipInfoEnd;
    ioEnrDom.giveTipInfos(tipInfoStart, tipInfoEnd);
    std :: vector< TipInfo >tipInfo = {tipInfoStart, tipInfoEnd};

    SpatialLocalizer *localizer = iDomain.giveSpatialLocalizer();

    for ( size_t tipIndex = 0; tipIndex < tipInfo.size(); tipIndex++ ) {
        // Construct circle points on an arc from -90 to 90 degrees
        double angle = -90.0 + mAngleInc;
        std :: vector< double >angles;
        while ( angle <= ( 90.0 - mAngleInc ) ) {
            angles.push_back(angle * M_PI / 180.0);
            angle += mAngleInc;
        }

        const FloatArray &xT    = tipInfo [ tipIndex ].mGlobalCoord;
        const FloatArray &t             = tipInfo [ tipIndex ].mTangDir;
        const FloatArray &n             = tipInfo [ tipIndex ].mNormalDir;

        // It is meaningless to propagate a tip that is not inside any element
        Element *el = localizer->giveElementContainingPoint(tipInfo [ tipIndex ].mGlobalCoord);
        if ( el != NULL ) {
            std :: vector< FloatArray >circPoints;

            for ( size_t i = 0; i < angles.size(); i++ ) {
                FloatArray tangent(2);
                tangent.zero();
                tangent.add(cos(angles [ i ]), t);
                tangent.add(sin(angles [ i ]), n);
                tangent.normalize();

                FloatArray x(xT);
                x.add(mRadius, tangent);
                circPoints.push_back(x);
            }



            std :: vector< double >sigTTArray, sigRTArray;

            // Loop over circle points
            for ( size_t pointIndex = 0; pointIndex < circPoints.size(); pointIndex++ ) {
                FloatArray stressVec;

                if ( mUseRadialBasisFunc ) {
                    // Interpolate stress with radial basis functions

                    // Choose a cut-off length l:
                    // take the distance between two nodes in the element containing the
                    // crack tip multiplied by a constant factor.
                    // ( This choice implies that we hope that the element has reasonable
                    // aspect ratio.)
                    const FloatArray &x1 = * ( el->giveDofManager(1)->giveCoordinates() );
                    const FloatArray &x2 = * ( el->giveDofManager(2)->giveCoordinates() );
                    const double l = 1.0 * x1.distance(x2);

                    // Use the octree to get all elements that have
                    // at least one Gauss point in a certain region around the tip.
                    const double searchRadius = 3.0 * l;
                    std :: set< int >elIndices;
                    localizer->giveAllElementsWithIpWithinBox(elIndices, circPoints [ pointIndex ], searchRadius);


                    // Loop over the elements and Gauss points obtained.
                    // Evaluate the interpolation.
                    FloatArray sumQiWiVi;
                    double sumWiVi = 0.0;
                    for ( int elIndex: elIndices ) {
                        Element *gpEl = iDomain.giveElement(elIndex);
                        IntegrationRule *iRule = gpEl->giveDefaultIntegrationRulePtr();

                        for ( GaussPoint *gp_i: *iRule ) {

                            ////////////////////////////////////////
                            // Compute global gp coordinates
                            FloatArray N;
                            FEInterpolation *interp = gpEl->giveInterpolation();
                            interp->evalN( N, * ( gp_i->giveCoordinates() ), FEIElementGeometryWrapper(gpEl) );


                            // Compute global coordinates of Gauss point
                            FloatArray globalCoord(2);
                            globalCoord.zero();

                            for ( int i = 1; i <= gpEl->giveNumberOfDofManagers(); i++ ) {
                                DofManager *dMan = gpEl->giveDofManager(i);
                                globalCoord.at(1) += N.at(i) * dMan->giveCoordinate(1);
                                globalCoord.at(2) += N.at(i) * dMan->giveCoordinate(2);
                            }


                            ////////////////////////////////////////
                            // Compute weight of kernel function

                            FloatArray tipToGP;
                            tipToGP.beDifferenceOf(globalCoord, xT);
                            bool inFrontOfCrack = true;
//.........这里部分代码省略.........

int
LoadBalancer :: unpackMigratingData(Domain *d, ProcessCommunicator &pc)
{
    // create temp space for dofManagers and elements
    // merging should be made by domain ?
    // maps of new dofmanagers and elements indexed by global number

    // we can put local dofManagers and elements into maps (should be done before unpacking)
    // int nproc=this->giveEngngModel()->giveNumberOfProcesses();
    int myrank = d->giveEngngModel()->giveRank();
    int iproc = pc.giveRank();
    int _mode, _globnum, _type;
    bool _newentry;
    classType _etype;
    IntArray _partitions, local_partitions;
    //LoadBalancer::DofManMode dmode;
    DofManager *dofman;
    DomainTransactionManager *dtm = d->giveTransactionManager();

    //  **************************************************
    //  Unpack migrating data to remote partition
    //  **************************************************

    if ( iproc == myrank ) {
        return 1;                // skip local partition
    }

    // query process communicator to use
    ProcessCommunicatorBuff *pcbuff = pc.giveProcessCommunicatorBuff();
    ProcessCommDataStream pcDataStream(pcbuff);

    pcbuff->unpackInt(_type);
    // unpack dofman data
    while ( _type != LOADBALANCER_END_DATA ) {
        _etype = ( classType ) _type;
        pcbuff->unpackInt(_mode);
        switch ( _mode ) {
        case LoadBalancer :: DM_Remote:
            // receiving new local dofManager
            pcbuff->unpackInt(_globnum);
            /*
             * _newentry = false;
             * if ( ( dofman = dtm->giveDofManager(_globnum) ) == NULL ) {
             *  // data not available -> create a new one
             *  _newentry = true;
             *  dofman = CreateUsrDefDofManagerOfType(_etype, 0, d);
             * }
             */
            _newentry = true;
            dofman = CreateUsrDefDofManagerOfType(_etype, 0, d);

            dofman->setGlobalNumber(_globnum);
            // unpack dofman state (this is the local dofman, not available on remote)
            dofman->restoreContext(& pcDataStream, CM_Definition | CM_DefinitionGlobal | CM_State | CM_UnknownDictState);
            // unpack list of new partitions
            pcbuff->unpackIntArray(_partitions);
            dofman->setPartitionList(& _partitions);
            dofman->setParallelMode(DofManager_local);
            // add transaction if new entry allocated; otherwise existing one has been modified via returned dofman
            if ( _newentry ) {
                dtm->addTransaction(DomainTransactionManager :: DTT_ADD, DomainTransactionManager :: DCT_DofManager, _globnum, dofman);
            }

            //dmanMap[_globnum] = dofman;
            break;

        case LoadBalancer :: DM_Shared:
            // receiving new shared dofManager, that was local on sending partition
            // should be received only once (from partition where was local)
            pcbuff->unpackInt(_globnum);
            /*
             * _newentry = false;
             * if ( ( dofman = dtm->giveDofManager(_globnum) ) == NULL ) {
             *  // data not available -> mode should be SharedUpdate
             *  _newentry = true;
             *  dofman = CreateUsrDefDofManagerOfType(_etype, 0, d);
             * }
             */
            _newentry = true;
            dofman = CreateUsrDefDofManagerOfType(_etype, 0, d);


            dofman->setGlobalNumber(_globnum);
            // unpack dofman state (this is the local dofman, not available on remote)
            dofman->restoreContext(& pcDataStream, CM_Definition | CM_DefinitionGlobal | CM_State | CM_UnknownDictState);
            // unpack list of new partitions
            pcbuff->unpackIntArray(_partitions);
            dofman->setPartitionList(& _partitions);
            dofman->setParallelMode(DofManager_shared);
#ifdef __VERBOSE_PARALLEL
            fprintf(stderr, "[%d] received Shared new dofman [%d]\n", myrank, _globnum);
#endif
            // add transaction if new entry allocated; otherwise existing one has been modified via returned dofman
            if ( _newentry ) {
                dtm->addTransaction(DomainTransactionManager :: DTT_ADD, DomainTransactionManager :: DCT_DofManager, _globnum, dofman);
            }

            //dmanMap[_globnum] = dofman;
            break;

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

//.........这里部分代码省略.........
        for(size_t i = 0; i < numTracEl; i++) {
            std::vector<FloatArray> arcPos;
            double xiS = 0.0, xiE = 0.0;
            iBC.giveTractionElArcPos(i, xiS, xiE);
            arcPos.push_back( FloatArray{xiS} );
            arcPos.push_back( FloatArray{xiE} );

            arcPosArray.push_back(arcPos);
        }

        std :: stringstream strArcPos;
        strArcPos << "ArcPosGnuplotTime" << time << ".dat";
        std :: string nameArcPos = strArcPos.str();

        WritePointsToGnuplot(nameArcPos, arcPosArray);


        // Traction (normal, tangent)
        std::vector< std::vector<FloatArray> > nodeTractionNTArray;
        for(size_t i = 0; i < numTracEl; i++) {

            std::vector<FloatArray> tractions;
            FloatArray tS, tE;

            iBC.giveTraction(i, tS, tE, VM_Total, tStep);
            FloatArray n,t;
            iBC.giveTractionElNormal(i, n, t);


            double tSn = tS.dotProduct(n,2);
            double tSt = tS.dotProduct(t,2);
            tractions.push_back( {tSn ,tSt} );

            double tEn = tE.dotProduct(n,2);
            double tEt = tE.dotProduct(t,2);
            tractions.push_back( {tEn, tEt} );
            nodeTractionNTArray.push_back(tractions);
        }

        std :: stringstream strTractionsNT;
        strTractionsNT << "TractionsNormalTangentGnuplotTime" << time << ".dat";
        std :: string nameTractionsNT = strTractionsNT.str();

        WritePointsToGnuplot(nameTractionsNT, nodeTractionNTArray);



        // Boundary points and displacements
        IntArray boundaries, bNodes;
        iBC.giveBoundaries(boundaries);

        std::vector< std::vector<FloatArray> > bndNodes;

        for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {

            Element *e = iBC.giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
            int boundary = boundaries.at(pos * 2);

            e->giveInterpolation()->boundaryGiveNodes(bNodes, boundary);

            std::vector<FloatArray> bndSegNodes;

            // Add the start and end nodes of the segment
            DofManager *startNode   = e->giveDofManager( bNodes[0] );
            FloatArray xS    = *(startNode->giveCoordinates());

            Dof *dSu = startNode->giveDofWithID(D_u);
            double dU = dSu->giveUnknown(VM_Total, tStep);
            xS.push_back(dU);

            Dof *dSv = startNode->giveDofWithID(D_v);
            double dV = dSv->giveUnknown(VM_Total, tStep);
            xS.push_back(dV);

            bndSegNodes.push_back(xS);

            DofManager *endNode     = e->giveDofManager( bNodes[1] );
            FloatArray xE    = *(endNode->giveCoordinates());

            Dof *dEu = endNode->giveDofWithID(D_u);
            dU = dEu->giveUnknown(VM_Total, tStep);
            xE.push_back(dU);

            Dof *dEv = endNode->giveDofWithID(D_v);
            dV = dEv->giveUnknown(VM_Total, tStep);
            xE.push_back(dV);

            bndSegNodes.push_back(xE);

            bndNodes.push_back(bndSegNodes);
        }

        std :: stringstream strBndNodes;
        strBndNodes << "BndNodesGnuplotTime" << time << ".dat";
        std :: string nameBndNodes = strBndNodes.str();

        WritePointsToGnuplot(nameBndNodes, bndNodes);

    }
}

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
//.........这里部分代码省略.........

int
LoadBalancer :: packMigratingData(Domain *d, ProcessCommunicator &pc)
{
    int myrank = d->giveEngngModel()->giveRank();
    int iproc = pc.giveRank();
    int idofman, ndofman;
    classType dtype;
    DofManager *dofman;
    LoadBalancer :: DofManMode dmode;

    //  **************************************************
    //  Pack migrating data to remote partition
    //  **************************************************

    // pack dofManagers
    if ( iproc == myrank ) {
        return 1;                // skip local partition
    }

    // query process communicator to use
    ProcessCommunicatorBuff *pcbuff = pc.giveProcessCommunicatorBuff();
    ProcessCommDataStream pcDataStream(pcbuff);
    // loop over dofManagers
    ndofman = d->giveNumberOfDofManagers();
    for ( idofman = 1; idofman <= ndofman; idofman++ ) {
        dofman = d->giveDofManager(idofman);
        dmode = this->giveDofManState(idofman);
        dtype = dofman->giveClassID();
        // sync data to remote partition
        // if dofman already present on remote partition then there is no need to sync
        //if ((this->giveDofManPartitions(idofman)->findFirstIndexOf(iproc))) {
        if ( ( this->giveDofManPartitions(idofman)->findFirstIndexOf(iproc) ) &&
            ( !dofman->givePartitionList()->findFirstIndexOf(iproc) ) ) {
            pcbuff->packInt(dtype);
            pcbuff->packInt(dmode);
            pcbuff->packInt( dofman->giveGlobalNumber() );

            // pack dofman state (this is the local dofman, not available on remote)
            /* this is a potential performance leak, sending shared dofman to a partition,
             * in which is already shared does not require to send context (is already there)
             * here for simplicity it is always send */
            dofman->saveContext(& pcDataStream, CM_Definition | CM_DefinitionGlobal | CM_State | CM_UnknownDictState);
            // send list of new partitions
            pcbuff->packIntArray( * ( this->giveDofManPartitions(idofman) ) );
        }
    }

    // pack end-of-dofman-section record
    pcbuff->packInt(LOADBALANCER_END_DATA);

    int ielem, nelem = d->giveNumberOfElements(), nsend = 0;

    Element *elem;

    for ( ielem = 1; ielem <= nelem; ielem++ ) { // begin loop over elements
        elem = d->giveElement(ielem);
        if ( ( elem->giveParallelMode() == Element_local ) &&
            ( this->giveElementPartition(ielem) == iproc ) ) {
            // pack local element (node numbers shuld be global ones!!!)
            // pack type
            pcbuff->packInt( elem->giveClassID() );
            // nodal numbers shuld be packed as global !!
            elem->saveContext(& pcDataStream, CM_Definition | CM_DefinitionGlobal | CM_State);
            nsend++;
        }
    } // end loop over elements

    // pack end-of-element-record
    pcbuff->packInt(LOADBALANCER_END_DATA);

    OOFEM_LOG_RELEVANT("[%d] LoadBalancer:: sending %d migrating elements to %d\n", myrank, nsend, iproc);

    return 1;
}

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.
//.........这里部分代码省略.........

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

    // Initiates the total displacement to zero.
    if ( tStep->isTheFirstStep() ) {
        Domain *d = this->giveDomain(1);
        for ( int i = 1; i <= d->giveNumberOfDofManagers(); i++ ) {
            DofManager *dofman = d->giveDofManager(i);
            for ( int j = 1; j <= dofman->giveNumberOfDofs(); j++ ) {
                dofman->giveDof(j)->updateUnknownsDictionary(tStep, VM_Total_Old, 0.);
                dofman->giveDof(j)->updateUnknownsDictionary(tStep, VM_Total, 0.);
                // This is actually redundant now;
                //dofman->giveDof(j)->updateUnknownsDictionary(tStep, VM_Incremental, 0.);
            }
        }

	int nbc = d->giveNumberOfBoundaryConditions();
	for ( int ibc = 1; ibc <= nbc; ++ibc ) {
	  GeneralBoundaryCondition *bc = d->giveBc(ibc);
	  ActiveBoundaryCondition *abc;

	  if ( ( abc = dynamic_cast< ActiveBoundaryCondition * >( bc ) ) ) {
	    int ndman = abc->giveNumberOfInternalDofManagers();
	    for ( int i = 1; i <= ndman; i++ ) {
	      DofManager *dofman = abc->giveInternalDofManager(i);
	      for ( int j = 1; j <= dofman->giveNumberOfDofs(); j++ ) {
                dofman->giveDof(j)->updateUnknownsDictionary(tStep, VM_Total_Old, 0.);
                dofman->giveDof(j)->updateUnknownsDictionary(tStep, VM_Total, 0.);
                // This is actually redundant now;
                //dofman->giveDof(j)->updateUnknownsDictionary(tStep, VM_Incremental, 0.);
	      }
	    }
	  }
	}
    }

    // Apply dirichlet b.c's on total values
    Domain *d = this->giveDomain(1);
    for ( int i = 1; i <= d->giveNumberOfDofManagers(); i++ ) {
        DofManager *dofman = d->giveDofManager(i);
        for ( int j = 1; j <= dofman->giveNumberOfDofs(); j++ ) {
            Dof *d = dofman->giveDof(j);
            double tot = d->giveUnknown(VM_Total_Old, tStep);
            if ( d->hasBc(tStep) ) {
                tot += d->giveBcValue(VM_Incremental, tStep);
            }

            d->updateUnknownsDictionary(tStep, VM_Total, tot);
        }
    }


#ifdef VERBOSE
    OOFEM_LOG_RELEVANT( "Solving [step number %8d, time %15e]\n", tStep->giveNumber(), tStep->giveTargetTime() );
#endif

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

    if (neq == 0) { // Allows for fully prescribed/empty problems.
        return;
    }

    incrementOfDisplacementVector.resize(neq);
    incrementOfDisplacementVector.zero();

#ifdef VERBOSE
    OOFEM_LOG_INFO("Assembling load\n");
#endif
    // Assembling the element part of load vector
    internalLoadVector.resize(neq);
    internalLoadVector.zero();
    this->assembleVector( internalLoadVector, tStep, EID_MomentumBalance, InternalForcesVector,
                          VM_Total, EModelDefaultEquationNumbering(), this->giveDomain(1) );

    loadVector.resize(neq);
    loadVector.zero();
    this->assembleVector( loadVector, tStep, EID_MomentumBalance, ExternalForcesVector,
                          VM_Total, EModelDefaultEquationNumbering(), this->giveDomain(1) );

    loadVector.subtract(internalLoadVector);

#ifdef VERBOSE
    OOFEM_LOG_INFO("Assembling stiffness matrix\n");
#endif
    if ( stiffnessMatrix ) {
        delete stiffnessMatrix;
    }

    stiffnessMatrix = classFactory.createSparseMtrx(sparseMtrxType);
    if ( stiffnessMatrix == NULL ) {
        _error("solveYourselfAt: sparse matrix creation failed");
    }

    stiffnessMatrix->buildInternalStructure( this, 1, EID_MomentumBalance, EModelDefaultEquationNumbering() );
    stiffnessMatrix->zero();
    this->assemble( stiffnessMatrix, tStep, EID_MomentumBalance, StiffnessMatrix,
                   EModelDefaultEquationNumbering(), this->giveDomain(1) );

#ifdef VERBOSE
//.........这里部分代码省略.........

bool
NRSolver :: checkConvergence(FloatArray &RT, FloatArray &F, FloatArray &rhs,  FloatArray &ddX, FloatArray &X,
                             double RRT, const FloatArray &internalForcesEBENorm,
                             int nite, bool &errorOutOfRange, TimeStep *tNow)
{
    double forceErr, dispErr;
    FloatArray dg_forceErr, dg_dispErr, dg_totalLoadLevel, dg_totalDisp;
    bool answer;
    EModelDefaultEquationNumbering dn;
 #ifdef __PARALLEL_MODE
  #ifdef __PETSC_MODULE
    PetscContext *parallel_context = engngModel->givePetscContext(this->domain->giveNumber());
    Natural2LocalOrdering *n2l = parallel_context->giveN2Lmap();
  #endif
 #endif

    /*
     * The force errors are (if possible) evaluated as relative errors.
     * If the norm of applied load vector is zero (one may load by temperature, etc)
     * then the norm of reaction forces is used in relative norm evaluation.
     *
     * Note: This is done only when all dofs are included (nccdg = 0). Not implemented if
     * multiple convergence criteria are used.
     *
     */

    answer = true;
    errorOutOfRange = false;

    if ( internalForcesEBENorm.giveSize() > 1 ) { // Special treatment when just one norm is given; No grouping
        int nccdg = this->domain->giveMaxDofID();
        // Keeps tracks of which dof IDs are actually in use;
        IntArray idsInUse(nccdg);
        idsInUse.zero();
        // zero error norms per group
        dg_forceErr.resize(nccdg); dg_forceErr.zero();
        dg_dispErr.resize(nccdg); dg_dispErr.zero();
        dg_totalLoadLevel.resize(nccdg); dg_totalLoadLevel.zero();
        dg_totalDisp.resize(nccdg); dg_totalDisp.zero();
        // loop over dof managers
        int ndofman = domain->giveNumberOfDofManagers();
        for ( int idofman = 1; idofman <= ndofman; idofman++ ) {
            DofManager *dofman = domain->giveDofManager(idofman);
 #if ( defined ( __PARALLEL_MODE ) && defined ( __PETSC_MODULE ) )
            if ( !parallel_context->isLocal(dofman) ) {
                continue;
            }

 #endif

            // loop over individual dofs
            int ndof = dofman->giveNumberOfDofs();
            for ( int idof = 1; idof <= ndof; idof++ ) {
                Dof *dof = dofman->giveDof(idof);
                if ( !dof->isPrimaryDof() ) continue;
                int eq = dof->giveEquationNumber(dn);
                int dofid = dof->giveDofID();
                if ( !eq ) continue;
 
                dg_forceErr.at(dofid) += rhs.at(eq) * rhs.at(eq);
                dg_dispErr.at(dofid) += ddX.at(eq) * ddX.at(eq);
                dg_totalLoadLevel.at(dofid) += RT.at(eq) * RT.at(eq);
                dg_totalDisp.at(dofid) += X.at(eq) * X.at(eq);
                idsInUse.at(dofid) = 1;
            } // end loop over DOFs
        } // end loop over dof managers

        // loop over elements and their DOFs
        int nelem = domain->giveNumberOfElements();
        for ( int ielem = 1; ielem <= nelem; ielem++ ) {
            Element *elem = domain->giveElement(ielem);
 #ifdef __PARALLEL_MODE
            if ( elem->giveParallelMode() != Element_local ) {
                continue;
            }

 #endif
            // loop over element internal Dofs
            for ( int idofman = 1; idofman <= elem->giveNumberOfInternalDofManagers(); idofman++) {
                DofManager *dofman = elem->giveInternalDofManager(idofman);
                int ndof = dofman->giveNumberOfDofs();
                // loop over individual dofs
                for ( int idof = 1; idof <= ndof; idof++ ) {
                    Dof *dof = dofman->giveDof(idof);
                    if ( !dof->isPrimaryDof() ) continue;
                    int eq = dof->giveEquationNumber(dn);
                    int dofid = dof->giveDofID();
                    
                    if ( !eq ) continue;
 #if ( defined ( __PARALLEL_MODE ) && defined ( __PETSC_MODULE ) )
                    if ( engngModel->isParallel() && !n2l->giveNewEq(eq) ) continue;
 #endif
                    dg_forceErr.at(dofid) += rhs.at(eq) * rhs.at(eq);
                    dg_dispErr.at(dofid) += ddX.at(eq) * ddX.at(eq);
                    dg_totalLoadLevel.at(dofid) += RT.at(eq) * RT.at(eq);
                    dg_totalDisp.at(dofid) += X.at(eq) * X.at(eq);
                    idsInUse.at(dofid) = 1;
                } // end loop over DOFs
            } // end loop over element internal dofmans
        } // end loop over elements
//.........这里部分代码省略.........

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
ProblemCommunicator :: setUpCommunicationMapsForElementCut(EngngModel *pm,
                                                           bool excludeSelfCommFlag)
{
    Domain *domain = pm->giveDomain(1);
    int nnodes = domain->giveNumberOfDofManagers();
    int i, j, partition;

    if ( this->mode == ProblemCommMode__ELEMENT_CUT ) {
        /*
         * Initially, each partition knows for which nodes a receive
         * is needed (and can therefore compute easily the recv map),
         * but does not know for which nodes it should send data to which
         * partition. Hence, the communication setup is performed by
         * broadcasting "send request" lists of nodes for which
         * a partition expects to receive data (ie. of those nodes
         * which the partition uses, but does not own) to all
         * collaborating processes. The "send request" list are
         * converted into send maps.
         */

        // receive maps can be build locally,
        // but send maps should be assembled from broadcasted lists (containing
        // expected receive nodes) of remote partitions.

        // first build local receive map
        IntArray domainNodeRecvCount(size);
        const IntArray *partitionList;
        DofManager *dofMan;
        //Element    *element;
        int domainRecvListSize = 0, domainRecvListPos = 0;
        //int nelems;
        int result = 1;

        for ( i = 1; i <= nnodes; i++ ) {
            partitionList = domain->giveDofManager(i)->givePartitionList();
            if ( domain->giveDofManager(i)->giveParallelMode() == DofManager_remote ) {
                // size of partitionList should be 1 <== only ine master
                for ( j = 1; j <= partitionList->giveSize(); j++ ) {
                    if ( !( excludeSelfCommFlag && ( this->rank == partitionList->at(j) ) ) ) {
                        domainRecvListSize++;
                        domainNodeRecvCount.at(partitionList->at(j) + 1)++;
                    }
                }
            }
        }

        // build maps simultaneously
        IntArray pos(size);
        IntArray **maps = new IntArray * [ size ];
        for ( i = 0; i < size; i++ ) {
            maps [ i ] = new IntArray( domainNodeRecvCount.at(i + 1) );
        }

        // allocate also domain receive list to be broadcasted
        IntArray domainRecvList(domainRecvListSize);

        if ( domainRecvListSize ) {
            for ( i = 1; i <= nnodes; i++ ) {
                // test if node is remote DofMan
                dofMan = domain->giveDofManager(i);
                if ( dofMan->giveParallelMode() == DofManager_remote ) {
                    domainRecvList.at(++domainRecvListPos) = dofMan->giveGlobalNumber();

                    partitionList = domain->giveDofManager(i)->givePartitionList();
                    // size of partitionList should be 1 <== only ine master
                    for ( j = 1; j <= partitionList->giveSize(); j++ ) {
                        if ( !( excludeSelfCommFlag && ( this->rank == partitionList->at(j) ) ) ) {
                            partition = partitionList->at(j);
                            maps [ partition ]->at( ++pos.at(partition + 1) ) = i;
                        }
                    }
                }
            }
        }

        // set up process recv communicator maps
        for ( i = 0; i < size; i++ ) {
            this->setProcessCommunicatorToRecvArry(this->giveProcessCommunicator(i), * maps [ i ]);
            //this->giveDomainCommunicator(i)->setToRecvArry (this->engngModel, *maps[i]);
        }

        // delete local maps
        for ( i = 0; i < size; i++ ) {
            delete maps [ i ];
        }

        delete maps;

        // to assemble send maps, we must analyze broadcasted remote domain send lists
        // and we must also broadcast our send list.

#ifdef __VERBOSE_PARALLEL
        VERBOSEPARALLEL_PRINT("ProblemCommunicator::setUpCommunicationMaps", "Element-cut broadcasting started", rank);
#endif


        StaticCommunicationBuffer commBuff(MPI_COMM_WORLD);
        IntArray remoteDomainRecvList;
        IntArray toSendMap;
//.........这里部分代码省略.........

/////////////////////////////////////////////////
// Help functions
void GnuplotExportModule::outputReactionForces(TimeStep *tStep)
{
    // Add sum of reaction forces to arrays
    // Compute sum of reaction forces for each BC number
    Domain *domain = emodel->giveDomain(1);
    StructuralEngngModel *seMod = dynamic_cast<StructuralEngngModel* >(emodel);
    if(seMod == NULL) {
        OOFEM_ERROR("failed to cast to StructuralEngngModel.");
    }

    IntArray ielemDofMask;
    FloatArray reactions;
    IntArray dofManMap, dofidMap, eqnMap;

    // test if solution step output is active
    if ( !testTimeStepOutput(tStep) ) {
        return;
    }

    // map contains corresponding dofmanager and dofs numbers corresponding to prescribed equations
    // sorted according to dofmanger number and as a minor crit. according to dof number
    // this is necessary for extractor, since the sorted output is expected
    seMod->buildReactionTable(dofManMap, dofidMap, eqnMap, tStep, 1);

    // compute reaction forces
    seMod->computeReaction(reactions, tStep, 1);

    // Find highest index of prescribed dofs
    int maxIndPresDof = 0;
    for ( int i = 1; i <= dofManMap.giveSize(); i++ ) {
        maxIndPresDof = std::max(maxIndPresDof, dofidMap.at(i));
    }

    int numBC = domain->giveNumberOfBoundaryConditions();

    while ( mReactionForceHistory.size() < size_t(numBC) ) {
        std::vector<FloatArray> emptyArray;
        mReactionForceHistory.push_back( emptyArray );
    }

    maxIndPresDof = domain->giveNumberOfSpatialDimensions();

    while ( mDispHist.size() < size_t(numBC) ) {
        std::vector<double> emptyArray;
        mDispHist.push_back( emptyArray );
    }

    for(int bcInd = 0; bcInd < numBC; bcInd++) {
        FloatArray fR(maxIndPresDof), disp(numBC);
        fR.zero();


        for ( int i = 1; i <= dofManMap.giveSize(); i++ ) {
            DofManager *dMan = domain->giveDofManager( dofManMap.at(i) );
            Dof *dof = dMan->giveDofWithID( dofidMap.at(i) );

            if ( dof->giveBcId() == bcInd+1 ) {
                fR.at( dofidMap.at(i) ) += reactions.at( eqnMap.at(i) );

                // Slightly dirty
                BoundaryCondition *bc = dynamic_cast<BoundaryCondition*> (domain->giveBc(bcInd+1));
                if ( bc != NULL ) {
                    disp.at(bcInd+1) = bc->give(dof, VM_Total, tStep->giveTargetTime());
                }
                ///@todo This function should be using the primaryfield instead of asking BCs directly. / Mikael
            }
        }

        mDispHist[bcInd].push_back(disp.at(bcInd+1));
        mReactionForceHistory[bcInd].push_back(fR);



        // X
        FILE * pFileX;
        char fileNameX[100];
        sprintf(fileNameX, "ReactionForceGnuplotBC%dX.dat", bcInd+1);
        pFileX = fopen ( fileNameX , "wb" );

        fprintf(pFileX, "#u Fx\n");
        for ( size_t j = 0; j < mDispHist[bcInd].size(); j++ ) {
            fprintf(pFileX, "%e %e\n", mDispHist[bcInd][j], mReactionForceHistory[bcInd][j].at(1) );
        }

        fclose(pFileX);

        // Y
        FILE * pFileY;
        char fileNameY[100];
        sprintf(fileNameY, "ReactionForceGnuplotBC%dY.dat", bcInd+1);
        pFileY = fopen ( fileNameY , "wb" );

        fprintf(pFileY, "#u Fx\n");
        for ( size_t j = 0; j < mDispHist[bcInd].size(); j++ ) {
            if( mReactionForceHistory[bcInd][j].giveSize() >= 2 ) {
                fprintf(pFileY, "%e %e\n", mDispHist[bcInd][j], mReactionForceHistory[bcInd][j].at(2) );
            }
        }
//.........这里部分代码省略.........