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

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 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 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 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
Quad1MindlinShell3D :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    // We need to overload this for practical reasons (this 3d shell has all 9 dofs, but the shell part only cares for the first 8)
    // This elements adds an additional stiffness for the so called drilling dofs, meaning we need to work with all 9 components.
    FloatMatrix b, d;
    FloatArray n, strain, stress;
    FloatArray shellUnknowns(20), drillUnknowns(4), unknowns;

    this->computeVectorOf(EID_MomentumBalance, VM_Total, tStep, unknowns);
    // Split this for practical reasons into normal shell dofs and drilling dofs
    for ( int i = 0; i < 4; ++i ) {
        shellUnknowns(0 + i*5) = unknowns(0 + i*6);
        shellUnknowns(1 + i*5) = unknowns(1 + i*6);
        shellUnknowns(2 + i*5) = unknowns(2 + i*6);
        shellUnknowns(3 + i*5) = unknowns(3 + i*6);
        shellUnknowns(4 + i*5) = unknowns(4 + i*6);
        drillUnknowns(i) = unknowns(5 + i*6);
    }

    FloatArray shellForces(20), drillMoment(4);
    shellForces.zero();
    drillMoment.zero();
    StructuralCrossSection *cs = this->giveStructuralCrossSection();
    double drillCoeff = cs->give(CS_DrillingStiffness);

    IntegrationRule *iRule = integrationRulesArray [ 0 ];
    for ( int i = 0; i < iRule->giveNumberOfIntegrationPoints(); i++ ) {
        GaussPoint *gp = iRule->getIntegrationPoint(i);
        this->computeBmatrixAt(gp, b);
        double dV = this->computeVolumeAround(gp);

        if ( useUpdatedGpRecord ) {
            stress = static_cast< StructuralMaterialStatus * >( this->giveMaterial()->giveStatus(gp) )->giveStressVector();
        } else {
            strain.beProductOf(b, shellUnknowns);
            cs->giveRealStress_Shell(stress, gp, strain, tStep);
        }
        shellForces.plusProduct(b, stress, dV);

        // Drilling stiffness is here for improved numerical properties
        if (drillCoeff > 0.) {
            this->interp.evalN(n, *gp->giveCoordinates(), FEIVoidCellGeometry());
            for ( int j = 0; j < 4; j++) {
                n(j) -= 0.25;
            }
            double dtheta = n.dotProduct(drillUnknowns);
            drillMoment.add(drillCoeff * dV * dtheta, n); ///@todo Decide on how to alpha should be defined.
        }
    }

    answer.resize(24);
    answer.zero();
    answer.assemble(shellForces, this->shellOrdering);

    if (drillCoeff > 0.) {
        answer.assemble(drillMoment, this->drillOrdering);
    }
}

void
Quad1MindlinShell3D :: computeBodyLoadVectorAt(FloatArray &answer, Load *forLoad, TimeStep *stepN, ValueModeType mode)
{
    // Only gravity load
    double dV, density;
    GaussPoint *gp;
    FloatArray forceX, forceY, forceZ, glob_gravity, gravity, n;

    if ( ( forLoad->giveBCGeoType() != BodyLoadBGT ) || ( forLoad->giveBCValType() != ForceLoadBVT ) ) {
        _error("computeBodyLoadVectorAt: unknown load type");
    }

    // note: force is assumed to be in global coordinate system.
    forLoad->computeComponentArrayAt(glob_gravity, stepN, mode);
    // Transform the load into the local c.s.
    gravity.beProductOf(this->lcsMatrix, glob_gravity); ///@todo Check potential transpose here.

    if ( gravity.giveSize() ) {
        IntegrationRule *ir = integrationRulesArray [ 0 ];
        for ( int i = 0; i < ir->giveNumberOfIntegrationPoints(); ++i) {
            gp = ir->getIntegrationPoint(i);

            this->interp.evalN(n, *gp->giveCoordinates(), FEIVoidCellGeometry());
            dV = this->computeVolumeAround(gp) * this->giveCrossSection()->give(CS_Thickness);
            density = this->giveMaterial()->give('d', gp);

            forceX.add(density * gravity.at(1) * dV, n);
            forceY.add(density * gravity.at(2) * dV, n);
            forceZ.add(density * gravity.at(3) * dV, n);
        }

        answer.resize(24);
        answer.zero();

        answer.at(1)  = forceX.at(1);
        answer.at(2)  = forceY.at(1);
        answer.at(3)  = forceZ.at(1);

        answer.at(7)  = forceX.at(2);
        answer.at(8)  = forceY.at(2);
        answer.at(9)  = forceZ.at(2);

        answer.at(13) = forceX.at(3);
        answer.at(14) = forceY.at(3);
        answer.at(15) = forceZ.at(3);

        answer.at(19) = forceX.at(4);
        answer.at(20) = forceY.at(4);
        answer.at(21) = forceZ.at(4);

    } else {
        answer.resize(0);
    }
}

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
Quad1MindlinShell3D :: computeSurfaceLoadVectorAt(FloatArray &answer, Load *load,
                                                int iSurf, TimeStep *tStep, ValueModeType mode)
{
    BoundaryLoad *surfLoad = static_cast< BoundaryLoad * >(load);
    if ( dynamic_cast< ConstantPressureLoad * >(surfLoad) ) { // Just checking the type of b.c.
        // EXPERIMENTAL CODE:
        IntegrationRule *iRule;
        FloatArray n, gcoords, pressure;

        answer.resize( 24 );
        answer.zero();

        //int approxOrder = surfLoad->giveApproxOrder() + this->giveApproxOrder();

        iRule = this->integrationRulesArray[ 0 ];
        for ( int i = 0; i < iRule->giveNumberOfIntegrationPoints(); i++ ) {
            GaussPoint *gp = iRule->getIntegrationPoint(i);
            double dV = this->computeVolumeAround(gp);
            this->interp.evalN(n, *gp->giveCoordinates(), FEIVoidCellGeometry());
            this->interp.local2global(gcoords, *gp->giveCoordinates(), FEIElementGeometryWrapper(this));
            surfLoad->computeValueAt(pressure, tStep, gcoords, mode);

            answer.at( 3) += n.at(1) * pressure.at(1) * dV;
            answer.at( 9) += n.at(2) * pressure.at(1) * dV;
            answer.at(15) += n.at(3) * pressure.at(1) * dV;
            answer.at(21) += n.at(4) * pressure.at(1) * dV;
        }
        // Second surface is the outside;
        if ( iSurf == 2 ) {
            answer.negated();
        }
    } else {
        OOFEM_ERROR("Quad1MindlinShell3D only supports constant pressure boundary load.");
    }
}

void Tr1Darcy :: computeEdgeBCSubVectorAt(FloatArray &answer, Load *load, int iEdge, TimeStep *tStep)
{
    /*
     * Given the load *load, return it's contribution.
     *
     */

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

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

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

        GaussIntegrationRule iRule(1, this, 1, 1);
        GaussPoint *gp;
        FloatArray N, loadValue, reducedAnswer;
        reducedAnswer.resize(3);
        reducedAnswer.zero();
        IntArray mask;

        iRule.setUpIntegrationPoints(_Line, numberOfEdgeIPs, _Unknown);

        for ( int i = 0; i < iRule.getNumberOfIntegrationPoints(); i++ ) {
            gp = iRule.getIntegrationPoint(i);
            FloatArray *lcoords = gp->giveCoordinates();
            this->interpolation_lin.edgeEvalN(N, *lcoords, FEIElementGeometryWrapper(this));
            double dV = this->computeEdgeVolumeAround(gp, iEdge);

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

            reducedAnswer.add(loadValue.at(1) * dV, N);
        }

        this->interpolation_lin.computeLocalEdgeMapping(mask, iEdge);
        answer.assemble(reducedAnswer, mask);
    }
}

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

void
Quad1MindlinShell3D :: computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    // We need to overload this for practical reasons (this 3d shell has all 9 dofs, but the shell part only cares for the first 8)
    // This elements adds an additional stiffness for the so called drilling dofs, meaning we need to work with all 9 components.
    FloatMatrix d, b, db;
    FloatArray n;

    FloatMatrix shellStiffness(20, 20), drillStiffness(4, 4);
    shellStiffness.zero();
    drillStiffness.zero();
    double drillCoeff = this->giveStructuralCrossSection()->give(CS_DrillingStiffness);

    IntegrationRule *iRule = integrationRulesArray [ 0 ];
    for ( int i = 0; i < iRule->giveNumberOfIntegrationPoints(); i++ ) {
        GaussPoint *gp = iRule->getIntegrationPoint(i);
        this->computeBmatrixAt(gp, b);
        double dV = this->computeVolumeAround(gp);

        this->computeConstitutiveMatrixAt(d, rMode, gp, tStep);

        db.beProductOf(d, b);
        shellStiffness.plusProductSymmUpper(b, db, dV);

        // Drilling stiffness is here for improved numerical properties
        if (drillCoeff > 0.) {
            this->interp.evalN(n, *gp->giveCoordinates(), FEIVoidCellGeometry());
            for ( int j = 0; j < 4; j++) {
                n(j) -= 0.25;
            }
            drillStiffness.plusDyadSymmUpper(n, drillCoeff * dV);
        }
    }
    shellStiffness.symmetrized();

    answer.resize(24, 24);
    answer.zero();
    answer.assemble(shellStiffness, this->shellOrdering);

    if (drillCoeff > 0.) {
        drillStiffness.symmetrized();
        answer.assemble(drillStiffness, this->drillOrdering);
    }
}

void
Quad1Mindlin :: computeBodyLoadVectorAt(FloatArray &answer, Load *forLoad, TimeStep *stepN, ValueModeType mode)
{
    // Only gravity load
    double dV, load;
    GaussPoint *gp;
    FloatArray force, gravity, n;

    if ( ( forLoad->giveBCGeoType() != BodyLoadBGT ) || ( forLoad->giveBCValType() != ForceLoadBVT ) ) {
        _error("computeBodyLoadVectorAt: unknown load type");
    }

    // note: force is assumed to be in global coordinate system.
    forLoad->computeComponentArrayAt(gravity, stepN, mode);

    force.resize(0);
    if ( gravity.giveSize() ) {
        IntegrationRule *ir = integrationRulesArray [ 0 ]; ///@todo Other/higher integration for lumped mass matrices perhaps?
        for ( int i = 0; i < ir->getNumberOfIntegrationPoints(); ++i) {
            gp = ir->getIntegrationPoint(i);

            this->interp_lin.evalN(n, *gp->giveCoordinates(), FEIElementGeometryWrapper(this));
            dV = this->computeVolumeAround(gp) * this->giveCrossSection()->give(CS_Thickness);
            load = this->giveMaterial()->give('d', gp) * gravity.at(3) * dV;

            force.add(load, n);
        }

        answer.resize(12);
        answer.zero();

        answer.at(1)  = force.at(1);
        answer.at(4)  = force.at(2);
        answer.at(7)  = force.at(3);
        answer.at(10) = force.at(4);

    } else {
        answer.resize(0);
    }
}

void
Lattice2d_mt :: updateInternalState(TimeStep *stepN)
// Updates the receiver at end of step.
{
    int i, j;
    IntegrationRule *iRule;
    FloatArray f, r;
    FloatMatrix n;
    TransportMaterial *mat = ( ( TransportMaterial * ) this->giveMaterial() );
    GaussPoint *gp;

    // force updating ip values
    for ( i = 0; i < numberOfIntegrationRules; i++ ) {
        iRule = integrationRulesArray [ i ];
        for ( j = 0; j < iRule->giveNumberOfIntegrationPoints(); j++ ) {
            gp = iRule->getIntegrationPoint(j);
            this->computeNmatrixAt( n, *gp->giveCoordinates() );
            this->computeVectorOf(EID_ConservationEquation, VM_Total, stepN, r);
            f.beProductOf(n, r);
            mat->updateInternalState(f, gp, stepN);
        }
    }
}

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

}

//.........这里部分代码省略.........
         *        }
         *
         *        // this fixes a bug in ELIXIR
         *        if ( fabs(s [ i ]) < 1.0e-6 ) {
         *            s [ i ] = 1.0e-6;
         *        }
         *    }
         *
         *    if ( context.getScalarAlgo() == SA_ZPROFILE ) {
         *        EASValsSetColor( context.getDeformedElementColor() );
         *        EASValsSetLineWidth(OOFEG_DEFORMED_GEOMETRY_WIDTH);
         *        tr =  CreateQuad3D(p);
         *        EGWithMaskChangeAttributes(WIDTH_MASK | COLOR_MASK | LAYER_MASK, tr);
         *    } else {
         *        context.updateFringeTableMinMax(s, 4);
         *        tr =  CreateQuadWD3D(p, s [ 0 ], s [ 1 ], s [ 2 ], s [ 3 ]);
         *        EGWithMaskChangeAttributes(LAYER_MASK, tr);
         *    }
         *
         *    EMAddGraphicsToModel(ESIModel(), tr);
         * }
         */
    } else if ( context.giveIntVarMode() == ISM_local ) {
        // ========== plot the local values (raw data) =====================
        if ( numberOfGaussPoints != 4 ) {
            return;
        }

        int ip;
        GaussPoint *gp;
        IntArray ind(4);
        FloatArray *gpCoords;
        WCRec pp [ 9 ];

        for ( i = 0; i < 8; i++ ) {
            if ( context.getInternalVarsDefGeoFlag() ) {
                // use deformed geometry
                defScale = context.getDefScale();
                pp [ i ].x = ( FPNum ) this->giveNode(i + 1)->giveUpdatedCoordinate(1, tStep, defScale);
                pp [ i ].y = ( FPNum ) this->giveNode(i + 1)->giveUpdatedCoordinate(2, tStep, defScale);
                pp [ i ].z = 0.;
            } else {
                pp [ i ].x = ( FPNum ) this->giveNode(i + 1)->giveCoordinate(1);
                pp [ i ].y = ( FPNum ) this->giveNode(i + 1)->giveCoordinate(2);
                pp [ i ].z = 0.;
            }
        }

        pp [ 8 ].x = 0.25 * ( pp [ 0 ].x + pp [ 1 ].x + pp [ 2 ].x + pp [ 3 ].x );
        pp [ 8 ].y = 0.25 * ( pp [ 0 ].y + pp [ 1 ].y + pp [ 2 ].y + pp [ 3 ].y );
        pp [ 8 ].z = 0.;

        for ( ip = 1; ip <= integrationRulesArray [ 0 ]->giveNumberOfIntegrationPoints(); ip++ ) {
            gp = integrationRulesArray [ 0 ]->getIntegrationPoint(ip - 1);
            gpCoords = gp->giveCoordinates();
            if ( ( gpCoords->at(1) > 0. ) && ( gpCoords->at(2) > 0. ) ) {
                ind.at(1) = 0;
                ind.at(2) = 4;
                ind.at(3) = 8;
                ind.at(4) = 7;
            } else if ( ( gpCoords->at(1) < 0. ) && ( gpCoords->at(2) > 0. ) ) {
                ind.at(1) = 4;
                ind.at(2) = 1;
                ind.at(3) = 5;
                ind.at(4) = 8;
            } else if ( ( gpCoords->at(1) < 0. ) && ( gpCoords->at(2) < 0. ) ) {
                ind.at(1) = 5;
                ind.at(2) = 2;
                ind.at(3) = 6;
                ind.at(4) = 8;
            } else {
                ind.at(1) = 6;
                ind.at(2) = 3;
                ind.at(3) = 7;
                ind.at(4) = 8;
            }

            if ( giveIPValue(v [ 0 ], gp, context.giveIntVarType(), tStep) == 0 ) {
                return;
            }

            indx = context.giveIntVarIndx();

            for ( i = 1; i <= 4; i++ ) {
                s [ i - 1 ] = v [ 0 ].at(indx);
            }

            for ( i = 0; i < 4; i++ ) {
                p [ i ].x = pp [ ind.at(i + 1) ].x;
                p [ i ].y = pp [ ind.at(i + 1) ].y;
                p [ i ].z = pp [ ind.at(i + 1) ].z;
            }

            context.updateFringeTableMinMax(s, 4);
            tr =  CreateQuadWD3D(p, s [ 0 ], s [ 1 ], s [ 2 ], s [ 3 ]);
            EGWithMaskChangeAttributes(LAYER_MASK, tr);
            EMAddGraphicsToModel(ESIModel(), tr);
        }
    }
}