C++ FloatMatrix类代码示例

void
SolutionbasedShapeFunction :: computeCorrectionFactors(modeStruct &myMode, IntArray *Dofs, double *am, double *ap)
{
    /*
     * *Compute c0, cp, cm, Bp, Bm, Ap and Am
     */

    double A0p = 0.0, App = 0.0, A0m = 0.0, Amm = 0.0, Bp = 0.0, Bm = 0.0, c0 = 0.0, cp = 0.0, cm = 0.0;

    EngngModel *m = myMode.myEngngModel;
    Set *mySet = m->giveDomain(1)->giveSet(externalSet);

    IntArray BoundaryList = mySet->giveBoundaryList();

    for ( int i = 0; i < BoundaryList.giveSize() / 2; i++ ) {
        int ElementID = BoundaryList(2 * i);
        int Boundary = BoundaryList(2 * i + 1);

        Element *thisElement = m->giveDomain(1)->giveElement(ElementID);
        FEInterpolation *geoInterpolation = thisElement->giveInterpolation();
        IntArray bnodes, zNodes, pNodes, mNodes;
        FloatMatrix nodeValues;

        geoInterpolation->boundaryGiveNodes(bnodes, Boundary);

        nodeValues.resize( this->dofs.giveSize(), bnodes.giveSize() );
        nodeValues.zero();

        // Change to global ID for bnodes and identify the intersection of bnodes and the zero boundary
        splitBoundaryNodeIDs(myMode, * thisElement, bnodes, pNodes, mNodes, zNodes, nodeValues);

        std :: unique_ptr< IntegrationRule >iRule(geoInterpolation->giveBoundaryIntegrationRule(order, Boundary));

        for ( GaussPoint *gp: *iRule ) {
            FloatArray *lcoords = gp->giveCoordinates();
            FloatArray gcoords, normal, N;
            FloatArray Phi;

            double detJ = fabs( geoInterpolation->boundaryGiveTransformationJacobian( Boundary, * lcoords, FEIElementGeometryWrapper(thisElement) ) ) * gp->giveWeight();

            geoInterpolation->boundaryEvalNormal( normal, Boundary, * lcoords, FEIElementGeometryWrapper(thisElement) );
            geoInterpolation->boundaryEvalN( N, Boundary, * lcoords, FEIElementGeometryWrapper(thisElement) );
            geoInterpolation->boundaryLocal2Global( gcoords, Boundary, * lcoords, FEIElementGeometryWrapper(thisElement) );

            FloatArray pPhi, mPhi, zPhi;
            pPhi.resize( Dofs->giveSize() );
            pPhi.zero();
            mPhi.resize( Dofs->giveSize() );
            mPhi.zero();
            zPhi.resize( Dofs->giveSize() );
            zPhi.zero();

            // Build phi (analytical averaging, not projected onto the mesh)
            computeBaseFunctionValueAt(Phi, gcoords, * Dofs, * myMode.myEngngModel);

            // Build zPhi for this DofID
            for ( int l = 1; l <= zNodes.giveSize(); l++ ) {
                int nodeID = zNodes.at(l);
                for ( int m = 1; m <= this->dofs.giveSize(); m++ ) {
                    zPhi.at(m) = zPhi.at(m) + N.at(nodeID) * nodeValues.at(m, nodeID);
                }
            }


            // Build pPhi for this DofID
            for ( int l = 1; l <= pNodes.giveSize(); l++ ) {
                int nodeID = pNodes.at(l);
                for ( int m = 1; m <= this->dofs.giveSize(); m++ ) {
                    pPhi.at(m) = pPhi.at(m) + N.at(nodeID) * nodeValues.at(m, nodeID);
                }
            }

            // Build mPhi for this DofID
            for ( int l = 1; l <= mNodes.giveSize(); l++ ) {
                int nodeID = mNodes.at(l);
                for ( int m = 1; m <= this->dofs.giveSize(); m++ ) {
                    mPhi.at(m) = mPhi.at(m) + N.at(nodeID) * nodeValues.at(m, nodeID);
                }
            }

            c0 = c0 + zPhi.dotProduct(normal, 3) * detJ;
            cp = cp + pPhi.dotProduct(normal, 3) * detJ;
            cm = cm + mPhi.dotProduct(normal, 3) * detJ;

            App = App + pPhi.dotProduct(pPhi, 3) * detJ;
            Amm = Amm + mPhi.dotProduct(mPhi, 3) * detJ;
            A0p = A0p + zPhi.dotProduct(pPhi, 3) * detJ;
            A0m = A0m + zPhi.dotProduct(mPhi, 3) * detJ;

            Bp = Bp + Phi.dotProduct(pPhi, 3) * detJ;
            Bm = Bm + Phi.dotProduct(mPhi, 3) * detJ;
        }
    }

    * am = -( A0m * cp * cp - Bm * cp * cp - A0p * cm * cp + App * c0 * cm + Bp * cm * cp ) / ( App * cm * cm + Amm * cp * cp );
    * ap = -( A0p * cm * cm - Bp * cm * cm - A0m * cm * cp + Amm * c0 * cp + Bm * cm * cp ) / ( App * cm * cm + Amm * cp * cp );
}

void
IntElLine2 :: computeNmatrixAt(GaussPoint *ip, FloatMatrix &answer)
{

    // Returns the modified N-matrix which multiplied with u give the spatial jump.
    FloatArray N;
    answer.resize(2, 12);
    answer.zero();

    if(linear) {
		interpLin.evalN( N, ip->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

		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(3);

		answer.at(1, 7) = answer.at(2, 8) = N.at(1);
		answer.at(1, 9) = answer.at(2, 10) = N.at(2);
//		answer.at(1, 11) = answer.at(2, 12) = N.at(3);
    }
    else {
		interp.evalN( N, ip->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

		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(3);

		answer.at(1, 7) = answer.at(2, 8) = N.at(1);
		answer.at(1, 9) = answer.at(2, 10) = N.at(2);
		answer.at(1, 11) = answer.at(2, 12) = N.at(3);
    }
}

void LineSurfaceTension :: computeTangent(FloatMatrix &answer, TimeStep *tStep)
{
#if 1
    ///@todo Not sure if it's a good idea to use this tangent.
    answer.resize(4,4);
    answer.zero();
#else
    domainType dt = this->giveDomain()->giveDomainType();
    int ndofs = this->computeNumberOfDofs(EID_MomentumBalance);
    Node *node1, *node2;
    double x1, x2, y1, y2, dx, dy, vx, vy, length, width, gamma_s;

    gamma_s = this->giveMaterial()->give('g',NULL);

    node1 = giveNode(1);
    node2 = giveNode(2);

    x1 = node1->giveCoordinate(1);
    x2 = node2->giveCoordinate(1);

    y1 = node1->giveCoordinate(2);
    y2 = node2->giveCoordinate(2);

    dx = x2-x1;
    dy = y2-y1;

    length = sqrt(dx*dx + dy*dy);

    vx = dx/length;
    vy = dy/length;

    FloatArray Ah(4);
    Ah.at(1) = -vx;
    Ah.at(2) = -vy;
    Ah.at(3) = vx;
    Ah.at(4) = vy;

    FloatMatrix NpTNp(4,4);
    NpTNp.zero();
    NpTNp.at(1,1) = 1;
    NpTNp.at(2,2) = 1;
    NpTNp.at(3,3) = 1;
    NpTNp.at(4,4) = 1;
    NpTNp.at(1,3) = -1;
    NpTNp.at(2,4) = -1;
    NpTNp.at(3,1) = -1;
    NpTNp.at(4,2) = -1;

    answer.resize(ndofs,ndofs);
    answer.zero();
    if (dt == _3dAxisymmMode) {
        OOFEM_WARNING("Not tested");
        FloatArray Bh(4);
        Bh.zero();
        Bh.at(1) = 1;
        Bh.at(3) = 1;

        // It was simpler to write this in index notation.
        // Also using 0-based, to reduce typing
        double rJinv = (x1+x2)/length;
        answer.zero();
        for (int i = 0; i < 4; i++) for (int j = 0; j < 4; j++) {
            answer(i,j) = M_PI*(Ah(i)*Bh(j) + Bh(i)*Ah(j)+ rJinv*(NpTNp(i,j) - Ah(i)*Ah(j)));
        }
    }
    else {
        width = 1;
        answer.beDyadicProductOf(Ah,Ah);
        answer.add(NpTNp);
        answer.times(width/length);
    }

    answer.times(gamma_s);
#endif
}

void SUPGTangentAssembler :: matrixFromElement(FloatMatrix &answer, Element &el, TimeStep *tStep) const
{
    SUPGElement *element = static_cast< SUPGElement * >( &el );
    IntArray vloc, ploc;
    FloatMatrix h;
    int size = element->computeNumberOfDofs();
    element->giveLocalVelocityDofMap(vloc);
    element->giveLocalPressureDofMap(ploc);
    answer.resize(size, size);
    answer.zero();

    element->computeAccelerationTerm_MB(h, tStep);
    answer.assemble(h, vloc);
    element->computeAdvectionDerivativeTerm_MB(h, tStep);
    h.times( alpha * tStep->giveTimeIncrement() );
    answer.assemble(h, vloc);
    element->computeDiffusionDerivativeTerm_MB(h, rmode, tStep);

    h.times( alpha * tStep->giveTimeIncrement() );
    answer.assemble(h, vloc);
    element->computePressureTerm_MB(h, tStep);
    answer.assemble(h, vloc, ploc);
    element->computeLSICStabilizationTerm_MB(h, tStep);
    h.times( alpha * tStep->giveTimeIncrement() * lscale / ( dscale * uscale * uscale ) );
    answer.assemble(h, vloc);
    element->computeBCLhsTerm_MB(h, tStep);
    if ( h.isNotEmpty() ) {
        h.times( alpha * tStep->giveTimeIncrement() );
        answer.assemble(h, vloc);
    }

    element->computeBCLhsPressureTerm_MB(h, tStep);
    if ( h.isNotEmpty() ) {
        answer.assemble(h, vloc, ploc);
    }

    // conservation eq part
    element->computeLinearAdvectionTerm_MC(h, tStep);
    h.times( alpha * tStep->giveTimeIncrement() * 1.0 / ( dscale * uscale ) );
    answer.assemble(h, ploc, vloc);
    element->computeAdvectionDerivativeTerm_MC(h, tStep);
    h.times( alpha * tStep->giveTimeIncrement() );
    answer.assemble(h, ploc, vloc);
    element->computeAccelerationTerm_MC(h, tStep);
    answer.assemble(h, ploc, vloc);
    element->computeBCLhsPressureTerm_MC(h, tStep);
    if ( h.isNotEmpty() ) {
        answer.assemble(h, ploc, vloc);
    }

    element->computeDiffusionDerivativeTerm_MC(h, tStep);
    h.times( alpha * tStep->giveTimeIncrement() );
    answer.assemble(h, ploc, vloc);
    element->computePressureTerm_MC(h, tStep);
    answer.assemble(h, ploc);
}

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
FiberedCrossSection :: give3dBeamStiffMtrx(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
//
// General strain fiber vector has one of the following forms:
// 1) strainVector3d {eps_x,eps_y,eps_z,gamma_yz,gamma_zx,gamma_xy}
//
// returned strain or stress vector has the form:
// 2) strainVectorShell {eps_x, gamma_xz, gamma_xy, \der{phi_x}{x}, kappa_y, kappa_z}
//
{
    FloatMatrix fiberMatrix;
    GaussPoint *fiberGp;
    double fiberThick, fiberWidth, fiberZCoord, fiberYCoord;
    double fiberZCoord2, fiberYCoord2, Ip = 0.0, A = 0.0, Ik, G = 0.0;

    // if (form != ReducedForm) error ("give3dShellMaterialStiffness : full form unsupported");

    answer.resize(6, 6);
    answer.zero();
    // perform integration over layers

    for ( int i = 1; i <= numberOfFibers; i++ ) {
        fiberGp = giveSlaveGaussPoint(gp, i - 1);
        this->giveFiberMaterialStiffnessMatrix(fiberMatrix, rMode, fiberGp, tStep);
        //
        // resolve current layer z-coordinate
        //
        fiberThick  = this->fiberThicks.at(i);
        fiberWidth  = this->fiberWidths.at(i);
        fiberZCoord = fiberZcoords.at(i);
        fiberYCoord = fiberYcoords.at(i);
        fiberYCoord2 = fiberYCoord * fiberYCoord;
        fiberZCoord2 = fiberZCoord * fiberZCoord;
        //
        // perform integration
        //
        // 1) membrane terms N, Qz, Qy
        answer.at(1, 1) += fiberMatrix.at(1, 1) * fiberWidth * fiberThick;

        answer.at(2, 2) += fiberMatrix.at(2, 2) * fiberWidth * fiberThick;

        answer.at(3, 3) += fiberMatrix.at(3, 3) * fiberWidth * fiberThick;

        // 2) bending terms mx, my, mz

        Ip             += fiberWidth * fiberThick * fiberZCoord2 + fiberWidth * fiberThick * fiberYCoord2;
        A              += fiberWidth * fiberThick;
        G               = fiberMatrix.at(2, 2) * fiberWidth * fiberThick;

        answer.at(5, 5) += fiberMatrix.at(1, 1) * fiberWidth * fiberThick * fiberZCoord2;
        answer.at(6, 6) += fiberMatrix.at(1, 1) * fiberWidth * fiberThick * fiberYCoord2;
    }

    ///@todo This must be wrong, it will use the last evaluated G (from the last fiber), outside the loop. FIXME!
    G /= A;
    Ik = A * A * A * A / ( 40.0 * Ip );
    answer.at(4, 4) = G * Ik;
}

void
LSpaceBB :: computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int li, int ui)
// Returns the [6x24] strain-displacement matrix {B} of the receiver, eva-
// luated at gp.
// B matrix  -  6 rows : epsilon-X, epsilon-Y, epsilon-Z, gamma-YZ, gamma-ZX, gamma-XY  :
{
    int i;
    FloatMatrix dnx, dnx0;
    FloatArray coord(3);

    answer.resize(6, 24);
    answer.zero();
    coord.zero();


    LSpace :: interpolation.evaldNdx( dnx, * gp->giveCoordinates(), FEIElementGeometryWrapper(this) );
    LSpace :: interpolation.evaldNdx( dnx0, coord, FEIElementGeometryWrapper(this) );

    // deviatoric part fully integrated, volumetric part in one point
    // here we follow BBar approach
    //
    // construct BB = B(gp) + Pv [B(0)-B(gp)]
    // where Pv is volumetric projection mtrx
    // B(gp) is original geometrical matrix evalueated at gp
    // B(0)  is geometrical matrix evalueated at centroid
    //
    // assemble Pv [B(0)-B(gp)]
    for ( i = 1; i <= 8; i++ ) {
        answer.at(1, 3 * i - 2) = answer.at(2, 3 * i - 2) = answer.at(3, 3 * i - 2) = ( dnx0.at(i, 1) - dnx.at(i, 1) ) / 3.0;
        answer.at(1, 3 * i - 1) = answer.at(2, 3 * i - 1) = answer.at(3, 3 * i - 1) = ( dnx0.at(i, 2) - dnx.at(i, 2) ) / 3.0;
        answer.at(1, 3 * i - 0) = answer.at(2, 3 * i - 0) = answer.at(3, 3 * i - 0) = ( dnx0.at(i, 3) - dnx.at(i, 3) ) / 3.0;
    }

    // add B(gp)
    for ( i = 1; i <= 8; i++ ) {
        answer.at(1, 3 * i - 2) += dnx.at(i, 1);
        answer.at(2, 3 * i - 1) += dnx.at(i, 2);
        answer.at(3, 3 * i - 0) += dnx.at(i, 3);
    }

    for ( i = 1; i <= 8; i++ ) {
        answer.at(4, 3 * i - 1) += dnx.at(i, 3);
        answer.at(4, 3 * i - 0) += dnx.at(i, 2);

        answer.at(5, 3 * i - 2) += dnx.at(i, 3);
        answer.at(5, 3 * i - 0) += dnx.at(i, 1);

        answer.at(6, 3 * i - 2) += dnx.at(i, 2);
        answer.at(6, 3 * i - 1) += dnx.at(i, 1);
    }
}

void FE2FluidMaterial :: giveStiffnessMatrices(FloatMatrix &dsdd, FloatArray &dsdp, FloatArray &dedd, double &dedp,
                                               MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    FE2FluidMaterialStatus *ms = static_cast< FE2FluidMaterialStatus * >( this->giveStatus(gp) );
    ms->computeTangents(tStep);
    if ( mode == TangentStiffness ) {
        dsdd = ms->giveDeviatoricTangent();
        dsdp = ms->giveDeviatoricPressureTangent();
        dedd = ms->giveVolumetricDeviatoricTangent();
        dedp = ms->giveVolumetricPressureTangent();
#if 0
        // Numerical ATS for debugging
        FloatMatrix numericalATS(6, 6);
        FloatArray dsig;
        FloatArray tempStrain(6);

        tempStrain.zero();
        FloatArray sig, strain, sigPert;
        double epspvol;
        computeDeviatoricStressVector(sig, epspvol, gp, tempStrain, 0., tStep);
        double h = 0.001; // Linear problem, size of this doesn't matter.
        for ( int k = 1; k <= 6; ++k ) {
            strain = tempStrain;
            strain.at(k) += h;
            double tmp = strain.at(1) + strain.at(2) + strain.at(3);
            strain.at(1) -= tmp/3.0;
            strain.at(2) -= tmp/3.0;
            strain.at(3) -= tmp/3.0;
            strain.printYourself();
            computeDeviatoricStressVector(sigPert, epspvol, gp, strain, 0., tStep);
            sigPert.printYourself();
            dsig.beDifferenceOf(sigPert, sig);
            numericalATS.setColumn(dsig, k);
        }
        numericalATS.times(1. / h);

        printf("Analytical deviatoric tangent = ");
        dsdd.printYourself();
        printf("Numerical deviatoric tangent = ");
        numericalATS.printYourself();
        numericalATS.subtract(dsdd);
        double norm = numericalATS.computeFrobeniusNorm();
        if ( norm > dsdd.computeFrobeniusNorm() * DEBUG_ERR && norm > 0.0 ) {
            OOFEM_ERROR("Error in deviatoric tangent");
        }
#endif
#if 0
        // Numerical ATS for debugging
        FloatArray strain(3);
        strain.zero();
        FloatArray sig, sigh;
        double epspvol, pressure = 0.0;
        double h = 1.00; // Linear problem, size of this doesn't matter.
        computeDeviatoricStressVector(sig, epspvol, gp, strain, pressure, tStep);
        computeDeviatoricStressVector(sigh, epspvol, gp, strain, pressure + h, tStep);

        FloatArray dsigh;
        dsigh.beDifferenceOf(sigh, sig);
        dsigh.times(1 / h);

        printf("Analytical deviatoric pressure tangent = ");
        dsdp.printYourself();
        printf("Numerical deviatoric pressure tangent = ");
        dsigh.printYourself();
        dsigh.subtract(dsdp);
        double norm = dsigh.computeNorm();
        if ( norm > dsdp.computeNorm() * DEBUG_ERR && norm > 0.0 ) {
            OOFEM_ERROR("Error in deviatoric pressure tangent");
        }
#endif
#if 0
        // Numerical ATS for debugging
        FloatArray tempStrain(3);
        tempStrain.zero();
        FloatArray sig, strain;
        double epspvol, epspvol11, epspvol22, epspvol12, pressure = 0.0;
        double h = 1.0; // Linear problem, size of this doesn't matter.

        computeDeviatoricStressVector(sig, epspvol, gp, tempStrain, pressure, tStep);
        strain = tempStrain;
        strain.at(1) += h;
        computeDeviatoricStressVector(sig, epspvol11, gp, strain, pressure, tStep);
        strain = tempStrain;
        strain.at(2) += h;
        computeDeviatoricStressVector(sig, epspvol22, gp, strain, pressure, tStep);
        strain = tempStrain;
        strain.at(3) += h;
        computeDeviatoricStressVector(sig, epspvol12, gp, strain, pressure, tStep);

        FloatArray dvol(3);
        dvol.at(1) = ( epspvol11 - epspvol ) / h;
        dvol.at(2) = ( epspvol22 - epspvol ) / h;
        dvol.at(3) = ( epspvol12 - epspvol ) / h;
        dvol.at(1) += 1.0;
        dvol.at(2) += 1.0;

        printf("Analytical volumetric deviatoric tangent = ");
        dedd.printYourself();
        printf("Numerical volumetric deviatoric tangent = ");
        dvol.printYourself();
//.........这里部分代码省略.........

void
Structural3DElement :: computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int li, int ui)
// Returns the [ 6 x (nno*3) ] strain-displacement matrix {B} of the receiver, eva-
// luated at gp.
// B matrix  -  6 rows : epsilon-X, epsilon-Y, epsilon-Z, gamma-YZ, gamma-ZX, gamma-XY  :
{
    FEInterpolation *interp = this->giveInterpolation();
    FloatMatrix dNdx; 
    interp->evaldNdx( dNdx, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
    
    answer.resize(6, dNdx.giveNumberOfRows() * 3);
    answer.zero();

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

        answer.at(5, 3 * i - 2) = answer.at(4, 3 * i - 1) = dNdx.at(i, 3);
        answer.at(6, 3 * i - 2) = answer.at(4, 3 * i - 0) = dNdx.at(i, 2);
        answer.at(6, 3 * i - 1) = answer.at(5, 3 * i - 0) = dNdx.at(i, 1);
    }
}

void PrescribedGradientBCWeak :: assembleGPContrib(SparseMtrx &answer, TimeStep *tStep,
                      CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s, TracSegArray &iEl, GaussPoint &iGP)
{

    SpatialLocalizer *localizer = domain->giveSpatialLocalizer();

    ///////////////
    // Gamma_plus
	FloatMatrix contrib;
	assembleTangentGPContributionNew(contrib, iEl, iGP, -1.0, iGP.giveGlobalCoordinates());

    // Compute vector of traction unknowns
    FloatArray tracUnknowns;
    iEl.mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), VM_Total, tStep);

    IntArray trac_rows;
    iEl.giveTractionLocationArray(trac_rows, type, r_s);


    FloatArray dispElLocCoord, closestPoint;
    Element *dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, iGP.giveGlobalCoordinates() );

    IntArray disp_cols;
    dispEl->giveLocationArray(disp_cols, c_s);

    answer.assemble(trac_rows, disp_cols, contrib);

    FloatMatrix contribT;
    contribT.beTranspositionOf(contrib);
    answer.assemble(disp_cols, trac_rows, contribT);



    ///////////////
    // Gamma_minus
    contrib.clear();
	FloatArray xMinus;
	this->giveMirroredPointOnGammaMinus(xMinus, iGP.giveGlobalCoordinates() );
	assembleTangentGPContributionNew(contrib, iEl, iGP, 1.0, xMinus);

    // Compute vector of traction unknowns
	tracUnknowns.clear();
	iEl.mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), VM_Total, tStep);

    trac_rows.clear();
    iEl.giveTractionLocationArray(trac_rows, type, r_s);


    dispElLocCoord.clear(); closestPoint.clear();
    dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, xMinus );

    disp_cols.clear();
    dispEl->giveLocationArray(disp_cols, c_s);

    answer.assemble(trac_rows, disp_cols, contrib);

    contribT.clear();
    contribT.beTranspositionOf(contrib);
    answer.assemble(disp_cols, trac_rows, contribT);


    // Assemble zeros on diagonal (required by PETSc solver)
	FloatMatrix KZero(1,1);
	KZero.zero();
    for( int i :  trac_rows) {
        answer.assemble(IntArray({i}), IntArray({i}), KZero);
    }
}

void PrescribedGradientBCWeak :: computeField(FloatArray &sigma, TimeStep *tStep)
{
	double Lx = mUC[0] - mLC[0];
	double Ly = mUC[1] - mLC[1];
	double dSize = Lx*Ly;
//    printf("dSize: %e\n", dSize);

    const int dim = domain->giveNumberOfSpatialDimensions();
    FloatMatrix stressMatrix(dim, dim);

    for( auto *el: mpTracElNew ) {

        for ( GaussPoint *gp: *(el->mIntRule.get()) ) {

            // For now, assume piecewise constant approx
            FloatArray Ntrac = FloatArray { 1.0 };

            // N-matrix
            FloatMatrix Nmat;
            Nmat.beNMatrixOf(Ntrac, dim);

        	// Fetch global coordinate x
        	const FloatArray &x = gp->giveGlobalCoordinates();

            // Compute vector of traction unknowns
            FloatArray tracUnknowns;
            el->mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), VM_Total, tStep);


            FloatArray traction;
            traction.beProductOf(Nmat, tracUnknowns);

            FloatArray tmp;
            giveBoundaryCoordVector(tmp, x);

            FloatMatrix contrib;
            contrib.beDyadicProductOf(traction, tmp);

            double detJ = 0.5 * el->giveLength();
            contrib.times( detJ * gp->giveWeight() );

            for ( int m = 0; m < dim; m++ ) {
                for ( int n = 0; n < dim; n++ ) {
                    stressMatrix(m, n) += contrib(m, n);
                }
            }

        }

    }

    if ( dim == 2 ) {
        sigma = {
            stressMatrix(0, 0), stressMatrix(1, 1), 0.0, 0.0, 0.0, 0.5*(stressMatrix(0, 1) + stressMatrix(1, 0))
        };
    } else   {
        sigma.beVectorForm(stressMatrix);
    }

    sigma.times(1.0 / dSize);

}

void
MisesMatGrad :: giveInternalLength(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    answer.resize(1, 1);
    answer.at(1, 1) = L;
}

void
StructuralFE2Material :: give3dMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    if ( useNumTangent ) {
        // Numerical tangent
        StructuralFE2MaterialStatus *status = static_cast<StructuralFE2MaterialStatus*>( this->giveStatus( gp ) );
        double h = 1.0e-9;

        const FloatArray &epsRed = status->giveTempStrainVector();
        FloatArray eps;
        StructuralMaterial::giveFullSymVectorForm(eps, epsRed, gp->giveMaterialMode() );


        int dim = eps.giveSize();
        answer.resize(dim, dim);
        answer.zero();

        FloatArray sig, sigPert, epsPert;

        for(int i = 1; i <= dim; i++) {
            // Add a small perturbation to the strain
            epsPert = eps;
            epsPert.at(i) += h;

            giveRealStressVector_3d(sigPert, gp, epsPert, tStep);
            answer.setColumn(sigPert, i);
        }

        giveRealStressVector_3d(sig, gp, eps, tStep);

        for(int i = 1; i <= dim; i++) {
            for(int j = 1; j <= dim; j++) {
                answer.at(j,i) -= sig.at(j);
                answer.at(j,i) /= h;
            }
        }

    } else {

        StructuralFE2MaterialStatus *ms = static_cast< StructuralFE2MaterialStatus * >( this->giveStatus(gp) );
        ms->computeTangent(tStep);
        const FloatMatrix &ans9 = ms->giveTangent();

        // Compute the (minor) symmetrized tangent:
        answer.resize(6, 6);
        for ( int i = 0; i < 6; ++i ) {
            for ( int j = 0; j < 6; ++j ) {
                answer(i, j) = ans9(i, j);
            }
        }
        for ( int i = 0; i < 6; ++i ) {
            for ( int j = 6; j < 9; ++j ) {
                answer(i, j-3) += ans9(i, j);
                answer(j-3, i) += ans9(j, i);
            }
        }
        for ( int i = 6; i < 9; ++i ) {
            for ( int j = 6; j < 9; ++j ) {
                answer(j-3, i-3) += ans9(j, i);
            }
        }
        for ( int i = 0; i < 6; ++i ) {
            for ( int j = 3; j < 6; ++j ) {
                answer(j, i) *= 0.5;
                answer(i, j) *= 0.5;
            }
        }
#if 0
        // Numerical ATS for debugging
        FloatMatrix numericalATS(6, 6);
        FloatArray dsig;
        // Note! We need a copy of the temp strain, since the pertubations might change it.
        FloatArray tempStrain = ms->giveTempStrainVector();

        FloatArray sig, strain, sigPert;
        giveRealStressVector_3d(sig, gp, tempStrain, tStep);
        double hh = 1e-6;
        for ( int k = 1; k <= 6; ++k ) {
            strain = tempStrain;
            strain.at(k) += hh;
            giveRealStressVector_3d(sigPert, gp, strain, tStep);
            dsig.beDifferenceOf(sigPert, sig);
            numericalATS.setColumn(dsig, k);
        }
        numericalATS.times(1. / hh);
        giveRealStressVector_3d(sig, gp, tempStrain, tStep); // Reset

        //answer.printYourself("Analytical deviatoric tangent");
        //numericalATS.printYourself("Numerical deviatoric tangent");

        numericalATS.subtract(answer);
        double norm = numericalATS.computeFrobeniusNorm();
        if ( norm > answer.computeFrobeniusNorm() * 1e-3 && norm > 0.0 ) {
            OOFEM_ERROR("Error in deviatoric tangent");
        }
#endif
    }
}

void
LIBeam3dNL :: computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    int i, j, k;
    double s1, s2;
    FloatMatrix d, x, xt(12, 6), dxt, sn, sm, sxd, y;
    FloatArray n(3), m(3), xd(3), TotalStressVector;
    IntegrationRule *iRule = integrationRulesArray [ giveDefaultIntegrationRule() ];
    GaussPoint *gp = iRule->getIntegrationPoint(0);

    answer.resize( this->computeNumberOfDofs(EID_MomentumBalance), this->computeNumberOfDofs(EID_MomentumBalance) );
    answer.zero();

    // linear part

    this->updateTempTriad(tStep); // update temp triad
    this->computeXMtrx(x, tStep);
    xt.zero();
    for ( i = 1; i <= 12; i++ ) {
        for ( j = 1; j <= 3; j++ ) {
            for ( k = 1; k <= 3; k++ ) {
                // compute x*Tbar, taking into account sparsity of Tbar
                xt.at(i, j)   += x.at(i, k) * tempTc.at(k, j);
                xt.at(i, j + 3) += x.at(i, k + 3) * tempTc.at(k, j);
            }
        }
    }

    this->computeConstitutiveMatrixAt(d, rMode, gp, tStep);
    dxt.beProductTOf(d, xt);
    answer.beProductOf(xt, dxt);
    answer.times(1. / this->l0);

    // geometric stiffness ks = ks1+ks2
    // ks1
    this->computeStressVector(TotalStressVector, gp, tStep);

    for ( i = 1; i <= 3; i++ ) {
        s1 = s2 = 0.0;
        for ( j = 1; j <= 3; j++ ) {
            s1 += tempTc.at(i, j) * TotalStressVector.at(j);
            s2 += tempTc.at(i, j) * TotalStressVector.at(j + 3);
        }

        n.at(i)   = s1;
        m.at(i)   = s2;
    }

    this->computeSMtrx(sn, n);
    this->computeSMtrx(sm, m);

    for ( i = 1; i <= 3; i++ ) {
        for ( j = 1; j <= 3; j++ ) {
            answer.at(i, j + 3)   += sn.at(i, j);
            answer.at(i, j + 9)   += sn.at(i, j);
            answer.at(i + 3, j + 3) += sm.at(i, j);
            answer.at(i + 3, j + 9) += sm.at(i, j);

            answer.at(i + 6, j + 3) -= sn.at(i, j);
            answer.at(i + 6, j + 9) -= sn.at(i, j);
            answer.at(i + 9, j + 3) -= sm.at(i, j);
            answer.at(i + 9, j + 9) -= sm.at(i, j);
        }
    }

    // ks2
    this->computeXdVector(xd, tStep);
    this->computeSMtrx(sxd, xd);

    y.beProductOf(sxd, sn);
    y.times(0.5);

    for ( i = 1; i <= 3; i++ ) {
        for ( j = 1; j <= 3; j++ ) {
            answer.at(i + 3, j)     -= sn.at(i, j);
            answer.at(i + 3, j + 3)   += y.at(i, j);
            answer.at(i + 3, j + 6)   += sn.at(i, j);
            answer.at(i + 3, j + 9)   += y.at(i, j);

            answer.at(i + 9, j)     -= sn.at(i, j);
            answer.at(i + 9, j + 3)   += y.at(i, j);
            answer.at(i + 9, j + 6)   += sn.at(i, j);
            answer.at(i + 9, j + 9)   += y.at(i, j);
        }
    }
}

void
Structural3DElement :: computeBHmatrixAt(GaussPoint *gp, FloatMatrix &answer)
// Returns the [ 9 x (nno * 3) ] displacement gradient matrix {BH} of the receiver,
// evaluated at gp.
// BH matrix  -  9 rows : du/dx, dv/dy, dw/dz, dv/dz, du/dz, du/dy, dw/dy, dw/dx, dv/dx
{
    FEInterpolation *interp = this->giveInterpolation();
    FloatMatrix dNdx; 
    interp->evaldNdx( dNdx, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
    
    answer.resize(9, dNdx.giveNumberOfRows() * 3);
    answer.zero();

    for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
        answer.at(1, 3 * i - 2) = dNdx.at(i, 1);     // du/dx
        answer.at(2, 3 * i - 1) = dNdx.at(i, 2);     // dv/dy
        answer.at(3, 3 * i - 0) = dNdx.at(i, 3);     // dw/dz
        answer.at(4, 3 * i - 1) = dNdx.at(i, 3);     // dv/dz
        answer.at(7, 3 * i - 0) = dNdx.at(i, 2);     // dw/dy
        answer.at(5, 3 * i - 2) = dNdx.at(i, 3);     // du/dz
        answer.at(8, 3 * i - 0) = dNdx.at(i, 1);     // dw/dx
        answer.at(6, 3 * i - 2) = dNdx.at(i, 2);     // du/dy
        answer.at(9, 3 * i - 1) = dNdx.at(i, 1);     // dv/dx
    }

}