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

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

bool
BeamElementErrorCheckingRule :: check(Domain *domain, TimeStep *tStep)
{
    // Rule doesn't apply yet.
    if ( tStep->giveNumber() != tstep ) {
        return true;
    }

    FloatArray val;
    Element *element = domain->giveGlobalElement(number);
    if ( !element ) {
        if ( domain->giveEngngModel()->isParallel() ) {
            return true;
        } else {
            OOFEM_WARNING("Element %d not found.", number);
            return false;
        }
    }
    if ( element->giveParallelMode() != Element_local ) {
        return true;
    }

    if (ist == BET_localEndDisplacement) {
        element->computeVectorOf(VM_Total, tStep, val);
    } else if (ist ==  BET_localEndForces) {
        if(Beam2d* b = dynamic_cast<Beam2d*>(element)) b->giveEndForcesVector(val, tStep);
        else if(Beam3d* b = dynamic_cast<Beam3d*>(element)) b->giveEndForcesVector(val, tStep);
        else {
            OOFEM_WARNING("Element %d has no beam interface.", number);
            return false;
        }
    }

    if ( component > val.giveSize() || component < 1 ) {
        OOFEM_WARNING("Check failed in: beam_element %d, ist %d, component %d:\n"
                      "Component not found!",
                      number, ist, component);
        val.printYourself();
        return false;
    }

    double elementValue = val.at(component);
    bool check = checkValue(elementValue);
    if ( !check ) {
        OOFEM_WARNING("Check failed in: tstep %d, beam_element %d, ist %d, component %d:\n"
                      "value is %.8e, but should be %.8e ( error is %e but tolerance is %e )",
                      tstep, number, ist, component,
                      elementValue, value, fabs(elementValue-value), tolerance );
        val.printYourself();
    }
    return check;
}

void PrescribedGradientBCWeak :: giveMirroredPointOnGammaMinus(FloatArray &oPosMinus, const FloatArray &iPosPlus) const
{
//#if 0
	if(mMirrorFunction == 0) {
		oPosMinus = iPosPlus;
		const double distTol = 1.0e-12;

//		if ( iPosPlus.distance(mUC) < distTol ) {
//			printf("iPosPlus: %.12e %.12e\n", iPosPlus [ 0 ], iPosPlus [ 1 ]);
//			OOFEM_ERROR("Unmappable point.")
//		}

		bool mappingPerformed = false;

		if ( iPosPlus [ 0 ] > mUC [ 0 ] - distTol ) {
			oPosMinus [ 0 ] = mLC [ 0 ];
			mappingPerformed = true;
			return;
		}

		if ( iPosPlus [ 1 ] > mUC [ 1 ] - distTol ) {
			oPosMinus [ 1 ] = mLC [ 1 ];
			mappingPerformed = true;
			return;
		}

		if ( !mappingPerformed ) {
			iPosPlus.printYourself();
			OOFEM_ERROR("Mapping failed.")
		}

void GnuplotExportModule::outputBoundaryCondition(PrescribedGradientBCNeumann &iBC, TimeStep *tStep)
{
    FloatArray stress;
    iBC.computeField(stress, tStep);

    printf("Mean stress computed in Gnuplot export module: "); stress.printYourself();

    double time = 0.0;

    TimeStep *ts = emodel->giveCurrentStep();
    if ( ts != NULL ) {
        time = ts->giveTargetTime();
    }

    int bcIndex = iBC.giveNumber();

    std :: stringstream strMeanStress;
    strMeanStress << "PrescribedGradientGnuplotMeanStress" << bcIndex << "Time" << time << ".dat";
    std :: string nameMeanStress = strMeanStress.str();
    std::vector<double> componentArray, stressArray;

    for(int i = 1; i <= stress.giveSize(); i++) {
        componentArray.push_back(i);
        stressArray.push_back(stress.at(i));
    }

    XFEMDebugTools::WriteArrayToGnuplot(nameMeanStress, componentArray, stressArray);

    // Homogenized strain
    
    FloatArray grad;
    iBC.giveGradientVoigt(grad);
    outputGradient(iBC.giveNumber(), *iBC.giveDomain(), grad, tStep);
}

bool
ElementErrorCheckingRule :: check(Domain *domain, TimeStep *tStep)
{
    // Rule doesn't apply yet.
    if ( tStep->giveNumber() != tstep ) {
        return true;
    }

    FloatArray ipval;
    Element *element = domain->giveGlobalElement(number);
    if ( !element ) {
        if ( domain->giveEngngModel()->isParallel() ) {
            return true;
        } else {
            OOFEM_WARNING("Element %d not found.", number);
            return false;
        }
    }
    if ( element->giveParallelMode() != Element_local ) {
        return true;
    }

    // note! GPs are numbered from 0 internally, but written with 1-index, inconsistent!
    GaussPoint *gp = element->giveIntegrationRule(irule)->getIntegrationPoint(gpnum-1);
    element->giveIPValue(ipval, gp, ist, tStep);

    if ( component > ipval.giveSize() || component < 1 ) {
        OOFEM_WARNING("Check failed in: element %d, gpnum %d, ist %d, component %d:\n"
                      "Component not found!",
                      number, gpnum, ist, component);
        ipval.printYourself();
        return false;
    }

    double elementValue = ipval.at(component);
    bool check = checkValue(elementValue);
    if ( !check ) {
        OOFEM_WARNING("Check failed in: tstep %d, element %d, gpnum %d, ist %d, component %d:\n"
                      "value is %.8e, but should be %.8e ( error is %e but tolerance is %e )",
                      tstep, number, gpnum, ist, component,
                      elementValue, value, fabs(elementValue-value), tolerance );
        ipval.printYourself();
    }
    return check;
}

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

void GnuplotExportModule::outputGradient(int bc, Domain &d, FloatArray &grad, TimeStep *tStep)
{
    // Homogenized strain
    double timeFactor = d.giveBc(bc)->giveTimeFunction()->evaluateAtTime(tStep->giveTargetTime());
    printf("timeFactor: %e\n", timeFactor );
    grad.times(timeFactor);
    printf("Mean grad computed in Gnuplot export module: "); grad.printYourself();

    double time = tStep->giveTargetTime();


    std :: stringstream strMeanGrad;
    strMeanGrad << "PrescribedGradientGnuplotMeanGrad" << bc << "Time" << time << ".dat";
    std :: string nameMeanGrad = strMeanGrad.str();
    std::vector<double> componentArrayGrad, gradArray;

    for(int i = 1; i <= grad.giveSize(); i++) {
        componentArrayGrad.push_back(i);
        gradArray.push_back(grad.at(i));
    }

    XFEMDebugTools::WriteArrayToGnuplot(nameMeanGrad, componentArrayGrad, gradArray);
}

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();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                            a.at(j, i) = aj + cg * ak;
                            b.at(j, i) = bj + cg * bk;
                            a.at(i, k) = ak + ca * aj;
                            b.at(i, k) = bk + ca * bj;
                        }
                    }
                }                           // label 190

                ak = a.at(k, k);
                bk = b.at(k, k);
                a.at(k, k) = ak + 2.0 *ca *a.at(j, k) + ca *ca *a.at(j, j);
                b.at(k, k) = bk + 2.0 *ca *b.at(j, k) + ca *ca *b.at(j, j);
                a.at(j, j) = a.at(j, j) + 2.0 *cg *a.at(j, k) + cg * cg * ak;
                b.at(j, j) = b.at(j, j) + 2.0 *cg *b.at(j, k) + cg * cg * bk;
                a.at(j, k) = 0.0;
                b.at(j, k) = 0.0;
                //
                // update the eigenvector matrix after each rotation
                //
                for ( int i = 1; i <= n; i++ ) {
                    xj = x.at(i, j);
                    xk = x.at(i, k);
                    x.at(i, j) = xj + cg * xk;
                    x.at(i, k) = xk + ca * xj;
                }                        // label 200
            }
        }                                // label 210

        //
        // update the eigenvalues after each sweep
        //
#ifdef DETAILED_REPORT
        OOFEM_LOG_DEBUG("GJacobi: a,b after sweep\n");
        a.printYourself();
        b.printYourself();
#endif
        for ( int i = 1; i <= n; i++ ) {
            // in original uncommented
            //      if ((a.at(i,i) <= 0.) || (b.at(i,i) <= 0.))
            //        error ("solveYourselfAt: Matrices are not positive definite");
            eigv.at(i) = a.at(i, i) / b.at(i, i);
        }                                          // label 220

# ifdef DETAILED_REPORT
        OOFEM_LOG_DEBUG("GJacobi: current eigenvalues are:\n");
        eigv.printYourself();
        OOFEM_LOG_DEBUG("GJacobi: current eigenvectors are:\n");
        x.printYourself();
# endif
        //
        // check for convergence
        //
        for ( int i = 1; i <= n; i++ ) {       // label 230
            double tol = rtol * d.at(i);
            double dif = ( eigv.at(i) - d.at(i) );
            if ( fabs(dif) > tol ) {
                goto label280;
            }
        }                                 // label 240

        //
        // check all off-diagonal elements to see if another sweep is required
        //
        eps = rtol * rtol;
        for ( int j = 1; j <= nr; j++ ) {
            int jj = j + 1;