C++ GaussPoint类代码示例

void
TrPlaneStrRot3d :: printOutputAt(FILE *file, TimeStep *tStep)
// Performs end-of-step operations.
{
    FloatArray v;

    fprintf( file, "element %d (%8d) :\n", this->giveLabel(), this->giveNumber() );

    for ( int i = 0; i < numberOfIntegrationRules; i++ ) {
        for ( int j = 0; j < integrationRulesArray [ i ]->giveNumberOfIntegrationPoints(); j++ ) {
            GaussPoint *gp = integrationRulesArray [ i ]->getIntegrationPoint(j);

            // gp   -> printOutputAt(file,stepN) ;

            fprintf( file, "  GP %2d.%-2d :", i + 1, gp->giveNumber() );

            this->giveIPValue(v, gp, IST_ShellStrainCurvatureTensor, tStep);
            fprintf(file, "  strains ");
            fprintf( file,
                    " % .4e % .4e % .4e % .4e % .4e % .4e % .4e % .4e % .4e % .4e % .4e % .4e ",
                    v.at(1), v.at(2), v.at(3),  2. * v.at(4), 2. * v.at(5), 2. * v.at(6),
                    v.at(7), v.at(8), v.at(9),  2. * v.at(10), 2. * v.at(11), 2. * v.at(12) );

            this->giveIPValue(v, gp, IST_ShellForceMomentumTensor, tStep);
            fprintf(file, "\n              stresses");
            fprintf( file,
                    " % .4e % .4e % .4e % .4e % .4e % .4e % .4e % .4e % .4e % .4e % .4e % .4e ",
                    v.at(1), v.at(2), v.at(3),  v.at(4), v.at(5), v.at(6),
                    v.at(7), v.at(8), v.at(9),  v.at(10), v.at(11), v.at(12) );

            fprintf(file, "\n");
        }
    }
}

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

double ShapeFunction::det_jacobian(const Element *el, const GaussPoint &g) const
{
  double result;
  if(sfcache[g.order()]->query_jac(el, g, result))
    return result;

  // don't be tempted to rewrite this in terms of
  // det_jacobian(Element*, MasterCoord&) because that one doesn't use
  // the precomputed values of the shape function derivatives!
#if DIM==2
  result = el->jacobian(0, 0, g) * el->jacobian(1, 1, g) -
    el->jacobian(0, 1, g) * el->jacobian(1, 0, g);
#elif DIM==3

  // typing out a closed form in code is messy for 3d
  double m[DIM][DIM]; 
  int ii, jj;
  for(ii=0; ii<DIM; ++ii) {
    for(jj=0; jj<DIM; ++jj) {
      m[ii][jj] = el->jacobian(ii,jj,g);
    }
  }
  result = vtkMath::Determinant3x3(m);

#endif

  sfcache[g.order()]->store_jac(el, g, result);
  return result;
}

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
LIBeam3dNL :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    GaussPoint *gp = this->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0);
    FloatArray nm(6), stress, strain;
    FloatMatrix x;
    double s1, s2;

    // update temp triad
    this->updateTempTriad(tStep);

    if ( useUpdatedGpRecord == 1 ) {
        stress = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() )->giveStrainVector();
    } else {
        this->computeStrainVector(strain, gp, tStep);
        this->computeStressVector(stress, strain, gp, tStep);
    }

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

        nm.at(i)   = s1;
        nm.at(i + 3) = s2;
    }

    this->computeXMtrx(x, tStep);
    answer.beProductOf(x, nm);
}

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

double
Lattice2d :: giveCrackWidth()
{
    GaussPoint *gp = this->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0);
    LatticeMaterialStatus *status = static_cast< LatticeMaterialStatus * >( gp->giveMaterialStatus() );
    double crackWidth = 0;
    crackWidth = status->giveCrackWidth();

    return crackWidth;
}

int
Lattice2d :: giveCrackFlag()
{
    GaussPoint *gp = this->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0);
    LatticeMaterialStatus *status = static_cast< LatticeMaterialStatus * >( gp->giveMaterialStatus() );
    int crackFlag = 0;
    crackFlag = status->giveCrackFlag();

    return crackFlag;
}

double
Lattice2d :: giveDeltaDissipation()
{
    GaussPoint *gp = this->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0);
    LatticeMaterialStatus *status = static_cast< LatticeMaterialStatus * >( gp->giveMaterialStatus() );
    double deltaDissipation = 0;
    deltaDissipation = status->giveDeltaDissipation();

    return deltaDissipation;
}

double
Lattice2d :: giveNormalStress()
{
    GaussPoint *gp = this->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0);
    LatticeMaterialStatus *status = static_cast< LatticeMaterialStatus * >( gp->giveMaterialStatus() );
    double normalStress = 0;
    normalStress = status->giveNormalStress();

    return normalStress;
}

double
Lattice2d_mt :: givePressure()
{
    LatticeTransportMaterialStatus *status;
    IntegrationRule *iRule = this->giveDefaultIntegrationRulePtr();
    GaussPoint *gp = iRule->getIntegrationPoint(0);

    status = static_cast< LatticeTransportMaterialStatus * >( gp->giveMaterialStatus() );

    return status->givePressure();
}

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
FiberedCrossSection :: giveGeneralizedStress_Beam3d(FloatArray &answer, GaussPoint *gp, const FloatArray &strain, TimeStep *tStep)
{
    double fiberThick, fiberWidth, fiberZCoord, fiberYCoord;
    FloatArray fiberStrain, reducedFiberStress;
    StructuralElement *element = static_cast< StructuralElement * >( gp->giveElement() );
    FiberedCrossSectionInterface *interface;

    if ( ( interface = static_cast< FiberedCrossSectionInterface * >( element->giveInterface(FiberedCrossSectionInterfaceType) ) ) == NULL ) {
        OOFEM_ERROR("element with no fiber support encountered");
    }

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

    for ( int i = 1; i <= numberOfFibers; i++ ) {
        GaussPoint *fiberGp = this->giveSlaveGaussPoint(gp, i - 1);
        StructuralMaterial *fiberMat = static_cast< StructuralMaterial * >( domain->giveMaterial( fiberMaterials.at(i) ) );
        // the question is whether this function should exist ?
        // if yes the element details will be hidden.
        // good idea also should be existence of element::GiveBmatrixOfLayer
        // and computing strains here - but first idea looks better
        // but treating of geometric non-linearities may become more complicated
        // another approach - use several functions with assumed kinematic constraints

        // resolve current layer z-coordinate
        fiberThick  = this->fiberThicks.at(i);
        fiberWidth  = this->fiberWidths.at(i);
        fiberYCoord = fiberGp->giveNaturalCoordinate(1);
        fiberZCoord = fiberGp->giveNaturalCoordinate(2);

        interface->FiberedCrossSectionInterface_computeStrainVectorInFiber(fiberStrain, strain, fiberGp, tStep);

        fiberMat->giveRealStressVector_Fiber(reducedFiberStress, fiberGp, fiberStrain, tStep);

        // perform integration
        // 1) membrane terms N, Qz, Qy
        answer.at(1) += reducedFiberStress.at(1) * fiberWidth * fiberThick;
        answer.at(2) += reducedFiberStress.at(2) * fiberWidth * fiberThick;
        answer.at(3) += reducedFiberStress.at(3) * fiberWidth * fiberThick;
        // 2) bending terms mx, my, mxy
        answer.at(4) += ( reducedFiberStress.at(2) * fiberWidth * fiberThick * fiberYCoord -
                          reducedFiberStress.at(3) * fiberWidth * fiberThick * fiberZCoord );
        answer.at(5) += reducedFiberStress.at(1) * fiberWidth * fiberThick * fiberZCoord;
        answer.at(6) -= reducedFiberStress.at(1) * fiberWidth * fiberThick * fiberYCoord;
    }

    // now we must update master gp ///@ todo simply chosen the first fiber material as master material /JB
    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >
                                       ( domain->giveMaterial( fiberMaterials.at(1) )->giveStatus(gp) );
    status->letTempStrainVectorBe(strain);
    status->letTempStressVectorBe(answer);
}