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

double FEI2dLineLin :: edgeEvalNormal(FloatArray &normal, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    normal.resize(2);
    normal.at(1) = cellgeo.giveVertexCoordinates(2)->at(xind) - cellgeo.giveVertexCoordinates(1)->at(xind);
    normal.at(2) = -(cellgeo.giveVertexCoordinates(2)->at(yind) - cellgeo.giveVertexCoordinates(1)->at(yind));
    return normal.normalize()*0.5;
}

void PolygonLine :: giveTangent(FloatArray &oTangent, const double &iArcPosition) const
{
    double L = computeLength();
    double xSegStart = 0.0, xSegEnd = 0.0;
    double xiSegStart = 0.0, xiSegEnd = 0.0;
    size_t numSeg = mVertices.size() - 1;
    const double xiTol = 1.0e-9;

    for ( size_t i = 0; i < numSeg; i++ ) {
        xSegEnd += mVertices [ i ].distance(mVertices [ i + 1 ]);

        xiSegStart      = xSegStart / L;
        xiSegEnd        = xSegEnd / L;

        if ( iArcPosition > xiSegStart-xiTol && iArcPosition < xiSegEnd+xiTol ) {
            // The given point is within the segment

            const FloatArray &p1 = mVertices [ i ];
            const FloatArray &p2 = mVertices [ i+1 ];

            oTangent = {p2(0) - p1(0), p2(1) - p1(1)};

            oTangent.normalize();

            return;
        }

    }

    OOFEM_ERROR("Arc position not found.")
}

double
FEI3dHexaLin :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    FloatArray a, b, dNdksi(4), dNdeta(4);
    double ksi, eta;
    IntArray snodes;
    
    this->computeLocalSurfaceMapping(snodes, isurf);

    ksi = lcoords.at(1);
    eta = lcoords.at(2);

    // No need to divide by 1/4, we'll normalize anyway;
    dNdksi.at(1) =  ( 1. + eta );
    dNdksi.at(2) = -( 1. + eta );
    dNdksi.at(3) = -( 1. - eta );
    dNdksi.at(4) =  ( 1. - eta );

    dNdeta.at(1) =  ( 1. + ksi );
    dNdeta.at(2) =  ( 1. - ksi );
    dNdeta.at(3) = -( 1. - ksi );
    dNdeta.at(4) = -( 1. + ksi );

    for (int i = 1; i <= 4; ++i) {
        a.add(dNdksi.at(i), *cellgeo.giveVertexCoordinates(snodes.at(i)));
        b.add(dNdeta.at(i), *cellgeo.giveVertexCoordinates(snodes.at(i)));
    }
    
    answer.beVectorProductOf(a, b);
    return answer.normalize()*0.0625;
}

void
DKTPlate :: computeMidPlaneNormal(FloatArray &answer, const GaussPoint *gp)
// returns normal vector to midPlane in GaussPoinr gp of receiver
{
    FloatArray u, v;
    u.beDifferenceOf( * this->giveNode(2)->giveCoordinates(), * this->giveNode(1)->giveCoordinates() );
    v.beDifferenceOf( * this->giveNode(3)->giveCoordinates(), * this->giveNode(1)->giveCoordinates() );

    answer.beVectorProductOf(u, v);
    answer.normalize();
}

void PrescribedGradientBCWeak :: giveTractionElNormal(size_t iElInd, FloatArray &oNormal, FloatArray &oTangent) const
{
	FloatArray xS, xE;
	giveTractionElCoord(iElInd, xS, xE);

    oTangent.beDifferenceOf(xE, xS);
    oTangent.normalize();

    oNormal = {
        oTangent [ 1 ], -oTangent [ 0 ]
    };
}

void
Lattice2d :: giveCrossSectionCoordinates(FloatArray &coords)
{
    double x1, y1, x2, y2;
    x1 = this->giveNode(1)->giveCoordinate(1);
    y1 = this->giveNode(1)->giveCoordinate(2);
    x2 = this->giveNode(2)->giveCoordinate(1);
    y2 = this->giveNode(2)->giveCoordinate(2);

    //Compute normal and shear direction
    FloatArray normalDirection;
    FloatArray shearDirection;
    normalDirection.resize(2);
    normalDirection.zero();
    shearDirection.resize(2);
    shearDirection.zero();
    normalDirection.at(1) = x2 - x1;
    normalDirection.at(2) = y2 - y1;
    normalDirection.normalize();
    if ( normalDirection.at(2) == 0. ) {
        shearDirection.at(1) = 0.;
        shearDirection.at(2) = 1.;
    } else {
        shearDirection.at(1) = 1.0;
        shearDirection.at(2) =
            -normalDirection.at(1) / normalDirection.at(2);
    }

    shearDirection.normalize();

    coords.resize(6);
    coords.at(1) = this->gpCoords.at(1) - shearDirection.at(1) * this->width / 2.;
    coords.at(2) = this->gpCoords.at(2) - shearDirection.at(2) * this->width / 2.;
    coords.at(3) = 0.;
    coords.at(4) = this->gpCoords.at(1) + shearDirection.at(1) * this->width / 2.;
    coords.at(5) = this->gpCoords.at(2) + shearDirection.at(2) * this->width / 2.;
    coords.at(6) = 0.;

    return;
}

double
FEI2dQuadLin :: edgeEvalNormal(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    int nodeA, nodeB;
    IntArray edgeNodes;
    this->computeLocalEdgeMapping(edgeNodes, iedge);
    nodeA = edgeNodes.at(1);
    nodeB = edgeNodes.at(2);

    answer.resize(2);
    answer.at(1) = -(cellgeo.giveVertexCoordinates(nodeB)->at(yind) - cellgeo.giveVertexCoordinates(nodeA)->at(yind) );
    answer.at(2) =  (cellgeo.giveVertexCoordinates(nodeB)->at(xind) - cellgeo.giveVertexCoordinates(nodeA)->at(xind) );
    return answer.normalize() * 0.5;
}

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

void
IntElLine1PhF :: computeTransformationMatrixAt(GaussPoint *gp, FloatMatrix &answer)
{
    // Transformation matrix to the local coordinate system
    FloatArray G;
    this->computeCovarBaseVectorAt(gp, G);
    G.normalize();

    answer.resize(2, 2);
    answer.at(1, 1) =  G.at(1);
    answer.at(2, 1) = -G.at(2);
    answer.at(1, 2) =  G.at(2);
    answer.at(2, 2) =  G.at(1);

}

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
Structural2DElement :: giveMaterialOrientationAt(FloatArray &x, FloatArray &y, const FloatArray &lcoords)
{
    if ( this->elemLocalCS.isNotEmpty() ) { // User specified orientation
        x = {
            elemLocalCS.at(1, 1), elemLocalCS.at(2, 1)
        };
        y = {
            -x(1), x(0)
        };
    } else {
        FloatMatrix jac;
        this->giveInterpolation()->giveJacobianMatrixAt( jac, lcoords, * this->giveCellGeometryWrapper() );
        x.beColumnOf(jac, 1); // This is {dx/dxi, dy/dxi, dz/dxi}
        x.normalize();
        y = {
            -x(1), x(0)
        };
    }
}

double FEI2dTrQuad :: edgeEvalNormal(FloatArray &normal, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    IntArray edgeNodes;
    this->computeLocalEdgeMapping(edgeNodes, iedge);
    double xi = lcoords(0);
    double dN1dxi = -0.5 + xi;
    double dN2dxi =  0.5 + xi;
    double dN3dxi = -2.0 * xi;

    normal.resize(2);

    normal.at(1) = dN1dxi * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(yind) +
    dN2dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(yind) +
    dN3dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(yind);

    normal.at(2) = -dN1dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(xind) +
    - dN2dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(xind) +
    - dN3dxi *cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(xind);

    return normal.normalize();
}

double
FEI3dHexaQuad :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    FloatArray a, b, dNdksi(8), dNdeta(8);
    double ksi, eta;
    IntArray snodes;
    
    this->computeLocalSurfaceMapping(snodes, isurf);

    ksi = lcoords.at(1);
    eta = lcoords.at(2);

    // No need to divide by 1/4, we'll normalize anyway;
    dNdksi.at(1) =  0.25 * ( 1. + eta ) * ( 2.0 * ksi + eta );
    dNdksi.at(2) = -0.25 * ( 1. + eta ) * ( -2.0 * ksi + eta );
    dNdksi.at(3) = -0.25 * ( 1. - eta ) * ( -2.0 * ksi - eta );
    dNdksi.at(4) =  0.25 * ( 1. - eta ) * ( 2.0 * ksi - eta );
    dNdksi.at(5) = -ksi * ( 1. + eta );
    dNdksi.at(6) = -0.5 * ( 1. - eta * eta );
    dNdksi.at(7) = -ksi * ( 1. - eta );
    dNdksi.at(8) =  0.5 * ( 1. - eta * eta );
    
    dNdeta.at(1) =  0.25 * ( 1. + ksi ) * ( 2.0 * eta + ksi );
    dNdeta.at(2) =  0.25 * ( 1. - ksi ) * ( 2.0 * eta - ksi );
    dNdeta.at(3) = -0.25 * ( 1. - ksi ) * ( -2.0 * eta - ksi );
    dNdeta.at(4) = -0.25 * ( 1. + ksi ) * ( -2.0 * eta + ksi );
    dNdeta.at(5) =  0.5 * ( 1. - ksi * ksi );
    dNdeta.at(6) = -eta * ( 1. - ksi );
    dNdeta.at(7) = -0.5 * ( 1. - ksi * ksi );
    dNdeta.at(8) = -eta * ( 1. + ksi );

    for ( int i = 1; i <= 8; ++i ) {
        a.add(dNdksi.at(i), *cellgeo.giveVertexCoordinates(snodes.at(i)));
        b.add(dNdeta.at(i), *cellgeo.giveVertexCoordinates(snodes.at(i)));
    }
    
    answer.beVectorProductOf(a, b);
    return answer.normalize();
}

void PLCrackPrescribedDir::propagateInterfaces(EnrichmentDomain &ioEnrDom) {
	printf("Entering PLCrackPrescribedDir::propagateInterfaces().\n");

	// Fetch crack tip data
	std::vector<TipInfo> tipInfo;
	ioEnrDom.giveTipInfos(tipInfo);

	int tipIndex = 1;
	FloatArray dir;
	double angleRad = mAngle*M_PI/180.0;
	dir.setValues(2, cos(angleRad), sin(angleRad));
	dir.normalize();


	std::vector<TipPropagation> tipPropagations;
	TipPropagation tipProp;
	tipProp.mTipIndex = tipIndex;
	tipProp.mPropagationDir = dir;
	tipProp.mPropagationLength = mIncrementLength;
	tipPropagations.push_back(tipProp);

	ioEnrDom.propagateTips(tipPropagations);
}

void
LEPlic :: doCellDLS(FloatArray &fvgrad, int ie, bool coord_upd, bool vof_temp_flag)
{
    int i, ineighbr, nneighbr;
    double fvi, fvk, wk, dx, dy;
    bool isBoundaryCell = false;
    LEPlicElementInterface *interface, *ineghbrInterface;
    FloatMatrix lhs(2, 2);
    FloatArray rhs(2), xi(2), xk(2);
    IntArray currCell(1), neighborList;
    ConnectivityTable *contable = domain->giveConnectivityTable();

    if ( ( interface = ( LEPlicElementInterface * ) ( domain->giveElement(ie)->giveInterface(LEPlicElementInterfaceType) ) ) ) {
        if ( vof_temp_flag ) {
            fvi = interface->giveTempVolumeFraction();
        } else {
            fvi = interface->giveVolumeFraction();
        }

        if ( ( fvi > 0. ) && ( fvi <= 1.0 ) ) {
            // potentially boundary cell

            if ( ( fvi > 0. ) && ( fvi < 1.0 ) ) {
                isBoundaryCell = true;
            }

            /* DLS (Differential least square reconstruction)
             *
             * In the DLS method, volume fraction Taylor series expansion of vf (volume fraction)
             * is formed from each reference cell volume fraction vf at element center x(i) to each
             * cell neighbor at point x(k). The sum (vf(i)-vf(k))^2 over all immediate neighbors
             * is then minimized inthe least square sense.
             */
            // get list of neighbours to current cell including current cell
            currCell.at(1) = ie;
            contable->giveElementNeighbourList(neighborList, currCell);
            // loop over neighbors to assemble normal equations
            nneighbr = neighborList.giveSize();
            interface->giveElementCenter(this, xi, coord_upd);
            lhs.zero();
            rhs.zero();
            for ( i = 1; i <= nneighbr; i++ ) {
                ineighbr = neighborList.at(i);
                if ( ineighbr == ie ) {
                    continue;         // skip itself
                }

                if ( ( ineghbrInterface =
                          ( LEPlicElementInterface * ) ( domain->giveElement(ineighbr)->giveInterface(LEPlicElementInterfaceType) ) ) ) {
                    if ( vof_temp_flag ) {
                        fvk = ineghbrInterface->giveTempVolumeFraction();
                    } else {
                        fvk = ineghbrInterface->giveVolumeFraction();
                    }

                    if ( fvk < 1.0 ) {
                        isBoundaryCell = true;
                    }

                    ineghbrInterface->giveElementCenter(this, xk, coord_upd);
                    wk = xk.distance(xi);
                    dx = ( xk.at(1) - xi.at(1) ) / wk;
                    dy = ( xk.at(2) - xi.at(2) ) / wk;
                    lhs.at(1, 1) += dx * dx;
                    lhs.at(1, 2) += dx * dy;
                    lhs.at(2, 2) += dy * dy;

                    rhs.at(1) += ( fvi - fvk ) * dx / wk;
                    rhs.at(2) += ( fvi - fvk ) * dy / wk;
                }
            }

            if ( isBoundaryCell ) {
                // symmetry
                lhs.at(2, 1) = lhs.at(1, 2);

                // solve normal equation for volume fraction gradient
                lhs.solveForRhs(rhs, fvgrad);

                // compute unit normal
                fvgrad.normalize();
                fvgrad.negated();
#ifdef __OOFEG
                /*
                 * EASValsSetLayer(OOFEG_DEBUG_LAYER);
                 * WCRec p[2];
                 * double tx = -fvgrad.at(2), ty=fvgrad.at(1);
                 * p[0].x=xi.at(1)-tx*0.1;
                 * p[0].y=xi.at(2)-ty*0.1;
                 * p[1].x=xi.at(1)+tx*0.1;
                 * p[1].y=xi.at(2)+ty*0.1;
                 * p[0].z = p[1].z = 0.0;
                 * GraphicObj *go = CreateLine3D(p);
                 * EGWithMaskChangeAttributes(LAYER_MASK, go);
                 * EMAddGraphicsToModel(ESIModel(), go);
                 * ESIEventLoop (YES, "Cell DLS finished; Press Ctrl-p to continue");
                 */
#endif
            } else {
                fvgrad.zero();
//.........这里部分代码省略.........