C++ DofManager::giveCoordinates方法代码示例

double
DelaunayTriangle :: giveEdgeLength(int nodeA, int nodeB)
{
    DofManager *dmanA = domain->giveDofManager( giveNode(nodeA) );
    DofManager *dmanB = domain->giveDofManager( giveNode(nodeB) );

    return dmanA->giveCoordinates()->distance( dmanB->giveCoordinates() );
}

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

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

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

double TransportGradientPeriodic :: giveUnknown(double val, ValueModeType mode, TimeStep *tStep, ActiveDof *dof)
{
    DofManager *master = this->domain->giveDofManager(this->slavemap[dof->giveDofManager()->giveNumber()]);
    FloatArray dx, g;
    dx.beDifferenceOf(* dof->giveDofManager()->giveCoordinates(), * master->giveCoordinates());
    this->grad->giveUnknownVector(g, this->grad_ids, mode, tStep);
    return val + g.dotProduct(dx);
}

void TransportGradientPeriodic :: computeDofTransformation(ActiveDof *dof, FloatArray &masterContribs)
{
    DofManager *master = this->domain->giveDofManager(this->slavemap[dof->giveDofManager()->giveNumber()]);
    FloatArray *coords = dof->giveDofManager()->giveCoordinates();
    FloatArray *masterCoords = master->giveCoordinates();

    FloatArray dx;
    dx.beDifferenceOf(* coords, * masterCoords );

    masterContribs.resize(dx.giveSize() + 1);

    masterContribs.at(1) = 1.; // Master dof is always weight 1.0
    for ( int i = 1; i <= dx.giveSize(); ++i ) {
        masterContribs.at(i+1) = dx.at(i);
    }
}

void
SolutionbasedShapeFunction :: copyDofManagersToSurfaceData(modeStruct *mode, IntArray nodeList, bool isPlus, bool isMinus, bool isZeroBoundary)
{
    for ( int i = 1; i <= nodeList.giveSize(); i++ ) {
        FloatArray values;

        IntArray DofIDs;
        DofManager *dman = mode->myEngngModel->giveDomain(1)->giveDofManager(nodeList.at(i));

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

/* <<<<<<< HEAD
=======
        for ( int j = 1; j <= this->dofs.giveSize(); j++ ) {
            SurfaceDataStruct *surfaceData = new(SurfaceDataStruct);
            Dof *d = dman->giveDofWithID( dofs.at(j) );
>>>>>>> 147f565295394adef603dae296a820af5f28d9cd
*/
        // Check that current node contains current DofID
        for (int j=1; j<=this->dofs.giveSize(); j++) {
            for (Dof *d: *dman ){ //int k=1; k<= dman->giveNumberOfDofs(); k++ ) {

                //Dof *d = dman->dofArray.at(k);// giveDof(k);

                if (d->giveDofID() == this->dofs.at(j)) {
                    SurfaceDataStruct *surfaceData = new(SurfaceDataStruct);
                    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);
                }
            }
        }
    }
}

void
SolutionbasedShapeFunction :: splitBoundaryNodeIDs(modeStruct &mode, Element &e, IntArray &bnodes, IntArray &pList, IntArray &mList, IntArray &zList, FloatMatrix &nodeValues)
{
    pList.clear();
    mList.clear();
    zList.clear();

    for ( int j = 1; j <= bnodes.giveSize(); j++ ) {
        DofManager *dman = e.giveDofManager( bnodes.at(j) );

        bool isZero = false;
        bool isPlus = false;
        bool isMinus = false;

        whichBoundary(* dman->giveCoordinates(), isPlus, isMinus, isZero);

        if ( isZero ) {
            zList.insertSorted(j);
        } else if ( isPlus ) {
            pList.insertSorted(j);
        } else if ( isMinus ) {
            mList.insertSorted(j);
        }

        // Find global DofManager and fetch nodal values
        for ( size_t k = 0; k < mode.SurfaceData.size(); k++ ) {
            if ( mode.SurfaceData.at(k)->DofMan == dman ) {
                int IndexOfDofIDItem = 0;
                for ( int l = 1; l <= dofs.giveSize(); l++ ) {
                    if ( dofs.at(l) == mode.SurfaceData.at(k)->DofID ) {
                        IndexOfDofIDItem = l;
                        break;
                    }
                }
                nodeValues.at(IndexOfDofIDItem, j) = mode.SurfaceData.at(k)->value;
            }
        }
    }
}

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

double
UserDefDirichletBC :: give(Dof *dof, ValueModeType mode, TimeStep *stepN)
{
    double factor = this->giveLoadTimeFunction()->evaluate(stepN, mode);
    DofManager *dMan = dof->giveDofManager();


    /*
     * The Python function takes two input arguments:
     * 	1) An array with node coordinates
     * 	2) The dof id
     */
    int numArgs = 3;

    // Create array with node coordinates
    int dim = dMan->giveCoordinates()->giveSize();
    PyObject *pArgArray = PyList_New(dim);

    PyObject *pArgs = PyTuple_New(numArgs);

    for (int i = 0; i < dim; i++) {
        PyList_SET_ITEM(pArgArray, i, PyFloat_FromDouble( dMan->giveCoordinate(i+1) ));
    }
    // PyTuple_SetItem takes over responsibility for objects passed
    // to it -> no DECREF
    PyTuple_SetItem(pArgs, 0, pArgArray);


    // Dof number
    PyObject *pValDofNum = PyLong_FromLong(dof->giveDofID());
    PyTuple_SetItem(pArgs, 1, pValDofNum);


    // Time
    PyObject *pTargetTime = PyFloat_FromDouble( stepN->giveTargetTime() );
    PyTuple_SetItem(pArgs, 2, pTargetTime);

    // Value returned from the Python function
    PyObject *pRetVal;

    if ( PyCallable_Check(mpFunc) ) {
        pRetVal = PyObject_CallObject(mpFunc, pArgs);
    } else {
        OOFEM_ERROR("UserDefDirichletBC :: give: Python function is not callable.");
    }

    // Get return value
    double retVal = 0.0;
    if ( pRetVal != NULL ) {
        retVal = PyFloat_AsDouble(pRetVal);
    }
    else {
        OOFEM_ERROR("UserDefDirichletBC :: give: Failed to fetch Python return value.");
    }

    // Decrement reference count on pointers
    Py_DECREF(pArgs);
    Py_DECREF(pRetVal);

    return retVal*factor;
}

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