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

int Line :: computeNumberOfIntersectionPoints(Element *element)
{

    int numIntersec = 0;
    const double relTol = 1.0e-3;
    const double LineLength = giveLength();
    const double absTol = relTol*std::max(LineLength, XfemTolerances::giveCharacteristicElementLength() );

    const int numEdges = element->giveInterpolation()->giveNumberOfEdges();

    for ( int edgeIndex = 1; edgeIndex <= numEdges; edgeIndex++ ) {
        IntArray bNodes;
        element->giveInterpolation()->boundaryGiveNodes(bNodes, edgeIndex);

        const int nsLoc = bNodes.at(1);
        const int neLoc = bNodes.at( bNodes.giveSize() );

        FloatArray xS = *(element->giveNode(nsLoc)->giveCoordinates() );
        xS.resizeWithValues(2);
        FloatArray xE = *(element->giveNode(neLoc)->giveCoordinates() );
        xE.resizeWithValues(2);

        const double dist = BasicGeometry :: computeLineDistance(xS, xE, mVertices[0], mVertices[1]);

        if(dist < absTol) {
            numIntersec++;
        }
    }

    return numIntersec;
}

int
StructuralInterfaceMaterial :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
{
    StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus(gp) );
    if ( type == IST_InterfaceJump ) {
        answer = status->giveJump();
        answer.resizeWithValues(3); // In case some model is not storing all components.
        return 1;
    } else if ( type == IST_InterfaceTraction ) {
        answer = status->giveTraction();
        answer.resizeWithValues(3);
        return 1;
    } else if ( type == IST_InterfaceFirstPKTraction ) {
        answer = status->giveFirstPKTraction();
        answer = status->giveTempFirstPKTraction();
        answer.resizeWithValues(3);
        return 1;
    } else if ( type == IST_DeformationGradientTensor ) {
        answer.beVectorForm( status->giveF() );
        return 1;
    } else if ( type == IST_InterfaceNormal ) {
        answer = status->giveNormal();
        return 1;
    } else {
        return Material :: giveIPValue(answer, gp, type, tStep);
    }
}

void DofManager :: giveUnknownVector(FloatArray &answer, const IntArray &dofIDArry,
                                     PrimaryField &field, ValueModeType mode, TimeStep *tStep, bool padding)
{
    answer.resize( dofIDArry.giveSize() );

    int k = 0;
    for ( auto &dofid: dofIDArry ) {
        auto pos = this->findDofWithDofId( ( DofIDItem ) dofid );
        if ( pos == this->end() ) {
            if ( padding ) {
                answer.at(++k) = 0.;
                continue;
            } else {
                continue;
            }
        }
        answer.at(++k) = (*pos)->giveUnknown(field, mode, tStep);
    }
    answer.resizeWithValues(k);

    // Transform to global c.s.
    FloatMatrix L2G;
    if ( this->computeL2GTransformation(L2G, dofIDArry) ) {
        answer.rotatedWith(L2G, 'n');
    }
}

double
Tetrah1_ht :: SpatialLocalizerI_giveDistanceFromParametricCenter(const FloatArray &coords)
{
    FloatArray lcoords(3), gcoords;
    double dist;
    int size, gsize;

    lcoords.zero();
    this->computeGlobalCoordinates(gcoords, lcoords);

    if ( ( size = coords.giveSize() ) < ( gsize = gcoords.giveSize() ) ) {
        OOFEM_ERROR("coordinates size mismatch");
    }

    if ( size == gsize ) {
        dist = coords.distance(gcoords);
    } else {
        FloatArray helpCoords = coords;

        helpCoords.resizeWithValues(gsize);
        dist = helpCoords.distance(gcoords);
    }

    return dist;
}

void
IntElPoint :: computeLocalSlipDir(FloatArray &normal)
{
    normal.resizeWithValues(3);
    if ( this->referenceNode ) {
        // normal
        normal.beDifferenceOf(*domain->giveNode(this->referenceNode)->giveCoordinates(), *this->giveNode(1)->giveCoordinates());
    } else {
        if ( normal.at(1) == 0 && normal.at(2) == 0 && normal.at(3) == 0 ) {
            OOFEM_ERROR("Normal is not defined (referenceNode=0,normal=(0,0,0))");
        }
    }

    normal.normalize();
}

int
TransportMaterial :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
// IST_Humidity must be overriden!
{
    TransportMaterialStatus *ms = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
    if ( type == IST_Temperature || type == IST_MassConcentration_1 || type == IST_Humidity ) {
        FloatArray vec = ms->giveField();
        answer = FloatArray{vec.at( ( type == IST_Temperature ) ? 1 : 2 ) };
        return 1;
    } else if ( type == IST_TemperatureFlow ) {
        TransportElement *transpElem = static_cast< TransportElement * >( gp->giveElement() );
        transpElem->computeFlow(answer, gp, tStep);
        return 1;
    } else if ( type == IST_Velocity ) { ///@todo Shouldn't be named velocity.. instead, "MassFlow" or something suitable like that.
        answer = ms->giveFlux();
        answer.resizeWithValues(3);
        return 1;
    } else if ( type == IST_PressureGradient ) {
        answer = ms->giveGradient();
        answer.resizeWithValues(3);
        return 1;
    } else if ( type == IST_Density ) {
        answer = FloatArray{ this->give('d', gp) };
        return 1;
    } else if ( type == IST_HeatCapacity ) {
        answer = FloatArray{ this->give('c', gp) };
        return 1;
    } else if ( type == IST_ThermalConductivityIsotropic ) {
        answer = FloatArray{ this->give('k', gp) };
        return 1;
    } else if ( type == IST_Maturity ) {
        answer = FloatArray{ ms->giveMaturity() };
        return 1;
    }
    return Material :: giveIPValue(answer, gp, type, tStep);
}

void GradDpElement :: computeNonlocalDegreesOfFreedom(FloatArray &answer, TimeStep *tStep)
{
    ///@todo Just use this instead! (should work, but I don't have any tests right now)
#if 0
    this->giveStructuralElement()->computeVectorOf({G_0}, VM_Total, tStep, answer);
#else
    StructuralElement *elem = this->giveStructuralElement();
    FloatArray u;
    answer.resize(nlSize);

    elem->computeVectorOf(VM_Total, tStep, u);
    u.resizeWithValues(totalSize);
    for  ( int i = 1; i <= nlSize; i++ ) {
        answer.at(i) = u.at( locK.at(i) );
    }
#endif
}

void
InterfaceElem1d :: computeLocalSlipDir(FloatArray &normal)
{
    normal.resizeWithValues(3);
    if ( this->referenceNode ) {
        // normal
        normal.at(1) = domain->giveNode(this->referenceNode)->giveCoordinate(1) - this->giveNode(1)->giveCoordinate(1);
        normal.at(2) = domain->giveNode(this->referenceNode)->giveCoordinate(2) - this->giveNode(1)->giveCoordinate(2);
        normal.at(3) = domain->giveNode(this->referenceNode)->giveCoordinate(3) - this->giveNode(1)->giveCoordinate(3);
    } else {
        if ( normal.at(1) == 0 && normal.at(2) == 0 && normal.at(3) == 0 ) {
            _error("computeLocalSlipDir: normal is not defined (referenceNode=0,normal=(0,0,0))");
        }
    }

    normal.normalize();
}

void PrescribedGradientBCWeak :: createTractionMesh(bool iEnforceCornerPeriodicity, int iNumSides)
{
	bool split_at_holes = true;

    const double l_s = mUC[0] - mLC[0];
    const double minPointDist = 1.0e-4*l_s;

	// Find holes intersecting the RVE boundary so that these can be excluded
	std::vector<FloatArray> holeCoordUnsorted, allCoordUnsorted;
	findHoleCoord(holeCoordUnsorted, allCoordUnsorted);


	// Add corner points
    holeCoordUnsorted.push_back( {mUC[0], mLC[1]} );
    allCoordUnsorted.push_back( {mUC[0], mLC[1]} );

    holeCoordUnsorted.push_back( {mUC[0], mUC[1]} );
    allCoordUnsorted.push_back( {mUC[0], mUC[1]} );

    holeCoordUnsorted.push_back( {mLC[0], mUC[1]} );
    allCoordUnsorted.push_back( {mLC[0], mUC[1]} );


    // Add crack-boundary intersections
    findCrackBndIntersecCoord(holeCoordUnsorted);



    // Add periodicity points
    findPeriodicityCoord(holeCoordUnsorted);


    // Sort arrays in terms of arc length along the RVE boundary
    std :: sort( holeCoordUnsorted.begin(), holeCoordUnsorted.end(), ArcPosSortFunction4( mLC, mUC, 1.0e-4 ) );
    std :: sort( allCoordUnsorted.begin(), allCoordUnsorted.end(), ArcPosSortFunction4( mLC, mUC, 1.0e-4 ) );


    // Remove points that are too close to each other
    removeClosePoints(holeCoordUnsorted, minPointDist);

    removeClosePoints(allCoordUnsorted, minPointDist);


	// Create two arrays of segments, where each array represents the coarsest possible traction
	// mesh on one side of the RVE
    ArcPosSortFunction4 arcPosFunc( mLC, mUC, 1.0e-4 );

    std :: vector< TracSegArray * > tracElNew0, tracElNew1;
    tracElNew0.push_back(new TracSegArray());
    tracElNew1.push_back(new TracSegArray());

    for(size_t i = 1; i < holeCoordUnsorted.size(); i++) {

    	FloatArray xS = holeCoordUnsorted[i-1];
    	xS.resizeWithValues(2);
    	FloatArray xE = holeCoordUnsorted[i];
    	xE.resizeWithValues(2);
    	const FloatArray xC = {0.5*(xS[0]+xE[0]), 0.5*(xS[1]+xE[1])};

    	if( arcPosFunc.calcArcPos(xC) < 2.*l_s) {
    		tracElNew0[0]->mInteriorSegments.push_back( Line(xS, xE) );
    	}
    	else {
    		tracElNew1[0]->mInteriorSegments.push_back( Line(xS, xE) );
    	}
    }

    // Remove segments located in holes
    removeSegOverHoles(*(tracElNew0[0]), 1.0e-4);
    removeSegOverHoles(*(tracElNew1[0]), 1.0e-4);

    // TODO: Refinement.


    if(split_at_holes) {
    	splitSegments(tracElNew0);
    	splitSegments(tracElNew1);
    }

	// Identify additional points that can be used to refine the traction mesh







    //////////////////////////////////////////////////
    // Create traction dofs
    int numNodes = domain->giveNumberOfDofManagers();
    int totNodesCreated = 0;


    // For each side (0 and 1), loop over over elements

    // We may always create the first node on the element.
    // For the linear approximation, it may need to be a slave node,
    // depending on which element it is.

    // For now, consider only piecewise constant approximations. Then,
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                residualStrength = E * e0 * exp(-1. * ( tempKappa - e0 ) / ef);
            } else {
                OOFEM_ERROR("Unknown softening type for cohesive crack model.");
            }
        }

        if ( status ) {
            status->setResidualTensileStrength(residualStrength);
        }

#endif
    } else {
        // damage grows
        tempKappa = kappa;

        FloatArray crackPlaneNormal(3);
        for ( int i = 1; i <= principalStress.giveSize(); i++ ) {
            crackPlaneNormal.at(i) = principalDir.at(i, 1);
        }
        this->initDamaged(tempKappa, crackPlaneNormal, gp, tStep);
        this->computeDamage(omega, tempKappa, gp);
    }

    answer.zero();

    if ( omega > 0. ) {
        tempNominalStress = tempEffectiveStress;

        if ( this->isotropic ) {
            // convert effective stress to nominal stress
            tempNominalStress.times(1. - omega);
            answer.add(tempNominalStress);
        } else {
            // stress transformation matrix
            FloatMatrix Tstress;

            // compute principal nominal stresses by multiplying effective stresses by damage
            for ( int i = 1; i <= principalStress.giveSize(); i++ ) {
                if ( principalStress.at(i) > 0. ) {
                    // convert principal effective stress to nominal stress
                    principalStress.at(i) *= ( 1. - omega );
                }
            }

            if ( mode == _PlaneStress ) {
                principalStress.resizeWithValues(3);
                givePlaneStressVectorTranformationMtrx(Tstress, principalDir, true);
            } else {
                principalStress.resizeWithValues(6);
                giveStressVectorTranformationMtrx(Tstress, principalDir, true);
            }

            principalStress.rotatedWith(Tstress, 'n');


            if ( mode == _PlaneStress ) { // plane stress
                answer.add(principalStress);
            } else if ( this->giveSizeOfVoigtSymVector(mode) != this->giveSizeOfVoigtSymVector(_3dMat) ) { // mode = _PlaneStrain or axial symmetry
                StressVector redFormStress(mode);
                redFormStress.convertFromFullForm(principalStress, mode);
                answer.add(redFormStress);
            } else { // 3D
                answer.add(principalStress);
            }
        }
    } else {
        answer.add(tempEffectiveStress);
    }

#ifdef supplementary_info

    if ( ( omega == 0. ) || ( sigma1 <= 0 ) ) {
        status->setCrackWidth(0.);
    } else {
        FloatArray principalStrains;
        this->computePrincipalValues(principalStrains, totalStrain, principal_strain);

        double crackWidth;
        //case when the strain localizes into narrow band and after a while the stresses relax

        double strainWithoutTemperShrink = principalStrains.at(1);
        strainWithoutTemperShrink -=  status->giveTempThermalStrain();
        strainWithoutTemperShrink -=  status->giveTempDryingShrinkageStrain();
        strainWithoutTemperShrink -=  status->giveTempAutogenousShrinkageStrain();


        crackWidth = status->giveCharLength() * omega * strainWithoutTemperShrink;

        status->setCrackWidth(crackWidth);
    }

#endif

    // update gp
    status->letTempStrainVectorBe(totalStrain);
    status->letTempStressVectorBe(answer);
    status->letTempViscoelasticStressVectorBe(tempEffectiveStress);
    status->setTempKappa(tempKappa);
    status->setTempDamage(omega);
}

void PolygonLine :: computeTangentialSignDist(double &oDist, const FloatArray &iPoint, double &oMinArcDist) const
{
	FloatArray point = iPoint;
	point.resizeWithValues(2);

    const int numSeg = this->giveNrVertices() - 1;

    double xi = 0.0, xiUnbounded = 0.0;

    if(numSeg == 0) {
    	FloatArray crackP1 = giveVertex ( 1 );
    	oDist = crackP1.distance(iPoint);
    	oMinArcDist = 0.0;
    	return;
    }

    if(numSeg == 1) {
    	FloatArray crackP1 = giveVertex ( 1 );
    	crackP1.resizeWithValues(2);
    	FloatArray crackP2 = giveVertex ( 2 );
    	crackP2.resizeWithValues(2);
        point.distance(crackP1, crackP2, xi, xiUnbounded);

        if( xiUnbounded < 0.0 ) {
            oDist = xiUnbounded*crackP1.distance(crackP2);
            oMinArcDist = 0.0;
            return;
        }

        if( xiUnbounded > 1.0 ) {
            oDist = -(xiUnbounded-1.0)*crackP1.distance(crackP2);
            oMinArcDist = 1.0;
            return;
        }

        const double L = computeLength();
        double distToStart  = xi*L;

        oDist = std::min(distToStart, (L - distToStart) );
        oMinArcDist = distToStart/L;
        return;
    }

    bool isBeforeStart = false, isAfterEnd = false;
    double distBeforeStart = 0.0, distAfterEnd = 0.0;

    ///////////////////////////////////////////////////////////////////
    // Check first segment
    FloatArray crackP1_start = giveVertex ( 1 );
    crackP1_start.resizeWithValues(2);
    FloatArray crackP2_start = giveVertex ( 2 );
    crackP2_start.resizeWithValues(2);
    const double distSeg_start = point.distance(crackP1_start, crackP2_start, xi, xiUnbounded);

    if( xiUnbounded < 0.0 ) {
        isBeforeStart = true;
        distBeforeStart = xiUnbounded*crackP1_start.distance(crackP2_start);
    }

    double arcPosPassed = crackP1_start.distance(crackP2_start);
    double distToStart  = xi*crackP1_start.distance(crackP2_start);

    double minGeomDist  = distSeg_start;




    ///////////////////////////////////////////////////////////////////
    // Check interior segments
    for ( int segId = 2; segId <= numSeg-1; segId++ ) {
    	FloatArray crackP1 = giveVertex ( segId      );
    	crackP1.resizeWithValues(2);
        FloatArray crackP2 = giveVertex ( segId+1    );
        crackP2.resizeWithValues(2);

        const double distSeg = point.distance(crackP1, crackP2, xi, xiUnbounded);

        if(distSeg < minGeomDist) {
            isBeforeStart = false;
            minGeomDist = distSeg;
            distToStart = arcPosPassed + xi*crackP1.distance(crackP2);
        }

        arcPosPassed += crackP1.distance(crackP2);

    }



    ///////////////////////////////////////////////////////////////////
    // Check last segment
    FloatArray crackP1_end = giveVertex ( numSeg );
    crackP1_end.resizeWithValues(2);
    FloatArray crackP2_end = giveVertex ( numSeg+1 );
    crackP2_end.resizeWithValues(2);
    const double distSeg_end = point.distance(crackP1_end, crackP2_end, xi, xiUnbounded);

    if(numSeg > 1) {
        if( xiUnbounded > 1.0 ) {
            arcPosPassed += xiUnbounded*crackP1_end.distance(crackP2_end);
//.........这里部分代码省略.........

bool Triangle :: pointIsInTriangle(const FloatArray &iP) const
{
    FloatArray P(iP);

    const double tol2 = 1.0e-18;

    // Compute triangle normal
    FloatArray p1p2;
    p1p2.beDifferenceOf(mVertices [ 1 ], mVertices [ 0 ]);

    FloatArray p1p3;
    p1p3.beDifferenceOf(mVertices [ 2 ], mVertices [ 0 ]);


    // Edge 1
    FloatArray t1;
    t1.beDifferenceOf(mVertices [ 1 ], mVertices [ 0 ]);
    if(t1.computeSquaredNorm() < tol2) {
        // The triangle is degenerated
        return false;
    }
    else {
        t1.normalize();
    }


    FloatArray a1;

    // Edge 2
    FloatArray t2;
    t2.beDifferenceOf(mVertices [ 2 ], mVertices [ 1 ]);
    if(t2.computeSquaredNorm() < tol2) {
        // The triangle is degenerated
        return false;
    }
    else {
        t2.normalize();
    }

    FloatArray a2;


    // Edge 3
    FloatArray t3;
    t3.beDifferenceOf(mVertices [ 0 ], mVertices [ 2 ]);
    if(t3.computeSquaredNorm() < tol2) {
        // The triangle is degenerated
        return false;
    }
    else {
        t3.normalize();
    }

    FloatArray a3;


    // Project point onto triangle plane
    FloatArray pProj = P;

    if( p1p2.giveSize() == 2 ) {
        // 2D

    	pProj.resizeWithValues(2);

        a1 = {-t1[1], t1[0]};
        a2 = {-t2[1], t2[0]};
        a3 = {-t3[1], t3[0]};
    }
    else {
        // 3D
        FloatArray N;

        N.beVectorProductOf(p1p2, p1p3);

        if(N.computeSquaredNorm() < tol2) {
            // The triangle is degenerated
            return false;
        }
        else {
            N.normalize();
        }

        // Compute normal distance from triangle to point
        FloatArray p1p;
        p1p.beDifferenceOf(P, mVertices [ 0 ]);
        double d = p1p.dotProduct(N);

        pProj.add(-d, N);


        a1.beVectorProductOf(N, t1);
//        if(a1.computeSquaredNorm() < tol2) {
//            // The triangle is degenerated
//            return false;
//        }
//        else {
//            a1.normalize();
//        }

        a2.beVectorProductOf(N, t2);
//.........这里部分代码省略.........