C++ GaussPoint::giveWeight方法代码示例

void Tr1Darcy :: computeStiffnessMatrix(FloatMatrix &answer, TimeStep *atTime)
{
    /*
     * Return Ke = integrate(B^T K B)
     */

    FloatMatrix B, BT, K, KB;
    FloatArray *lcoords;
    GaussPoint *gp;

    TransportMaterial *mat = ( TransportMaterial * ) this->domain->giveMaterial(this->material);

    IntegrationRule *iRule = integrationRulesArray [ 0 ];

    answer.resize(3, 3);
    answer.zero();

    for ( int i = 0; i < iRule->getNumberOfIntegrationPoints(); i++ ) {
        gp = iRule->getIntegrationPoint(i);
        lcoords = gp->giveCoordinates();

        double detJ = this->interpolation_lin.giveTransformationJacobian( * lcoords, FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evaldNdx( BT, * lcoords, FEIElementGeometryWrapper(this) );
        
        mat->giveCharacteristicMatrix(K, FullForm, TangentStiffness, gp, atTime);

        B.beTranspositionOf(BT);
        KB.beProductOf(K, B);
        answer.plusProductUnsym(B, KB, detJ * gp->giveWeight() ); // Symmetric part is just a single value, not worth it.
    }
}

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

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

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


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

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

    Nmatrix.resize(2,12);

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

        gp = iRule->getIntegrationPoint(i);

        lcoords = gp->giveCoordinates();

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

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

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

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

    }

}

void Tr1Darcy :: computeInternalForcesVector(FloatArray &answer, TimeStep *atTime)
{
    FloatArray *lcoords, w, a, gradP, I;
    FloatMatrix BT;
    GaussPoint *gp;

    TransportMaterial *mat = ( TransportMaterial * ) this->domain->giveMaterial(this->material);
    IntegrationRule *iRule = integrationRulesArray [ 0 ];

    this->computeVectorOf(EID_ConservationEquation, VM_Total, atTime, a);

    answer.resize(3);
    answer.zero();

    for ( int i = 0; i < iRule->getNumberOfIntegrationPoints(); i++ ) {
        gp = iRule->getIntegrationPoint(i);
        lcoords = gp->giveCoordinates();

        double detJ = this->interpolation_lin.giveTransformationJacobian( * lcoords, FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evaldNdx( BT, * lcoords, FEIElementGeometryWrapper(this) );

        gradP.beTProductOf(BT, a);

        mat->giveFluxVector(w, gp, gradP, atTime);

        I.beProductOf(BT, w);
        answer.add(- gp->giveWeight() * detJ, I);
    }
}

void Tr21Stokes :: computeBodyLoadVectorAt(FloatArray &answer, Load *load, TimeStep *tStep)
{
    IntegrationRule *iRule = this->integrationRulesArray [ 0 ];
    GaussPoint *gp;
    FloatArray N, gVector, *lcoords, temparray(15);
    double dA, detJ, rho;

    load->computeComponentArrayAt(gVector, tStep, VM_Total);
    temparray.zero();
    if ( gVector.giveSize() ) {
        for ( int k = 0; k < iRule->getNumberOfIntegrationPoints(); k++ ) {
            gp = iRule->getIntegrationPoint(k);
            lcoords = gp->giveCoordinates();

            rho = this->giveMaterial()->giveCharacteristicValue(MRM_Density, gp, tStep);
            detJ = fabs( this->interpolation_quad.giveTransformationJacobian(* lcoords, FEIElementGeometryWrapper(this)) );
            dA = detJ * gp->giveWeight();

            this->interpolation_quad.evalN(N, * lcoords, FEIElementGeometryWrapper(this));
            for ( int j = 0; j < 6; j++ ) {
                temparray(2 * j)     += N(j) * rho * gVector(0) * dA;
                temparray(2 * j + 1) += N(j) * rho * gVector(1) * dA;
            }
        }
    }

    answer.resize(15);
    answer.zero();
    answer.assemble( temparray, this->ordering );
}

void PrescribedGradientBCWeak :: assembleTangentGPContributionNew(FloatMatrix &oTangent, TracSegArray &iEl, GaussPoint &iGP, const double &iScaleFactor, const FloatArray &iBndCoord)
{
    int dim = domain->giveNumberOfSpatialDimensions();
    double detJ = 0.5 * iEl.giveLength();

    //////////////////////////////////
    // Compute traction N-matrix
    // For now, assume piecewise constant approx
    FloatArray Ntrac = FloatArray { 1.0 };
    FloatMatrix NtracMat;
    NtracMat.beNMatrixOf(Ntrac, dim);


    //////////////////////////////////
    // Compute displacement N-matrix
    // Identify the displacement element
    // we are currently standing in
    // and compute local coordinates on
    // the displacement element
    SpatialLocalizer *localizer = domain->giveSpatialLocalizer();
    FloatArray dispElLocCoord, closestPoint;
    Element *dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, iBndCoord );


    // Compute basis functions
    XfemElementInterface *xfemElInt = dynamic_cast< XfemElementInterface * >( dispEl );
    FloatMatrix NdispMat;

    if ( xfemElInt != NULL && domain->hasXfemManager() ) {
        // If the element is an XFEM element, we use the XfemElementInterface to compute the N-matrix
        // of the enriched element.
        xfemElInt->XfemElementInterface_createEnrNmatrixAt(NdispMat, dispElLocCoord, * dispEl, false);
    } else   {
        // Otherwise, use the usual N-matrix.
        const int numNodes = dispEl->giveNumberOfDofManagers();
        FloatArray N(numNodes);

        const int dim = dispEl->giveSpatialDimension();

        NdispMat.resize(dim, dim * numNodes);
        NdispMat.zero();
        dispEl->giveInterpolation()->evalN( N, dispElLocCoord, FEIElementGeometryWrapper(dispEl) );

        NdispMat.beNMatrixOf(N, dim);
    }

    FloatMatrix contrib;
    contrib.beTProductOf(NtracMat, NdispMat);
    contrib.times( iScaleFactor * detJ * iGP.giveWeight() );

    oTangent = contrib;
}

void Line2SurfaceTension :: computeLoadVector(FloatArray &answer, ValueModeType mode, TimeStep *tStep)
{
    ///@todo Support axisymm.
    //domainType dt = this->giveDomain()->giveDomainType();
    IntegrationRule *iRule = this->integrationRulesArray [ 0 ];
    double t = 1, gamma_s;
    ///@todo Should i use this? Not used in FM module (but perhaps it should?) / Mikael.
    //t = this->giveDomain()->giveCrossSection(1)->give(CS_Thickness);
    gamma_s = this->giveMaterial()->give('g', NULL);

    FloatMatrix xy(2, 3);
    Node *node;
    for ( int i = 1; i <= 3; i++ ) {
        node = giveNode(i);
        xy.at(1, i) = node->giveCoordinate(1);
        xy.at(2, i) = node->giveCoordinate(2);
    }

    FloatArray A;
    FloatArray dNdxi(3);
    FloatArray es(2); // tangent vector to curve
    FloatMatrix BJ(2, 6);
    BJ.zero();

    answer.resize(6);
    answer.zero();

    for ( int k = 0; k < iRule->getNumberOfIntegrationPoints(); k++ ) {
        GaussPoint *gp = iRule->getIntegrationPoint(k);
        //interpolation.evaldNdx(dN, domain, dofManArray, * gp->giveCoordinates(), 0.0);
        double xi = gp->giveCoordinate(1);

        // Some simplifications can be performed, since the mapping J is a scalar.
        dNdxi.at(1) = -0.5 + xi;
        dNdxi.at(2) =  0.5 + xi;
        dNdxi.at(3) = -2.0 * xi;

        es.beProductOf(xy, dNdxi);
        double J = es.computeNorm();
        es.times(1 / J); //es.normalize();

        // dNds = dNdxi/J
        // B.at(1,1) = dNds.at(1); and so on.

        BJ.at(1, 1) = BJ.at(2, 2) = dNdxi.at(1);
        BJ.at(1, 3) = BJ.at(2, 4) = dNdxi.at(2);
        BJ.at(1, 5) = BJ.at(2, 6) = dNdxi.at(3);

        A.beTProductOf(BJ, es);
        answer.add( - gamma_s * t * gp->giveWeight(), A); // Note! Negative sign!
    }
}

void Tr21Stokes :: computeInternalForcesVector(FloatArray &answer, TimeStep *tStep)
{
    IntegrationRule *iRule = integrationRulesArray [ 0 ];
    FluidDynamicMaterial *mat = ( FluidDynamicMaterial * ) this->domain->giveMaterial(this->material);
    FloatArray a_pressure, a_velocity, devStress, epsp, BTs, Nh, dNv(12);
    double r_vol, pressure;
    FloatMatrix dN, B(3, 12);
    B.zero();
    
    this->computeVectorOf(EID_MomentumBalance, VM_Total, tStep, a_velocity);
    this->computeVectorOf(EID_ConservationEquation, VM_Total, tStep, a_pressure);
    
    FloatArray momentum(12), conservation(3);
    momentum.zero();
    conservation.zero();
    GaussPoint *gp;

    for ( int i = 0; i < iRule->getNumberOfIntegrationPoints(); i++ ) {
        gp = iRule->getIntegrationPoint(i);
        FloatArray *lcoords = gp->giveCoordinates();

        double detJ = fabs(this->interpolation_quad.giveTransformationJacobian(* lcoords, FEIElementGeometryWrapper(this)));
        this->interpolation_quad.evaldNdx(dN, * lcoords, FEIElementGeometryWrapper(this));
        this->interpolation_lin.evalN(Nh, * lcoords, FEIElementGeometryWrapper(this));
        double dA = detJ * gp->giveWeight();

        for ( int j = 0, k = 0; j < 6; j++, k += 2 ) {
            dNv(k)     = B(0, k)     = B(2, k + 1) = dN(j, 0);
            dNv(k + 1) = B(1, k + 1) = B(2, k)     = dN(j, 1);
        }

        pressure = Nh.dotProduct(a_pressure);
        epsp.beProductOf(B, a_velocity);

        mat->computeDeviatoricStressVector(devStress, r_vol, gp, epsp, pressure, tStep);
        BTs.beTProductOf(B, devStress);

        momentum.add(dA, BTs);
        momentum.add(-pressure*dA, dNv);
        conservation.add(r_vol*dA, Nh);
    }

    FloatArray temp(15);
    temp.zero();
    temp.addSubVector(momentum, 1);
    temp.addSubVector(conservation, 13);

    answer.resize(15);
    answer.zero();
    answer.assemble(temp, this->ordering);
}

void Tr21Stokes :: computeEdgeBCSubVectorAt(FloatArray &answer, Load *load, int iEdge, TimeStep *tStep)
{
    answer.resize(15);
    answer.zero();

    if ( load->giveType() == TransmissionBC ) { // Neumann boundary conditions (traction)
        BoundaryLoad *boundaryLoad = ( BoundaryLoad * ) load;

        int numberOfEdgeIPs = ( int ) ceil( ( boundaryLoad->giveApproxOrder() + 1. ) / 2. ) * 2;

        GaussIntegrationRule iRule(1, this, 1, 1);
        GaussPoint *gp;
        FloatArray N, t, f(6);
        IntArray edge_mapping;

        f.zero();
        iRule.setUpIntegrationPoints(_Line, numberOfEdgeIPs, _Unknown);

        for ( int i = 0; i < iRule.getNumberOfIntegrationPoints(); i++ ) {
            gp = iRule.getIntegrationPoint(i);
            FloatArray *lcoords = gp->giveCoordinates();

            this->interpolation_quad.edgeEvalN(N, * lcoords, FEIElementGeometryWrapper(this));
            double detJ = fabs(this->interpolation_quad.edgeGiveTransformationJacobian(iEdge, * lcoords, FEIElementGeometryWrapper(this)));
            double dS = gp->giveWeight() * detJ;

            if ( boundaryLoad->giveFormulationType() == BoundaryLoad :: BL_EntityFormulation ) { // Edge load in xi-eta system
                boundaryLoad->computeValueAt(t, tStep, * lcoords, VM_Total);
            } else   { // Edge load in x-y system
                FloatArray gcoords;
                this->interpolation_quad.edgeLocal2global(gcoords, iEdge, * lcoords, FEIElementGeometryWrapper(this));
                boundaryLoad->computeValueAt(t, tStep, gcoords, VM_Total);
            }

            // Reshape the vector
            for ( int j = 0; j < 3; j++ ) {
                f(2 * j)     += N(j) * t(0) * dS;
                f(2 * j + 1) += N(j) * t(1) * dS;
            }
        }

        answer.assemble(f, this->edge_ordering [ iEdge - 1 ]);
    } else   {
        OOFEM_ERROR("Tr21Stokes :: computeEdgeBCSubVectorAt - Strange boundary condition type");
    }
}

int
PatchIntegrationRule :: SetUpPointsOnTriangle(int nPoints, MaterialMode mode, GaussPoint ***arry)
{
    numberOfIntegrationPoints = GaussIntegrationRule :: SetUpPointsOnTriangle(nPoints, mode, arry);
    firstLocalStrainIndx = 1;
    lastLocalStrainIndx = 3;
    // convert patch coordinates into element based, update weights accordingly
    for ( int j = 0; j <  numberOfIntegrationPoints; j++ ) {
        GaussPoint *gp = ( * arry ) [ j ];
        patch->convertGPIntoParental(gp); // convert coordinates into parental
        Element *elg = ( Element * ) patch->giveParent();
        double parentArea = elg->computeArea();
        Triangle *tr = ( Triangle * ) patch;
        gp->setWeight(8.0 * gp->giveWeight() * tr->getArea() / parentArea); // update integration weight
    }

    return numberOfIntegrationPoints;
}

void Tr21Stokes :: giveElementFMatrix(FloatMatrix &answer)
{

    IntegrationRule *iRule = integrationRulesArray [ 0 ];
    GaussPoint *gp;
    double detJ;
    FloatArray N, N2, *lcoords;
    IntArray col;
    FloatMatrix temp;

    N2.resize(6);   N2.zero();
    col.resize(2);  col.at(1)=1;  col.at(2)=2;

    temp.resize(15,2);
    temp.zero();

    for (int i=0; i<iRule->getNumberOfIntegrationPoints(); i++) {
        gp = iRule->getIntegrationPoint(i);
        lcoords = gp->giveCoordinates();

        this->interpolation_quad.evalN(N, *lcoords, FEIElementGeometryWrapper(this));
        detJ = this->interpolation_quad.giveTransformationJacobian(*lcoords, FEIElementGeometryWrapper(this));
        N.times(gp->giveWeight()*detJ);
        //N.printYourself();
        N2.add(N);
    }

    for (int i=1; i<=6; i++) {
        temp.at(i*2-1, 1)=N2.at(i);
        temp.at(i*2, 2)=N2.at(i);
    }

    answer.resize(17,2);
    answer.zero();
    answer.assemble(temp, this->ordering, col);

}

void
tet21ghostsolid :: computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{

    IntegrationRule *iRule = integrationRulesArray [ giveDefaultIntegrationRule() ];
#ifdef __FM_MODULE
    FluidDynamicMaterial *fluidMaterial = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
#endif

    FloatMatrix Kf, G, Kx, D, B, Ed, EdB, dNx;
    FloatArray Nlin, dNv;

    for (int j = 0; j<iRule->giveNumberOfIntegrationPoints(); j++) {
        GaussPoint *gp = iRule->getIntegrationPoint(j);

        double detJ = fabs( ( this->interpolation.giveTransformationJacobian( * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) ) ) );
        double weight = gp->giveWeight();

        this->interpolation.evaldNdx( dNx, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evalN( Nlin, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

        dNv.resize(30); // dNv = [dN1/dx dN1/dy dN1/dz dN2/dx dN2/dy dN2/dz ... dN10/dz]

        for (int k = 0; k<dNx.giveNumberOfRows(); k++) {
            dNv.at(k*3+1) = dNx.at(k+1,1);
            dNv.at(k*3+2) = dNx.at(k+1,2);
            dNv.at(k*3+3) = dNx.at(k+1,3);
        }

        if (nlGeometry == 0) {

            this->computeBmatrixAt(gp, B);

            // Fluid part
            gp->setMaterialMode(_3dFlow);
#ifdef __FM_MODULE
            fluidMaterial->giveDeviatoricStiffnessMatrix(Ed, TangentStiffness, gp, tStep);
#else
            OOFEM_ERROR("Fluid module missing\n");
#endif
            gp->setMaterialMode(_3dMat);

            EdB.beProductOf(Ed, B);
            Kf.plusProductSymmUpper(B, EdB, detJ*weight);

            // Ghost solid part
            EdB.beProductOf(Dghost, B);
            Kx.plusProductSymmUpper(B, EdB, detJ*weight);

            // Incompressibility part
            G.plusDyadUnsym(dNv, Nlin, -detJ*weight);

        } else {
            OOFEM_ERROR ("No support for large deformations yet!");
        }

    }

    FloatMatrix GT;

    GT.beTranspositionOf(G);
    //GTdeltat.beTranspositionOf(G);
    //GTdeltat.times(deltat);
    Kf.symmetrized();
    Kx.symmetrized();
    //    Kf.printYourself();
    //    G.printYourself();
    //    GT.printYourself();
    //    Kx.printYourself();

    answer.resize(64, 64);
    answer.zero();
#define USEUNCOUPLED 0

#if USEUNCOUPLED == 1
    // Totaly uncoupled
    answer.assemble(Kf, momentum_ordering, momentum_ordering);
    answer.assemble(G, momentum_ordering, conservation_ordering);
    answer.assemble(GT, conservation_ordering, momentum_ordering);
    answer.assemble(Kx, ghostdisplacement_ordering, ghostdisplacement_ordering);
#else
    answer.assemble(Kf, ghostdisplacement_ordering, ghostdisplacement_ordering);
    answer.assemble(Kf, ghostdisplacement_ordering, momentum_ordering);
    answer.assemble(G, ghostdisplacement_ordering, conservation_ordering);
    answer.assemble(GT, conservation_ordering, ghostdisplacement_ordering);
    answer.assemble(GT, conservation_ordering, momentum_ordering);
    answer.assemble(Kx, momentum_ordering, ghostdisplacement_ordering);
#endif

    //answer.printYourself();

}

void
tet21ghostsolid :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{

    IntegrationRule *iRule = integrationRulesArray [ giveDefaultIntegrationRule() ];
#ifdef __FM_MODULE
    FluidDynamicMaterial *fluidMaterial = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
#endif

    FloatMatrix Kf, G, Kx, B, Ed, dNx;
    FloatArray Strain, Stress, Nlin, dNv, a, aVelocity, aPressure, aGhostDisplacement, fluidStress, epsf;
    FloatArray momentum, conservation, auxstress;
    double pressure, epsvol;

    this->computeVectorOf( VM_Total, tStep, a);
    if (!tStep->isTheFirstStep()) {
        // a.printYourself();
    }

    aVelocity.beSubArrayOf(a, momentum_ordering);
    aPressure.beSubArrayOf(a, conservation_ordering);
    aGhostDisplacement.beSubArrayOf(a, ghostdisplacement_ordering);

    for (int j = 0; j<iRule->giveNumberOfIntegrationPoints(); j++) {
        GaussPoint *gp = iRule->getIntegrationPoint(j);

        double detJ = fabs( ( this->interpolation.giveTransformationJacobian( * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) ) ) );
        double weight = gp->giveWeight();

        this->interpolation.evaldNdx( dNx, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evalN( Nlin, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

        dNv.resize(30);
        for (int k = 0; k<dNx.giveNumberOfColumns(); k++) {
            dNv.at(k*3+1) = dNx.at(1,k+1);
            dNv.at(k*3+2) = dNx.at(2,k+1);
            dNv.at(k*3+3) = dNx.at(3,k+1);
        }

        if (nlGeometry == 0) {

            this->computeBmatrixAt(gp, B);
            epsf.beProductOf(B, aVelocity);
            pressure = Nlin.dotProduct(aPressure);

            // Fluid part
            gp->setMaterialMode(_3dFlow);
#ifdef __FM_MODULE
            fluidMaterial->computeDeviatoricStressVector(fluidStress, epsvol, gp, epsf, pressure, tStep);
#else
            OOFEM_ERROR("Missing FM module");
#endif
            gp->setMaterialMode(_3dMat);

            momentum.plusProduct(B, fluidStress, detJ*weight);
            momentum.add(-pressure * detJ * weight, dNv);
            conservation.add(epsvol * detJ * weight, Nlin);

            // Ghost solid part
            Strain.beProductOf(B, aGhostDisplacement);
            Stress.beProductOf(Dghost, Strain);
            auxstress.plusProduct(B, Stress, detJ * weight);

        } else {
            OOFEM_ERROR("No support for large deformations yet!");
        }

    }
    answer.resize(64);
    answer.zero();

#if USEUNCOUPLED == 1
    // Totaly uncoupled
    answer.assemble(momentum, momentum_ordering);
    answer.assemble(conservation, conservation_ordering);
    answer.assemble(auxstress, ghostdisplacement_ordering);
#else
    answer.assemble(momentum, ghostdisplacement_ordering);
    answer.assemble(conservation, conservation_ordering);
    answer.assemble(auxstress, momentum_ordering);
#endif

    // Test linear

    /*    if (this->giveNumber() == 364) {
        FloatMatrix K;
        FloatArray ans;
        this->computeStiffnessMatrix(K, TangentStiffness, tStep);
        ans.beProductOf(K, a);
        ans.printYourself();
        answer.printYourself();
    } */


}

void
tet21ghostsolid :: computeLoadVector(FloatArray &answer, Load *load, CharType type, ValueModeType mode, TimeStep *tStep)
{

    answer.resize(64);
    answer.zero();

    if ( type != ExternalForcesVector ) {
        answer.resize(0);
        return;
    }

#ifdef __FM_MODULE
    FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
#endif
    IntegrationRule *iRule = this->integrationRulesArray [ 0 ];
    FloatArray N, gVector, temparray(30), dNv, u, inc, u_prev, vload;
    FloatMatrix dNx, G;

    load->computeComponentArrayAt(gVector, tStep, VM_Total);
    temparray.zero();

    vload.resize(4);
    vload.zero();

    this->giveDisplacementsIncrementData(u_prev, u, inc, tStep);

    for ( int k = 0; k < iRule->giveNumberOfIntegrationPoints(); k++ ) {
        GaussPoint *gp = iRule->getIntegrationPoint(k);
        FloatArray *lcoords = gp->giveNaturalCoordinates();

        double detJ = fabs( this->interpolation.giveTransformationJacobian( * lcoords, FEIElementGeometryWrapper(this) ) );
        double dA = detJ * gp->giveWeight();

        // Body load
        if ( gVector.giveSize() ) {
#ifdef __FM_MODULE
            double rho = mat->give('d', gp);
#else
            OOFEM_ERROR("Missing FM module");
            double rho = 1.0;
#endif
            this->interpolation.evalN( N, * lcoords, FEIElementGeometryWrapper(this) );

            for ( int j = 0; j < N.giveSize(); j++ ) {
                temparray(3 * j + 0) += N(j) * rho * gVector(0) * dA;
                temparray(3 * j + 1) += N(j) * rho * gVector(1) * dA;
                temparray(3 * j + 2) += N(j) * rho * gVector(2) * dA;
            }
        }

        // "load" from previous step
        this->interpolation.evaldNdx( dNx, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evalN( N, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

        dNv.resize(30);
        for (int k = 0; k<dNx.giveNumberOfColumns(); k++) {
            dNv.at(k*3+1) = dNx.at(1,k+1);
            dNv.at(k*3+2) = dNx.at(2,k+1);
            dNv.at(k*3+3) = dNx.at(3,k+1);
        }

        G.plusDyadUnsym(N, dNv, dA);
        FloatMatrix GT;
        GT.beTranspositionOf(G);

        vload.plusProduct(GT, u_prev, -0.0 );
        //vload.printYourself();

    }

#if USEUNCOUPLED == 1
    // Totaly uncoupled
    answer.assemble(temparray, this->momentum_ordering);
#else
    answer.assemble(temparray, this->ghostdisplacement_ordering);
    answer.assemble(vload, this->conservation_ordering);
#endif



    // answer.printYourself();
}

//.........这里部分代码省略.........
            }
#ifdef __VTK_MODULE
            cellVarsArray->SetTuple3(ielem-1,  mtrx.at(1, pos), mtrx.at(2,pos), mtrx.at(3,pos) );
#else
            fprintf( stream, "%f %f %f  ", mtrx.at(1, pos), mtrx.at(2, pos), mtrx.at(3, pos) );
#endif
        }

#ifdef __VTK_MODULE
        stream->GetCellData()->SetActiveVectors(__InternalStateTypeToString(type));
        stream->GetCellData()->SetVectors(cellVarsArray);
#endif
        break;

    default:
        bool reshape = false;
        InternalStateValueType vt = giveInternalStateValueType(type);
        if ( vt == ISVT_SCALAR ) {
            ncomponents = 1;
        } else if ( vt == ISVT_VECTOR ) {
            ncomponents = 3;
        } else {
            ncomponents = 9;
            reshape = true;
        }
#ifdef __VTK_MODULE
        cellVarsArray->SetNumberOfComponents(ncomponents);
        cellVarsArray->SetNumberOfTuples(nelem);
#else
        fprintf( stream, "<DataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"%d\" format=\"ascii\">\n", __InternalStateTypeToString(type), ncomponents );
#endif

        IntArray redIndx;
        for ( ielem = 1; ielem <= nelem; ielem++ ) {
            elem = d->giveElement(ielem);
            if ( (( region > 0 ) && ( this->smoother->giveElementVirtualRegionNumber(ielem) != region ))
                    || this->isElementComposite(elem) || !elem-> isActivated(tStep) ) { // composite cells exported individually
                continue;
            }
#ifdef __PARALLEL_MODE
            if ( elem->giveParallelMode() != Element_local ) {
                continue;
            }
#endif
            gptot = 0;
            answer.resize(0);
            iRule = elem->giveDefaultIntegrationRulePtr();
            if (iRule) {
                MaterialMode mmode = _Unknown;
                for (int i = 0; i < iRule->getNumberOfIntegrationPoints(); ++i) {
                    gp = iRule->getIntegrationPoint(i);
                    mmode = gp->giveMaterialMode();
                    elem->giveIPValue(temp, gp, type, tStep);
                    gptot += gp->giveWeight();
                    answer.add(gp->giveWeight(), temp);
                }
                answer.times(1./gptot);
                elem->giveMaterial()->giveIntVarCompFullIndx(redIndx, type, mmode);
            }
            // Reshape the Voigt vectors to include all components (duplicated if necessary, VTK insists on 9 components for tensors.)
            if ( reshape && answer.giveSize() != 9) { // If it has 9 components, then it is assumed to be proper already.
                FloatArray tmp = answer;
                this->makeFullForm(answer, tmp, vt, redIndx);
            } else if ( vt == ISVT_VECTOR && answer.giveSize() < 3) {
                answer.setValues(3,
                                 answer.giveSize() > 1 ? answer.at(1) : 0.0,
                                 answer.giveSize() > 2 ? answer.at(2) : 0.0,
                                 0.0);
            } else if ( ncomponents != answer.giveSize() ) { // Trying to gracefully handle bad cases, just output zeros.
                answer.resize(ncomponents);
                answer.zero();
            }
            for (int i = 1; i <= ncomponents; ++i) {
#ifdef __VTK_MODULE
                cellVarsArray->SetComponent(ielem-1, i-1, answer.at(i));
#else
                fprintf( stream, "%e ", answer.at(i) );
#endif
            }
#ifndef __VTK_MODULE
            fprintf( stream, "\n" );
#endif
        }
#ifdef __VTK_MODULE
        if (ncomponents == 1) {
            stream->GetCellData()->SetActiveScalars(__InternalStateTypeToString(type));
            stream->GetCellData()->SetScalars(cellVarsArray);
        } else if (ncomponents == 3) {
            stream->GetCellData()->SetActiveVectors(__InternalStateTypeToString(type));
            stream->GetCellData()->SetVectors(cellVarsArray);
        } else {
            stream->GetCellData()->SetActiveTensors(__InternalStateTypeToString(type));
            stream->GetCellData()->SetTensors(cellVarsArray);
        }
#endif
    }
#ifndef __VTK_MODULE
    fprintf(stream, "</DataArray>\n");
#endif
}

void Tr21Stokes :: computeStiffnessMatrix(FloatMatrix &answer, TimeStep *tStep)
{
    // Note: Working with the components; [K, G+Dp; G^T+Dv^T, C] . [v,p]
    FluidDynamicMaterial *mat = ( FluidDynamicMaterial * ) this->domain->giveMaterial(this->material);
    IntegrationRule *iRule = this->integrationRulesArray [ 0 ];
    GaussPoint *gp;
    FloatMatrix B(3, 12), EdB, K(12,12), G, Dp, DvT, C, Ed, dN;
    FloatArray *lcoords, dN_V(12), Nlin, Ep, Cd, tmpA, tmpB;
    double Cp;

    K.zero();
    G.zero();

    for ( int i = 0; i < iRule->getNumberOfIntegrationPoints(); i++ ) {
        // Compute Gauss point and determinant at current element
        gp = iRule->getIntegrationPoint(i);
        lcoords = gp->giveCoordinates();

        double detJ = fabs(this->interpolation_quad.giveTransformationJacobian(* lcoords, FEIElementGeometryWrapper(this)));
        double dA = detJ * gp->giveWeight();

        this->interpolation_quad.evaldNdx(dN, * lcoords, FEIElementGeometryWrapper(this));
        this->interpolation_lin.evalN(Nlin, * lcoords, FEIElementGeometryWrapper(this));
        for ( int j = 0, k = 0; j < 6; j++, k += 2 ) {
            dN_V(k)     = B(0, k)     = B(2, k + 1) = dN(j, 0);
            dN_V(k + 1) = B(1, k + 1) = B(2, k)     = dN(j, 1);
        }

        // Computing the internal forces should have been done first.
        mat->giveDeviatoricStiffnessMatrix(Ed, TangentStiffness, gp, tStep); // dsigma_dev/deps_dev
        mat->giveDeviatoricPressureStiffness(Ep, TangentStiffness, gp, tStep); // dsigma_dev/dp
        mat->giveVolumetricDeviatoricStiffness(Cd, TangentStiffness, gp, tStep); // deps_vol/deps_dev
        mat->giveVolumetricPressureStiffness(Cp, TangentStiffness, gp, tStep); // deps_vol/dp

        EdB.beProductOf(Ed,B);
        K.plusProductSymmUpper(B, EdB, dA);
        G.plusDyadUnsym(dN_V, Nlin, -dA);
        C.plusDyadSymmUpper(Nlin, Nlin, Cp*dA);

        tmpA.beTProductOf(B, Ep);
        Dp.plusDyadUnsym(tmpA, Nlin, dA);

        tmpB.beTProductOf(B, Cd);
        DvT.plusDyadUnsym(Nlin, tmpB, dA);
    }

    K.symmetrized();
    C.symmetrized();
    FloatMatrix GTDvT, GDp;
    GTDvT.beTranspositionOf(G);
    GTDvT.add(DvT);
    GDp = G;
    GDp.add(Dp);

    FloatMatrix temp(15, 15);
    temp.setSubMatrix(K, 1, 1);
    temp.setSubMatrix(GTDvT, 13, 1);
    temp.setSubMatrix(GDp, 1, 13);
    temp.setSubMatrix(C, 13, 13);

    answer.resize(15, 15);
    answer.zero();
    answer.assemble(temp, this->ordering);
}