C++ FEInterpolation::evalN方法代码示例

void
AxisymElement :: computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int li, int ui)
//
// Returns the [ 6 x (nno*2) ] strain-displacement matrix {B} of the receiver,
// evaluated at gp.
// (epsilon_x,epsilon_y,...,Gamma_xy) = B . r
// r = ( u1,v1,u2,v2,u3,v3,u4,v4)
{

    FEInterpolation *interp = this->giveInterpolation();
     
    FloatArray N;
    interp->evalN( N, gp->giveNaturalCoordinates(), *this->giveCellGeometryWrapper() );
    double r = 0.0;
    for ( int i = 1; i <= this->giveNumberOfDofManagers(); i++ ) {
        double x = this->giveNode(i)->giveCoordinate(1);
        r += x * N.at(i);
    } 
    
    FloatMatrix dNdx;
    interp->evaldNdx( dNdx, gp->giveNaturalCoordinates(), *this->giveCellGeometryWrapper() );
    answer.resize(6, dNdx.giveNumberOfRows() * 2);
    answer.zero();

    for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
        answer.at(1, i * 2 - 1) = dNdx.at(i, 1);
        answer.at(2, i * 2 - 0) = dNdx.at(i, 2);
        answer.at(3, i * 2 - 1) = N.at(i) / r;
        answer.at(6, 2 * i - 1) = dNdx.at(i, 2);
        answer.at(6, 2 * i - 0) = dNdx.at(i, 1);
    }
  
}

bool EnrichmentItem :: tipIsTouchingEI(const TipInfo &iTipInfo)
{
    double tol = 1.0e-9;
    SpatialLocalizer *localizer = giveDomain()->giveSpatialLocalizer();

    Element *tipEl = localizer->giveElementContainingPoint(iTipInfo.mGlobalCoord);
    if ( tipEl != NULL ) {
        // Check if the candidate tip is located on the current crack
        FloatArray N;
        FloatArray locCoord;
        tipEl->computeLocalCoordinates(locCoord, iTipInfo.mGlobalCoord);
        FEInterpolation *interp = tipEl->giveInterpolation();
        interp->evalN( N, locCoord, FEIElementGeometryWrapper(tipEl) );

        double normalSignDist;
        evalLevelSetNormal( normalSignDist, iTipInfo.mGlobalCoord, N, tipEl->giveDofManArray() );

        double tangSignDist;
        evalLevelSetTangential( tangSignDist, iTipInfo.mGlobalCoord, N, tipEl->giveDofManArray() );

        if ( fabs(normalSignDist) < tol && tangSignDist > tol ) {
            return true;
        }
    }

    return false;
}

int
DofManValueField :: evaluateAt(FloatArray &answer, const FloatArray &coords, ValueModeType mode, TimeStep *tStep)
{
    int result = 0; // assume ok
    FloatArray lc, n;

    // request element containing target point
    Element *elem = this->domain->giveSpatialLocalizer()->giveElementContainingPoint(coords);
    if ( elem ) { // ok element containing target point found
        FEInterpolation *interp = elem->giveInterpolation();
        if ( interp ) {
            // map target point to element local coordinates
            if ( interp->global2local( lc, coords, FEIElementGeometryWrapper(elem) ) ) {
                // evaluate interpolation functions at target point
                interp->evalN( n, lc, FEIElementGeometryWrapper(elem) );
                // loop over element nodes
                for ( int i = 1; i <= n.giveSize(); i++ ) {
                    // multiply nodal value by value of corresponding shape function and add this to answer
                    answer.add( n.at(i), this->dmanvallist[elem->giveDofManagerNumber(i)-1] );
                }
            } else { // mapping from global to local coordinates failed
                result = 1; // failed
            }
        } else {  // element without interpolation
            result = 1; // failed
        }
    } else { // no element containing given point found
        result = 1; // failed
    }
    return result;
}

void PlaneStressStructuralElementEvaluator :: computeNMatrixAt(FloatMatrix &answer, GaussPoint *gp)
{
    FloatArray N;
    FEInterpolation *interp = gp->giveElement()->giveInterpolation();
    interp->evalN( N, gp->giveNaturalCoordinates(), FEIIGAElementGeometryWrapper( gp->giveElement(), gp->giveIntegrationRule()->giveKnotSpan() ) );
    answer.beNMatrixOf(N, 2);
}

void
L4Axisymm :: computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int li, int ui)
{
    // Returns the [ 6 x (nno*2) ] strain-displacement matrix {B} of the receiver,
    // evaluated at gp. Uses reduced integration.
    // (epsilon_x,epsilon_y,...,Gamma_xy) = B . r
    // r = ( u1,v1,u2,v2,u3,v3,u4,v4)

    FloatArray N, NRed, redCoord;
    if ( numberOfFiAndShGaussPoints == 1 ) { // Reduced integration
        redCoord  = {0.0, 0.0}; // eval in centroid
    } else {
        redCoord = gp->giveNaturalCoordinates();
    }


    FEInterpolation *interp = this->giveInterpolation();
        

    interp->evalN( N, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
    interp->evalN( NRed, redCoord, FEIElementGeometryWrapper(this) );
    
    // Evaluate radius at center
    double r = 0.0;
    for ( int i = 1; i <= this->giveNumberOfDofManagers(); i++ ) {
        double x = this->giveNode(i)->giveCoordinate(1);
        r += x * NRed.at(i);
    } 
    
    FloatMatrix dNdx, dNdxRed;
    interp->evaldNdx( dNdx, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
    interp->evaldNdx( dNdxRed, redCoord, FEIElementGeometryWrapper(this) );
    answer.resize(6, dNdx.giveNumberOfRows() * 2);
    answer.zero();

    for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
        answer.at(1, i * 2 - 1) = dNdx.at(i, 1);
        answer.at(2, i * 2 - 0) = dNdx.at(i, 2);
        answer.at(3, i * 2 - 1) = NRed.at(i) / r;
        answer.at(6, 2 * i - 1) = dNdxRed.at(i, 2);
        answer.at(6, 2 * i - 0) = dNdxRed.at(i, 1);
    }
}

/* 3D Space Elements */
void Space3dStructuralElementEvaluator :: computeNMatrixAt(FloatMatrix &answer, GaussPoint *gp)
{
    FloatArray N;
    Element *element = this->giveElement();
    FEInterpolation *interp = element->giveInterpolation();

    interp->evalN( N, gp->giveNaturalCoordinates(), FEIIGAElementGeometryWrapper( element, gp->giveIntegrationRule()->giveKnotSpan() ) );

    answer.beNMatrixOf(N, 3);
}

void
IntElLine1PhF :: computeNmatrixAt(GaussPoint *ip, FloatMatrix &answer)
{
    // Returns the modified N-matrix which multiplied with u give the spatial jump.

    FloatArray N;
    FEInterpolation *interp = this->giveInterpolation();
    interp->evalN( N, ip->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

    answer.resize(2, 8);
    answer.zero();
    answer.at(1, 1) = answer.at(2, 2) = -N.at(1);
    answer.at(1, 3) = answer.at(2, 4) = -N.at(2);

    answer.at(1, 5) = answer.at(2, 6) = N.at(1);
    answer.at(1, 7) = answer.at(2, 8) = N.at(2);
}

bool Inclusion :: isMaterialModified(GaussPoint &iGP, Element &iEl, CrossSection * &opCS) const
{
    // Check if the point is located inside the inclusion

    FloatArray N;
    FEInterpolation *interp = iEl.giveInterpolation();
    interp->evalN( N, * iGP.giveNaturalCoordinates(), FEIElementGeometryWrapper(& iEl) );

    const IntArray &elNodes = iEl.giveDofManArray();

    double levelSetGP = 0.0;
    interpLevelSet(levelSetGP, N, elNodes);

    if ( levelSetGP < 0.0 ) {
        opCS = mpCrossSection;
        return true;
    }

    return false;
}

int
SmoothedNodalInternalVariableField :: evaluateAt(FloatArray &answer, FloatArray &coords, ValueModeType mode, TimeStep *tStep)
{
    int result = 0; // assume ok
    FloatArray lc, n;
    const FloatArray *nodalValue;

    // use whole domain recovery
    // create a new set containing all elements
    Set elemSet(0, this->domain);
    elemSet.addAllElements();
    this->smoother->recoverValues(elemSet, istType, tStep);
    // request element containing target point
    Element *elem = this->domain->giveSpatialLocalizer()->giveElementContainingPoint(coords);
    if ( elem ) { // ok element containing target point found
        FEInterpolation *interp = elem->giveInterpolation();
        if ( interp ) {
            // map target point to element local coordinates
            if ( interp->global2local( lc, coords, FEIElementGeometryWrapper(elem) ) ) {
                // evaluate interpolation functions at target point
                interp->evalN( n, lc, FEIElementGeometryWrapper(elem) );
                // loop over element nodes
                for ( int i = 1; i <= n.giveSize(); i++ ) {
                    // request nodal value
                    this->smoother->giveNodalVector( nodalValue, elem->giveDofManagerNumber(i) );
                    // multiply nodal value by value of corresponding shape function and add this to answer
                    answer.add(n.at(i), * nodalValue);
                }
            } else { // mapping from global to local coordinates failed
                result = 1; // failed
            }
        } else {  // element without interpolation
            result = 1; // failed
        }
    } else { // no element containing given point found
        result = 1; // failed
    }
    return result;
}

void
ZZErrorEstimatorInterface :: ZZErrorEstimatorI_computeElementContributions(double &eNorm, double &sNorm,
                                                                           ZZErrorEstimator :: NormType norm,
                                                                           InternalStateType type,
                                                                           TimeStep *tStep)
{
    int nDofMans;
    FEInterpolation *interpol = element->giveInterpolation();
    const FloatArray *recoveredStress;
    FloatArray sig, lsig, diff, ldiff, n;
    FloatMatrix nodalRecoveredStreses;

    nDofMans = element->giveNumberOfDofManagers();
    // assemble nodal recovered stresses
    for ( int i = 1; i <= element->giveNumberOfNodes(); i++ ) {
        element->giveDomain()->giveSmoother()->giveNodalVector( recoveredStress,
                                                            element->giveDofManager(i)->giveNumber() );
        if ( i == 1 ) {
            nodalRecoveredStreses.resize( nDofMans, recoveredStress->giveSize() );
        }
        for ( int j = 1; j <= recoveredStress->giveSize(); j++ ) {
            nodalRecoveredStreses.at(i, j) = recoveredStress->at(j);
        }
    }
    /* Note: The recovered stresses should be in global coordinate system. This is important for shells, for example, to make
     * sure that forces and moments in the same directions are averaged. For elements where local and global coordina systems
     * are the same this does not matter.
     */

    eNorm = sNorm = 0.0;

    // compute  the e-norm and s-norm
    if ( norm == ZZErrorEstimator :: L2Norm ) {
        for ( GaussPoint *gp: *this->ZZErrorEstimatorI_giveIntegrationRule() ) {
            double dV = element->computeVolumeAround(gp);
            interpol->evalN( n, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(element) );

            diff.beTProductOf(nodalRecoveredStreses, n);

            element->giveIPValue(sig, gp, type, tStep);
            /* the internal stress difference is in global coordinate system */
            diff.subtract(sig);

            eNorm += diff.computeSquaredNorm() * dV;
            sNorm += sig.computeSquaredNorm() * dV;
        }
    } else if ( norm == ZZErrorEstimator :: EnergyNorm ) {
        FloatArray help, ldiff_reduced, lsig_reduced;
        FloatMatrix D, DInv;
        StructuralElement *selem = static_cast< StructuralElement * >(element);

        for ( GaussPoint *gp: *this->ZZErrorEstimatorI_giveIntegrationRule() ) {
            double dV = element->computeVolumeAround(gp);
            interpol->evalN( n, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(element) );
            selem->computeConstitutiveMatrixAt(D, TangentStiffness, gp, tStep);
            DInv.beInverseOf(D);

            diff.beTProductOf(nodalRecoveredStreses, n);

            element->giveIPValue(sig, gp, type, tStep); // returns full value now
            diff.subtract(sig);
            /* the internal stress difference is in global coordinate system */
            /* needs to be transformed into local system to compute associated energy */
            this->ZZErrorEstimatorI_computeLocalStress(ldiff, diff);
            StructuralMaterial :: giveReducedSymVectorForm( ldiff_reduced, ldiff, gp->giveMaterialMode() );

            help.beProductOf(DInv, ldiff_reduced);
            eNorm += ldiff_reduced.dotProduct(help) * dV;
            this->ZZErrorEstimatorI_computeLocalStress(lsig, sig);
            StructuralMaterial :: giveReducedSymVectorForm( lsig_reduced, lsig, gp->giveMaterialMode() );
            help.beProductOf(DInv, lsig_reduced);
            sNorm += lsig_reduced.dotProduct(help) * dV;
        }
    } else {
        OOFEM_ERROR("unsupported norm type");
    }

    eNorm = sqrt(eNorm);
    sNorm = sqrt(sNorm);
}

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

void
PrescribedMean :: assemble(SparseMtrx &answer, TimeStep *tStep, CharType type,
                           const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s)
{

    if ( type != TangentStiffnessMatrix && type != StiffnessMatrix ) {
        return;
    }

    computeDomainSize();

    IntArray c_loc, r_loc;
    lambdaDman->giveLocationArray(lambdaIDs, r_loc, r_s);
    lambdaDman->giveLocationArray(lambdaIDs, c_loc, c_s);

    for ( int i=1; i<=elements.giveSize(); i++ ) {
        int elementID = elements.at(i);
        Element *thisElement = this->giveDomain()->giveElement(elementID);
        FEInterpolation *interpolator = thisElement->giveInterpolation(DofIDItem(dofid));

        IntegrationRule *iRule = (elementEdges) ? (interpolator->giveBoundaryIntegrationRule(3, sides.at(i))) :
                                 (interpolator->giveIntegrationRule(3));

        for ( GaussPoint * gp: * iRule ) {
            FloatArray lcoords = gp->giveNaturalCoordinates();
            FloatArray N; //, a;
            FloatMatrix temp, tempT;
            double detJ = 0.0;
            IntArray boundaryNodes, dofids= {(DofIDItem) this->dofid}, r_Sideloc, c_Sideloc;

            if (elementEdges) {
                // Compute boundary integral
                interpolator->boundaryGiveNodes( boundaryNodes, sides.at(i) );
                interpolator->boundaryEvalN(N, sides.at(i), lcoords, FEIElementGeometryWrapper(thisElement));
                detJ = fabs ( interpolator->boundaryGiveTransformationJacobian(sides.at(i), lcoords, FEIElementGeometryWrapper(thisElement)) );
                // Retrieve locations for dofs on boundary
                thisElement->giveBoundaryLocationArray(r_Sideloc, boundaryNodes, dofids, r_s);
                thisElement->giveBoundaryLocationArray(c_Sideloc, boundaryNodes, dofids, c_s);
            } else {
                interpolator->evalN(N, lcoords, FEIElementGeometryWrapper(thisElement));
                detJ = fabs ( interpolator->giveTransformationJacobian(lcoords, FEIElementGeometryWrapper(thisElement) ) );
                IntArray DofIDStemp, rloc, cloc;

                thisElement->giveLocationArray(rloc, r_s, &DofIDStemp);
                thisElement->giveLocationArray(cloc, c_s, &DofIDStemp);

                r_Sideloc.clear();
                c_Sideloc.clear();
                for (int j=1; j<=DofIDStemp.giveSize(); j++) {
                    if (DofIDStemp.at(j)==dofids.at(1)) {
                        r_Sideloc.followedBy({rloc.at(j)});
                        c_Sideloc.followedBy({cloc.at(j)});
                    }
                }
            }

            // delta p part:
            temp = N*detJ*gp->giveWeight()*(1.0/domainSize);
            tempT.beTranspositionOf(temp);

            answer.assemble(r_Sideloc, c_loc, temp);
            answer.assemble(r_loc, c_Sideloc, tempT);
        }

        delete iRule;

    }

}

void
PrescribedMean :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep,
        CharType type, ValueModeType mode,
        const UnknownNumberingScheme &s, FloatArray *eNorm)
{
    computeDomainSize();

    // Fetch unknowns of this boundary condition
    IntArray lambdaLoc;
    FloatArray lambda;
    lambdaDman->giveUnknownVector(lambda, lambdaIDs, mode, tStep);
    lambdaDman->giveLocationArray(lambdaIDs, lambdaLoc, s);

    for ( int i=1; i<=elements.giveSize(); i++ ) {
        int elementID = elements.at(i);
        Element *thisElement = this->giveDomain()->giveElement(elementID);
        FEInterpolation *interpolator = thisElement->giveInterpolation(DofIDItem(dofid));

        IntegrationRule *iRule = (elementEdges) ? (interpolator->giveBoundaryIntegrationRule(3, sides.at(i))) :
                                 (interpolator->giveIntegrationRule(3));

        for ( GaussPoint * gp: * iRule ) {
            FloatArray lcoords = gp->giveNaturalCoordinates();
            FloatArray a, N, pressureEqns, lambdaEqns;
            IntArray boundaryNodes, dofids= {(DofIDItem) this->dofid}, locationArray;
            double detJ=0.0;

            if (elementEdges) {
                // Compute integral
                interpolator->boundaryGiveNodes( boundaryNodes, sides.at(i) );
                thisElement->computeBoundaryVectorOf(boundaryNodes, dofids, VM_Total, tStep, a);
                interpolator->boundaryEvalN(N, sides.at(i), lcoords, FEIElementGeometryWrapper(thisElement));
                detJ = fabs ( interpolator->boundaryGiveTransformationJacobian(sides.at(i), lcoords, FEIElementGeometryWrapper(thisElement)) );

                // Retrieve locations for dofs with dofids
                thisElement->giveBoundaryLocationArray(locationArray, boundaryNodes, dofids, s);
            } else {
                thisElement->computeVectorOf(dofids, VM_Total, tStep, a);
                interpolator->evalN(N, lcoords, FEIElementGeometryWrapper(thisElement));
                detJ = fabs ( interpolator->giveTransformationJacobian(lcoords, FEIElementGeometryWrapper(thisElement)));

                IntArray DofIDStemp, loc;

                thisElement->giveLocationArray(loc, s, &DofIDStemp);

                locationArray.clear();
                for (int j=1; j<=DofIDStemp.giveSize(); j++) {
                    if (DofIDStemp.at(j)==dofids.at(1)) {
                        locationArray.followedBy({loc.at(j)});
                    }
                }
            }

            // delta p part:
            pressureEqns = N*detJ*gp->giveWeight()*lambda.at(1)*(1.0/domainSize);

            // delta lambda part
            lambdaEqns.resize(1);
            lambdaEqns.at(1) = N.dotProduct(a);
            lambdaEqns.times(detJ*gp->giveWeight()*1.0/domainSize);
            lambdaEqns.at(1) = lambdaEqns.at(1);


            // delta p part
            answer.assemble(pressureEqns, locationArray);

            // delta lambda part
            answer.assemble(lambdaEqns, lambdaLoc);
        }
        delete iRule;
    }

}

void HangingNode :: postInitialize()
{
    Node :: postInitialize();

    Element *e;
    FEInterpolation *fei;
    FloatArray lcoords, masterContribution;

#ifdef __OOFEG
    if ( initialized ) {
        return;
    }
    initialized = true;
#endif

    // First check element and interpolation
    if ( masterElement == -1 ) { // Then we find it by taking the closest (probably containing element)
        FloatArray closest;
        SpatialLocalizer *sp = this->domain->giveSpatialLocalizer();
        sp->init();
        // Closest point or containing point? It should be contained, but with numerical errors it might be slightly outside
        // so the closest point is more robust.
        if ( !( e = sp->giveElementClosestToPoint(lcoords, closest, coordinates, this->masterRegion) ) ) {
            OOFEM_ERROR("Couldn't find closest element (automatically).");
        }
        this->masterElement = e->giveNumber();
    } else if ( !( e = this->giveDomain()->giveElement(this->masterElement) ) ) {
        OOFEM_ERROR("Requested element %d doesn't exist.", this->masterElement);
    }
    if ( !( fei = e->giveInterpolation() ) ) {
        OOFEM_ERROR("Requested element %d doesn't have a interpolator.", this->masterElement);
    }

    if ( lcoords.giveSize() == 0 ) { // we don't need to do this again if the spatial localizer was used.
        fei->global2local( lcoords, coordinates, FEIElementGeometryWrapper(e) );
    }

    // Initialize slave dofs (inside check of consistency of receiver and master dof)
    const IntArray &masterNodes = e->giveDofManArray();
    for ( Dof *dof: *this ) {
        SlaveDof *sdof = dynamic_cast< SlaveDof * >(dof);
        if ( sdof ) {
            DofIDItem id = sdof->giveDofID();
            fei = e->giveInterpolation(id);
            if ( !fei ) {
                OOFEM_ERROR("Requested interpolation for dof id %d doesn't exist in element %d.",
                             id, this->masterElement);
            }
#if 0 // This won't work (yet), as it requires some more general FEI classes, or something similar.
            if ( fei->hasMultiField() ) {
                FloatMatrix multiContribution;
                IntArray masterDofIDs, masterNodesDup, dofids;
                fei->evalMultiN(multiContribution, dofids, lcoords, FEIElementGeometryWrapper(e), 0.0);
                masterContribution.flatten(multiContribution);
                masterDofIDs.clear();
                for ( int i = 0; i <= multiContribution.giveNumberOfColumns(); ++i ) {
                    masterDofIDs.followedBy(dofids);
                    masterNodesDup.followedBy(masterNodes);
                }
                sdof->initialize(masterNodesDup, & masterDofIDs, masterContribution);
            } else { }
#else
            // Note: There can be more masterNodes than masterContributions, since all the
            // FEI classes are based on that the first nodes correspond to the simpler/linear interpolation.
            // If this assumption is changed in FEIElementGeometryWrapper + friends,
            // masterNode will also need to be modified for each dof accordingly.
            fei->evalN( masterContribution, lcoords, FEIElementGeometryWrapper(e) );
            sdof->initialize(masterNodes, IntArray(), masterContribution);
#endif
        }
    }
}

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