C++ FloatArray::computeNorm方法代码示例

void FE2FluidMaterial :: giveDeviatoricPressureStiffness(FloatArray &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    FE2FluidMaterialStatus *ms = static_cast<FE2FluidMaterialStatus*> (this->giveStatus(gp));
    ms->computeTangents(tStep);
    if ( mode == TangentStiffness ) {
        answer = ms->giveDeviatoricPressureTangent();
#ifdef DEBUG_TANGENT
        // 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 = "); answer.printYourself();
        printf("Numerical deviatoric pressure tangent = "); dsigh.printYourself();
        dsigh.subtract(answer);
        double norm = dsigh.computeNorm();
        if (norm > answer.computeNorm()*DEBUG_ERR && norm > 0.0) {
            OOFEM_ERROR("Error in deviatoric pressure tangent");
        }
#endif
    } else {
        OOFEM_ERROR("Mode not implemented");
    }
}

void
LIBeam3dNL ::  computeRotMtrx(FloatMatrix &answer, FloatArray &psi)
{
    FloatMatrix S(3, 3), SS(3, 3);
    double psiSize;

    if ( psi.giveSize() != 3 ) {
        _error("computeSMtrx: psi param size mismatch");
    }

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

    psiSize = psi.computeNorm();
    answer.at(1, 1) = answer.at(2, 2) = answer.at(3, 3) = 1.;

    if ( psiSize <= 1.e-40 ) {
        return;
    }

    this->computeSMtrx(S, psi);
    SS.beProductOf(S, S);
    S.times(sin(psiSize) / psiSize);
    SS.times( ( 1. - cos(psiSize) ) / ( psiSize * psiSize ) );

    answer.add(S);
    answer.add(SS);
}

double
ScalarErrorIndicator :: giveElementError(EE_ErrorType type, Element *elem, TimeStep *tStep)
{
    FloatArray val;
    int result = 1;
    double sval, maxVal = 0.0;

    if ( type != indicatorET ) {
        return 0.0;
    }

    if ( this->skipRegion( elem->giveRegionNumber() ) ) {
        return 0.0;
    }

    for ( GaussPoint *gp: *elem->giveDefaultIntegrationRulePtr() ) {
        result = elem->giveIPValue(val, gp, varType, tStep);
        if ( result ) {
            sval = val.computeNorm();
            if ( gp->giveNumber() == 1 ) {
                maxVal = sval;
            } else {
                maxVal = max(maxVal, sval);
            }
        }
    }

    return maxVal;
}

void 
IntMatCoulombContact :: computeEngTraction(double &normalStress, FloatArray &shearStress,
                                         FloatArray &tempShearStressShift, double normalJump, const FloatArray &shearJump )
{

    double maxShearStress = 0.0;
    double shift = -this->kn * this->stiffCoeff * normalClearance;

    if ( normalJump + normalClearance <= 0. ) {
        normalStress = this->kn * ( normalJump + normalClearance ) + shift; //in compression and after the clearance gap closed
        maxShearStress = fabs( normalStress ) * this->frictCoeff;
    } else {
        normalStress = this->kn * this->stiffCoeff * ( normalJump + normalClearance ) + shift;
        maxShearStress = 0.;
    }

    if ( shearJump.giveSize() != 0 ) {

        shearStress = this->kn * shearJump - tempShearStressShift;
        double dp = shearStress.computeNorm();
        double eps = 1.0e-15; // small number
        if ( dp > maxShearStress ) {
            shearStress.times( maxShearStress / ( dp + eps ) );
        }
        tempShearStressShift = this->kn * shearJump - shearStress;

    } else { // 1d -> no shear stresses
        return;
    }
}

double Line :: computeTangentialDistanceToEnd(FloatArray *point)
{
    FloatArray projection;
    this->computeProjection(projection);
    FloatArray tmp;
    tmp.beDifferenceOf(* point, mVertices [ 1 ]);
    return tmp.dotProduct(projection) / projection.computeNorm();
}

double LineSurfaceTension :: SpatialLocalizerI_giveDistanceFromParametricCenter(const FloatArray &coords)
{
    FloatArray c;
    c = *this->giveNode(1)->giveCoordinates();
    c.add(*this->giveNode(2)->giveCoordinates());
    c.times(0.5);
    c.subtract(coords);
    return c.computeNorm();
}

double
FEI3dTrQuad :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    FloatArray G1, G2; // local curvilinear base vectors
    this->surfaceEvalBaseVectorsAt(G1, G2, lcoords, cellgeo);
    answer.beVectorProductOf(G1, G2);
    double J = answer.computeNorm();
    answer.times(1 / J);
    return J;
}

void
NonlinearMassTransferMaterial :: giveFluxVector(FloatArray &answer, GaussPoint *gp, const FloatArray &grad, const FloatArray &field, TimeStep *tStep)
{
    TransportMaterialStatus *ms = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );

    double gradPNorm = grad.computeNorm();
    answer.beScaled( -(1. + C * pow(gradPNorm, alpha)), grad);

    ms->setTempGradient(grad);
    ms->setTempField(field);
    ms->setTempFlux(answer);
}

void Line :: transformIntoPolar(FloatArray *point, FloatArray &answer)
{
    FloatArray xp;
    FloatMatrix Qt;
    FloatArray help;
    this->computeTransformationMatrix(Qt);
    help.beDifferenceOf(* point, mVertices [ 1 ]);
    xp.beProductOf(Qt, help);
    answer.resize(2);
    answer.at(1) = xp.computeNorm();
    answer.at(2) = atan2( xp.at(2), xp.at(1) );
}

void
NonlinearMassTransferMaterial :: giveFluxVector(FloatArray &answer, GaussPoint *gp, const FloatArray &eps, TimeStep *tStep)
{
    AnisotropicMassTransferMaterialStatus *thisMaterialStatus;
    thisMaterialStatus = ( ( AnisotropicMassTransferMaterialStatus * ) this->giveStatus(gp) );

    thisMaterialStatus->setPressureGradient(eps);

    double gradPNorm = eps.computeNorm();

    answer = eps;
    answer.times( 1 + C * pow(gradPNorm, alpha) );
    answer.negated();

    thisMaterialStatus->setSeepageValocity(answer);
}

double
ParallelContext :: localNorm(const FloatArray &src)
{
#ifdef __PARALLEL_MODE
    if ( emodel->isParallel() ) {
        double norm2 = 0.0, norm2_tot;
        int size = src.giveSize();
        Natural2LocalOrdering *n2l = this->giveN2Lmap();
        for ( int i = 0; i < size; i++ ) {
            if ( n2l->giveNewEq(i + 1) ) {
                norm2 += src(i) * src(i);
            }
        }
        MPI_Allreduce( & norm2, & norm2_tot, 1, MPI_DOUBLE, MPI_SUM, this->emodel->giveParallelComm() );
        return sqrt(norm2_tot);
    }
#endif
    return src.computeNorm();
}

void
NonlinearMassTransferMaterial :: giveCharacteristicMatrix(FloatMatrix &answer,
                                                          MatResponseForm form,
                                                          MatResponseMode mode,
                                                          GaussPoint *gp,
                                                          TimeStep *atTime)
{
    MaterialMode mMode = gp->giveMaterialMode();
    AnisotropicMassTransferMaterialStatus *status = ( ( AnisotropicMassTransferMaterialStatus * ) this->giveStatus(gp) );
    FloatArray eps = status->giveGradP();
    double gradPNorm;
    FloatMatrix t1, t2;

    gradPNorm = eps.computeNorm();

    t1.beDyadicProductOf(eps, eps);
    if ( gradPNorm != 0.0 ) {
        t1.times( C * alpha * pow(gradPNorm, alpha - 2) );
    }

    switch  ( mMode ) {
    case _1dHeat:
        t2.resize(1, 1);
        t2.at(1, 1) = 1;
        break;
    case _2dHeat:
        t2.resize(2, 2);
        t2.at(1, 1) = t2.at(2, 2) = 1;
        break;
    case _3dHeat:
        t2.resize(3, 3);
        t2.at(1, 1) = t2.at(2, 2) = t2.at(3, 3) = 1;
        break;
    default:
        _error2( "giveCharacteristicMatrix : unknown mode (%s)", __MaterialModeToString(mMode) );
    }
    
    answer.beEmptyMtrx();
    answer.add(t1);
    answer.add(1 + C * pow(gradPNorm, alpha), t2);

}

void FE2FluidMaterial :: giveVolumetricDeviatoricStiffness(FloatArray &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    FE2FluidMaterialStatus *ms = static_cast<FE2FluidMaterialStatus*> (this->giveStatus(gp));
    ms->computeTangents(tStep);
    if ( mode == TangentStiffness ) {
        answer = ms->giveVolumetricDeviatoricTangent();
#ifdef DEBUG_TANGENT
        // 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 = "); answer.printYourself();
        printf("Numerical volumetric deviatoric tangent = "); dvol.printYourself();
        dvol.subtract(answer);
        double norm = dvol.computeNorm();
        if (norm > answer.computeNorm()*DEBUG_ERR && norm > 0.0) {
            OOFEM_ERROR("Error in volumetric deviatoric tangent");
        }
#endif
    } else {
        OOFEM_ERROR("Mode not implemented");
    }
}

double
FEI3dTetQuad :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    IntArray snodes(3);
    FloatArray a,b;
    this->computeLocalSurfaceMapping(snodes, isurf);

    double l1, l2, l3;
    l1 = lcoords.at(1);
    l2 = lcoords.at(2);
    l3 = 1.0 - l1 - l2;

    FloatArray dNdxi(6), dNdeta(6);

    dNdxi(0) = 4.0 * l1 - 1.0;
    dNdxi(1) = 0.0;
    dNdxi(2) = -1.0 * ( 4.0 * l3 - 1.0 );
    dNdxi(3) = 4.0 * l2;
    dNdxi(4) = -4.0 * l2;
    dNdxi(5) = 4.0 * l3 - 4.0 * l1;

    dNdeta(0) = 0.0;
    dNdeta(1) = 4.0 * l2 - 1.0;
    dNdeta(2) = -1.0 * ( 4.0 * l3 - 1.0 );
    dNdeta(3) = 4.0 * l1;
    dNdeta(4) = 4.0 * l3 - 4.0 * l2;
    dNdeta(5) = -4.0 * l1;

    for (int i = 0; i < 6; ++i) {
        a.add(dNdxi(i),  *cellgeo.giveVertexCoordinates(snodes(i)));
        b.add(dNdeta(i), *cellgeo.giveVertexCoordinates(snodes(i)));
    }
    answer.beVectorProductOf(a, b);
    double J = answer.computeNorm();
    answer.times(1/J);
    return J;
}

//.........这里部分代码省略.........
    //create previous solution from IC or from previous tStep
    if ( tStep->isTheFirstStep() ) {
        if ( !stepWhenIcApply ) {
            stepWhenIcApply.reset( new TimeStep( *tStep->givePreviousStep() ) );
        }
        this->applyIC(stepWhenIcApply.get()); //insert solution to hash=1(previous), if changes in equation numbering
    }

    double dTTau = tStep->giveTimeIncrement();
    double Tau = tStep->giveTargetTime() - ( 1. - alpha ) * tStep->giveTimeIncrement();
    //Time step in which material laws are taken into account
    TimeStep TauStep(tStep->giveNumber(), this, tStep->giveMetaStepNumber(), Tau, dTTau, tStep->giveSolutionStateCounter() + 1);

    //Predictor
    FloatArray *solutionVector;
    UnknownsField->advanceSolution(tStep);
    solutionVector = UnknownsField->giveSolutionVector(tStep);

    //Initialize and give solutionVector from previous solution
    if ( changingProblemSize ) {
        if ( !tStep->isTheFirstStep() ) {
            //copy recent solution to previous position, copy from hash=0 to hash=1(previous)
            copyUnknownsInDictionary( VM_Total, tStep, tStep->givePreviousStep() );
        }

        UnknownsField->initialize( VM_Total, tStep->givePreviousStep(), *solutionVector, EModelDefaultEquationNumbering() );
    } else {
        //copy previous solution vector to actual
        *solutionVector = *UnknownsField->giveSolutionVector( tStep->givePreviousStep() );
    }

    this->updateInternalState(& TauStep); //insert to hash=0(current), if changes in equation numbering

    FloatArray solutionVectorIncrement(neq);
    int nite = 0;

    OOFEM_LOG_INFO("Time            Iter       ResidNorm       IncrNorm\n__________________________________________________________\n");


    do {
        nite++;

        // Corrector
#ifdef VERBOSE
        // printf("\nAssembling conductivity and capacity matrices");
#endif

        if ( ( nite == 1 ) || ( NR_Mode == nrsolverFullNRM ) || ( ( NR_Mode == nrsolverAccelNRM ) && ( nite % MANRMSteps == 0 ) ) ) {
            conductivityMatrix->zero();
            //Assembling left hand side - start with conductivity matrix
            this->assemble( *conductivityMatrix, & TauStep, IntSourceLHSMatrix,
                           EModelDefaultEquationNumbering(), this->giveDomain(1) );
            conductivityMatrix->times(alpha);
            //Add capacity matrix
            this->assemble( *conductivityMatrix, & TauStep, NSTP_MidpointLhs,
                           EModelDefaultEquationNumbering(), this->giveDomain(1) );
        }

        rhs.resize(neq);
        rhs.zero();
        //edge or surface load on element
        this->assembleVectorFromElements( rhs, & TauStep, ElementBCTransportVector, VM_Total,
                                         EModelDefaultEquationNumbering(), this->giveDomain(1) );
        //add internal source vector on elements
        this->assembleVectorFromElements( rhs, & TauStep, ElementInternalSourceVector, VM_Total,
                                         EModelDefaultEquationNumbering(), this->giveDomain(1) );
        //add nodal load
        this->assembleVectorFromDofManagers( rhs, & TauStep, ExternalForcesVector, VM_Total,
                                            EModelDefaultEquationNumbering(), this->giveDomain(1) );

        // subtract the rhs part depending on previous solution
        assembleAlgorithmicPartOfRhs(rhs, EModelDefaultEquationNumbering(), & TauStep);
        // set-up numerical model
        this->giveNumericalMethod( this->giveCurrentMetaStep() );

        // call numerical model to solve arised problem
#ifdef VERBOSE
        //OOFEM_LOG_INFO("Solving ...\n");
#endif

        // compute norm of residuals from balance equations
        solutionErr = rhs.computeNorm();

        linSolver->solve(*conductivityMatrix, rhs, solutionVectorIncrement);
        solutionVector->add(solutionVectorIncrement);
        this->updateInternalState(& TauStep); //insert to hash=0(current), if changes in equation numbering
        // compute error in the solutionvector increment
        incrementErr = solutionVectorIncrement.computeNorm();

        // update solution state counter
        TauStep.incrementStateCounter();
        tStep->incrementStateCounter();

        OOFEM_LOG_INFO("%-15e %-10d %-15e %-15e\n", tStep->giveTargetTime(), nite, solutionErr, incrementErr);

        if ( nite >= nsmax ) {
            OOFEM_ERROR("convergence not reached after %d iterations", nsmax);
        }
    } while ( ( fabs(solutionErr) > rtol ) || ( fabs(incrementErr) > rtol ) );
}