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

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

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

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

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

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
        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);
//        if(a2.computeSquaredNorm() < tol2) {
//            // The triangle is degenerated
//            return false;
//.........这里部分代码省略.........